Blocks depth auto selection, DSM resolution based on GSD estimate

Former-commit-id: d0a9374e9d

modules/odm_dem2mesh/src/main.cpp

@@ -53,12 +53,13 @@ cmdLineParameter< char* >
     InputFile( "inputFile" ) ,
     OutputFile( "outputFile" );
 cmdLineParameter< int >
-    MaxVertexCount( "maxVertexCount" );
+    MaxVertexCount( "maxVertexCount" ) ,
+    MaxTileLength( "maxTileLength" );
 cmdLineReadable
     Verbose( "verbose" );
 
 cmdLineReadable* params[] = {
-    &InputFile , &OutputFile , &MaxVertexCount , &Verbose ,
+    &InputFile , &OutputFile , &MaxVertexCount , &MaxTileLength, &Verbose ,
     NULL
 };
 
@@ -67,6 +68,7 @@ void help(char *ex){
               << "\t -" << InputFile.name << " <input DSM raster>" << std::endl
               << "\t -" << OutputFile.name << " <output PLY mesh>" << std::endl
               << "\t [-" << MaxVertexCount.name << " <target number vertices> (Default: 100000)]" << std::endl
+              << "\t [-" << MaxTileLength.name << " <max length of a tile. Smaller values take longer to process but reduce memory usage by splitting the meshing process into tiles.> (Default: 1000)]" << std::endl
               << "\t [-" << Verbose.name << "]" << std::endl;
     exit(EXIT_FAILURE);
 }
@@ -282,6 +284,7 @@ int main(int argc, char **argv) {
     cmdLineParse( argc-1 , &argv[1] , params );
     if( !InputFile.set || !OutputFile.set ) help(argv[0]);
     if ( !MaxVertexCount.set ) MaxVertexCount.value = 100000;
+    if ( !MaxTileLength.set ) MaxTileLength.value = 1000;
 
     logWriter.verbose = Verbose.set;
     logWriter.outputFile = "odm_dem2mesh.txt";
@@ -308,14 +311,24 @@ int main(int argc, char **argv) {
 
         GDALRasterBand *band = dataset->GetRasterBand(1);
 
-        int qtreeLevels = 1;
-        subdivisions = (int)pow(2, qtreeLevels);
-        int numBlocks = subdivisions * subdivisions;
-        blockSizeX = arr_width / subdivisions;
-        blockSizeY = arr_height / subdivisions;
+        int qtreeLevels = 0;
+        int numBlocks = 1;
+        while(true){
+            subdivisions = (int)pow(2, qtreeLevels);
+            numBlocks = subdivisions * subdivisions;
+            blockSizeX = arr_width / subdivisions;
+            blockSizeY = arr_height / subdivisions;
+            if (blockSizeX > MaxTileLength.value || blockSizeY > MaxTileLength.value){
+                qtreeLevels++;
+            }else{
+                break;
+            }
+        }
+
         int blockXPad = 0; // Blocks > 0 need to re-add a column for seamless meshing
         int blockYPad = 0; // Blocks > 0 need to re-add a row for seamless meshing
 
+        logWriter("Blocks depth: %d\n", qtreeLevels);
         logWriter("Splitting area in %d\n", numBlocks);
         logWriter("Block size is %d, %d\n", blockSizeX, blockSizeY);
 

opendm/gsd.py

@@ -4,6 +4,22 @@ import numpy as np
 from repoze.lru import lru_cache
 from opendm import log
 
+def rounded_gsd(reconstruction_json, default_value=None, ndigits=0, ignore_gsd=False):
+    """
+    :param reconstruction_json path to OpenSfM's reconstruction.json
+    :return GSD value rounded. If GSD cannot be computed, or ignore_gsd is set, it returns a default value.
+    """
+    if ignore_gsd:
+        return default_value
+
+    gsd = opensfm_reconstruction_average_gsd(reconstruction_json)
+
+    if gsd is not None:
+        return round(gsd, ndigits)
+    else:
+        return default_value
+
+        
 def image_scale_factor(target_resolution, reconstruction_json, gsd_error_estimate = 0.1):
     """
     :param target_resolution resolution the user wants have in cm / pixel

scripts/odm_meshing.py

@@ -87,9 +87,15 @@ class ODMeshingCell(ecto.Cell):
               # a larger radius interolation --> less holes
               if args.fast_orthophoto:
                   dsm_radius *= 2
+
+              dsm_multiplier = max(1.0, gsd.rounded_gsd(tree.opensfm_reconstruction, default_value=4, ndigits=3, ignore_gsd=args.ignore_gsd))
               
-              # A good DSM size is 1/2 of the target orthophoto resolution
-              dsm_resolution = ortho_resolution * 2
+              # A good DSM size depends on the flight altitude.
+              # Flights at low altitude need more details (higher resolution) 
+              # Flights at higher altitude benefit from smoother surfaces (lower resolution)
+              dsm_resolution = ortho_resolution * dsm_multiplier
+
+              log.ODM_DEBUG('ODM 2.5D DSM resolution: %s' % dsm_resolution)
 
               mesh.create_25dmesh(infile, tree.odm_25dmesh,
                     dsm_radius=dsm_radius,