Truecolor FigSpec FS60-CL

Un piccolo script per avere un truecolor dal sensore iperspettrale FigSpex F60-CL

L'unica elaborazione e' uno stretch delle bande dei colori RGB 

 

Red band:   idx 123 -> 649.58 nm
Green band: idx 75 -> 549.73 nm
Blue band:  idx 28 -> 450.00 nm

#!/usr/bin/env python3

"""

Extract a true-color (RGB) composite from a large ENVI BIL hyperspectral cube

without loading the whole file into memory.

Usage:

python3 make_true_color.py figspec.hdr figspec.img true_color.png

If your data file has a different extension (.dat, .raw, .bin, or no

extension at all), just pass its actual path as the second argument.

"""

import sys

import os

import re

import argparse

import numpy as np

from PIL import Image

def parse_hdr(hdr_path):

txt = open(hdr_path, "r", errors="ignore").read()

def get_int(key):

m = re.search(rf"{key}\s*=\s*(\d+)", txt)

return int(m.group(1))

samples = get_int("samples")

lines = get_int("lines")

bands = get_int("bands")

dtype_code = get_int("data type")

m = re.search(r"interleave\s*=\s*(\w+)", txt, re.IGNORECASE)

interleave = m.group(1).lower()

m = re.search(r"byte order\s*=\s*(\d+)", txt)

byte_order = int(m.group(1)) if m else 0

m = re.search(r"wavelength\s*=\s*\{([^}]+)\}", txt, re.IGNORECASE)

wavelengths = [float(x) for x in m.group(1).split(",")] if m else None

# ENVI data type codes -> numpy dtype

dtype_map = {

1: np.uint8,

2: np.int16,

3: np.int32,

4: np.float32,

5: np.float64,

6: np.complex64,

9: np.complex128,

12: np.uint16,

13: np.uint32,

14: np.int64,

15: np.uint64,

}

np_dtype = np.dtype(dtype_map[dtype_code])

if byte_order == 1:

np_dtype = np_dtype.newbyteorder(">")

else:

np_dtype = np_dtype.newbyteorder("<")

return {

"samples": samples,

"lines": lines,

"bands": bands,

"dtype": np_dtype,

"interleave": interleave,

"wavelengths": wavelengths,

}

def closest_band(wavelengths, target):

return min(range(len(wavelengths)), key=lambda i: abs(wavelengths[i] - target))

def stretch(band, low_pct=2, high_pct=98):

"""Percentile contrast stretch to 0-255, robust to outliers/no-data."""

band = band.astype(np.float32)

valid = band[np.isfinite(band)]

lo, hi = np.percentile(valid, [low_pct, high_pct])

if hi <= lo:

hi = lo + 1

out = np.clip((band - lo) / (hi - lo), 0, 1)

return (out * 255).astype(np.uint8)

def main():

ap = argparse.ArgumentParser()

ap.add_argument("hdr_path")

ap.add_argument("data_path")

ap.add_argument("out_path", nargs="?", default="true_color.png")

ap.add_argument("--offset", type=int, default=None,

help="Byte offset into data file where the cube starts "

"(overrides auto-detection / header offset)")

ap.add_argument("--lines", type=int, default=None,

help="Override number of lines (in case header value is wrong)")

args = ap.parse_args()

hdr_path = args.hdr_path

data_path = args.data_path

out_path = args.out_path

meta = parse_hdr(hdr_path)

samples, lines, bands = meta["samples"], meta["lines"], meta["bands"]

dtype = meta["dtype"]

interleave = meta["interleave"]

wavelengths = meta["wavelengths"]

if args.lines is not None:

lines = args.lines

itemsize = dtype.itemsize

per_line_bytes = bands * samples * itemsize

file_size = os.path.getsize(data_path)

if args.offset is not None:

offset = args.offset

else:

expected = per_line_bytes * lines

if file_size == expected:

offset = 0

elif file_size > expected:

# extra trailing bytes; ignore them, assume no leading offset

offset = 0

else:

# file smaller than expected: leftover bytes at the START are

# likely a header the .hdr file didn't account for.

remainder = file_size % per_line_bytes

offset = remainder

lines = file_size // per_line_bytes

print(f"NOTE: file size ({file_size}) doesn't match header dims. "

f"Auto-detected a {offset}-byte leading offset and adjusted "

f"lines to {lines}. Override with --offset/--lines if wrong.")

print(f"Cube: {lines} lines x {bands} bands x {samples} samples, "

f"dtype={dtype}, interleave={interleave}, offset={offset} bytes")

if wavelengths:

r_idx = closest_band(wavelengths, 650)

g_idx = closest_band(wavelengths, 550)

b_idx = closest_band(wavelengths, 450)

print(f"Red band: idx {r_idx} -> {wavelengths[r_idx]:.2f} nm")

print(f"Green band: idx {g_idx} -> {wavelengths[g_idx]:.2f} nm")

print(f"Blue band: idx {b_idx} -> {wavelengths[b_idx]:.2f} nm")

else:

raise ValueError("No wavelength info in header; specify band indices manually.")

# Memory-map the file so we only touch the pages we actually read.

if interleave == "bil":

# file layout: [line][band][sample]

cube = np.memmap(data_path, dtype=dtype, mode="r", offset=offset,

shape=(lines, bands, samples))

red = np.array(cube[:, r_idx, :])

green = np.array(cube[:, g_idx, :])

blue = np.array(cube[:, b_idx, :])

elif interleave == "bsq":

# file layout: [band][line][sample]

cube = np.memmap(data_path, dtype=dtype, mode="r", offset=offset,

shape=(bands, lines, samples))

red = np.array(cube[r_idx, :, :])

green = np.array(cube[g_idx, :, :])

blue = np.array(cube[b_idx, :, :])

elif interleave == "bip":

# file layout: [line][sample][band]

cube = np.memmap(data_path, dtype=dtype, mode="r", offset=offset,

shape=(lines, samples, bands))

red = np.array(cube[:, :, r_idx])

green = np.array(cube[:, :, g_idx])

blue = np.array(cube[:, :, b_idx])

else:

raise ValueError(f"Unsupported interleave: {interleave}")

print("Bands extracted, stretching contrast...")

r8 = stretch(red)

g8 = stretch(green)

b8 = stretch(blue)

rgb = np.dstack([r8, g8, b8])

img = Image.fromarray(rgb, mode="RGB")

img.save(out_path)

print(f"Saved true-color image to {out_path} ({img.size[0]}x{img.size[1]})")

if __name__ == "__main__":

main()