Disable seam leveling for multispectral data, and add backstops to flag and fix any negative reflectances that slip by

opendm/orthophoto.py

@@ -19,6 +19,38 @@
 from osgeo import ogr
 
 
+def clamp_negative_pixels(raster_path):
+    """
+    Checks a raster (georeferenced or not, any number of bands) for negative
+    pixel values. If any are found, they are zeroed in place on disk and a
+    warning is logged.
+
+    Uses rasterio (rather than osgeo.gdal's Band.ReadAsArray/WriteArray) since
+    the latter lazily imports osgeo.gdal_array, a separate compiled extension
+    that can be ABI-incompatible with the installed numpy in some environments.
+    :return the total number of negative pixels found (and zeroed) across all bands
+    """
+    if not io.file_exists(raster_path):
+        return 0
+
+    total_negative = 0
+    try:
+        with rasterio.open(raster_path, 'r+') as ds:
+            arr = ds.read()
+            negative = arr < 0
+            total_negative = int(np.count_nonzero(negative))
+            if total_negative > 0:
+                arr[negative] = 0
+                ds.write(arr)
+    except Exception as e:
+        log.WARNING("Cannot open %s to check for negative pixel values: %s" % (raster_path, str(e)))
+        return 0
+
+    if total_negative > 0:
+        log.WARNING("%s: found and zeroed %s negative pixel values" % (raster_path, total_negative))
+
+    return total_negative
+
 def get_orthophoto_vars(args):
     return {
         'TILED': 'NO' if args.orthophoto_no_tiled else 'YES',

stages/mvstex.py

@@ -1,10 +1,11 @@
-import os, shutil
+import os, shutil, glob
 
 from opendm import log
 from opendm import io
 from opendm import system
 from opendm import context
 from opendm import types
+from opendm import orthophoto
 from opendm.multispectral import get_primary_band_name
 from opendm.photo import find_largest_photo_dim
 from opendm.objpacker import obj_pack
@@ -80,11 +81,18 @@ def add_run(nvm_file, primary=True, band=None):
 
                 # Format arguments to fit Mvs-Texturing app
                 skipGlobalSeamLeveling = ""
+                skipLocalSeamLeveling = ""
                 keepUnseenFaces = ""
                 nadir = ""
 
-                if args.texturing_skip_global_seam_leveling:
+                # for multispectral data, force disabling of global seam leveling to avoid negative texture and reflectance values,
+                # or if user requested disabling for RGB/thermal data
+                if reconstruction.multi_camera or args.texturing_skip_global_seam_leveling:
                     skipGlobalSeamLeveling = "--skip_global_seam_leveling"
+                # for multispectral data, force disabling of local seam leveling to avoid negative texture and reflectance values
+                if reconstruction.multi_camera:
+                    skipLocalSeamLeveling = "--skip_local_seam_leveling"
+                # TODO: allow user-requested --skip_local_seam_leveling for RGB/thermal data?
                 if args.texturing_keep_unseen_faces:
                     keepUnseenFaces = "--keep_unseen_faces"
                 if (r['nadir']):
@@ -98,6 +106,7 @@ def add_run(nvm_file, primary=True, band=None):
                     'dataTerm': 'gmi',
                     'outlierRemovalType': 'gauss_clamping',
                     'skipGlobalSeamLeveling': skipGlobalSeamLeveling,
+                    'skipLocalSeamLeveling': skipLocalSeamLeveling,
                     'keepUnseenFaces': keepUnseenFaces,
                     'toneMapping': 'none',
                     'nadirMode': nadir,
@@ -121,12 +130,20 @@ def add_run(nvm_file, primary=True, band=None):
                         '-t {toneMapping} '
                         '{intermediate} '
                         '{skipGlobalSeamLeveling} '
+                        '{skipLocalSeamLeveling} '
                         '{keepUnseenFaces} '
                         '{nadirMode} '
                         '{labelingFile} '
                         '{numThreads} '
                         '{maxTextureSize} '.format(**kwargs))
 
+                # For multispectral data, check the resulting texture atlases for negative pixel
+                # values (reflectances should be >= 0). Warn and zero if any found. Typically, disabling global and local
+                # seam leveling avoids negative pixel values, but force an explicit clamp in case something slips by
+                if reconstruction.multi_camera:
+                    for texfile in glob.glob(kwargs['out_dir'] + "*.tif"):
+                        orthophoto.clamp_negative_pixels(texfile)
+
                 if r['primary'] and (not r['nadir'] or args.skip_3dmodel):
                     # Single material?
                     if args.texturing_single_material:

stages/odm_orthophoto.py

@@ -109,6 +109,13 @@ def process(self, args, outputs):
                        '-outputCornerFile "{corners}" {bands} {depth_idx} {inpaint} '
                        '{utm_offsets} {a_srs} {vars} {gdal_configs} '.format(**kwargs), env_vars={'OMP_NUM_THREADS': args.max_concurrency})
 
+            # Multispectral reflectance values should never be negative, so catch, warn and zero
+            # any that slipped through before any further cutline/feathering/tiling steps
+            # operate on this raster. Typically, these will have been avoided in previous stages. However, if they
+            # show up in this stage, then we should figure out why and fix the problem (not just the symptom)
+            if reconstruction.multi_camera:
+                orthophoto.clamp_negative_pixels(kwargs['ortho'])
+
             # Create georeferenced GeoTiff
             if reconstruction.is_georeferenced():
                 bounds_file_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.gpkg')