diff --git a/opendm/report/dsm_gradient.png b/opendm/report/dsm_gradient.png new file mode 100644 index 00000000..c6fc4f2e Binary files /dev/null and b/opendm/report/dsm_gradient.png differ diff --git a/opendm/tiles/tiler.py b/opendm/tiles/tiler.py index d7801ab5..59dc185a 100644 --- a/opendm/tiles/tiler.py +++ b/opendm/tiles/tiler.py @@ -13,17 +13,32 @@ def generate_orthophoto_tiles(geotiff, output_dir, max_concurrency): except Exception as e: log.ODM_WARNING("Cannot generate orthophoto tiles: %s" % str(e)) -def generate_dem_tiles(geotiff, output_dir, max_concurrency): +def generate_colored_hillshade(geotiff): relief_file = os.path.join(os.path.dirname(__file__), "color_relief.txt") hsv_merge_script = os.path.join(os.path.dirname(__file__), "hsv_merge.py") colored_dem = io.related_file_path(geotiff, postfix="color") hillshade_dem = io.related_file_path(geotiff, postfix="hillshade") colored_hillshade_dem = io.related_file_path(geotiff, postfix="colored_hillshade") - try: + outputs = [colored_dem, hillshade_dem, colored_hillshade_dem] + + # Cleanup previous + for f in outputs: + if os.path.isfile(f): + os.remove(f) + system.run('gdaldem color-relief "%s" "%s" "%s" -alpha -co ALPHA=YES' % (geotiff, relief_file, colored_dem)) system.run('gdaldem hillshade "%s" "%s" -z 1.0 -s 1.0 -az 315.0 -alt 45.0' % (geotiff, hillshade_dem)) system.run('python3 "%s" "%s" "%s" "%s"' % (hsv_merge_script, colored_dem, hillshade_dem, colored_hillshade_dem)) + + return outputs + except Exception as e: + log.ODM_WARNING("Cannot generate colored hillshade: %s" % str(e)) + return (None, None, None) + +def generate_dem_tiles(geotiff, output_dir, max_concurrency): + try: + colored_dem, hillshade_dem, colored_hillshade_dem = generate_colored_hillshade(geotiff) generate_tiles(colored_hillshade_dem, output_dir, max_concurrency) # Cleanup diff --git a/opendm/utils.py b/opendm/utils.py index 4edd0860..27ca6db1 100644 --- a/opendm/utils.py +++ b/opendm/utils.py @@ -1,5 +1,6 @@ from opendm import log from opendm.photo import find_largest_photo_dim +from osgeo import gdal def get_depthmap_resolution(args, photos): if 'depthmap_resolution_is_set' in args: @@ -23,3 +24,19 @@ def get_depthmap_resolution(args, photos): else: log.ODM_WARNING("Cannot compute max image dimensions, going with default depthmap_resolution of 640") return 640 # Sensible default + +def get_raster_stats(geotiff): + stats = [] + gtif = gdal.Open(geotiff) + for b in range(gtif.RasterCount): + srcband = gtif.GetRasterBand(b + 1) + s = srcband.GetStatistics(True, True) + + stats.append({ + 'min': s[0], + 'max': s[1], + 'mean': s[2], + 'stddev': s[3] + }) + + return stats \ No newline at end of file diff --git a/stages/odm_report.py b/stages/odm_report.py index 861b0014..6b22fbfe 100644 --- a/stages/odm_report.py +++ b/stages/odm_report.py @@ -11,18 +11,21 @@ from opendm.shots import get_geojson_shots_from_opensfm from opendm.osfm import OSFMContext from opendm import gsd from opendm.point_cloud import export_summary_json - +from opendm.cropper import Cropper +from opendm.orthophoto import get_orthophoto_vars, get_max_memory +from opendm.tiles.tiler import generate_colored_hillshade +from opendm.utils import get_raster_stats def hms(seconds): h = seconds // 3600 m = seconds % 3600 // 60 s = seconds % 3600 % 60 if h > 0: - return '{}h:{}m:{}s'.format(h, m, s) + return '{}h:{}m:{}s'.format(h, m, round(s, 0)) elif m > 0: - return '{}m:{}s'.format(m, s) + return '{}m:{}s'.format(m, round(s, 0)) else: - return '{}s'.format(s) + return '{}s'.format(round(s, 0)) def generate_point_cloud_stats(input_point_cloud, pc_summary_file): @@ -96,9 +99,10 @@ class ODMReport(types.ODM_Stage): odm_stats['point_cloud_statistics']['dense'] = not args.fast_orthophoto # Add runtime stats + total_time = (system.now_raw() - outputs['start_time']).total_seconds() odm_stats['odm_processing_statistics'] = { - 'total_time': (system.now_raw() - outputs['start_time']).total_seconds(), - 'total_time_human': hms(121), + 'total_time': total_time, + 'total_time_human': hms(total_time), 'average_gsd': gsd.opensfm_reconstruction_average_gsd(octx.recon_file(), use_all_shots=reconstruction.has_gcp()), } @@ -115,7 +119,7 @@ class ODMReport(types.ODM_Stage): if odm_stats.get('point_cloud_statistics') and point_cloud_file and views_dimension: bounds = odm_stats['point_cloud_statistics'].get('summary', {}).get('bounds') if bounds: - diagram_target_size = 1600 # pixels + image_target_size = 1400 # pixels osfm_stats_dir = os.path.join(tree.opensfm, "stats") diagram_tiff = os.path.join(osfm_stats_dir, "overlap.tif") diagram_png = os.path.join(osfm_stats_dir, "overlap.png") @@ -123,7 +127,7 @@ class ODMReport(types.ODM_Stage): width = bounds.get('maxx') - bounds.get('minx') height = bounds.get('maxy') - bounds.get('miny') max_dim = max(width, height) - resolution = float(max_dim) / float(diagram_target_size) + resolution = float(max_dim) / float(image_target_size) radius = resolution * math.sqrt(2) # Larger radius for sparse point cloud diagram @@ -141,10 +145,35 @@ class ODMReport(types.ODM_Stage): resolution, views_dimension, radius)) report_assets = os.path.abspath(os.path.join(os.path.dirname(__file__), "../opendm/report")) overlap_color_map = os.path.join(report_assets, "overlap_color_map.txt") + + bounds_file_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.gpkg') + if args.crop > 0 and os.path.isfile(bounds_file_path): + Cropper.crop(bounds_file_path, diagram_tiff, get_orthophoto_vars(args), keep_original=False) + system.run("gdaldem color-relief \"{}\" \"{}\" \"{}\" -of PNG -alpha".format(diagram_tiff, overlap_color_map, diagram_png)) - # Copy legend - shutil.copy(os.path.join(report_assets, "overlap_diagram_legend.png"), os.path.join(osfm_stats_dir, "overlap_diagram_legend.png")) + # Copy assets + for asset in ["overlap_diagram_legend.png", "dsm_gradient.png"]: + shutil.copy(os.path.join(report_assets, asset), os.path.join(osfm_stats_dir, asset)) + + # Generate previews of ortho/dsm + if os.path.isfile(tree.odm_orthophoto_tif): + osfm_ortho = os.path.join(osfm_stats_dir, "ortho.png") + system.run("gdal_translate -outsize {} 0 -of png \"{}\" \"{}\" --config GDAL_CACHEMAX {}%".format(image_target_size, tree.odm_orthophoto_tif, osfm_ortho, get_max_memory())) + + dsm_file = tree.path("odm_dem", "dsm.tif") + if os.path.isfile(dsm_file): + log.ODM_INFO("Computing raster stats for %s" % dsm_file) + dsm_stats = get_raster_stats(dsm_file) + if len(dsm_stats) > 0: + odm_stats['dsm_statistics'] = dsm_stats[0] + + osfm_dsm = os.path.join(osfm_stats_dir, "dsm.png") + colored_dem, hillshade_dem, colored_hillshade_dem = generate_colored_hillshade(dsm_file) + system.run("gdal_translate -outsize {} 0 -of png \"{}\" \"{}\" --config GDAL_CACHEMAX {}%".format(image_target_size, colored_hillshade_dem, osfm_dsm, get_max_memory())) + for f in [colored_dem, hillshade_dem, colored_hillshade_dem]: + if os.path.isfile(f): + os.remove(f) else: log.ODM_WARNING("Cannot generate overlap diagram, cannot compute point cloud bounds") else: diff --git a/stages/openmvs.py b/stages/openmvs.py index 0edb5d35..43eedd68 100644 --- a/stages/openmvs.py +++ b/stages/openmvs.py @@ -48,7 +48,7 @@ class ODMOpenMVSStage(types.ODM_Stage): "--min-resolution %s" % depthmap_resolution, "--max-resolution %s" % int(outputs['undist_image_max_size']), "--max-threads %s" % args.max_concurrency, - "--number-views-fuse 2", + "--number-views-fuse 3", '-w "%s"' % depthmaps_dir, "-v 0", ]