Added support of radiometric calibration for DJI Mavic 2 Enterprize Advanced. Added methods for adding FLIR sensor into reconstruction

pull/1391/head
markFieldman 2021-12-29 17:11:09 +02:00
rodzic 3957278c2e
commit 27e6116977
7 zmienionych plików z 475 dodań i 4 usunięć

Wyświetl plik

@ -617,6 +617,9 @@ class ODM_Photo:
self.gps_xy_stddev = self.gps_z_stddev = dop
def is_thermal(self):
#Added for support M2EA camera sensor
if(self.camera_make == "DJI"):
return self.camera_model == "MAVIC2-ENTERPRISE-ADVANCED" and self.width == 640 and self.height == 512
return self.band_name.upper() in ["LWIR"] # TODO: more?
def camera_id(self):

Wyświetl plik

@ -1,4 +1,5 @@
from opendm import log
from opendm.thermal_tools import dji_unpack
import cv2
def resize_to_match(image, match_photo = None):
@ -17,16 +18,16 @@ def resize_to_match(image, match_photo = None):
interpolation=cv2.INTER_LANCZOS4)
return image
def dn_to_temperature(photo, image):
def dn_to_temperature(photo, image, dataset_tree):
"""
Convert Digital Number values to temperature (C) values
:param photo ODM_Photo
:param image numpy array containing image data
:param resize_to_photo ODM_Photo that photo should be resized to (to match its dimensions)
:param dataset_tree path to original source image to read data using PIL for DJI thermal photos
:return numpy array with temperature (C) image values
"""
image = image.astype("float32")
# Handle thermal bands
if photo.is_thermal():
@ -34,12 +35,18 @@ def dn_to_temperature(photo, image):
# The following will work for MicaSense Altum cameras
# but not necessarily for others
if photo.camera_make == "MicaSense" and photo.camera_model == "Altum":
image = image.astype("float32")
image -= (273.15 * 100.0) # Convert Kelvin to Celsius
image *= 0.01
return image
elif photo.camera_make == "DJI" and photo.camera_model == "MAVIC2-ENTERPRISE-ADVANCED":
image = dji_unpack.extract_temperatures_dji(photo, image, dataset_tree)
image = image.astype("float32")
return image
else:
log.ODM_WARNING("Unsupported camera [%s %s], thermal band will have digital numbers." % (photo.camera_make, photo.camera_model))
else:
image = image.astype("float32")
log.ODM_WARNING("Tried to radiometrically calibrate a non-thermal image with temperature values (%s)" % photo.filename)
return image

Wyświetl plik

@ -0,0 +1,51 @@
from PIL import Image
import numpy as np
from opendm import system
from opendm import log
from opendm.thermal_tools.thermal_utils import sensor_vals_to_temp
def extract_temperatures_dji(photo, image, dataset_tree):
"""Extracts the DJI-encoded thermal image as 2D floating-point numpy array with temperatures in degC.
The raw sensor values are obtained using the sample binaries provided in the official Thermal SDK by DJI.
The executable file is run and generates a 16 bit unsigned RAW image with Little Endian byte order.
Link to DJI Forum post: https://forum.dji.com/forum.php?mod=redirect&goto=findpost&ptid=230321&pid=2389016
"""
# Harcoded metadata for mean of values
# This is added to support the possibility of extracting RJPEG from DJI M2EA
meta = {
"Emissivity": 0.95,
"ObjectDistance": 50, #This is mean value of flights for better results. Need to be changed later, or improved by bypassing options from task broker
"AtmosphericTemperature": 20,
"ReflectedApparentTemperature": 30,
"IRWindowTemperature": 20,
"IRWindowTransmission": 1,
"RelativeHumidity": 40,
"PlanckR1": 21106.77,
"PlanckB": 1501,
"PlanckF": 1,
"PlanckO": -7340,
"PlanckR2": 0.012545258,
}
if photo.camera_model == "MAVIC2-ENTERPRISE-ADVANCED":
# Adding support for MAVIC2-ENTERPRISE-ADVANCED Camera images
im = Image.open(f"{dataset_tree}/{photo.filename}")
# concatenate APP3 chunks of data
a = im.applist[3][1]
for i in range(4, 14):
a += im.applist[i][1]
# create image from bytes
try:
img = Image.frombytes("I;16L", (640, 512), a)
except ValueError as e:
log.ODM_ERROR(e)
log.ODM_ERROR(photo.filename)
else:
log.ODM_DEBUG("Only DJI M2EA currently supported, please wait for new updates")
return image
# Extract raw sensor values from generated image into numpy array
raw_sensor_np = np.array(img)
## extracting the temperatures from thermal images
thermal_np = sensor_vals_to_temp(raw_sensor_np, **meta)
return thermal_np

Wyświetl plik

@ -0,0 +1,271 @@
"""
THIS IS WIP, DON'T USE THIS FILE, IT IS HERE FOR FURTHER IMPROVEMENT
Tools for extracting thermal data from FLIR images.
Derived from https://bitbucket.org/nimmerwoner/flyr/src/master/
"""
import os
from io import BufferedIOBase, BytesIO
from typing import BinaryIO, Dict, Optional, Tuple, Union
import numpy as np
from PIL import Image
# Constants
SEGMENT_SEP = b"\xff"
APP1_MARKER = b"\xe1"
MAGIC_FLIR_DEF = b"FLIR\x00"
CHUNK_APP1_BYTES_COUNT = len(APP1_MARKER)
CHUNK_LENGTH_BYTES_COUNT = 2
CHUNK_MAGIC_BYTES_COUNT = len(MAGIC_FLIR_DEF)
CHUNK_SKIP_BYTES_COUNT = 1
CHUNK_NUM_BYTES_COUNT = 1
CHUNK_TOT_BYTES_COUNT = 1
CHUNK_PARTIAL_METADATA_LENGTH = CHUNK_APP1_BYTES_COUNT + CHUNK_LENGTH_BYTES_COUNT + CHUNK_MAGIC_BYTES_COUNT
CHUNK_METADATA_LENGTH = (
CHUNK_PARTIAL_METADATA_LENGTH + CHUNK_SKIP_BYTES_COUNT + CHUNK_NUM_BYTES_COUNT + CHUNK_TOT_BYTES_COUNT
)
def unpack(path_or_stream: Union[str, BinaryIO]) -> np.ndarray:
"""Unpacks the FLIR image, meaning that it will return the thermal data embedded in the image.
Parameters
----------
path_or_stream : Union[str, BinaryIO]
Either a path (string) to a FLIR file, or a byte stream such as
BytesIO or file opened as `open(file_path, "rb")`.
Returns
-------
FlyrThermogram
When successful, a FlyrThermogram object containing thermogram data.
"""
if isinstance(path_or_stream, str) and os.path.isfile(path_or_stream):
with open(path_or_stream, "rb") as flirh:
return unpack(flirh)
elif isinstance(path_or_stream, BufferedIOBase):
stream = path_or_stream
flir_app1_stream = extract_flir_app1(stream)
flir_records = parse_flir_app1(flir_app1_stream)
raw_np = parse_thermal(flir_app1_stream, flir_records)
return raw_np
else:
raise ValueError("Incorrect input")
def extract_flir_app1(stream: BinaryIO) -> BinaryIO:
"""Extracts the FLIR APP1 bytes.
Parameters
---------
stream : BinaryIO
A full bytes stream of a JPEG file, expected to be a FLIR file.
Raises
------
ValueError
When the file is invalid in one the next ways, a
ValueError is thrown.
* File is not a JPEG
* A FLIR chunk number occurs more than once
* The total chunks count is inconsistent over multiple chunks
* No APP1 segments are successfully parsed
Returns
-------
BinaryIO
A bytes stream of the APP1 FLIR segments
"""
# Check JPEG-ness
_ = stream.read(2)
chunks_count: Optional[int] = None
chunks: Dict[int, bytes] = {}
while True:
b = stream.read(1)
if b == b"":
break
if b != SEGMENT_SEP:
continue
parsed_chunk = parse_flir_chunk(stream, chunks_count)
if not parsed_chunk:
continue
chunks_count, chunk_num, chunk = parsed_chunk
chunk_exists = chunks.get(chunk_num, None) is not None
if chunk_exists:
raise ValueError("Invalid FLIR: duplicate chunk number")
chunks[chunk_num] = chunk
# Encountered all chunks, break out of loop to process found metadata
if chunk_num == chunks_count:
break
if chunks_count is None:
raise ValueError("Invalid FLIR: no metadata encountered")
flir_app1_bytes = b""
for chunk_num in range(chunks_count + 1):
flir_app1_bytes += chunks[chunk_num]
flir_app1_stream = BytesIO(flir_app1_bytes)
flir_app1_stream.seek(0)
return flir_app1_stream
def parse_flir_chunk(stream: BinaryIO, chunks_count: Optional[int]) -> Optional[Tuple[int, int, bytes]]:
"""Parse flir chunk."""
# Parse the chunk header. Headers are as follows (definition with example):
#
# \xff\xe1<length: 2 bytes>FLIR\x00\x01<chunk nr: 1 byte><chunk count: 1 byte>
# \xff\xe1\xff\xfeFLIR\x00\x01\x01\x0b
#
# Meaning: Exif APP1, 65534 long, FLIR chunk 1 out of 12
marker = stream.read(CHUNK_APP1_BYTES_COUNT)
length_bytes = stream.read(CHUNK_LENGTH_BYTES_COUNT)
length = int.from_bytes(length_bytes, "big")
length -= CHUNK_METADATA_LENGTH
magic_flir = stream.read(CHUNK_MAGIC_BYTES_COUNT)
if not (marker == APP1_MARKER and magic_flir == MAGIC_FLIR_DEF):
# Seek back to just after byte b and continue searching for chunks
stream.seek(-len(marker) - len(length_bytes) - len(magic_flir), 1)
return None
stream.seek(1, 1) # skip 1 byte, unsure what it is for
chunk_num = int.from_bytes(stream.read(CHUNK_NUM_BYTES_COUNT), "big")
chunks_tot = int.from_bytes(stream.read(CHUNK_TOT_BYTES_COUNT), "big")
# Remember total chunks to verify metadata consistency
if chunks_count is None:
chunks_count = chunks_tot
if ( # Check whether chunk metadata is consistent
chunks_tot is None or chunk_num < 0 or chunk_num > chunks_tot or chunks_tot != chunks_count
):
raise ValueError(f"Invalid FLIR: inconsistent total chunks, should be 0 or greater, but is {chunks_tot}")
return chunks_tot, chunk_num, stream.read(length + 1)
def parse_thermal(stream: BinaryIO, records: Dict[int, Tuple[int, int, int, int]]) -> np.ndarray:
"""Parse thermal."""
RECORD_IDX_RAW_DATA = 1
raw_data_md = records[RECORD_IDX_RAW_DATA]
_, _, raw_data = parse_raw_data(stream, raw_data_md)
return raw_data
def parse_flir_app1(stream: BinaryIO) -> Dict[int, Tuple[int, int, int, int]]:
"""Parse flir app1."""
# 0x00 - string[4] file format ID = "FFF\0"
# 0x04 - string[16] file creator: seen "\0","MTX IR\0","CAMCTRL\0"
# 0x14 - int32u file format version = 100
# 0x18 - int32u offset to record directory
# 0x1c - int32u number of entries in record directory
# 0x20 - int32u next free index ID = 2
# 0x24 - int16u swap pattern = 0 (?)
# 0x28 - int16u[7] spares
# 0x34 - int32u[2] reserved
# 0x3c - int32u checksum
# 1. Read 0x40 bytes and verify that its contents equals AFF\0 or FFF\0
_ = stream.read(4)
# 2. Read FLIR record directory metadata (ref 3)
stream.seek(16, 1)
_ = int.from_bytes(stream.read(4), "big")
record_dir_offset = int.from_bytes(stream.read(4), "big")
record_dir_entries_count = int.from_bytes(stream.read(4), "big")
stream.seek(28, 1)
_ = int.from_bytes(stream.read(4), "big")
# 3. Read record directory (which is a FLIR record entry repeated
# `record_dir_entries_count` times)
stream.seek(record_dir_offset)
record_dir_stream = BytesIO(stream.read(32 * record_dir_entries_count))
# First parse the record metadata
record_details: Dict[int, Tuple[int, int, int, int]] = {}
for record_nr in range(record_dir_entries_count):
record_dir_stream.seek(0)
details = parse_flir_record_metadata(stream, record_nr)
if details:
record_details[details[1]] = details
# Then parse the actual records
# for (entry_idx, type, offset, length) in record_details:
# parse_record = record_parsers[type]
# stream.seek(offset)
# record = BytesIO(stream.read(length + 36)) # + 36 needed to find end
# parse_record(record, offset, length)
return record_details
def parse_flir_record_metadata(stream: BinaryIO, record_nr: int) -> Optional[Tuple[int, int, int, int]]:
"""Parse flir record metadata."""
# FLIR record entry (ref 3):
# 0x00 - int16u record type
# 0x02 - int16u record subtype: RawData 1=BE, 2=LE, 3=PNG; 1 for other record types
# 0x04 - int32u record version: seen 0x64,0x66,0x67,0x68,0x6f,0x104
# 0x08 - int32u index id = 1
# 0x0c - int32u record offset from start of FLIR data
# 0x10 - int32u record length
# 0x14 - int32u parent = 0 (?)
# 0x18 - int32u object number = 0 (?)
# 0x1c - int32u checksum: 0 for no checksum
entry = 32 * record_nr
stream.seek(entry)
record_type = int.from_bytes(stream.read(2), "big")
if record_type < 1:
return None
_ = int.from_bytes(stream.read(2), "big")
_ = int.from_bytes(stream.read(4), "big")
_ = int.from_bytes(stream.read(4), "big")
record_offset = int.from_bytes(stream.read(4), "big")
record_length = int.from_bytes(stream.read(4), "big")
_ = int.from_bytes(stream.read(4), "big")
_ = int.from_bytes(stream.read(4), "big")
_ = int.from_bytes(stream.read(4), "big")
return (entry, record_type, record_offset, record_length)
def parse_raw_data(stream: BinaryIO, metadata: Tuple[int, int, int, int]):
"""Parse raw data."""
(_, _, offset, length) = metadata
stream.seek(offset)
stream.seek(2, 1)
width = int.from_bytes(stream.read(2), "little")
height = int.from_bytes(stream.read(2), "little")
stream.seek(offset + 32)
# Read the bytes with the raw thermal data and decode using PIL
thermal_bytes = stream.read(length)
thermal_stream = BytesIO(thermal_bytes)
thermal_img = Image.open(thermal_stream)
thermal_np = np.array(thermal_img)
# Check shape
if thermal_np.shape != (height, width):
msg = "Invalid FLIR: metadata's width and height don't match thermal data's actual width\
and height ({} vs ({}, {})"
msg = msg.format(thermal_np.shape, height, width)
raise ValueError(msg)
# FLIR PNG data is in the wrong byte order, fix that
fix_byte_order = np.vectorize(lambda x: (x >> 8) + ((x & 0x00FF) << 8))
thermal_np = fix_byte_order(thermal_np)
return width, height, thermal_np

Wyświetl plik

@ -0,0 +1,139 @@
"""Thermal Image manipulation utilities."""
"""Based on https://github.com/detecttechnologies/thermal_base"""
import numpy as np
def sensor_vals_to_temp(
raw,
Emissivity=1.0,
ObjectDistance=1,
AtmosphericTemperature=20,
ReflectedApparentTemperature=20,
IRWindowTemperature=20,
IRWindowTransmission=1,
RelativeHumidity=50,
PlanckR1=21106.77,
PlanckB=1501,
PlanckF=1,
PlanckO=-7340,
PlanckR2=0.012545258,
**kwargs,):
"""Convert raw values from the thermographic sensor sensor to temperatures in °C. Tested for Flir and DJI cams."""
# this calculation has been ported to python from https://github.com/gtatters/Thermimage/blob/master/R/raw2temp.R
# a detailed explanation of what is going on here can be found there
# constants
ATA1 = 0.006569
ATA2 = 0.01262
ATB1 = -0.002276
ATB2 = -0.00667
ATX = 1.9
# transmission through window (calibrated)
emiss_wind = 1 - IRWindowTransmission
refl_wind = 0
# transmission through the air
h2o = (RelativeHumidity / 100) * np.exp(
1.5587
+ 0.06939 * (AtmosphericTemperature)
- 0.00027816 * (AtmosphericTemperature) ** 2
+ 0.00000068455 * (AtmosphericTemperature) ** 3
)
tau1 = ATX * np.exp(-np.sqrt(ObjectDistance / 2) * (ATA1 + ATB1 * np.sqrt(h2o))) + (1 - ATX) * np.exp(
-np.sqrt(ObjectDistance / 2) * (ATA2 + ATB2 * np.sqrt(h2o))
)
tau2 = ATX * np.exp(-np.sqrt(ObjectDistance / 2) * (ATA1 + ATB1 * np.sqrt(h2o))) + (1 - ATX) * np.exp(
-np.sqrt(ObjectDistance / 2) * (ATA2 + ATB2 * np.sqrt(h2o))
)
# radiance from the environment
raw_refl1 = PlanckR1 / (PlanckR2 * (np.exp(PlanckB / (ReflectedApparentTemperature + 273.15)) - PlanckF)) - PlanckO
# Reflected component
raw_refl1_attn = (1 - Emissivity) / Emissivity * raw_refl1
# Emission from atmosphere 1
raw_atm1 = (
PlanckR1 / (PlanckR2 * (np.exp(PlanckB / (AtmosphericTemperature + 273.15)) - PlanckF)) - PlanckO
)
# attenuation for atmospheric 1 emission
raw_atm1_attn = (1 - tau1) / Emissivity / tau1 * raw_atm1
# Emission from window due to its own temp
raw_wind = (
PlanckR1 / (PlanckR2 * (np.exp(PlanckB / (IRWindowTemperature + 273.15)) - PlanckF)) - PlanckO
)
# Componen due to window emissivity
raw_wind_attn = (
emiss_wind / Emissivity / tau1 / IRWindowTransmission * raw_wind
)
# Reflection from window due to external objects
raw_refl2 = (
PlanckR1 / (PlanckR2 * (np.exp(PlanckB / (ReflectedApparentTemperature + 273.15)) - PlanckF)) - PlanckO
)
# component due to window reflectivity
raw_refl2_attn = (
refl_wind / Emissivity / tau1 / IRWindowTransmission * raw_refl2
)
# Emission from atmosphere 2
raw_atm2 = (
PlanckR1 / (PlanckR2 * (np.exp(PlanckB / (AtmosphericTemperature + 273.15)) - PlanckF)) - PlanckO
)
# attenuation for atmospheric 2 emission
raw_atm2_attn = (
(1 - tau2) / Emissivity / tau1 / IRWindowTransmission / tau2 * raw_atm2
)
raw_obj = (
raw / Emissivity / tau1 / IRWindowTransmission / tau2
- raw_atm1_attn
- raw_atm2_attn
- raw_wind_attn
- raw_refl1_attn
- raw_refl2_attn
)
val_to_log = PlanckR1 / (PlanckR2 * (raw_obj + PlanckO)) + PlanckF
if any(val_to_log.ravel() < 0):
raise Exception("Image seems to be corrupted")
# temperature from radiance
return PlanckB / np.log(val_to_log) - 273.15
def parse_from_exif_str(temp_str):
"""String to float parser."""
# we assume degrees celsius for temperature, metres for length
if isinstance(temp_str, str):
return float(temp_str.split()[0])
return float(temp_str)
def normalize_temp_matrix(thermal_np):
"""Normalize a temperature matrix to the 0-255 uint8 image range."""
num = thermal_np - np.amin(thermal_np)
den = np.amax(thermal_np) - np.amin(thermal_np)
thermal_np = num / den
return thermal_np
def clip_temp_to_roi(thermal_np, thermal_roi_values):
"""
Given an RoI within a temperature matrix, this function clips the temperature values in the entire thermal.
Image temperature values above and below the max/min temperatures within the RoI are clipped to said max/min.
Args:
thermal_np (np.ndarray): Floating point array containing the temperature matrix.
thermal_roi_values (np.ndarray / list): Any iterable containing the temperature values within the RoI.
Returns:
np.ndarray: The clipped temperature matrix.
"""
maximum = np.amax(thermal_roi_values)
minimum = np.amin(thermal_roi_values)
thermal_np[thermal_np > maximum] = maximum
thermal_np[thermal_np < minimum] = minimum
return thermal_np
def scale_with_roi(thermal_np, thermal_roi_values):
"""Alias for clip_temp_to_roi, to be deprecated in the future."""
return clip_temp_to_roi(thermal_np, thermal_roi_values)

Wyświetl plik

@ -113,7 +113,7 @@ class ODMOpenSfMStage(types.ODM_Stage):
def radiometric_calibrate(shot_id, image):
photo = reconstruction.get_photo(shot_id)
if photo.is_thermal():
return thermal.dn_to_temperature(photo, image)
return thermal.dn_to_temperature(photo, image, tree.dataset_raw)
else:
return multispectral.dn_to_reflectance(photo, image, use_sun_sensor=args.radiometric_calibration=="camera+sun")