From 94c8eaecfd43eb0d211f3efd107ed334eab1463f Mon Sep 17 00:00:00 2001
From: Piero Toffanin <pt@masseranolabs.com>
Date: Thu, 5 Mar 2020 15:39:16 +0000
Subject: [PATCH] Radiance calculation

Former-commit-id: 5baf67e81a0cfaac618cc3f1707d81144cccdc80
---
 opendm/multispectral.py | 85 +++++++++++++++++++++++++++++++++++++++--
 opendm/photo.py         | 39 ++++++++++++++++++-
 2 files changed, 120 insertions(+), 4 deletions(-)

diff --git a/opendm/multispectral.py b/opendm/multispectral.py
index 2483e2f3..3d5ddbc1 100644
--- a/opendm/multispectral.py
+++ b/opendm/multispectral.py
@@ -1,3 +1,5 @@
+# Loosely based on https://github.com/micasense/imageprocessing/blob/master/micasense/utils.py
+
 def dn_to_radiance(photo, image):
     """
     Convert Digital Number values to Radiance values
@@ -5,6 +7,83 @@ def dn_to_radiance(photo, image):
     :param image numpy array containing image data
     :return numpy array with radiance image values
     """
-    # TODO!
-    print("HI")
-    return image
\ No newline at end of file
+
+    a1, a2, a3 = photo.get_radiometric_calibration()
+    dark_level = photo.get_dark_level()
+
+    exposure_time = photo.exposure_time
+    gain = photo.get_gain()
+
+    V, x, y = vignette_map(photo)
+    if x is None:
+        x, y = np.meshgrid(np.arange(photo.width), np.arange(photo.height))
+
+    if dark_level is not None:
+        image -= dark_level
+
+    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
+
+    bps = photo.bits_per_sample
+    if bps:
+        bit_depth_max = float(2 ** bps)
+    else:
+        # Infer from array dtype
+        info = np.iinfo(image.dtype)
+        bit_depth_max = info.max - info.min
+    
+    image = image.astype(float)
+
+    if gain is not None and exposure_time is not None:
+        image /= (gain * exposure_time)
+    
+    if a1 is not None:
+        image *= a1
+    
+    image /= bit_depth_max
+    
+    return image
+
+def vignette_map(photo):
+    x_vc, y_vc = photo.get_vignetting_center()
+    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()
+        polynomial.append(1.0)
+        vignette_poly = np.array(polynomial)
+
+        # perform vignette correction
+        # get coordinate grid across image
+        x, y = np.meshgrid(np.arange(photo.width), np.arange(photo.height))
+
+        # meshgrid returns transposed arrays
+        x = x.T
+        y = y.T
+
+        # compute matrix of distances from image center
+        r = np.hypot((x - x_vc), (y - y_vc))
+
+        # compute the vignette polynomial for each distance - we divide by the polynomial so that the
+        # corrected image is image_corrected = image_original * vignetteCorrection
+
+        vignette = 1.0 / np.polyval(vignette_poly, r)
+        return vignette, x, y
+    
+    return None, None, None
diff --git a/opendm/photo.py b/opendm/photo.py
index b421c307..f907784d 100644
--- a/opendm/photo.py
+++ b/opendm/photo.py
@@ -3,6 +3,7 @@ import logging
 import re
 
 import exifread
+import numpy as np
 from six import string_types
 
 import log
@@ -76,9 +77,11 @@ class ODM_Photo:
                 if 'GPS GPSLongitude' in tags and 'GPS GPSLongitudeRef' in tags:
                     self.longitude = self.dms_to_decimal(tags['GPS GPSLongitude'], tags['GPS GPSLongitudeRef'])
                 
-                # if 'BlackLevel' in tags:
                 if 'Image Tag 0xC61A' in tags:
                     self.black_level = self.list_values(tags['Image Tag 0xC61A'])
+                elif 'BlackLevel' in tags:
+                    self.black_level = self.list_values(tags['BlackLevel'])
+                
                 if 'EXIF ExposureTime' in tags:
                     self.exposure_time = self.float_value(tags['EXIF ExposureTime'])
                 if 'EXIF ISOSpeed' in tags:
@@ -126,6 +129,7 @@ class ODM_Photo:
             # print(self.bits_per_sample)
             # print(self.vignetting_center)
             # print(self.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
@@ -202,3 +206,36 @@ class ODM_Photo:
     def list_values(self, tag):
         return " ".join(map(str, tag.values))
 
+    def get_radiometric_calibration(self):
+        if self.radiometric_calibration:
+            parts = self.radiometric_calibration.split(" ")
+            if len(parts) == 3:
+                return list(map(float, parts))
+
+        return [None, None, None]                
+    
+    def get_dark_level(self):
+        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):
+        if self.iso_speed:
+            return self.iso_speed / 100.0
+        return None
+
+    def get_vignetting_center(self):
+        if self.vignetting_center:
+            parts = self.vignetting_center.split(" ")
+            if len(parts) == 2:
+                return list(map(float, parts))
+        return [None, None]
+
+    def get_vignetting_polynomial(self):
+        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