Reflectance calculations

Former-commit-id: d9089e48dd
pull/1161/head
Piero Toffanin 2020-03-05 16:34:16 -05:00
rodzic 94c8eaecfd
commit 0e8a325996
6 zmienionych plików z 256 dodań i 13 usunięć

Wyświetl plik

@ -163,11 +163,11 @@ def config():
parser.add_argument('--radiometric-calibration',
metavar='<string>',
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: '
'Can be set to one of: [none, camera, camera+sun]. Default: '
'%(default)s'))

145
opendm/dls.py 100644
Wyświetl plik

@ -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

Wyświetl plik

@ -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()
@ -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

Wyświetl plik

@ -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,6 +95,23 @@ 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))
@ -120,6 +144,22 @@ class ODM_Photo:
'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)
# print(self.radiometric_calibration)
@ -128,17 +168,20 @@ 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:
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):
@ -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

Wyświetl plik

@ -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

Wyświetl plik

@ -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)