From 27e61169774c30565f1425060d5cee316de700ef Mon Sep 17 00:00:00 2001 From: markFieldman Date: Wed, 29 Dec 2021 17:11:09 +0200 Subject: [PATCH] Added support of radiometric calibration for DJI Mavic 2 Enterprize Advanced. Added methods for adding FLIR sensor into reconstruction --- opendm/photo.py | 3 + opendm/thermal.py | 13 +- opendm/thermal_tools/__init__.py | 0 opendm/thermal_tools/dji_unpack.py | 51 +++++ opendm/thermal_tools/flir_unpack.py | 271 ++++++++++++++++++++++++++ opendm/thermal_tools/thermal_utils.py | 139 +++++++++++++ stages/run_opensfm.py | 2 +- 7 files changed, 475 insertions(+), 4 deletions(-) create mode 100644 opendm/thermal_tools/__init__.py create mode 100644 opendm/thermal_tools/dji_unpack.py create mode 100644 opendm/thermal_tools/flir_unpack.py create mode 100644 opendm/thermal_tools/thermal_utils.py diff --git a/opendm/photo.py b/opendm/photo.py index 46404e32..afbb65f7 100644 --- a/opendm/photo.py +++ b/opendm/photo.py @@ -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): diff --git a/opendm/thermal.py b/opendm/thermal.py index 5a5267d5..c64ad223 100644 --- a/opendm/thermal.py +++ b/opendm/thermal.py @@ -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 diff --git a/opendm/thermal_tools/__init__.py b/opendm/thermal_tools/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/opendm/thermal_tools/dji_unpack.py b/opendm/thermal_tools/dji_unpack.py new file mode 100644 index 00000000..8ade4c7e --- /dev/null +++ b/opendm/thermal_tools/dji_unpack.py @@ -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 \ No newline at end of file diff --git a/opendm/thermal_tools/flir_unpack.py b/opendm/thermal_tools/flir_unpack.py new file mode 100644 index 00000000..903d34d9 --- /dev/null +++ b/opendm/thermal_tools/flir_unpack.py @@ -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\xe1FLIR\x00\x01 + # \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 \ No newline at end of file diff --git a/opendm/thermal_tools/thermal_utils.py b/opendm/thermal_tools/thermal_utils.py new file mode 100644 index 00000000..6dbfdf5f --- /dev/null +++ b/opendm/thermal_tools/thermal_utils.py @@ -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) \ No newline at end of file diff --git a/stages/run_opensfm.py b/stages/run_opensfm.py index 37bdd639..b41a9c2c 100644 --- a/stages/run_opensfm.py +++ b/stages/run_opensfm.py @@ -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")