From 0e8a325996f88ecacffad099dcf7e166ba0fe8f4 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 5 Mar 2020 16:34:16 -0500 Subject: [PATCH] Reflectance calculations Former-commit-id: d9089e48dddbcc2749483ee8a78c8bda3a2de7b0 --- opendm/config.py | 6 +- opendm/dls.py | 145 ++++++++++++++++++++++++++++++++++++++++ opendm/multispectral.py | 56 +++++++++++++++- opendm/photo.py | 59 ++++++++++++++-- requirements.txt | 1 + stages/run_opensfm.py | 2 +- 6 files changed, 256 insertions(+), 13 deletions(-) create mode 100644 opendm/dls.py diff --git a/opendm/config.py b/opendm/config.py index 16aa062a..b2607254 100644 --- a/opendm/config.py +++ b/opendm/config.py @@ -163,11 +163,11 @@ def config(): parser.add_argument('--radiometric-calibration', metavar='', default='none', - choices=['none', 'camera', 'sun', 'sunangle'], + choices=['none', 'camera', 'camera+sun'], help=('Set the radiometric calibration to perform on images. ' 'When processing multispectral images you should set this option ' - 'to obtain reflectance values (otherwise you will get digital number (DN) values). ' - 'Can be set to one of: [none, camera, sun, sunangle]. Default: ' + 'to obtain reflectance values (otherwise you will get digital number (DN) values). ' + 'Can be set to one of: [none, camera, camera+sun]. Default: ' '%(default)s')) diff --git a/opendm/dls.py b/opendm/dls.py new file mode 100644 index 00000000..7e0bf980 --- /dev/null +++ b/opendm/dls.py @@ -0,0 +1,145 @@ +""" +MicaSense Downwelling Light Sensor Utilities + +Copyright 2017 MicaSense, Inc. + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in the +Software without restriction, including without limitation the rights to use, +copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the +Software, and to permit persons to whom the Software is furnished to do so, +subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +""" + + +import numpy as np +# for DLS correction, we need the sun position at the time the image was taken +# this can be computed using the pysolar package (ver 0.6) +# https://pypi.python.org/pypi/Pysolar/0.6 +# we import multiple times with checking here because the case of Pysolar is +# different depending on the python version :( +import imp + +havePysolar = False + +try: + import pysolar.solar as pysolar + havePysolar = True +except ImportError: + try: + import Pysolar.solar as pysolar + havePysolar = True + except ImportError: + import pysolar.solar as pysolar + havePysolar = True +finally: + if not havePysolar: + print("Unable to import pysolar") + +def fresnel(phi): + return __multilayer_transmission(phi, n=[1.000277,1.6,1.38]) + +# define functions to compute the DLS-Sun angle: +def __fresnel_transmission(phi, n1=1.000277, n2=1.38, polarization=[.5, .5]): + """compute fresnel transmission between media with refractive indices n1 and n2""" + # computes the reflection and transmittance + # for incidence angles phi for transition from medium + # with refractive index n1 to n2 + # teflon e.g. n2=1.38 + # polycarbonate n2=1.6 + # polarization=[.5,.5] - unpolarized light + # polarization=[1.,0] - s-polarized light - perpendicular to plane of incidence + # polarization=[0,1.] - p-polarized light - parallel to plane of incidence + f1 = np.cos(phi) + f2 = np.sqrt(1-(n1/n2*np.sin(phi))**2) + Rs = ((n1*f1-n2*f2)/(n1*f1+n2*f2))**2 + Rp = ((n1*f2-n2*f1)/(n1*f2+n2*f1))**2 + T = 1.-polarization[0]*Rs-polarization[1]*Rp + if T > 1: T= 0. + if T < 0: T = 0. + if np.isnan(T): T = 0. + return T + +def __multilayer_transmission(phi, n, polarization=[.5, .5]): + T = 1.0 + phi_eff = np.copy(phi) + for i in range(0,len(n)-1): + n1 = n[i] + n2 = n[i+1] + phi_eff = np.arcsin(np.sin(phi_eff)/n1) + T *= __fresnel_transmission(phi_eff, n1, n2, polarization=polarization) + return T + +# get the position of the sun in North-East-Down (NED) coordinate system +def ned_from_pysolar(sunAzimuth, sunAltitude): + """Convert pysolar coordinates to NED coordinates.""" + elements = ( + np.cos(sunAzimuth) * np.cos(sunAltitude), + np.sin(sunAzimuth) * np.cos(sunAltitude), + -np.sin(sunAltitude), + ) + return np.array(elements).transpose() + +# get the sensor orientation in North-East-Down coordinates +# pose is a yaw/pitch/roll tuple of angles measured for the DLS +# ori is the 3D orientation vector of the DLS in body coordinates (typically [0,0,-1]) +def get_orientation(pose, ori): + """Generate an orientation vector from yaw/pitch/roll angles in radians.""" + yaw, pitch, roll = pose + c1 = np.cos(-yaw) + s1 = np.sin(-yaw) + c2 = np.cos(-pitch) + s2 = np.sin(-pitch) + c3 = np.cos(-roll) + s3 = np.sin(-roll) + Ryaw = np.array([[c1, s1, 0], [-s1, c1, 0], [0, 0, 1]]) + Rpitch = np.array([[c2, 0, -s2], [0, 1, 0], [s2, 0, c2]]) + Rroll = np.array([[1, 0, 0], [0, c3, s3], [0, -s3, c3]]) + R = np.dot(Ryaw, np.dot(Rpitch, Rroll)) + n = np.dot(R, ori) + return n + +# from the current position (lat,lon,alt) tuple +# and time (UTC), as well as the sensor orientation (yaw,pitch,roll) tuple +# compute a sensor sun angle - this is needed as the actual sun irradiance +# (for clear skies) is related to the measured irradiance by: + +# I_measured = I_direct * cos (sun_sensor_angle) + I_diffuse +# For clear sky, I_direct/I_diffuse ~ 6 and we can simplify this to +# I_measured = I_direct * (cos (sun_sensor_angle) + 1/6) + +def compute_sun_angle( + position, + pose, + utc_datetime, + sensor_orientation, +): + """ compute the sun angle using pysolar functions""" + altitude = 0 + azimuth = 0 + import warnings + with warnings.catch_warnings(): # Ignore pysolar leap seconds offset warning + warnings.simplefilter("ignore") + try: + altitude = pysolar.get_altitude(position[0], position[1], utc_datetime) + azimuth = pysolar.get_azimuth(position[0], position[1], utc_datetime) + except AttributeError: # catch 0.6 version of pysolar required for python 2.7 support + altitude = pysolar.GetAltitude(position[0], position[1], utc_datetime) + azimuth = 180-pysolar.GetAzimuth(position[0], position[1], utc_datetime) + sunAltitude = np.radians(np.array(altitude)) + sunAzimuth = np.radians(np.array(azimuth)) + sunAzimuth = sunAzimuth % (2 * np.pi ) #wrap range 0 to 2*pi + nSun = ned_from_pysolar(sunAzimuth, sunAltitude) + nSensor = np.array(get_orientation(pose, sensor_orientation)) + angle = np.arccos(np.dot(nSun, nSensor)) + return nSun, nSensor, angle, sunAltitude, sunAzimuth diff --git a/opendm/multispectral.py b/opendm/multispectral.py index 3d5ddbc1..782e9616 100644 --- a/opendm/multispectral.py +++ b/opendm/multispectral.py @@ -1,3 +1,5 @@ +from opendm import dls +import math # Loosely based on https://github.com/micasense/imageprocessing/blob/master/micasense/utils.py def dn_to_radiance(photo, image): @@ -7,7 +9,13 @@ def dn_to_radiance(photo, image): :param image numpy array containing image data :return numpy array with radiance image values """ - + # Handle thermal bands (experimental) + if photo.band_name == 'LWIR': + image -= (273.15 * 100.0) # Convert Kelvin to Celsius + image = image.astype(float) * 0.01 + return image + + # All others a1, a2, a3 = photo.get_radiometric_calibration() dark_level = photo.get_dark_level() @@ -56,7 +64,7 @@ def dn_to_radiance(photo, image): image *= a1 image /= bit_depth_max - + return image def vignette_map(photo): @@ -87,3 +95,47 @@ def vignette_map(photo): return vignette, x, y return None, None, None + +def dn_to_reflectance(photo, image, use_sun_sensor=True): + radiance = dn_to_radiance(photo, image) + irradiance_scale = compute_irradiance_scale_factor(photo, use_sun_sensor=use_sun_sensor) + return radiance * irradiance_scale + +def compute_irradiance_scale_factor(photo, use_sun_sensor=True): + # Thermal? + if photo.band_name == "LWIR": + return 1.0 + + if photo.irradiance is not None: + horizontal_irradiance = photo.irradiance + return math.pi / horizontal_irradiance + + if use_sun_sensor: + # Estimate it + dls_orientation_vector = np.array([0,0,-1]) + 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, + dls_orientation_vector) + + angular_correction = dls.fresnel(sun_sensor_angle) + + # TODO: support for direct and scattered irradiance + + direct_to_diffuse_ratio = 6.0 # Assumption + spectral_irradiance = photo.sun_sensor # TODO: support for XMP:SpectralIrradiance + + percent_diffuse = 1.0 / direct_to_diffuse_ratio + sensor_irradiance = spectral_irradiance / angular_correction + + # find direct irradiance in the plane normal to the sun + untilted_direct_irr = sensor_irradiance / (percent_diffuse + np.cos(sun_sensor_angle)) + direct_irradiance = untilted_direct_irr + scattered_irradiance = untilted_direct_irr*percent_diffuse + + # compute irradiance on the ground using the solar altitude angle + horizontal_irradiance = direct_irradiance * np.sin(solar_elevation) + scattered_irradiance + return math.pi / horizontal_irradiance + + return 1.0 \ No newline at end of file diff --git a/opendm/photo.py b/opendm/photo.py index f907784d..b791e02c 100644 --- a/opendm/photo.py +++ b/opendm/photo.py @@ -5,6 +5,8 @@ import re import exifread import numpy as np from six import string_types +from datetime import datetime, timedelta +import pytz import log import system @@ -33,17 +35,22 @@ class ODM_Photo: self.band_name = 'RGB' self.band_index = 0 - # Multi-spectral fields (not all cameras implement all these values) + # Multi-spectral fields self.radiometric_calibration = None self.black_level = None - # Capture info (most cameras implement these) + # Capture info self.exposure_time = None self.iso_speed = None self.bits_per_sample = None self.vignetting_center = None self.vignetting_polynomial = None + self.irradiance = None + self.sun_sensor = None + self.utc_time = None + # self.center_wavelength = None + # self.bandwidth = None # parse values from metadata self.parse_exif_values(path_file) @@ -88,7 +95,24 @@ class ODM_Photo: self.iso_speed = self.int_value(tags['EXIF ISOSpeed']) if 'Image BitsPerSample' in tags: self.bits_per_sample = self.int_value(tags['Image BitsPerSample']) - + if 'EXIF DateTimeOriginal' in tags: + str_time = tags['EXIF DateTimeOriginal'].values.encode('utf8') + utc_time = datetime.strptime(str_time, "%Y:%m:%d %H:%M:%S") + subsec = 0 + if 'EXIF SubSecTime' in tags: + subsec = self.int_value(tags['EXIF SubSecTime']) + negative = 1.0 + if subsec < 0: + negative = -1.0 + subsec *= -1.0 + subsec = float('0.{}'.format(int(subsec))) + subsec *= negative + ms = subsec * 1e3 + utc_time += timedelta(milliseconds = ms) + timezone = pytz.timezone('UTC') + epoch = timezone.localize(datetime.utcfromtimestamp(0)) + self.utc_time = (timezone.localize(utc_time) - epoch).total_seconds() * 1000.0 + except IndexError as e: log.ODM_WARNING("Cannot read EXIF tags for %s: %s" % (_path_file, e.message)) @@ -119,6 +143,22 @@ class ODM_Photo: 'Camera:VignettingPolynomial', 'Sentera:VignettingPolynomial', ]) + + self.set_attr_from_xmp_tag('irradiance', tags, [ + 'Camera:Irradiance' + ], float) + + self.set_attr_from_xmp_tag('sun_sensor', tags, [ + 'Camera:SunSensor' + ], float) + + # self.set_attr_from_xmp_tag('center_wavelength', tags, [ + # 'Camera:CentralWavelength' + # ], float) + + # self.set_attr_from_xmp_tag('bandwidth', tags, [ + # 'Camera:WavelengthFWHM' + # ], float) # print(self.band_name) # print(self.band_index) @@ -128,18 +168,21 @@ class ODM_Photo: # print(self.iso_speed) # print(self.bits_per_sample) # print(self.vignetting_center) - # print(self.vignetting_polynomial) + # print(self.sun_sensor) # exit(1) self.width, self.height = get_image_size.get_image_size(_path_file) # Sanitize band name since we use it in folder paths self.band_name = re.sub('[^A-Za-z0-9]+', '', self.band_name) - def set_attr_from_xmp_tag(self, attr, xmp_tags, tags): + def set_attr_from_xmp_tag(self, attr, xmp_tags, tags, cast=None): v = self.get_xmp_tag(xmp_tags, tags) if v is not None: - setattr(self, attr, v) - + if cast is None: + setattr(self, attr, v) + else: + setattr(self, attr, cast(v)) + def get_xmp_tag(self, xmp_tags, tags): if isinstance(tags, str): tags = [tags] @@ -153,6 +196,8 @@ class ODM_Photo: elif isinstance(t, dict): items = t.get('rdf:Seq', {}).get('rdf:li', {}) if items: + if isinstance(items, string_types): + return items return " ".join(items) elif isinstance(t, int) or isinstance(t, float): return t diff --git a/requirements.txt b/requirements.txt index ed110215..238ecac0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -13,6 +13,7 @@ networkx==2.2 scipy==1.2.1 numpy==1.15.4 pyproj==2.2.2 +Pysolar==0.6 psutil==5.6.3 joblib==0.13.2 Fiona==1.8.9.post2 diff --git a/stages/run_opensfm.py b/stages/run_opensfm.py index a65537a1..553e3ea5 100644 --- a/stages/run_opensfm.py +++ b/stages/run_opensfm.py @@ -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_radiance(photo, image) + return dn_to_reflectance(photo, image, use_sun_sensor=args.radiometric_calibration=="camera+sun") octx.convert_and_undistort(self.rerun(), radiometric_calibrate)