Merge pull request #1610 from pierotofy/305

Point spacing estimation

SuperBuild/cmake/External-FPCFilter.cmake

@@ -8,7 +8,7 @@ ExternalProject_Add(${_proj_name}
   #--Download step--------------
   DOWNLOAD_DIR      ${SB_DOWNLOAD_DIR}
   GIT_REPOSITORY    https://github.com/OpenDroneMap/FPCFilter
-  GIT_TAG           main
+  GIT_TAG           305
   #--Update/Patch step----------
   UPDATE_COMMAND    ""
   #--Configure step-------------

opendm/dem/commands.py

@@ -371,3 +371,11 @@ def window_filter_2d(arr, nodata, window, kernel_size, filter):
     win_arr[nodata_locs] = nodata
     win_arr = win_arr[window[0] - expanded_window[0] : window[2] - expanded_window[0], window[1] - expanded_window[1] : window[3] - expanded_window[1]]
     return win_arr
+
+
+def get_dem_radius_steps(stats_file, steps, resolution, multiplier = 1.0):
+    radius_steps = [point_cloud.get_spacing(stats_file, resolution) * multiplier]
+    for _ in range(steps - 1):
+        radius_steps.append(radius_steps[-1] * math.sqrt(2))
+    
+    return radius_steps

opendm/mesh.py

@@ -5,10 +5,11 @@ from opendm import system
 from opendm import log
 from opendm import context
 from opendm import concurrency
+from opendm import point_cloud
 from scipy import signal
 import numpy as np
 
-def create_25dmesh(inPointCloud, outMesh, dsm_radius=0.07, dsm_resolution=0.05, depth=8, samples=1, maxVertexCount=100000, available_cores=None, method='gridded', smooth_dsm=True):
+def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=0.05, depth=8, samples=1, maxVertexCount=100000, available_cores=None, method='gridded', smooth_dsm=True):
     # Create DSM from point cloud
 
     # Create temporary directory
@@ -19,17 +20,13 @@ def create_25dmesh(inPointCloud, outMesh, dsm_radius=0.07, dsm_resolution=0.05,
     os.mkdir(tmp_directory)
     log.ODM_INFO('Created temporary directory: %s' % tmp_directory)
 
-    radius_steps = [dsm_radius]
-    for _ in range(2):
-        radius_steps.append(radius_steps[-1] * math.sqrt(2)) # sqrt(2) is arbitrary
-
     log.ODM_INFO('Creating DSM for 2.5D mesh')
 
     commands.create_dem(
             inPointCloud,
             'mesh_dsm',
             output_type='max',
-            radiuses=list(map(str, radius_steps)),
+            radiuses=radius_steps,
             gapfill=True,
             outdir=tmp_directory,
             resolution=dsm_resolution,

opendm/point_cloud.py

@@ -71,7 +71,7 @@ def split(input_point_cloud, outdir, filename_template, capacity, dims=None):
     return [os.path.join(outdir, f) for f in os.listdir(outdir)]
 
 
-def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=16, sample_radius=0, boundary=None, max_concurrency=1):
+def filter(input_point_cloud, output_point_cloud, output_stats, standard_deviation=2.5, sample_radius=0, boundary=None, max_concurrency=1):
     """
     Filters a point cloud
     """
@@ -89,10 +89,11 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
         log.ODM_INFO("Sampling points around a %sm radius" % sample_radius)
         args.append('--radius %s' % sample_radius)
 
-    if standard_deviation > 0 and meank > 0:
+    meank = 16
-        log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation))
+    log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation))
-        args.append('--meank %s' % meank)
+    args.append('--meank %s' % meank)
-        args.append('--std %s' % standard_deviation)
+    args.append('--std %s' % standard_deviation)
+    args.append('--stats "%s"' % output_stats)
     
     if boundary is not None:
         log.ODM_INFO("Boundary {}".format(boundary))
@@ -107,6 +108,26 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
     if not os.path.exists(output_point_cloud):
         log.ODM_WARNING("{} not found, filtering has failed.".format(output_point_cloud))
 
+
+def get_spacing(stats_file, resolution_fallback=5.0):
+    def fallback():
+        log.ODM_WARNING("Cannot read %s, falling back to resolution estimate" % stats_file)
+        return (resolution_fallback / 100.0) / 2.0
+
+    if not os.path.isfile(stats_file):
+        return fallback()
+    
+    with open(stats_file, 'r') as f:
+        j = json.loads(f.read())
+        if "spacing" in j:
+            d = j["spacing"]
+            if d > 0:
+                return round(d, 3)
+            else:
+                return fallback()
+        else:
+            return fallback()
+
 def export_info_json(pointcloud_path, info_file_path):
     system.run('pdal info --dimensions "X,Y,Z" "{0}" > "{1}"'.format(pointcloud_path, info_file_path))
 

opendm/types.py

@@ -290,6 +290,7 @@ class ODM_Tree(object):
 
         # filter points
         self.filtered_point_cloud = os.path.join(self.odm_filterpoints, "point_cloud.ply")
+        self.filtered_point_cloud_stats = os.path.join(self.odm_filterpoints, "point_cloud_stats.json")
 
         # odm_meshing
         self.odm_mesh = os.path.join(self.odm_meshing, 'odm_mesh.ply')

stages/odm_dem.py

@@ -29,12 +29,8 @@ class ODMDEMStage(types.ODM_Stage):
             ignore_resolution = True
             pseudo_georeference = True
 
-        # It is probably not reasonable to have accurate DEMs a the same resolution as the source photos, so reduce it
-        # by a factor!
-        gsd_scaling = 2.0
-
         resolution = gsd.cap_resolution(args.dem_resolution, tree.opensfm_reconstruction, 
-                                        gsd_scaling=gsd_scaling,
+                                        gsd_scaling=1.0,
                                         ignore_gsd=args.ignore_gsd,
                                         ignore_resolution=ignore_resolution and args.ignore_gsd,
                                         has_gcp=reconstruction.has_gcp())
@@ -88,9 +84,7 @@ class ODMDEMStage(types.ODM_Stage):
                 if args.dsm or (args.dtm and args.dem_euclidean_map): products.append('dsm')
                 if args.dtm: products.append('dtm')
 
-                radius_steps = [(resolution / 100.0) / 2.0]
+                radius_steps = commands.get_dem_radius_steps(tree.filtered_point_cloud_stats, args.dem_gapfill_steps, resolution)
-                for _ in range(args.dem_gapfill_steps - 1):
-                    radius_steps.append(radius_steps[-1] * math.sqrt(2)) # sqrt(2) is arbitrary, maybe there's a better value?
 
                 for product in products:
                     commands.create_dem(

stages/odm_filterpoints.py

@@ -49,7 +49,7 @@ class ODMFilterPoints(types.ODM_Stage):
                 else:
                     log.ODM_WARNING("Not a georeferenced reconstruction, will ignore --auto-boundary")
                     
-            point_cloud.filter(inputPointCloud, tree.filtered_point_cloud, 
+            point_cloud.filter(inputPointCloud, tree.filtered_point_cloud, tree.filtered_point_cloud_stats,
                                 standard_deviation=args.pc_filter, 
                                 sample_radius=args.pc_sample,
                                 boundary=boundary_offset(outputs.get('boundary'), reconstruction.get_proj_offset()),

stages/odm_meshing.py

@@ -7,6 +7,7 @@ from opendm import context
 from opendm import mesh
 from opendm import gsd
 from opendm import types
+from opendm.dem import commands
 
 class ODMeshingStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -40,28 +41,18 @@ class ODMeshingStage(types.ODM_Stage):
           if not io.file_exists(tree.odm_25dmesh) or self.rerun():
 
               log.ODM_INFO('Writing ODM 2.5D Mesh file in: %s' % tree.odm_25dmesh)
-              ortho_resolution = gsd.cap_resolution(args.orthophoto_resolution, tree.opensfm_reconstruction, 
-                                                    ignore_gsd=args.ignore_gsd,
-                                                    ignore_resolution=(not reconstruction.is_georeferenced()) and args.ignore_gsd,
-                                                    has_gcp=reconstruction.has_gcp()) / 100.0 
-              
-              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 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
-              
-              dsm_radius = dsm_resolution * math.sqrt(2)
 
-              if args.fast_orthophoto:
+              multiplier = math.pi / 2.0
-                  dsm_radius *= 2
+              radius_steps = commands.get_dem_radius_steps(tree.filtered_point_cloud_stats, 3, args.orthophoto_resolution, multiplier=multiplier)
-                  dsm_resolution *= 8
+              dsm_resolution = radius_steps[0] / multiplier
 
               log.ODM_INFO('ODM 2.5D DSM resolution: %s' % dsm_resolution)
               
+              if args.fast_orthophoto:
+                  dsm_resolution *= 8.0
+
               mesh.create_25dmesh(tree.filtered_point_cloud, tree.odm_25dmesh,
-                    dsm_radius=dsm_radius,
+                    radius_steps,
                     dsm_resolution=dsm_resolution, 
                     depth=self.params.get('oct_tree'),
                     maxVertexCount=self.params.get('max_vertex'),