diff --git a/opendm/config.py b/opendm/config.py index 388108d64..8ac4b7cfd 100755 --- a/opendm/config.py +++ b/opendm/config.py @@ -299,15 +299,26 @@ def config(argv=None, parser=None): metavar='', 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='', + 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='', action=StoreValue, diff --git a/opendm/multispectral.py b/opendm/multispectral.py index ab1d56c58..6f3a55176 100644 --- a/opendm/multispectral.py +++ b/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) diff --git a/opendm/panel.py b/opendm/panel.py new file mode 100644 index 000000000..9e937d3da --- /dev/null +++ b/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 diff --git a/opendm/photo.py b/opendm/photo.py index eaa3433cb..e8a24e71d 100644 --- a/opendm/photo.py +++ b/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: diff --git a/stages/dataset.py b/stages/dataset.py index 924b8c0a8..966ed4a05 100644 --- a/stages/dataset.py +++ b/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: diff --git a/stages/run_opensfm.py b/stages/run_opensfm.py index 850c6ce83..be4da1f9f 100644 --- a/stages/run_opensfm.py +++ b/stages/run_opensfm.py @@ -2,6 +2,7 @@ import os import shutil import glob +import json from opendm import log from opendm import io @@ -99,6 +100,24 @@ def cleanup_disk_space(): largest_photo = None undistort_pipeline = [] + # Load per-band panel irradiance computed during the dataset stage + # (see stages/dataset.py). Only used for --radiometric-calibration camera+panel. + panel_irradiance_by_band = None + if args.radiometric_calibration == "camera+panel": + panels_file = os.path.join(tree.root_path, 'panels.json') + if os.path.exists(panels_file): + try: + with open(panels_file) as f: + panel_irradiance_by_band = json.load(f).get('irradiance', None) + except Exception as e: + log.ODM_WARNING("Cannot read panel irradiance database: %s" % str(e)) + + if panel_irradiance_by_band: + log.ODM_INFO("Using panel-derived irradiance for band(s): %s" % ", ".join(panel_irradiance_by_band.keys())) + else: + log.ODM_WARNING("camera+panel selected but no panel irradiance is available; " + "reflectance will fall back to stored/DLS irradiance.") + def undistort_callback(shot_id, image): for func in undistort_pipeline: image = func(shot_id, image) @@ -116,7 +135,12 @@ def radiometric_calibrate(shot_id, image): if photo.is_thermal(): return thermal.dn_to_temperature(photo, image, tree.dataset_raw) else: - return multispectral.dn_to_reflectance(photo, image, use_sun_sensor=args.radiometric_calibration=="camera+sun") + panel_irradiance = None + if panel_irradiance_by_band is not None: + panel_irradiance = panel_irradiance_by_band.get(photo.band_name) + return multispectral.dn_to_reflectance(photo, image, + use_sun_sensor=args.radiometric_calibration=="camera+sun", + panel_irradiance=panel_irradiance) def align_to_primary_band(shot_id, image):