diff --git a/SuperBuild/cmake/External-FPCFilter.cmake b/SuperBuild/cmake/External-FPCFilter.cmake index 78f5c398..7de6f668 100644 --- a/SuperBuild/cmake/External-FPCFilter.cmake +++ b/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------------- diff --git a/opendm/dem/commands.py b/opendm/dem/commands.py index 5cdfda10..4a948965 100755 --- a/opendm/dem/commands.py +++ b/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 \ No newline at end of file diff --git a/opendm/mesh.py b/opendm/mesh.py index 6708ab55..da393b79 100644 --- a/opendm/mesh.py +++ b/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, diff --git a/opendm/point_cloud.py b/opendm/point_cloud.py index 8d763594..f178c08c 100644 --- a/opendm/point_cloud.py +++ b/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: - log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation)) - args.append('--meank %s' % meank) - args.append('--std %s' % standard_deviation) + meank = 16 + log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation)) + args.append('--meank %s' % meank) + 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)) diff --git a/opendm/types.py b/opendm/types.py index a6110126..950e1d99 100644 --- a/opendm/types.py +++ b/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') diff --git a/stages/odm_dem.py b/stages/odm_dem.py index 67603591..cc8128cb 100755 --- a/stages/odm_dem.py +++ b/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] - 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? + radius_steps = commands.get_dem_radius_steps(tree.filtered_point_cloud_stats, args.dem_gapfill_steps, resolution) for product in products: commands.create_dem( diff --git a/stages/odm_filterpoints.py b/stages/odm_filterpoints.py index b988aa90..73db2655 100644 --- a/stages/odm_filterpoints.py +++ b/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()), diff --git a/stages/odm_meshing.py b/stages/odm_meshing.py index 72231edf..bfef268a 100644 --- a/stages/odm_meshing.py +++ b/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: - dsm_radius *= 2 - dsm_resolution *= 8 + multiplier = math.pi / 2.0 + radius_steps = commands.get_dem_radius_steps(tree.filtered_point_cloud_stats, 3, args.orthophoto_resolution, multiplier=multiplier) + 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'),