Merge pull request #1760 from pierotofy/fastcut

Skip feathered raster generation when possible

VERSION

@@ -1 +1 @@
-3.3.4
+3.5.1

configure.sh

@@ -127,6 +127,9 @@ installreqs() {
     installdepsfromsnapcraft build openmvs
     
     set -e
+
+    # edt requires numpy to build
+    pip install --ignore-installed numpy==1.23.1
     pip install --ignore-installed -r requirements.txt
     #if [ ! -z "$GPU_INSTALL" ]; then
     #fi

opendm/config.py

@@ -60,14 +60,14 @@ rerun_stages = {
     'orthophoto_no_tiled': 'odm_orthophoto',
     'orthophoto_png': 'odm_orthophoto',
     'orthophoto_resolution': 'odm_orthophoto',
-    'pc_classify': 'odm_dem',
+    'pc_classify': 'odm_georeferencing',
     'pc_copc': 'odm_georeferencing',
     'pc_csv': 'odm_georeferencing',
     'pc_ept': 'odm_georeferencing',
     'pc_filter': 'openmvs',
     'pc_las': 'odm_georeferencing',
     'pc_quality': 'opensfm',
-    'pc_rectify': 'odm_dem',
+    'pc_rectify': 'odm_georeferencing',
     'pc_sample': 'odm_filterpoints',
     'pc_skip_geometric': 'openmvs',
     'primary_band': 'dataset',
@@ -217,7 +217,7 @@ def config(argv=None, parser=None):
     parser.add_argument('--feature-type',
                         metavar='<string>',
                         action=StoreValue,
-                        default='sift',
+                        default='dspsift',
                         choices=['akaze', 'dspsift', 'hahog', 'orb', 'sift'],
                         help=('Choose the algorithm for extracting keypoints and computing descriptors. '
                             'Can be one of: %(choices)s. Default: '
@@ -485,7 +485,7 @@ def config(argv=None, parser=None):
                         metavar='<positive float>',
                         action=StoreValue,
                         type=float,
-                        default=2.5,
+                        default=5,
                         help='Filters the point cloud by removing points that deviate more than N standard deviations from the local mean. Set to 0 to disable filtering. '
                              'Default: %(default)s')
     
@@ -838,7 +838,7 @@ def config(argv=None, parser=None):
                         type=float,
                         action=StoreValue,
                         metavar='<positive float>',
-                        default=10,
+                        default=3,
                         help='Set a value in meters for the GPS Dilution of Precision (DOP) '
                         'information for all images. If your images are tagged '
                         'with high precision GPS information (RTK), this value will be automatically '

opendm/geo.py

@@ -12,6 +12,9 @@ class GeoFile:
 
         with open(self.geo_path, 'r') as f:
             contents = f.read().strip()
+        
+        # Strip eventual BOM characters
+        contents = contents.replace('\ufeff', '')
 
         lines = list(map(str.strip, contents.split('\n')))
         if lines:

opendm/gltf.py

@@ -279,9 +279,10 @@ def obj2glb(input_obj, output_glb, rtc=(None, None), draco_compression=True, _in
     )
 
     gltf.extensionsRequired = ['KHR_materials_unlit']
+    gltf.extensionsUsed = ['KHR_materials_unlit']
 
     if rtc != (None, None) and len(rtc) >= 2:
-        gltf.extensionsUsed = ['CESIUM_RTC', 'KHR_materials_unlit']
+        gltf.extensionsUsed.append('CESIUM_RTC')
         gltf.extensions = {
             'CESIUM_RTC': {
                 'center': [float(rtc[0]), float(rtc[1]), 0.0]

opendm/mesh.py

@@ -187,7 +187,7 @@ def screened_poisson_reconstruction(inPointCloud, outMesh, depth = 8, samples =
             if threads < 1:
                 break
             else:
-                log.ODM_WARNING("PoissonRecon failed with %s threads, let's retry with %s..." % (threads, threads // 2))
+                log.ODM_WARNING("PoissonRecon failed with %s threads, let's retry with %s..." % (threads * 2, threads))
 
 
     # Cleanup and reduce vertex count if necessary

opendm/osfm.py

@@ -64,7 +64,6 @@ class OSFMContext:
                             "Check that the images have enough overlap, "
                             "that there are enough recognizable features "
                             "and that the images are in focus. "
-                            "You could also try to increase the --min-num-features parameter."
                             "The program will now exit.")
 
         if rolling_shutter_correct:
@@ -794,3 +793,12 @@ def get_all_submodel_paths(submodels_path, *all_paths):
                 result.append([os.path.join(submodels_path, f, ap) for ap in all_paths])
 
     return result
+
+def is_submodel(opensfm_root):
+    # A bit hackish, but works without introducing additional markers / flags
+    # Look at the path of the opensfm directory and see if "submodel_" is part of it
+    parts = os.path.abspath(opensfm_root).split(os.path.sep)
+
+    return (len(parts) >= 2 and parts[-2][:9] == "submodel_") or \
+           os.path.isfile(os.path.join(opensfm_root, "split_merge_stop_at_reconstruction.txt")) or \
+           os.path.isfile(os.path.join(opensfm_root, "features", "empty"))

opendm/photo.py

@@ -430,7 +430,8 @@ class ODM_Photo:
                         camera_projection = camera_projection.lower()
 
                         # Parrot Sequoia's "fisheye" model maps to "fisheye_opencv"
-                        if camera_projection == "fisheye" and self.camera_make.lower() == "parrot" and self.camera_model.lower() == "sequoia":
+                        # or better yet, replace all fisheye with fisheye_opencv, but wait to change API signature
+                        if camera_projection == "fisheye":
                             camera_projection = "fisheye_opencv"
 
                         if camera_projection in projections:

opendm/point_cloud.py

@@ -9,6 +9,8 @@ from opendm.concurrency import parallel_map
 from opendm.utils import double_quote
 from opendm.boundary import as_polygon, as_geojson
 from opendm.dem.pdal import run_pipeline
+from opendm.opc import classify
+from opendm.dem import commands
 
 def ply_info(input_ply):
     if not os.path.exists(input_ply):
@@ -274,6 +276,32 @@ def merge_ply(input_point_cloud_files, output_file, dims=None):
     system.run(' '.join(cmd))
 
 def post_point_cloud_steps(args, tree, rerun=False):
+    # Classify and rectify before generating derivate files
+    if args.pc_classify:
+        pc_classify_marker = os.path.join(tree.odm_georeferencing, 'pc_classify_done.txt')
+
+        if not io.file_exists(pc_classify_marker) or rerun:
+            log.ODM_INFO("Classifying {} using Simple Morphological Filter (1/2)".format(tree.odm_georeferencing_model_laz))
+            commands.classify(tree.odm_georeferencing_model_laz,
+                                args.smrf_scalar, 
+                                args.smrf_slope, 
+                                args.smrf_threshold, 
+                                args.smrf_window
+                            )
+
+            log.ODM_INFO("Classifying {} using OpenPointClass (2/2)".format(tree.odm_georeferencing_model_laz))
+            classify(tree.odm_georeferencing_model_laz, args.max_concurrency)
+
+            with open(pc_classify_marker, 'w') as f:
+                f.write('Classify: smrf\n')
+                f.write('Scalar: {}\n'.format(args.smrf_scalar))
+                f.write('Slope: {}\n'.format(args.smrf_slope))
+                f.write('Threshold: {}\n'.format(args.smrf_threshold))
+                f.write('Window: {}\n'.format(args.smrf_window))
+    
+    if args.pc_rectify:
+        commands.rectify(tree.odm_georeferencing_model_laz)
+
     # XYZ point cloud output
     if args.pc_csv:
         log.ODM_INFO("Creating CSV file (XYZ format)")

opendm/rollingshutter.py

@@ -19,6 +19,7 @@ RS_DATABASE = {
     
     'dji fc220': 64, # DJI Mavic Pro (Platinum)
     'hasselblad l1d-20c': lambda p: 47 if p.get_capture_megapixels() < 17 else 56, # DJI Mavic 2 Pro (at 16:10 => 16.8MP 47ms, at 3:2 => 19.9MP 56ms. 4:3 has 17.7MP with same image height as 3:2 which can be concluded as same sensor readout)
+    'hasselblad l2d-20c': 16.6, # DJI Mavic 3 (not enterprise version)
 
     'dji fc3582': lambda p: 26 if p.get_capture_megapixels() < 48 else 60, # DJI Mini 3 pro (at 48MP readout is 60ms, at 12MP it's 26ms) 
 

opendm/types.py

@@ -13,6 +13,7 @@ from opendm import log
 from opendm import io
 from opendm import system
 from opendm import context
+from opendm import multispectral
 
 from opendm.progress import progressbc
 from opendm.photo import ODM_Photo
@@ -45,19 +46,51 @@ class ODM_Reconstruction(object):
             band_photos[p.band_name].append(p)
             
         bands_count = len(band_photos)
-        if bands_count >= 2 and bands_count <= 8:
+        
+        # Band name with the minimum number of photos
+        max_band_name = None
+        max_photos = -1
+        for band_name in band_photos:
+            if len(band_photos[band_name]) > max_photos:
+                max_band_name = band_name
+                max_photos = len(band_photos[band_name])
+
+        if bands_count >= 2 and bands_count <= 10:
             # Validate that all bands have the same number of images,
             # otherwise this is not a multi-camera setup
-            img_per_band = len(band_photos[p.band_name])
-            for band in band_photos:
-                if len(band_photos[band]) != img_per_band:
-                    log.ODM_ERROR("Multi-camera setup detected, but band \"%s\" (identified from \"%s\") has only %s images (instead of %s), perhaps images are missing or are corrupted. Please include all necessary files to process all bands and try again." % (band, band_photos[band][0].filename, len(band_photos[band]), img_per_band))
-                    raise RuntimeError("Invalid multi-camera images")
+            img_per_band = len(band_photos[max_band_name])
             
             mc = []
             for band_name in band_indexes:
                 mc.append({'name': band_name, 'photos': band_photos[band_name]})
             
+            filter_missing = False
+            for band in band_photos:
+                if len(band_photos[band]) < img_per_band:
+                    log.ODM_WARNING("Multi-camera setup detected, but band \"%s\" (identified from \"%s\") has only %s images (instead of %s), perhaps images are missing or are corrupted." % (band, band_photos[band][0].filename, len(band_photos[band]), len(band_photos[max_band_name])))
+                    filter_missing = True
+            
+            if filter_missing:
+                # Calculate files to ignore
+                _, p2s = multispectral.compute_band_maps(mc, max_band_name)
+
+                max_files_per_band = 0
+
+                for filename in p2s:
+                    max_files_per_band = max(max_files_per_band, len(p2s[filename]))
+
+                for filename in p2s:
+                    if len(p2s[filename]) < max_files_per_band:
+                        photos_to_remove = p2s[filename] + [p for p in self.photos if p.filename == filename]
+                        for photo in photos_to_remove:
+                            log.ODM_WARNING("Excluding %s" % photo.filename)
+
+                            self.photos = [p for p in self.photos if p != photo]
+                            for i in range(len(mc)):
+                                mc[i]['photos'] = [p for p in mc[i]['photos'] if p != photo]
+                
+                log.ODM_INFO("New image count: %s" % len(self.photos))
+
             # We enforce a normalized band order for all bands that we can identify
             # and rely on the manufacturer's band_indexes as a fallback for all others
             normalized_band_order = {
@@ -94,7 +127,7 @@ class ODM_Reconstruction(object):
 
             for c, d in enumerate(mc):
                 log.ODM_INFO(f"Band {c + 1}: {d['name']}")
-
+            
             return mc
 
         return None

opendm/video/srtparser.py

@@ -54,8 +54,10 @@ class SrtFileParser:
         if not self.gps_data:
             for d in self.data:
                 lat, lon, alt = d.get('latitude'), d.get('longitude'), d.get('altitude')
+                if alt is None:
+                    alt = 0
                 tm = d.get('start')
-
+                
                 if lat is not None and lon is not None:
                     if self.ll_to_utm is None:
                         self.ll_to_utm, self.utm_to_ll = location.utm_transformers_from_ll(lon, lat)
@@ -127,6 +129,20 @@ class SrtFileParser:
         # 00:00:35,000 --> 00:00:36,000
         # F/6.3, SS 60, ISO 100, EV 0, RTK (120.083799, 30.213635, 28), HOME (120.084146, 30.214243, 103.55m), D 75.36m, H 76.19m, H.S 0.30m/s, V.S 0.00m/s, F.PRY (-5.3°, 2.1°, 28.3°), G.PRY (-40.0°, 0.0°, 28.2°)
 
+        # DJI Unknown Model #1
+        # 1
+        # 00:00:00,000 --> 00:00:00,033
+        # <font size="28">SrtCnt : 1, DiffTime : 33ms
+        # 2024-01-18 10:23:26.397
+        # [iso : 150] [shutter : 1/5000.0] [fnum : 170] [ev : 0] [ct : 5023] [color_md : default] [focal_len : 240] [dzoom_ratio: 10000, delta:0],[latitude: -22.724555] [longitude: -47.602414] [rel_alt: 0.300 abs_alt: 549.679] </font>
+
+        # DJI Mavic 2 Zoom
+        # 1
+        # 00:00:00,000 --> 00:00:00,041
+        # <font size="36">FrameCnt : 1, DiffTime : 41ms
+        # 2023-07-15 11:55:16,320,933
+        # [iso : 100] [shutter : 1/400.0] [fnum : 280] [ev : 0] [ct : 5818] [color_md : default] [focal_len : 240] [latitude : 0.000000] [longtitude : 0.000000] [altitude: 0.000000] </font>
+
         with open(self.filename, 'r') as f:
 
             iso = None
@@ -197,12 +213,14 @@ class SrtFileParser:
 
                 latitude = match_single([
                     ("latitude: ([\d\.\-]+)", lambda v: float(v) if v != 0 else None),
+                    ("latitude : ([\d\.\-]+)", lambda v: float(v) if v != 0 else None),
                     ("GPS \([\d\.\-]+,? ([\d\.\-]+),? [\d\.\-]+\)", lambda v: float(v) if v != 0 else None),
                     ("RTK \([-+]?\d+\.\d+, (-?\d+\.\d+), -?\d+\)", lambda v: float(v) if v != 0 else None),
                 ], line)
 
                 longitude = match_single([
                     ("longitude: ([\d\.\-]+)", lambda v: float(v) if v != 0 else None),
+                    ("longtitude : ([\d\.\-]+)", lambda v: float(v) if v != 0 else None),
                     ("GPS \(([\d\.\-]+),? [\d\.\-]+,? [\d\.\-]+\)", lambda v: float(v) if v != 0 else None),
                     ("RTK \((-?\d+\.\d+), [-+]?\d+\.\d+, -?\d+\)", lambda v: float(v) if v != 0 else None),
                 ], line)
@@ -211,4 +229,5 @@ class SrtFileParser:
                     ("altitude: ([\d\.\-]+)", lambda v: float(v) if v != 0 else None),
                     ("GPS \([\d\.\-]+,? [\d\.\-]+,? ([\d\.\-]+)\)", lambda v: float(v) if v != 0 else None),
                     ("RTK \([-+]?\d+\.\d+, [-+]?\d+\.\d+, (-?\d+)\)", lambda v: float(v) if v != 0 else None),
+                    ("abs_alt: ([\d\.\-]+)", lambda v: float(v) if v != 0 else None),
                 ], line)

stages/mvstex.py

@@ -110,7 +110,7 @@ class ODMMvsTexStage(types.ODM_Stage):
 
                 mvs_tmp_dir = os.path.join(r['out_dir'], 'tmp')
 
-                # Make sure tmp directory is empty
+                # mvstex creates a tmp directory, so make sure it is empty
                 if io.dir_exists(mvs_tmp_dir):
                     log.ODM_INFO("Removing old tmp directory {}".format(mvs_tmp_dir))
                     shutil.rmtree(mvs_tmp_dir)

stages/odm_dem.py

@@ -12,7 +12,6 @@ from opendm.cropper import Cropper
 from opendm import pseudogeo
 from opendm.tiles.tiler import generate_dem_tiles
 from opendm.cogeo import convert_to_cogeo
-from opendm.opc import classify
 
 class ODMDEMStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -35,7 +34,6 @@ class ODMDEMStage(types.ODM_Stage):
                                         ignore_resolution=ignore_resolution and args.ignore_gsd,
                                         has_gcp=reconstruction.has_gcp())
 
-        log.ODM_INFO('Classify: ' + str(args.pc_classify))
         log.ODM_INFO('Create DSM: ' + str(args.dsm))
         log.ODM_INFO('Create DTM: ' + str(args.dtm))
         log.ODM_INFO('DEM input file {0} found: {1}'.format(dem_input, str(pc_model_found)))
@@ -45,34 +43,9 @@ class ODMDEMStage(types.ODM_Stage):
         if not io.dir_exists(odm_dem_root):
             system.mkdir_p(odm_dem_root)
 
-        if args.pc_classify and pc_model_found:
-            pc_classify_marker = os.path.join(odm_dem_root, 'pc_classify_done.txt')
-
-            if not io.file_exists(pc_classify_marker) or self.rerun():
-                log.ODM_INFO("Classifying {} using Simple Morphological Filter (1/2)".format(dem_input))
-                commands.classify(dem_input,
-                                  args.smrf_scalar, 
-                                  args.smrf_slope, 
-                                  args.smrf_threshold, 
-                                  args.smrf_window
-                                )
-
-                log.ODM_INFO("Classifying {} using OpenPointClass (2/2)".format(dem_input))
-                classify(dem_input, args.max_concurrency)
-
-                with open(pc_classify_marker, 'w') as f:
-                    f.write('Classify: smrf\n')
-                    f.write('Scalar: {}\n'.format(args.smrf_scalar))
-                    f.write('Slope: {}\n'.format(args.smrf_slope))
-                    f.write('Threshold: {}\n'.format(args.smrf_threshold))
-                    f.write('Window: {}\n'.format(args.smrf_window))
-            
         progress = 20
         self.update_progress(progress)
 
-        if args.pc_rectify:
-            commands.rectify(dem_input)
-
         # Do we need to process anything here?
         if (args.dsm or args.dtm) and pc_model_found:
             dsm_output_filename = os.path.join(odm_dem_root, 'dsm.tif')
@@ -120,7 +93,7 @@ class ODMDEMStage(types.ODM_Stage):
                     if args.cog:
                         convert_to_cogeo(dem_geotiff_path, max_workers=args.max_concurrency)
 
-                    progress += 30
+                    progress += 40
                     self.update_progress(progress)
             else:
                 log.ODM_WARNING('Found existing outputs in: %s' % odm_dem_root)

stages/odm_orthophoto.py

@@ -7,6 +7,7 @@ from opendm import context
 from opendm import types
 from opendm import gsd
 from opendm import orthophoto
+from opendm.osfm import is_submodel
 from opendm.concurrency import get_max_memory_mb
 from opendm.cutline import compute_cutline
 from opendm.utils import double_quote
@@ -114,6 +115,7 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
                     
                 # Cutline computation, before cropping
                 # We want to use the full orthophoto, not the cropped one.
+                submodel_run = is_submodel(tree.opensfm)
                 if args.orthophoto_cutline:
                     cutline_file = os.path.join(tree.odm_orthophoto, "cutline.gpkg")
 
@@ -122,15 +124,18 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
                                     cutline_file,
                                     args.max_concurrency,
                                     scale=0.25)
-
-                    orthophoto.compute_mask_raster(tree.odm_orthophoto_tif, cutline_file, 
-                                           os.path.join(tree.odm_orthophoto, "odm_orthophoto_cut.tif"),
-                                           blend_distance=20, only_max_coords_feature=True)
+                    
+                    if submodel_run:
+                        orthophoto.compute_mask_raster(tree.odm_orthophoto_tif, cutline_file, 
+                                            os.path.join(tree.odm_orthophoto, "odm_orthophoto_cut.tif"),
+                                            blend_distance=20, only_max_coords_feature=True)
+                    else:
+                        log.ODM_INFO("Not a submodel run, skipping mask raster generation")
 
                 orthophoto.post_orthophoto_steps(args, bounds_file_path, tree.odm_orthophoto_tif, tree.orthophoto_tiles, resolution)
 
                 # Generate feathered orthophoto also
-                if args.orthophoto_cutline:
+                if args.orthophoto_cutline and submodel_run:
                     orthophoto.feather_raster(tree.odm_orthophoto_tif, 
                             os.path.join(tree.odm_orthophoto, "odm_orthophoto_feathered.tif"),
                             blend_distance=20

stages/openmvs.py

@@ -94,7 +94,7 @@ class ODMOpenMVSStage(types.ODM_Stage):
                 extra_config.append("--ignore-mask-label 0")
 
             with open(densify_ini_file, 'w+') as f:
-                f.write("Optimize = 7\n")
+                f.write("Optimize = 7\nMin Views Filter = 1\n")
 
             def run_densify():
                 system.run('"%s" "%s" %s' % (context.omvs_densify_path, 
@@ -110,7 +110,7 @@ class ODMOpenMVSStage(types.ODM_Stage):
                     log.ODM_WARNING("OpenMVS failed with GPU, is your graphics card driver up to date? Falling back to CPU.")
                     gpu_config = ["--cuda-device -2"]
                     run_densify()
-                elif (e.errorCode == 137 or e.errorCode == 3221226505) and not pc_tile:
+                elif (e.errorCode == 137 or e.errorCode == 143 or e.errorCode == 3221226505) and not pc_tile:
                     log.ODM_WARNING("OpenMVS ran out of memory, we're going to turn on tiling to see if we can process this.")
                     pc_tile = True
                     config.append("--fusion-mode 1")
@@ -127,7 +127,7 @@ class ODMOpenMVSStage(types.ODM_Stage):
 
                 subscene_densify_ini_file = os.path.join(tree.openmvs, 'subscene-config.ini')
                 with open(subscene_densify_ini_file, 'w+') as f:
-                    f.write("Optimize = 0\nEstimation Geometric Iters = 0\n")
+                    f.write("Optimize = 0\nEstimation Geometric Iters = 0\nMin Views Filter = 1\n")
 
                 config = [
                     "--sub-scene-area 660000", # 8000
@@ -223,7 +223,7 @@ class ODMOpenMVSStage(types.ODM_Stage):
                         try:
                             system.run('"%s" %s' % (context.omvs_densify_path, ' '.join(config + gpu_config + extra_config)))
                         except system.SubprocessException as e:
-                            if e.errorCode == 137 or e.errorCode == 3221226505:
+                            if e.errorCode == 137 or e.errorCode == 143 or e.errorCode == 3221226505:
                                 log.ODM_WARNING("OpenMVS filtering ran out of memory, visibility checks will be skipped.")
                                 skip_filtering()
                             else: