better support for radiometric panels

opendm/config.py

@@ -299,15 +299,26 @@ def config(argv=None, parser=None):
             metavar='<string>',
             action=StoreValue,
             default='none',
-            choices=['none', 'camera', 'camera+sun'],
+            choices=['none', 'camera', 'camera+sun', 'camera+panel'],
             help=('Set the radiometric calibration to perform on images. '
                 'When processing multispectral and thermal images you should set this option '
                 'to obtain reflectance/temperature values (otherwise you will get digital number values). '
                 '[camera] applies black level, vignetting, row gradient gain/exposure compensation (if appropriate EXIF tags are found) and computes absolute temperature values. '
                 '[camera+sun] is experimental, applies all the corrections of [camera], plus compensates for spectral radiance registered via a downwelling light sensor (DLS) taking in consideration the angle of the sun. '
+                '[camera+panel] applies all the corrections of [camera], plus derives per-band irradiance from in-field calibration reflectance panel captures (MicaSense and compatible). Panel captures are auto-detected from image metadata. '
                 'Can be one of: %(choices)s. Default: '
                 '%(default)s'))
 
+    parser.add_argument('--panel-reflectance',
+            metavar='<string>',
+            action=StoreValue,
+            default=None,
+            help=('Override the calibration panel reflectance (albedo, 0-1) used by '
+                '[--radiometric-calibration camera+panel]. Useful when the panel metadata '
+                'is missing or unreadable. Provide either a single value applied to all bands (e.g. 0.49), '
+                'or per-band values as a JSON object (e.g. \'{"Red": 0.49, "Green": 0.49}\'). Default: '
+                '%(default)s'))
+
     parser.add_argument('--max-concurrency',
                         metavar='<positive integer>',
                         action=StoreValue,

opendm/multispectral.py

@@ -114,26 +114,30 @@ def vignette_map(photo):
 
     return None, None, None
 
-def dn_to_reflectance(photo, image, use_sun_sensor=True):
+def dn_to_reflectance(photo, image, use_sun_sensor=True, panel_irradiance=None):
     radiance = dn_to_radiance(photo, image)
-    irradiance = compute_irradiance(photo, use_sun_sensor=use_sun_sensor)
+    irradiance = compute_irradiance(photo, use_sun_sensor=use_sun_sensor, panel_irradiance=panel_irradiance)
     reflectance = radiance * math.pi / irradiance
     reflectance[reflectance < 0.0] = 0.0
     reflectance[reflectance > 1.0] = 1.0
     return reflectance
 
-def compute_irradiance(photo, use_sun_sensor=True):
+def compute_irradiance(photo, use_sun_sensor=True, panel_irradiance=None):
     # Thermal (this should never happen, but just in case..)
     if photo.is_thermal():
         return 1.0
 
+    # Calibration panel-derived irradiance takes precedence when available.
+    # (Must be checked before the stored horizontal irradiance below, otherwise
+    # it would never be reached for DLS cameras that always write that tag.)
+    if panel_irradiance is not None:
+        return panel_irradiance
+
     # Some cameras (Micasense, DJI) store the value (nice! just return)
     hirradiance = photo.get_horizontal_irradiance()
     if hirradiance is not None:
         return hirradiance
 
-    # TODO: support for calibration panels
-
     if use_sun_sensor and photo.get_sun_sensor():
         # Estimate it
         dls_orientation_vector = np.array([0,0,-1])
@@ -166,6 +170,146 @@ def compute_irradiance(photo, use_sun_sensor=True):
 
     return 1.0
 
+
+def _panel_irradiance_for_photo(photo, images_path, panel_reflectance, saturation_limit):
+    """
+    Compute the panel-derived irradiance for a single panel capture / band.
+    Returns a float irradiance or None if it cannot be determined.
+    """
+    from opendm import panel as panel_utils
+
+    image_file = os.path.join(images_path, photo.filename)
+    raw = imread(image_file, unchanged=True, anydepth=True)
+    if raw is None:
+        log.ODM_WARNING("Cannot read panel image %s" % photo.filename)
+        return None
+    if len(raw.shape) == 2:
+        raw = raw[:, :, np.newaxis]
+
+    # Determine panel reflectance (albedo). User override wins over camera metadata.
+    albedo = None
+    if isinstance(panel_reflectance, dict):
+        albedo = panel_reflectance.get(photo.band_name)
+        if albedo is None:
+            albedo = panel_reflectance.get(photo.band_name.lower())
+    elif panel_reflectance is not None:
+        albedo = float(panel_reflectance)
+    if albedo is None:
+        albedo = photo.get_panel_albedo()
+
+    # Determine the panel active-area polygon from metadata
+    region = photo.get_panel_region()
+
+    if region is None:
+        log.ODM_WARNING("No panel region found for %s (band %s); skipping this panel" % (photo.filename, photo.band_name))
+        return None
+    if albedo is None:
+        log.ODM_WARNING("No panel reflectance available for %s (band %s). "
+                        "Supply --panel-reflectance to override. Skipping this panel" % (photo.filename, photo.band_name))
+        return None
+
+    # Saturation guard, computed on the raw digital numbers
+    bit_depth_max = photo.get_bit_depth_max()
+    sat_threshold = (0.99 * bit_depth_max) if bit_depth_max else None
+    _, _, npix, sat_frac = panel_utils.region_stats(raw[:, :, 0], region, saturation_threshold=sat_threshold)
+    if npix == 0:
+        log.ODM_WARNING("Panel region for %s contains no pixels; skipping this panel" % photo.filename)
+        return None
+    if sat_frac > saturation_limit:
+        log.ODM_WARNING("Panel region for %s is %.0f%% saturated; skipping this panel" % (photo.filename, sat_frac * 100.0))
+        return None
+
+    radiance = dn_to_radiance(photo, raw)
+    mean_radiance, _, _, _ = panel_utils.region_stats(radiance[:, :, 0], region)
+    if mean_radiance is None or albedo == 0:
+        return None
+
+    # Inverse of the reflectance equation: irradiance = radiance * pi / reflectance
+    return mean_radiance * math.pi / albedo
+
+
+def parse_panel_reflectance(value):
+    """
+    Parse the --panel-reflectance CLI argument into a usable override.
+
+    Accepts None, a single float (applied to all bands), or a JSON object
+    mapping band name -> reflectance. Returns None, float, or dict.
+    """
+    if value is None:
+        return None
+
+    if isinstance(value, (int, float, dict)):
+        return value
+
+    value = str(value).strip()
+    if value == "":
+        return None
+
+    # Try JSON object (per-band reflectance) first
+    if value.startswith("{"):
+        try:
+            import json
+            parsed = json.loads(value)
+            if isinstance(parsed, dict):
+                return {k: float(v) for k, v in parsed.items()}
+        except Exception as e:
+            log.ODM_WARNING("Cannot parse --panel-reflectance JSON (%s): %s" % (value, str(e)))
+            return None
+
+    # Otherwise a single float applied to all bands
+    try:
+        return float(value)
+    except ValueError:
+        log.ODM_WARNING("Cannot parse --panel-reflectance value: %s" % value)
+        return None
+
+
+def compute_irradiance_from_panels(panel_photos, images_path, panel_reflectance=None, max_concurrency=1, saturation_limit=0.1):
+    """
+    Compute flight-level, per-band irradiance from calibration panel captures.
+
+    :param panel_photos: list of ODM_Photo that are calibration panel captures
+    :param images_path: directory containing the raw images
+    :param panel_reflectance: optional override of panel albedo; either a single
+        float applied to all bands, or a dict mapping band name -> albedo
+    :param max_concurrency: reserved for future parallelism
+    :param saturation_limit: reject a panel if more than this fraction of the
+        panel region pixels are saturated
+    :return: dict mapping band name -> irradiance (only for bands successfully derived)
+    """
+    by_band = {}
+    for p in panel_photos:
+        by_band.setdefault(p.band_name, []).append(p)
+
+    irradiance_by_band = {}
+
+    for band_name, photos in by_band.items():
+        # Thermal/LWIR bands have no reflectance panel
+        if photos[0].is_thermal():
+            continue
+
+        values = []
+        for photo in photos:
+            try:
+                irr = _panel_irradiance_for_photo(photo, images_path, panel_reflectance, saturation_limit)
+            except Exception as e:
+                log.ODM_WARNING("Cannot compute panel irradiance for %s: %s" % (photo.filename, str(e)))
+                irr = None
+            if irr is not None and irr > 0:
+                values.append(irr)
+
+        if values:
+            # Robust aggregate across all panel captures for this band
+            irradiance_by_band[band_name] = float(np.median(values))
+            log.ODM_INFO("Panel irradiance for band %s: %s (from %s panel capture(s))" % (
+                band_name, irradiance_by_band[band_name], len(values)))
+        else:
+            log.ODM_WARNING("Could not derive panel irradiance for band %s; "
+                            "it will fall back to stored/DLS irradiance." % band_name)
+
+    return irradiance_by_band
+
+
 def get_photos_by_band(multi_camera, user_band_name):
     band_name = get_primary_band_name(multi_camera, user_band_name)
 

opendm/panel.py

@@ -0,0 +1,33 @@
+import numpy as np
+from skimage import measure
+
+
+def region_stats(image, region, saturation_threshold=None):
+    """
+    Compute statistics for an image over a polygonal region.
+
+    :param image: 2D numpy array (single band)
+    :param region: list/array of (x, y) coordinate tuples describing the polygon
+    :param saturation_threshold: optional value above which a pixel is considered saturated
+    :return: (mean, std, num_pixels, saturated_fraction)
+    """
+    region = np.asarray(region, dtype=float)
+
+    # skimage uses (row, col) ordering, image coordinates are (x, y)
+    rev_pts = np.fliplr(region)
+    h, w = image.shape[:2]
+    mask = measure.grid_points_in_poly((h, w), rev_pts)
+
+    num_pixels = int(mask.sum())
+    if num_pixels == 0:
+        return None, None, 0, 0.0
+
+    pixels = image[mask]
+    mean_value = float(pixels.mean())
+    stdev = float(pixels.std())
+
+    saturated_fraction = 0.0
+    if saturation_threshold is not None:
+        saturated_fraction = float((pixels > saturation_threshold).sum()) / num_pixels
+
+    return mean_value, stdev, num_pixels, saturated_fraction

opendm/photo.py

@@ -126,6 +126,13 @@ def __init__(self, path_file):
         self.gain = None
         self.gain_adjustment = None
 
+        # Calibration reflectance panel fields
+        # (MicaSense and compatible cameras tag in-field panel captures)
+        self.calibration_picture = None # 2 == camera flagged this as a calibration picture
+        self.panel_albedo = None # active panel surface albedo (0-1) as computed by the camera
+        self.panel_region = None # "x1,y1,x2,y2,..." image coordinates of the panel active area
+        self.panel_serial = None # panel serial number (e.g. RP06-1234567-SC)
+
         # Capture info
         self.exposure_time = None
         self.iso_speed = None
@@ -370,6 +377,28 @@ def parse_exif_values(self, _path_file):
                         '@drone-dji:SensorGainAdjustment'
                     ], float)
 
+                    # Calibration reflectance panel tags (MicaSense and compatible).
+                    # Parsed natively from XMP (no exiftool dependency).
+                    self.set_attr_from_xmp_tag('calibration_picture', xtags, [
+                        '@Camera:CalibrationPicture',
+                        'Camera:CalibrationPicture',
+                    ], int)
+
+                    self.set_attr_from_xmp_tag('panel_albedo', xtags, [
+                        '@Camera:Albedo',
+                        'Camera:Albedo',
+                    ], float)
+
+                    self.set_attr_from_xmp_tag('panel_region', xtags, [
+                        '@Camera:ReflectArea',
+                        'Camera:ReflectArea',
+                    ])
+
+                    self.set_attr_from_xmp_tag('panel_serial', xtags, [
+                        '@Camera:PanelSerial',
+                        'Camera:PanelSerial',
+                    ])
+
                     # Camera make / model for some cameras is stored in the XMP
                     if self.camera_make == '':
                         self.set_attr_from_xmp_tag('camera_make', xtags, [
@@ -772,6 +801,46 @@ def get_capture_id(self):
 
         return self.get_utc_time()
 
+    def get_panel_albedo(self):
+        # Active panel surface albedo (0-1) as computed by the camera
+        if self.panel_albedo is not None:
+            try:
+                return float(self.panel_albedo)
+            except (ValueError, TypeError):
+                return None
+        return None
+
+    def get_panel_region(self):
+        # A list of (x, y) image coordinate tuples describing the panel active area
+        if self.panel_region is None:
+            return None
+        try:
+            coords = [int(round(float(item))) for item in re.split(r"[,\s]+", str(self.panel_region).strip()) if item != ""]
+        except (ValueError, TypeError):
+            return None
+        if len(coords) < 8 or len(coords) % 2 != 0:
+            return None
+        return list(zip(coords[0::2], coords[1::2]))
+
+    def get_panel_serial(self):
+        if self.panel_serial:
+            return str(self.panel_serial)
+        return None
+
+    def is_calibration_picture(self):
+        # True if the camera flagged this frame as a calibration (panel) picture.
+        # This is a looser check than is_panel_image(): the camera may flag the
+        # frame without having extracted the full panel region/albedo metadata.
+        return self.calibration_picture == 2
+
+    def is_panel_image(self):
+        # True if this is an auto-detected calibration panel image with all the
+        # metadata required to compute irradiance.
+        return self.is_calibration_picture() and \
+            self.get_panel_albedo() is not None and \
+            self.get_panel_region() is not None and \
+            self.get_panel_serial() is not None
+
     def get_gps_dop(self):
         val = -9999
         if self.gps_xy_stddev is not None:

stages/dataset.py

@@ -102,6 +102,7 @@ def find_mask(photo_path, masks):
 
         # check if we rerun cell or not
         images_database_file = os.path.join(tree.root_path, 'images.json')
+        panels_database_file = os.path.join(tree.root_path, 'panels.json')
         if not io.file_exists(images_database_file) or self.rerun():
             if not os.path.exists(images_dir):
                 raise system.ExitException("There are no images in %s! Make sure that your project path and dataset name is correct. The current is set to: %s" % (images_dir, args.project_path))
@@ -289,6 +290,45 @@ def parallel_bg_filter(item):
 
                 # End bg removal
 
+                # Detect calibration reflectance panels (MicaSense and compatible).
+                # Panel captures are radiometric ground-truth frames (shot on the
+                # ground), not aerial survey images, so we compute their per-band
+                # irradiance once and then exclude them from the reconstruction.
+                if args.radiometric_calibration == "camera+panel":
+                    from opendm import multispectral
+
+                    panel_photos = [p for p in photos if p.is_calibration_picture()]
+                    if len(panel_photos) > 0:
+                        log.ODM_INFO("Found %s calibration panel image(s)" % len(panel_photos))
+
+                        panel_reflectance = multispectral.parse_panel_reflectance(args.panel_reflectance)
+                        irradiance = multispectral.compute_irradiance_from_panels(
+                            panel_photos, images_dir,
+                            panel_reflectance=panel_reflectance,
+                            max_concurrency=args.max_concurrency)
+
+                        try:
+                            with open(panels_database_file, 'w') as f:
+                                f.write(json.dumps({
+                                    'irradiance': irradiance,
+                                    'files': [p.filename for p in panel_photos],
+                                }))
+                            log.ODM_INFO("Wrote panel irradiance database: %s" % panels_database_file)
+                        except Exception as e:
+                            log.ODM_WARNING("Cannot write panel database: %s" % str(e))
+
+                        # Exclude panel captures from the reconstruction
+                        panel_filenames = set(p.filename for p in panel_photos)
+                        photos = [p for p in photos if p.filename not in panel_filenames]
+                        log.ODM_INFO("Excluded %s panel image(s) from reconstruction; %s image(s) remain" % (
+                            len(panel_filenames), len(photos)))
+
+                        if len(photos) == 0:
+                            raise system.ExitException("All images were detected as calibration panels. Nothing left to reconstruct.")
+                    else:
+                        log.ODM_WARNING("--radiometric-calibration camera+panel was set, but no calibration panel "
+                                        "images were detected. Reflectance will fall back to stored/DLS irradiance.")
+
                 # Save image database for faster restart
                 save_images_database(photos, images_database_file)
             else: