From 68c606a81c6ae2f6a6046bb4f292ef5f46222aac Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 5 Dec 2019 15:28:06 -0500 Subject: [PATCH] Refactored odm_orthophoto to use GDAL for TIFF writing Former-commit-id: f484cc6771adb983a4f5acbed04358554f4fcabd --- modules/odm_orthophoto/CMakeLists.txt | 4 +- modules/odm_orthophoto/src/OdmOrthoPhoto.cpp | 117 ++++++++++++++----- modules/odm_orthophoto/src/OdmOrthoPhoto.hpp | 12 +- opendm/orthophoto.py | 14 ++- opendm/types.py | 48 ++++++++ stages/odm_orthophoto.py | 2 +- tests/test_remote.py | 3 +- 7 files changed, 161 insertions(+), 39 deletions(-) diff --git a/modules/odm_orthophoto/CMakeLists.txt b/modules/odm_orthophoto/CMakeLists.txt index 6494abbb..4d8e7004 100644 --- a/modules/odm_orthophoto/CMakeLists.txt +++ b/modules/odm_orthophoto/CMakeLists.txt @@ -11,6 +11,8 @@ add_definitions(-Wall -Wextra) # Find pcl at the location specified by PCL_DIR find_package(VTK 6.0 REQUIRED) find_package(PCL 1.8 HINTS "${PCL_DIR}/share/pcl-1.8" REQUIRED) +find_package(GDAL REQUIRED) +include_directories(${GDAL_INCLUDE_DIR}) # Find OpenCV at the default location find_package(OpenCV HINTS "${OPENCV_DIR}" REQUIRED) @@ -31,4 +33,4 @@ aux_source_directory("./src" SRC_LIST) # Add exectuteable add_executable(${PROJECT_NAME} ${SRC_LIST}) -target_link_libraries(odm_orthophoto ${PCL_COMMON_LIBRARIES} ${PCL_IO_LIBRARIES} ${PCL_SURFACE_LIBRARIES} ${OpenCV_LIBS}) \ No newline at end of file +target_link_libraries(odm_orthophoto ${PCL_COMMON_LIBRARIES} ${PCL_IO_LIBRARIES} ${PCL_SURFACE_LIBRARIES} ${OpenCV_LIBS} ${GDAL_LIBRARY}) diff --git a/modules/odm_orthophoto/src/OdmOrthoPhoto.cpp b/modules/odm_orthophoto/src/OdmOrthoPhoto.cpp index ee13056c..8cd892ab 100644 --- a/modules/odm_orthophoto/src/OdmOrthoPhoto.cpp +++ b/modules/odm_orthophoto/src/OdmOrthoPhoto.cpp @@ -306,6 +306,50 @@ void OdmOrthoPhoto::printHelp() log_.setIsPrintingInCout(false); } +void OdmOrthoPhoto::saveTIFF(const std::string &filename, GDALDataType dataType){ + GDALAllRegister(); + GDALDriverH hDriver = GDALGetDriverByName( "GTiff" ); + if (!hDriver){ + std::cerr << "Cannot initialize GeoTIFF driver. Check your GDAL installation." << std::endl; + exit(1); + } + char **papszOptions = NULL; + GDALDatasetH hDstDS = GDALCreate( hDriver, filename.c_str(), width, height, static_cast(bands.size()), dataType, papszOptions ); + GDALRasterBandH hBand; + + for (size_t i = 0; i < bands.size(); i++){ + hBand = GDALGetRasterBand( hDstDS, static_cast(i) + 1 ); + if (GDALRasterIO( hBand, GF_Write, 0, 0, width, height, + bands[i], width, height, dataType, 0, 0 ) != CE_None){ + std::cerr << "Cannot write TIFF to " << filename << std::endl; + exit(1); + } +// for (int j = 0; j < height; j++){ +// GDALRasterIO( hBand, GF_Write, 0, j, width, 1, +// bands[i][j], width, 1, dataType, 0, 0 ); +// } + } + GDALClose( hDstDS ); +} + +template +inline T maxRange(){ + return static_cast(pow(2, sizeof(T) * 8) - 1); +} + +template +void OdmOrthoPhoto::initBands(const cv::Mat &texture){ + // Channels + alpha + for (int i = 0; i < texture.channels() + 1; i++){ + size_t pixelCount = static_cast(width * height); + T *arr = new T[pixelCount]; + for (size_t j = 0; j < pixelCount; j++){ + arr[j] = i == 4 ? 0 : maxRange(); // Set alpha at 0 + } + bands.push_back(static_cast(arr)); + } +} + void OdmOrthoPhoto::createOrthoPhoto() { if(inputFile_.empty()) @@ -333,7 +377,7 @@ void OdmOrthoPhoto::createOrthoPhoto() log_ << "\tSpecified -inputGeoRefFile, but no boundary points. The georeference system will be ignored.\n"; } - log_ << "Reading mesh file...\n"; + log_ << "Reading mesh file... " << inputFile_ << "\n"; // The textured mesh. pcl::TextureMesh mesh; loadObjFile(inputFile_, mesh); @@ -375,24 +419,24 @@ void OdmOrthoPhoto::createOrthoPhoto() log_ << "Ortho photo area : " << xDiff*yDiff << "m2\n"; // The resolution necessary to fit the area with the given resolution. - int rowRes = static_cast(std::ceil(resolution_*yDiff)); - int colRes = static_cast(std::ceil(resolution_*xDiff)); - log_ << "Ortho photo resolution, width x height : " << colRes << "x" << rowRes << '\n'; + height = static_cast(std::ceil(resolution_*yDiff)); + width = static_cast(std::ceil(resolution_*xDiff)); + log_ << "Ortho photo resolution, width x height : " << width << "x" << height << '\n'; // Check size of photo. - if(0 >= rowRes*colRes) + if(0 >= height*width) { - if(0 >= rowRes) + if(0 >= height) { - log_ << "Warning: ortho photo has zero area, height = " << rowRes << ". Forcing height = 1.\n"; - rowRes = 1; + log_ << "Warning: ortho photo has zero area, height = " << height << ". Forcing height = 1.\n"; + height = 1; } - if(0 >= colRes) + if(0 >= width) { - log_ << "Warning: ortho photo has zero area, width = " << colRes << ". Forcing width = 1.\n"; - colRes = 1; + log_ << "Warning: ortho photo has zero area, width = " << width << ". Forcing width = 1.\n"; + width = 1; } - log_ << "New ortho photo resolution, width x height : " << colRes << "x" << rowRes << '\n'; + log_ << "New ortho photo resolution, width x height : " << width << "x" << height << '\n'; } // Contains the vertices of the mesh. @@ -493,20 +537,22 @@ void OdmOrthoPhoto::createOrthoPhoto() if (textureDepth == -1){ try{ textureDepth = texture.depth(); - log_ << "Texture depth is " << textureDepth << "bit\n"; + log_ << "Texture channels: " << texture.channels() << "\n"; if (textureDepth == CV_8U){ - photo_ = cv::Mat::zeros(rowRes, colRes, CV_8UC4) + cv::Scalar(255, 255, 255, 0); + log_ << "Texture depth: 8bit\n"; + initBands(texture); }else if (textureDepth == CV_16U){ - photo_ = cv::Mat::zeros(rowRes, colRes, CV_16UC4) + cv::Scalar(65535, 65535, 65535, 0); + log_ << "Texture depth: 16bit\n"; + initBands(texture); }else{ std::cerr << "Unsupported bit depth value: " << textureDepth; exit(1); } - depth_ = cv::Mat::zeros(rowRes, colRes, CV_32F) - std::numeric_limits::infinity(); + depth_ = cv::Mat::zeros(height, width, CV_32F) - std::numeric_limits::infinity(); }catch(const cv::Exception &e){ - std::cerr << "Couldn't allocate enough memory to render the orthophoto (" << colRes << "x" << rowRes << " cells = " << ((long long)colRes * (long long)rowRes * 4) << " bytes). Try to increase the --orthophoto-resolution parameter to a larger integer or add more RAM.\n"; + std::cerr << "Couldn't allocate enough memory to render the orthophoto (" << width << "x" << height << " cells = " << ((long long)width * (long long)height * 4) << " bytes). Try to increase the --orthophoto-resolution parameter to a larger integer or add more RAM.\n"; exit(1); } } @@ -531,9 +577,9 @@ void OdmOrthoPhoto::createOrthoPhoto() // ... and draw it into the ortho photo. if (textureDepth == CV_8U){ - drawTexturedTriangle(texture, polygon, meshCloud, uvs, faceIndex+faceOff); + drawTexturedTriangle(texture, polygon, meshCloud, uvs, faceIndex+faceOff); }else if (textureDepth == CV_16U){ - drawTexturedTriangle(texture, polygon, meshCloud, uvs, faceIndex+faceOff); + drawTexturedTriangle(texture, polygon, meshCloud, uvs, faceIndex+faceOff); } } faceOff += faces.size(); @@ -543,7 +589,17 @@ void OdmOrthoPhoto::createOrthoPhoto() log_ << '\n'; log_ << "Writing ortho photo to " << outputFile_ << "\n"; - cv::imwrite(outputFile_, photo_); + + if (textureDepth == CV_8U){ + saveTIFF(outputFile_, GDT_Byte); + }else if (textureDepth == CV_16U){ + saveTIFF(outputFile_, GDT_UInt16); + }else{ + std::cerr << "Unsupported bit depth value: " << textureDepth; + exit(1); + } + +// cv::imwrite(outputFile_, photo_); if (!outputCornerFile_.empty()) { @@ -801,7 +857,7 @@ void OdmOrthoPhoto::drawTexturedTriangle(const cv::Mat &texture, const pcl::Vert // Check bounding box overlap. int xMin = static_cast(std::min(std::min(v1x, v2x), v3x)); - if(xMin > photo_.cols) + if(xMin > width) { return; // Completely outside to the right. } @@ -811,7 +867,7 @@ void OdmOrthoPhoto::drawTexturedTriangle(const cv::Mat &texture, const pcl::Vert return; // Completely outside to the left. } int yMin = static_cast(std::min(std::min(v1y, v2y), v3y)); - if(yMin > photo_.rows) + if(yMin > height) { return; // Completely outside to the top. } @@ -915,7 +971,7 @@ void OdmOrthoPhoto::drawTexturedTriangle(const cv::Mat &texture, const pcl::Vert // The first pixel row for the bottom part of the triangle. int rqStart = std::max(static_cast(std::floor(topR+0.5f)), 0); // The last pixel row for the top part of the triangle. - int rqEnd = std::min(static_cast(std::floor(midR+0.5f)), photo_.rows); + int rqEnd = std::min(static_cast(std::floor(midR+0.5f)), height); // Traverse along row from top to middle. for(int rq = rqStart; rq < rqEnd; ++rq) @@ -927,7 +983,7 @@ void OdmOrthoPhoto::drawTexturedTriangle(const cv::Mat &texture, const pcl::Vert // The first pixel column for the current row. int cqStart = std::max(static_cast(std::floor(0.5f+std::min(ctm, ctb))), 0); // The last pixel column for the current row. - int cqEnd = std::min(static_cast(std::floor(0.5f+std::max(ctm, ctb))), photo_.cols); + int cqEnd = std::min(static_cast(std::floor(0.5f+std::max(ctm, ctb))), width); for(int cq = cqStart; cq < cqEnd; ++cq) { @@ -973,7 +1029,7 @@ void OdmOrthoPhoto::drawTexturedTriangle(const cv::Mat &texture, const pcl::Vert // The first pixel row for the bottom part of the triangle. int rqStart = std::max(static_cast(std::floor(midR+0.5f)), 0); // The last pixel row for the bottom part of the triangle. - int rqEnd = std::min(static_cast(std::floor(botR+0.5f)), photo_.rows); + int rqEnd = std::min(static_cast(std::floor(botR+0.5f)), height); // Traverse along row from middle to bottom. for(int rq = rqStart; rq < rqEnd; ++rq) @@ -985,7 +1041,7 @@ void OdmOrthoPhoto::drawTexturedTriangle(const cv::Mat &texture, const pcl::Vert // The first pixel column for the current row. int cqStart = std::max(static_cast(std::floor(0.5f+std::min(cmb, ctb))), 0); // The last pixel column for the current row. - int cqEnd = std::min(static_cast(std::floor(0.5f+std::max(cmb, ctb))), photo_.cols); + int cqEnd = std::min(static_cast(std::floor(0.5f+std::max(cmb, ctb))), width); for(int cq = cqStart; cq < cqEnd; ++cq) { @@ -1071,10 +1127,11 @@ void OdmOrthoPhoto::renderPixel(int row, int col, float s, float t, const cv::Ma b += static_cast(bl[0]) * dr * dt; b += static_cast(br[0]) * dl * dt; - photo_.at >(row,col) = cv::Vec(static_cast(b), - static_cast(g), - static_cast(r), - static_cast(255)); // Alpha should always be in the 255 range + size_t idx = static_cast(row * width + col); + static_cast(bands[0])[idx] = static_cast(r); + static_cast(bands[1])[idx] = static_cast(g); + static_cast(bands[2])[idx] = static_cast(b); + static_cast(bands[3])[idx] = static_cast(255); // Alpha should always be in the 255 range } void OdmOrthoPhoto::getBarycentricCoordinates(pcl::PointXYZ v1, pcl::PointXYZ v2, pcl::PointXYZ v3, float x, float y, float &l1, float &l2, float &l3) const diff --git a/modules/odm_orthophoto/src/OdmOrthoPhoto.hpp b/modules/odm_orthophoto/src/OdmOrthoPhoto.hpp index 35e782b2..6e9f07aa 100644 --- a/modules/odm_orthophoto/src/OdmOrthoPhoto.hpp +++ b/modules/odm_orthophoto/src/OdmOrthoPhoto.hpp @@ -20,6 +20,10 @@ // OpenCV #include +// GDAL +#include "gdal_priv.h" +#include "cpl_conv.h" // for CPLMalloc() + // Logger #include "Logger.hpp" @@ -76,6 +80,7 @@ public: int run(int argc, char* argv[]); private: + int width, height; /*! * \brief parseArguments Parses command line arguments. @@ -123,6 +128,10 @@ private: * @return */ Eigen::Transform readTransform(std::string transformFile_) const; + + template + void initBands(const cv::Mat &texture); + void saveTIFF(const std::string &filename, GDALDataType dataType); /*! * \brief Renders a triangle into the ortho photo. @@ -222,7 +231,8 @@ private: Eigen::Vector2f boundaryPoint3_; /**< The third boundary point for the ortho photo, in local coordinates. */ Eigen::Vector2f boundaryPoint4_; /**< The fourth boundary point for the ortho photo, in local coordinates. */ - cv::Mat photo_; /**< The ortho photo as an OpenCV matrix, CV_8UC3. */ + std::vector bands; + cv::Mat depth_; /**< The depth of the ortho photo as an OpenCV matrix, CV_32F. */ bool multiMaterial_; /**< True if the mesh has multiple materials. **/ diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py index 26d442d8..16a3841e 100644 --- a/opendm/orthophoto.py +++ b/opendm/orthophoto.py @@ -82,11 +82,15 @@ def compute_mask_raster(input_raster, vector_mask, output_raster, blend_distance out_image, out_transform = mask(rast, shapes, nodata=0) if blend_distance > 0: - alpha_band = out_image[3] - dist_t = ndimage.distance_transform_edt(alpha_band) - dist_t[dist_t <= blend_distance] /= blend_distance - dist_t[dist_t > blend_distance] = 1 - np.multiply(alpha_band, dist_t, out=alpha_band, casting="unsafe") + if out_image.shape[0] >= 4: + # rast_mask = rast.dataset_mask() + rast_mask = out_image[-1] + dist_t = ndimage.distance_transform_edt(rast_mask) + dist_t[dist_t <= blend_distance] /= blend_distance + dist_t[dist_t > blend_distance] = 1 + np.multiply(rast_mask, dist_t, out=rast_mask, casting="unsafe") + else: + log.ODM_WARNING("%s does not have an alpha band, cannot blend cutline!" % input_raster) with rasterio.open(output_raster, 'w', **rast.profile) as dst: dst.write(out_image) diff --git a/opendm/types.py b/opendm/types.py index ef8f739f..88c15fb3 100644 --- a/opendm/types.py +++ b/opendm/types.py @@ -93,6 +93,54 @@ class ODM_Reconstruction(object): self.photos = photos self.georef = None self.gcp = None + self.multi_camera = self.detect_multi_camera() + + def detect_multi_camera(self): + """ + Looks at the reconstruction photos and determines if this + is a single or multi-camera setup. + """ + supported_ext_re = r"\.(" + "|".join([e[1:] for e in context.supported_extensions]) + ")$" + + # Match filename_1.tif, filename_2.tif, filename_3.tif, ... + # + multi_camera_patterns = { + 'MicaSense RedEdge-M': r'^IMG_\d+_(?P\d{1})', + 'Parrot Sequoia': r'^IMG_\d+_\d+_\d+_(?P[A-Z]{3})', + } + + for cam, regex in multi_camera_patterns.items(): + pattern = re.compile(regex + supported_ext_re, re.IGNORECASE) + mc = {} + + for p in self.photos: + matches = re.match(pattern, p.filename) + if matches: + band = matches.group("band") + if not band in mc: + mc[band] = [] + mc[band].append(p) + + # We support between 2 and 6 bands + # If we matched more or less bands, we probably just + # found filename patterns that do not match a multi-camera setup + bands_count = len(mc) + if bands_count >= 2 and bands_count <= 6: + + # Validate that all bands have the same number of images, + # otherwise this is not a multi-camera setup + img_per_band = len(mc[band]) + valid = True + for band in mc: + if len(mc[band]) != img_per_band: + log.ODM_WARNING("This might be a multi-camera setup, but band \"%s\" (identified from \"%s\") has only %s images (instead of %s), perhaps images are missing or are corrupted." % (band, mc[band][0].filename, len(mc[band]), img_per_band)) + valid = False + break + + if valid: + return mc + + return None def is_georeferenced(self): return self.georef is not None diff --git a/stages/odm_orthophoto.py b/stages/odm_orthophoto.py index 9321632d..9bb574d2 100644 --- a/stages/odm_orthophoto.py +++ b/stages/odm_orthophoto.py @@ -97,7 +97,7 @@ class ODMOrthoPhotoStage(types.ODM_Stage): system.run('gdal_translate -a_ullr {ulx} {uly} {lrx} {lry} ' '{vars} ' - '-b 1 -b 2 -b 3 -mask 4 ' + # '-b 1 -b 2 -b 3 -mask 4 ' '-a_srs \"{proj}\" ' '--config GDAL_CACHEMAX {max_memory}% ' '--config GDAL_TIFF_INTERNAL_MASK YES ' diff --git a/tests/test_remote.py b/tests/test_remote.py index 8fa22783..5483bd0d 100644 --- a/tests/test_remote.py +++ b/tests/test_remote.py @@ -35,10 +35,11 @@ class TestRemote(unittest.TestCase): self.queue_num = queue_num self.uuid = 'xxxxx-xxxxx-xxxxx-xxxxx-xxxx' + str(queue_num) - def info(self): + def info(self, with_output=None): class StatusMock: status = TaskStatus.RUNNING if self.running else TaskStatus.QUEUED processing_time = 1 + output = "test output" return StatusMock() def remove(self):