decode_hir/decode_hir.py
2025-03-13 21:24:44 +01:00

204 lines
5.3 KiB
Python
Executable file

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
# Decoder for the .HIR thermal image format of the Seek Shot camera.
# Copyright (C) 2024 Thomas Kolb
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import sys
import struct
import numpy as np
import matplotlib.pyplot as pp
from PIL import Image
def map_degree_celsius(ir_raw):
#MIN_CELSIUS = 24
#MAX_CELSIUS = 50
#MIN_RAW = 4101
#MAX_RAW = 5864
#return MIN_CELSIUS + (ir_raw - MIN_RAW) * (MAX_CELSIUS - MIN_CELSIUS) / (MAX_RAW - MIN_RAW)
return ir_raw / 64 - 40
def ycbcr2rgb(y, cb, cr):
KB = 0.0722
KR = 0.2126
KG = 1 - KB - KR
ys = (np.astype(y, np.float32) - 16) / (256 - 32)
pb = (np.astype(cb, np.float32) - 128) / (256 - 32)
pr = (np.astype(cr, np.float32) - 128) / (256 - 32)
r = ys + (2 - 2*KR) * pr
g = ys - KB/KG * (2 - 2*KB) * pb - KR/KG * (2 - 2*KR) * pr
b = ys + (2 - 2*KB) * pb
offset = 16/256
scale = 256/224
r = (r - offset) * scale
g = (g - offset) * scale
b = (b - offset) * scale
r[r > 1.0] = 1.0
g[g > 1.0] = 1.0
b[b > 1.0] = 1.0
r[r < 0] = 0
g[g < 0] = 0
b[b < 0] = 0
gamma = 1.0
r = r**gamma
g = g**gamma
b = b**gamma
return 255*r, 255*g, 255*b
hirdata = open(sys.argv[1], 'rb').read()
print("Finding IR data offset...")
readpos = 2
while True:
block_type, block_len = struct.unpack(">HH", hirdata[readpos:readpos+4])
print(f"Block type: {block_type:04X}, block length: {block_len:6d} / 0x{block_len:04X} byte")
if block_type == 0xFFDA:
break
else:
readpos += block_len + 2
print(f"Found JPEG image data at offset 0x{readpos:08X}")
print("Scanning for JPEG end marker...")
word = 0x0000
while word != 0xFFD9:
while hirdata[readpos] != 0xFF:
readpos += 1
word, = struct.unpack(">H", hirdata[readpos:readpos+2])
readpos += 2
print(f"End of JPEG data is at 0x{readpos:08X}")
METADATA_OFFSET = (readpos & 0xFFFFF000) + 0x1000
IR_DATA_OFFSET = METADATA_OFFSET + 0x1000
IR_SIZE_X = 206
IR_SIZE_Y = 156
IR_NPIXELS = IR_SIZE_X * IR_SIZE_Y
VIS_BW_DATA_OFFSET = IR_DATA_OFFSET + 0x10000
VIS_BW_SIZE_X = 1440
VIS_BW_SIZE_Y = 1080
VIS_BW_NPIXELS = VIS_BW_SIZE_X * VIS_BW_SIZE_Y
VIS_CBCR_DATA_OFFSET = VIS_BW_DATA_OFFSET + VIS_BW_NPIXELS
VIS_CBCR_SIZE_X = 1440
VIS_CBCR_SIZE_X_SINGLE = VIS_CBCR_SIZE_X // 2
VIS_CBCR_SIZE_Y = 540
VIS_CBCR_NPIXELS = VIS_CBCR_SIZE_X * VIS_CBCR_SIZE_Y
print(f"Metadata is at 0x{METADATA_OFFSET:08X}")
print(f"Thermal data is at 0x{IR_DATA_OFFSET:08X}")
print(f"Visual black/white data is at 0x{VIS_BW_DATA_OFFSET:08X}")
print(f"Visual color (cbcr) data is at 0x{VIS_CBCR_DATA_OFFSET:08X}")
# process IR data
irdata_raw = np.frombuffer(hirdata[IR_DATA_OFFSET : IR_DATA_OFFSET + IR_NPIXELS * 2], dtype = np.int16)
irdata = map_degree_celsius(np.astype(irdata_raw, np.float32))
irdata_2d = np.reshape(irdata, (IR_SIZE_Y, IR_SIZE_X))
# process black/white image (Y channel)
vis_bw_data_raw = np.frombuffer(hirdata[VIS_BW_DATA_OFFSET : VIS_BW_DATA_OFFSET + VIS_BW_NPIXELS], dtype = np.uint8)
vis_bw_data_2d = np.reshape(vis_bw_data_raw, (VIS_BW_SIZE_Y, VIS_BW_SIZE_X))
# deinterleave subsampled color images (U/V channels)
vis_uv_data_raw = np.frombuffer(hirdata[VIS_CBCR_DATA_OFFSET : VIS_CBCR_DATA_OFFSET + VIS_CBCR_NPIXELS], dtype = np.uint8)
vis_cb_2d = np.reshape(vis_uv_data_raw[0::2], (VIS_CBCR_SIZE_Y, VIS_CBCR_SIZE_X_SINGLE))
vis_cr_2d = np.reshape(vis_uv_data_raw[1::2], (VIS_CBCR_SIZE_Y, VIS_CBCR_SIZE_X_SINGLE))
# combine images into RGB
# scale U and V up to match the black/white size
vis_cb_2d_scaled = np.asarray(Image.fromarray(vis_cb_2d, mode='L').resize( (VIS_BW_SIZE_X, VIS_BW_SIZE_Y) ))
vis_cr_2d_scaled = np.asarray(Image.fromarray(vis_cr_2d, mode='L').resize( (VIS_BW_SIZE_X, VIS_BW_SIZE_Y) ))
y = vis_bw_data_2d
cb = vis_cb_2d_scaled
cr = vis_cr_2d_scaled
print(y[0, 0], cb[0, 0], cr[0, 0])
r, g, b = ycbcr2rgb(y, cb, cr)
r = np.astype(r, np.uint8)
g = np.astype(g, np.uint8)
b = np.astype(b, np.uint8)
#r[r > 255] = 255
#g[g > 255] = 255
#b[b > 255] = 255
rgb = np.empty((VIS_BW_SIZE_Y, VIS_BW_SIZE_X, 3), dtype=np.uint8)
rgb[:,:,0] = r
rgb[:,:,1] = g
rgb[:,:,2] = b
print(rgb.shape)
vis_im = Image.fromarray(rgb, mode="RGB")
vis_im.save('/tmp/test.png')
Image.fromarray(y, "L").save('/tmp/y.png')
Image.fromarray(cb, "L").save('/tmp/cb.png')
Image.fromarray(cr, "L").save('/tmp/cr.png')
Image.fromarray(r, "L").save('/tmp/r.png')
Image.fromarray(g, "L").save('/tmp/g.png')
Image.fromarray(b, "L").save('/tmp/b.png')
print(vis_im.getpixel((0, 0)))
print(vis_im.getpixel((1, 0)))
# plotting
pp.figure()
pp.imshow(irdata_2d, cmap="magma")
pp.colorbar()
#pp.figure()
#pp.imshow(y)
#pp.colorbar()
#
#fig, axes = pp.subplots(1, 2)
#axes[0].imshow(cb)
#axes[1].imshow(cr)
pp.figure()
pp.imshow(rgb)
pp.show()