From a61b5b308a5692ef0bce12d977ab4ff15e506c16 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Fri, 6 Mar 2020 13:04:34 -0500 Subject: [PATCH] More radiometric calibration Former-commit-id: f51b204205571107c3c6084bda1f9b5ee23fd2f7 --- opendm/multispectral.py | 74 +++++++++++++++++++++++++---------------- opendm/photo.py | 45 ++++++++++++++++++++----- stages/run_opensfm.py | 62 +++++++++++++++++----------------- 3 files changed, 113 insertions(+), 68 deletions(-) diff --git a/opendm/multispectral.py b/opendm/multispectral.py index 782e9616..2e40f1a5 100644 --- a/opendm/multispectral.py +++ b/opendm/multispectral.py @@ -1,5 +1,9 @@ from opendm import dls import math +import numpy as np +from opendm import log + +import cv2 # TODO REMOVE # Loosely based on https://github.com/micasense/imageprocessing/blob/master/micasense/utils.py def dn_to_radiance(photo, image): @@ -9,10 +13,13 @@ def dn_to_radiance(photo, image): :param image numpy array containing image data :return numpy array with radiance image values """ + + image = image.astype(float) + # Handle thermal bands (experimental) if photo.band_name == 'LWIR': image -= (273.15 * 100.0) # Convert Kelvin to Celsius - image = image.astype(float) * 0.01 + image *= 0.01 return image # All others @@ -21,6 +28,14 @@ def dn_to_radiance(photo, image): exposure_time = photo.exposure_time gain = photo.get_gain() + photometric_exp = photo.get_photometric_exposure() + + if a1 is None and photometric_exp is None: + log.ODM_WARNING("Cannot perform radiometric calibration, no FNumber/Exposure Time or Radiometric Calibration EXIF tags found in %s. Using Digital Number." % photo.filename) + return image + + if a1 is None and photometric_exp is not None: + a1 = photometric_exp V, x, y = vignette_map(photo) if x is None: @@ -28,25 +43,9 @@ def dn_to_radiance(photo, image): if dark_level is not None: image -= dark_level + print("Adjusted black") - if V is not None: - # vignette correction - image *= V - - if exposure_time and a2 is not None and a3 is not None: - # row gradient correction - R = 1.0 / (1.0 + a2 * y / exposure_time - a3 * y) - image *= R - - # Floor any negative radiances to zero (can happend due to noise around blackLevel) - if dark_level is not None: - image[image < 0] = 0 - - # apply the radiometric calibration - i.e. scale by the gain-exposure product and - # multiply with the radiometric calibration coefficient - # need to normalize by 2^16 for 16 bit images - # because coefficients are scaled to work with input values of max 1.0 - + # Normalize DN to 0 - 1.0 bps = photo.bits_per_sample if bps: bit_depth_max = float(2 ** bps) @@ -54,16 +53,34 @@ def dn_to_radiance(photo, image): # Infer from array dtype info = np.iinfo(image.dtype) bit_depth_max = info.max - info.min + image /= bit_depth_max + + if V is not None: + # vignette correction + image *= V + print("Adjusted vignette") + + if exposure_time and a2 is not None and a3 is not None: + # row gradient correction + R = 1.0 / (1.0 + a2 * y / exposure_time - a3 * y) + image *= R + + cv2.imwrite("/datasets/mica/R.tif", R) + + print("Row gradient") - image = image.astype(float) + # Floor any negative radiances to zero (can happend due to noise around blackLevel) + if dark_level is not None: + image[image < 0] = 0 + + # apply the radiometric calibration - i.e. scale by the gain-exposure product and + # multiply with the radiometric calibration coefficient if gain is not None and exposure_time is not None: image /= (gain * exposure_time) + print("Gain adjustment") - if a1 is not None: - image *= a1 - - image /= bit_depth_max + image *= a1 return image @@ -72,8 +89,7 @@ def vignette_map(photo): polynomial = photo.get_vignetting_polynomial() if x_vc and polynomial: - # reverse list and append 1., so that we can call with numpy polyval - polynomial.reverse() + # append 1., so that we can call with numpy polyval polynomial.append(1.0) vignette_poly = np.array(polynomial) @@ -82,8 +98,8 @@ def vignette_map(photo): x, y = np.meshgrid(np.arange(photo.width), np.arange(photo.height)) # meshgrid returns transposed arrays - x = x.T - y = y.T + # x = x.T + # y = y.T # compute matrix of distances from image center r = np.hypot((x - x_vc), (y - y_vc)) @@ -116,7 +132,7 @@ def compute_irradiance_scale_factor(photo, use_sun_sensor=True): sun_vector_ned, sensor_vector_ned, sun_sensor_angle, \ solar_elevation, solar_azimuth = dls.compute_sun_angle([photo.latitude, photo.longitude], (0,0,0), # TODO: add support for sun sensor pose - photo.utc_time, + photo.get_utc_time(), dls_orientation_vector) angular_correction = dls.fresnel(sun_sensor_angle) diff --git a/opendm/photo.py b/opendm/photo.py index b791e02c..ad984c4d 100644 --- a/opendm/photo.py +++ b/opendm/photo.py @@ -36,6 +36,7 @@ class ODM_Photo: self.band_index = 0 # Multi-spectral fields + self.fnumber = None self.radiometric_calibration = None self.black_level = None @@ -91,8 +92,18 @@ class ODM_Photo: if 'EXIF ExposureTime' in tags: self.exposure_time = self.float_value(tags['EXIF ExposureTime']) + + if 'EXIF FNumber' in tags: + self.fnumber = self.float_value(tags['EXIF FNumber']) + if 'EXIF ISOSpeed' in tags: self.iso_speed = self.int_value(tags['EXIF ISOSpeed']) + elif 'EXIF PhotographicSensitivity' in tags: + self.iso_speed = self.int_value(tags['EXIF PhotographicSensitivity']) + elif 'EXIF ISOSpeedRatings' in tags: + self.iso_speed = self.int_value(tags['EXIF ISOSpeedRatings']) + + if 'Image BitsPerSample' in tags: self.bits_per_sample = self.int_value(tags['Image BitsPerSample']) if 'EXIF DateTimeOriginal' in tags: @@ -169,7 +180,8 @@ class ODM_Photo: # print(self.bits_per_sample) # print(self.vignetting_center) # print(self.sun_sensor) - # exit(1) + #print(self.get_vignetting_polynomial()) + #exit(1) self.width, self.height = get_image_size.get_image_size(_path_file) # Sanitize band name since we use it in folder paths @@ -233,7 +245,10 @@ class ODM_Photo: ) def float_values(self, tag): - return map(lambda v: float(v.num) / float(v.den), tag.values) + if isinstance(tag.values, list): + return map(lambda v: float(v.num) / float(v.den), tag.values) + else: + return [float(tag.values.num) / float(tag.values.den)] def float_value(self, tag): v = self.float_values(tag) @@ -241,7 +256,10 @@ class ODM_Photo: return v[0] def int_values(self, tag): - return map(int, tag.values) + if isinstance(tag.values, list): + return map(int, tag.values) + else: + return [int(tag.values)] def int_value(self, tag): v = self.int_values(tag) @@ -263,12 +281,11 @@ class ODM_Photo: if self.black_level: levels = np.array([float(v) for v in self.black_level.split(" ")]) return levels.mean() - return None def get_gain(self): + #(gain = ISO/100) if self.iso_speed: return self.iso_speed / 100.0 - return None def get_vignetting_center(self): if self.vignetting_center: @@ -281,6 +298,18 @@ class ODM_Photo: if self.vignetting_polynomial: parts = self.vignetting_polynomial.split(" ") if len(parts) > 0: - return list(map(float, parts)) - - return None \ No newline at end of file + coeffs = list(map(float, parts)) + + # Different camera vendors seem to use different ordering for the coefficients + if self.camera_make != "Sentera": + coeffs.reverse() + return coeffs + + def get_utc_time(self): + if self.utc_time: + return datetime.utcfromtimestamp(self.utc_time / 1000) + + def get_photometric_exposure(self): + # H ~= (exposure_time) / (f_number^2) + if self.fnumber is not None and self.exposure_time > 0: + return self.exposure_time / (self.fnumber * self.fnumber) \ No newline at end of file diff --git a/stages/run_opensfm.py b/stages/run_opensfm.py index 553e3ea5..1694b42c 100644 --- a/stages/run_opensfm.py +++ b/stages/run_opensfm.py @@ -9,7 +9,7 @@ from opendm import gsd from opendm import point_cloud from opendm import types from opendm.osfm import OSFMContext -from opendm.multispectral import dn_to_radiance +from opendm import multispectral class ODMOpenSfMStage(types.ODM_Stage): def process(self, args, outputs): @@ -22,40 +22,40 @@ class ODMOpenSfMStage(types.ODM_Stage): exit(1) octx = OSFMContext(tree.opensfm) - # octx.setup(args, tree.dataset_raw, photos, reconstruction=reconstruction, rerun=self.rerun()) - # octx.extract_metadata(self.rerun()) - # self.update_progress(20) - # octx.feature_matching(self.rerun()) - # self.update_progress(30) - # octx.reconstruct(self.rerun()) - # octx.extract_cameras(tree.path("cameras.json"), self.rerun()) - # self.update_progress(70) + octx.setup(args, tree.dataset_raw, photos, reconstruction=reconstruction, rerun=self.rerun()) + octx.extract_metadata(self.rerun()) + self.update_progress(20) + octx.feature_matching(self.rerun()) + self.update_progress(30) + octx.reconstruct(self.rerun()) + octx.extract_cameras(tree.path("cameras.json"), self.rerun()) + self.update_progress(70) - # # If we find a special flag file for split/merge we stop right here - # if os.path.exists(octx.path("split_merge_stop_at_reconstruction.txt")): - # log.ODM_INFO("Stopping OpenSfM early because we found: %s" % octx.path("split_merge_stop_at_reconstruction.txt")) - # self.next_stage = None - # return + # If we find a special flag file for split/merge we stop right here + if os.path.exists(octx.path("split_merge_stop_at_reconstruction.txt")): + log.ODM_INFO("Stopping OpenSfM early because we found: %s" % octx.path("split_merge_stop_at_reconstruction.txt")) + self.next_stage = None + return - # if args.fast_orthophoto: - # output_file = octx.path('reconstruction.ply') - # elif args.use_opensfm_dense: - # output_file = tree.opensfm_model - # else: - # output_file = tree.opensfm_reconstruction + if args.fast_orthophoto: + output_file = octx.path('reconstruction.ply') + elif args.use_opensfm_dense: + output_file = tree.opensfm_model + else: + output_file = tree.opensfm_reconstruction - # updated_config_flag_file = octx.path('updated_config.txt') + updated_config_flag_file = octx.path('updated_config.txt') - # # Make sure it's capped by the depthmap-resolution arg, - # # since the undistorted images are used for MVS - # outputs['undist_image_max_size'] = max( - # gsd.image_max_size(photos, args.orthophoto_resolution, tree.opensfm_reconstruction, ignore_gsd=args.ignore_gsd, has_gcp=reconstruction.has_gcp()), - # args.depthmap_resolution - # ) + # Make sure it's capped by the depthmap-resolution arg, + # since the undistorted images are used for MVS + outputs['undist_image_max_size'] = max( + gsd.image_max_size(photos, args.orthophoto_resolution, tree.opensfm_reconstruction, ignore_gsd=args.ignore_gsd, has_gcp=reconstruction.has_gcp()), + args.depthmap_resolution + ) - # if not io.file_exists(updated_config_flag_file) or self.rerun(): - # octx.update_config({'undistorted_image_max_size': outputs['undist_image_max_size']}) - # octx.touch(updated_config_flag_file) + if not io.file_exists(updated_config_flag_file) or self.rerun(): + octx.update_config({'undistorted_image_max_size': outputs['undist_image_max_size']}) + octx.touch(updated_config_flag_file) # These will be used for texturing / MVS if args.radiometric_calibration == "none": @@ -63,7 +63,7 @@ class ODMOpenSfMStage(types.ODM_Stage): else: def radiometric_calibrate(shot_id, image): photo = reconstruction.get_photo(shot_id) - return dn_to_reflectance(photo, image, use_sun_sensor=args.radiometric_calibration=="camera+sun") + return multispectral.dn_to_reflectance(photo, image, use_sun_sensor=args.radiometric_calibration=="camera+sun") octx.convert_and_undistort(self.rerun(), radiometric_calibrate)