Crop to boundary in filterpoints

opendm/boundary.py

@@ -36,6 +36,20 @@ def load_boundary(boundary_json, reproject_to_proj4=None):
         
         return coords
 
+def boundary_offset(boundary, reconstruction_offset):
+    if boundary is None or reconstruction_offset is None:
+        return boundary
+    
+    res = []
+    dims = len(boundary[0])
+    for c in boundary:
+        if dims == 2:
+            res.append((c[0] - reconstruction_offset[0], c[1] - reconstruction_offset[1]))
+        else:
+            res.append((c[0] - reconstruction_offset[0], c[1] - reconstruction_offset[1], c[2]))
+    
+    return res
+
 def as_polygon(boundary):
     return "POLYGON((" + ", ".join([" ".join(map(str, c)) for c in boundary]) + "))"
 

opendm/mesh.py

@@ -126,7 +126,7 @@ def dem_to_mesh_gridded(inGeotiff, outMesh, maxVertexCount, verbose=False, maxCo
 
     system.run('"{reconstructmesh}" -i "{infile}" '
          '-o "{outfile}" '
-         '--remove-spikes 0 --remove-spurious 20 --smooth 0 '
+         '--remove-spikes 0 --remove-spurious 0 --smooth 0 '
          '--target-face-num {max_faces} '.format(**cleanupArgs))
 
     # Delete intermediate results

opendm/point_cloud.py

@@ -7,6 +7,7 @@ from opendm import entwine
 from opendm import io
 from opendm.concurrency import parallel_map
 from opendm.utils import double_quote
+from opendm.boundary import as_polygon
 
 def ply_info(input_ply):
     if not os.path.exists(input_ply):
@@ -38,7 +39,8 @@ def ply_info(input_ply):
     return {
         'has_normals': has_normals,
         'vertex_count': vertex_count,
-        'has_views': has_views
+        'has_views': has_views,
+        'header_lines': i + 1
     }
 
 
@@ -68,7 +70,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, verbose=False, max_concurrency=1):
+def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=16, sample_radius=0, boundary=None, verbose=False, max_concurrency=1):
     """
     Filters a point cloud
     """
@@ -86,6 +88,10 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
         log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation))
         filters.append('outlier')
         filters.append('range')
+    
+    if boundary is not None:
+        log.ODM_INFO("Boundary {}".format(boundary))
+        filters.append('crop')
 
     info = ply_info(input_point_cloud)
     dims = "x=float,y=float,z=float,"
@@ -116,7 +122,8 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
             filter(pcs['path'], io.related_file_path(pcs['path'], postfix="_filtered"), 
                         standard_deviation=standard_deviation, 
                         meank=meank, 
-                        sample_radius=sample_radius, 
+                        sample_radius=sample_radius,
+                        boundary=boundary,
                         verbose=verbose,
                         max_concurrency=1)
         # Filter
@@ -159,6 +166,9 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
             # Remove outliers
             cmd += "--filters.range.limits=\"Classification![7:7]\" "
 
+        if 'crop' in filters:
+            cmd += "--filters.crop.polygon=\"%s\"" % as_polygon(boundary)
+
         system.run(cmd)
 
     if not os.path.exists(output_point_cloud):

opendm/types.py

@@ -158,6 +158,12 @@ class ODM_Reconstruction(object):
     def get_proj_srs(self):
         if self.is_georeferenced():
             return self.georef.proj4()
+    
+    def get_proj_offset(self):
+        if self.is_georeferenced():
+            return (self.georef.utm_east_offset, self.georef.utm_north_offset)
+        else:
+            return (None, None)
 
     def get_photo(self, filename):
         for p in self.photos:

stages/dataset.py

@@ -160,4 +160,5 @@ class ODMLoadDatasetStage(types.ODM_Stage):
             if reconstruction.is_georeferenced():
                 outputs['boundary'] = boundary.load_boundary(args.boundary, reconstruction.get_proj_srs())
             else:
+                args.boundary = None
                 log.ODM_WARNING("Reconstruction is not georeferenced, but boundary file provided (will ignore boundary file)")

stages/odm_filterpoints.py

@@ -6,6 +6,7 @@ from opendm import system
 from opendm import context
 from opendm import point_cloud
 from opendm import types
+from opendm.boundary import boundary_offset
 
 class ODMFilterPoints(types.ODM_Stage):
     def process(self, args, outputs):
@@ -24,6 +25,7 @@ class ODMFilterPoints(types.ODM_Stage):
             point_cloud.filter(inputPointCloud, tree.filtered_point_cloud, 
                                 standard_deviation=args.pc_filter, 
                                 sample_radius=args.pc_sample,
+                                boundary=boundary_offset(outputs.get('boundary'), reconstruction.get_proj_offset()),
                                 verbose=args.verbose,
                                 max_concurrency=args.max_concurrency)
             

stages/odm_georeferencing.py

@@ -128,12 +128,6 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                         '--writers.las.vlrs="{\\\"filename\\\": \\\"%s\\\", \\\"user_id\\\": \\\"ODM_GCP\\\", \\\"description\\\": \\\"Ground Control Points (GML)\\\"}"' % gcp_gml_export_file.replace(os.sep, "/")
                     ]
                 
-                if 'boundary' in outputs:
-                    stages.append("crop")
-                    params += [
-                        '--filters.crop.polygon="%s"' % as_polygon(outputs['boundary'])
-                    ]
-                
                 system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
 
                 self.update_progress(50)