More radiometric calibration

Former-commit-id: f51b204205

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:
+    # Normalize DN to 0 - 1.0
-        # 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
-
     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 *= a1
-    
-    image /= bit_depth_max
 
     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
+        # append 1., so that we can call with numpy polyval
-        polynomial.reverse()
         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
+        # x = x.T
-        y = y.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)

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))
+                coeffs = list(map(float, parts))
-        
+
-        return None
+                # 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)

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.setup(args, tree.dataset_raw, photos, reconstruction=reconstruction, rerun=self.rerun())
-        # octx.extract_metadata(self.rerun())
+        octx.extract_metadata(self.rerun())
-        # self.update_progress(20)
+        self.update_progress(20)
-        # octx.feature_matching(self.rerun())
+        octx.feature_matching(self.rerun())
-        # self.update_progress(30)
+        self.update_progress(30)
-        # octx.reconstruct(self.rerun())
+        octx.reconstruct(self.rerun())
-        # octx.extract_cameras(tree.path("cameras.json"), self.rerun())
+        octx.extract_cameras(tree.path("cameras.json"), self.rerun())
-        # self.update_progress(70)
+        self.update_progress(70)
 
-        # # If we find a special flag file for split/merge we stop right here
+        # 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")):
+        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"))
+            log.ODM_INFO("Stopping OpenSfM early because we found: %s" % octx.path("split_merge_stop_at_reconstruction.txt"))
-        #     self.next_stage = None
+            self.next_stage = None
-        #     return
+            return
 
-        # if args.fast_orthophoto:
+        if args.fast_orthophoto:
-        #     output_file = octx.path('reconstruction.ply')
+            output_file = octx.path('reconstruction.ply')
-        # elif args.use_opensfm_dense:
+        elif args.use_opensfm_dense:
-        #     output_file = tree.opensfm_model
+            output_file = tree.opensfm_model
-        # else:
+        else:
-        #     output_file = tree.opensfm_reconstruction
+            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,
+        # Make sure it's capped by the depthmap-resolution arg,
-        # # since the undistorted images are used for MVS
+        # since the undistorted images are used for MVS
-        # outputs['undist_image_max_size'] = max(
+        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()),
+            gsd.image_max_size(photos, args.orthophoto_resolution, tree.opensfm_reconstruction, ignore_gsd=args.ignore_gsd, has_gcp=reconstruction.has_gcp()),
-        #     args.depthmap_resolution
+            args.depthmap_resolution
-        # )
+        )
 
-        # if not io.file_exists(updated_config_flag_file) or self.rerun():
+        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.update_config({'undistorted_image_max_size': outputs['undist_image_max_size']})
-        #     octx.touch(updated_config_flag_file)
+            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)