Empirical line on image reference FigSpec

Per ovviare al problema della calibrazione visto che nel precedente post (calibrazione da manuale del sensore) ho estratto la firma in radianza di 4 punti del pannello di bianco di riferimento direttamente dall'immagine (ho avuto l'accortezza di inserirlo nella scena) usando Snap ed estrando 4 spettri in csv

 

qui i risultati sono decimente migliori...e' corretto che il pannello di bianco di riferimento sia a 0.85 perche' e' calibrato all'80%

Calibrazione di fabbrica del riferimento di bianco

 

Confrontando con gli assorbimenti atmosferici

O₂ A-band~758-770 (minimo a ~760-762nm) Ossigeno Forte

H₂O~810-840 Vapore acqueo Debole-moderata

H₂O (complesso 940nm)~890-990 (minimo tipico ~935-945, ma variabile)   

Confrontando con gli spettri misurati  

 

python3 make_reflectance_cube.py /media/disco/Arpat/FigSpec_260522_171651/figspec.hdr  /media/disco/Arpat/FigSpec_260522_171651/figspec.spe     --black  /media/disco/Arpat/FigSpec_260522_171651/figspec.figspecblack     --white-csv bianco.csv  --panel  ./Reflectance_Calibration_Cloth.figspecref     --out figspec_reflectance --from-line 33000 --to-line 40000

 

 

#!/usr/bin/env python3

"""

Computes a full calibrated reflectance cube (all bands) from an ENVI BIL

hyperspectral cube using dark-current and white-reference correction files

plus the actual spectral reflectance curve of the calibration panel, and

saves the result as a standard ENVI file (float32, BIL).

Usage:

python3 make_reflectance_cube.py figspec.hdr figspec.spe \

--black figspec.figspecblack \

--white figspec.figspecwhite \

--panel Reflectance_Calibration_Cloth.figspecref \

--out reflectance_cube

Produces:

reflectance_cube.hdr

reflectance_cube.img (float32, BIL, same samples/lines/bands as input)

Reflectance is computed per band, per spatial column as:

reflectance = (raw - dark) / (white - dark + eps) * panel_reflectance(wavelength)

Processes the cube in chunks of lines to keep memory usage low regardless

of file size.

"""

import os

import re

import argparse

import numpy as np

from scipy.signal import savgol_filter

def auto_offset_and_lines(path, bands, samples, itemsize, forced_lines=None,

forced_offset=None):

per_line_bytes = bands * samples * itemsize

file_size = os.path.getsize(path)

if forced_offset is not None and forced_lines is not None:

return forced_offset, forced_lines

expected_lines = forced_lines if forced_lines is not None else file_size // per_line_bytes

expected = per_line_bytes * expected_lines

if file_size == expected:

return 0, expected_lines

remainder = file_size % per_line_bytes

lines = file_size // per_line_bytes

print(f" [{os.path.basename(path)}] size={file_size} -> "

f"auto offset={remainder} bytes, lines={lines}")

return remainder, lines

def parse_correction_file(path, n_bands):

"""Parse a .figspecblack/.figspecwhite correction file -> (bands, spatial_bins)."""

data = open(path, "rb").read()

text = data.decode("latin1")

matches = list(re.finditer(r"[0-9a-fA-F]{200,}", text))

if not matches:

raise ValueError(f"No hex-encoded correction block found in {path}")

m = matches[0]

hexstr = data[m.start():m.end()].decode("ascii")

raw = bytes.fromhex(hexstr)

arr = np.frombuffer(raw, dtype="<f4")

if arr.size % n_bands != 0:

raise ValueError(f"{path}: {arr.size} floats doesn't divide evenly by {n_bands} bands")

spatial_bins = arr.size // n_bands

return arr.reshape(n_bands, spatial_bins)

def parse_white_csv(path, n_bands, wavelengths, wl_tolerance=0.1):

"""Parse a CSV of raw radiance spectra extracted directly from the image at the

calibration panel's location (e.g. several 'Pin N' columns from ENVI Spectrum View).

Averages all numeric columns together and returns a (bands, 1) map (spatially

uniform - no across-track/vignetting correction, since it's just point samples).

"""

with open(path, "r", errors="ignore") as f:

header = f.readline().strip().split("\t")

if len(header) == 1: # try comma if tab didn't split

header = header[0].split(",")

sep = ","

else:

sep = "\t"

rows = []

for line in f:

line = line.strip()

if not line:

continue

parts = line.split(sep)

rows.append([float(p) for p in parts])

rows = np.array(rows) # (bands, 1+n_pins) : col0 = wavelength

csv_wl = rows[:, 0]

csv_vals = rows[:, 1:] # (bands, n_pins)

if len(csv_wl) != n_bands:

raise ValueError(f"{path}: has {len(csv_wl)} rows but cube has {n_bands} bands")

max_diff = np.max(np.abs(csv_wl - wavelengths))

if max_diff > wl_tolerance:

print(f" WARNING: wavelength mismatch between {path} and cube header "

f"(max diff {max_diff:.3f} nm) - check band alignment!")

avg_spectrum = csv_vals.mean(axis=1) # (bands,)

n_pins = csv_vals.shape[1]

print(f" Parsed {path}: {n_pins} replicate spectra, averaged -> single (bands,) profile")

return avg_spectrum.reshape(n_bands, 1).astype(np.float32) # (bands, spatial_bins=1)

def parse_panel_reflectance(path):

"""Parse the calibration cloth CSV -> sorted (wavelengths_nm, reflectance_fraction) arrays."""

data = open(path, "rb").read().decode("latin1")

rows = []

for line in data.splitlines()[1:]:

line = line.strip().strip(",")

if not line:

continue

parts = [p.strip() for p in line.split(",") if p.strip() != ""]

if len(parts) < 2:

continue

try:

wl = float(parts[0])

refl_pct = float(parts[1])

rows.append((wl, refl_pct))

except ValueError:

continue

rows.sort()

wl = np.array([r[0] for r in rows])

refl = np.array([r[1] for r in rows]) / 100.0

return wl, refl

def upsample_spatial_2d(profile_2d, target_len):

"""Upsample a (bands, spatial_bins) map to (bands, target_len) via nearest-neighbor repeat."""

bins = profile_2d.shape[1]

factor = target_len / bins

if float(factor).is_integer():

return np.repeat(profile_2d, int(factor), axis=1)

x_old = np.linspace(0, 1, bins)

x_new = np.linspace(0, 1, target_len)

out = np.empty((profile_2d.shape[0], target_len), dtype=np.float32)

for b in range(profile_2d.shape[0]):

out[b] = np.interp(x_new, x_old, profile_2d[b])

return out

def write_envi_hdr(out_hdr_path, samples, lines, bands, wavelengths, wl_units,

interleave="bil", description="Calibrated reflectance cube"):

with open(out_hdr_path, "w") as f:

f.write("ENVI\n")

f.write(f"description = {{{description}}}\n")

f.write(f"samples = {samples}\n")

f.write(f"lines = {lines}\n")

f.write(f"bands = {bands}\n")

f.write("header offset = 0\n")

f.write("file type = ENVI Standard\n")

f.write("data type = 4\n") # float32

f.write(f"interleave = {interleave}\n")

f.write("byte order = 0\n")

f.write(f"wavelength units = {wl_units}\n")

wl_str = ", ".join(f"{w:.4f}" for w in wavelengths)

f.write(f"wavelength = {{{wl_str}}}\n")

f.write("reflectance scale factor = 1.0\n")

if __name__ == "__main__":

main()

 

Script radianza riflettanza per FigSpec F60-C

Seguendo i due post precedenti (1 e 2) di seguito e' riportato lo script per convertire l'immagine da DN a radianza e riflettanza usando empirical line e il pannello di bianco di riferimento (85%) in dotazione con lo strumento 

Per lanciare lo script si puo' usare 

python3 make_reflectance_cube.py figspec.hdr figspec.spe \
    --black figspec.figspecblack \
    --white figspec.figspecwhite \
    --panel Reflectance_Calibration_Cloth.figspecref \
    --out figspec_reflectance

 visto che le strisciate possono essere anche molto lunghe si puo' selezionare un intervallo di righe

 

python3 make_reflectance_cube.py figspec.hdr figspec.spe \
    --black figspec.figspecblack \
    --white figspec.figspecwhite \
    --panel Reflectance_Calibration_Cloth.figspecref \
    --out figspec_reflectance_5k_10k \
    --from-line 5000 --to-line 10000 

E qui iniziano i problemi

Se si mettono a confronto la riflettanza (prima immagine) e la radianza (seconda immagine) si vede che la riflettanza e' decisamente piu' rumorosa  ed ha valori superiori ad 1 (dovrebbero essere sempre inferiori ad 1)

 

 se poi vado ad estrarre lo spettro sul pannello di riferimento lo spettro non e' piatto (il dato e' clippato ma e' un artefatto del programma, se si vede l'immagine successiva si deve che i DN sono intorno a 3000-3100 quando il limite superiore e' 4095)

Temo che questo derivi dal metodo con cui lo strumento viene calibrato...a terra viene effettuato una immagine di dark current con l'otturatore chiuso ed una immagine di bianco tenendo in mano il drone ed inquadrando il bianco di riferimento..in pratica il bianco non e' ricavato direttamente dalla scena e possono passare anche diversi minuti tra il bianco a terra e la ripresa di tutta la scena...il che vuol dire che le condizioni di illuminazione possono essere cambiate
 

#!/usr/bin/env python3

"""

Computes a full calibrated reflectance cube (all bands) from an ENVI BIL

hyperspectral cube using dark-current and white-reference correction files

plus the actual spectral reflectance curve of the calibration panel, and

saves the result as a standard ENVI file (float32, BIL).

Usage:

python3 make_reflectance_cube.py figspec.hdr figspec.spe \

--black figspec.figspecblack \

--white figspec.figspecwhite \

--panel Reflectance_Calibration_Cloth.figspecref \

--out reflectance_cube

Produces:

reflectance_cube.hdr

reflectance_cube.img (float32, BIL, same samples/lines/bands as input)

Reflectance is computed per band, per spatial column as:

reflectance = (raw - dark) / (white - dark + eps) * panel_reflectance(wavelength)

Processes the cube in chunks of lines to keep memory usage low regardless

of file size.

"""

import os

import re

import argparse

import numpy as np

from scipy.signal import savgol_filter

def auto_offset_and_lines(path, bands, samples, itemsize, forced_lines=None,

forced_offset=None):

per_line_bytes = bands * samples * itemsize

file_size = os.path.getsize(path)

if forced_offset is not None and forced_lines is not None:

return forced_offset, forced_lines

expected_lines = forced_lines if forced_lines is not None else file_size // per_line_bytes

expected = per_line_bytes * expected_lines

if file_size == expected:

return 0, expected_lines

remainder = file_size % per_line_bytes

lines = file_size // per_line_bytes

print(f" [{os.path.basename(path)}] size={file_size} -> "

f"auto offset={remainder} bytes, lines={lines}")

return remainder, lines

def parse_correction_file(path, n_bands):

"""Parse a .figspecblack/.figspecwhite correction file -> (bands, spatial_bins)."""

data = open(path, "rb").read()

text = data.decode("latin1")

matches = list(re.finditer(r"[0-9a-fA-F]{200,}", text))

if not matches:

raise ValueError(f"No hex-encoded correction block found in {path}")

m = matches[0]

hexstr = data[m.start():m.end()].decode("ascii")

raw = bytes.fromhex(hexstr)

arr = np.frombuffer(raw, dtype="<f4")

if arr.size % n_bands != 0:

raise ValueError(f"{path}: {arr.size} floats doesn't divide evenly by {n_bands} bands")

spatial_bins = arr.size // n_bands

return arr.reshape(n_bands, spatial_bins)

def parse_panel_reflectance(path):

"""Parse the calibration cloth CSV -> sorted (wavelengths_nm, reflectance_fraction) arrays."""

data = open(path, "rb").read().decode("latin1")

rows = []

for line in data.splitlines()[1:]:

line = line.strip().strip(",")

if not line:

continue

parts = [p.strip() for p in line.split(",") if p.strip() != ""]

if len(parts) < 2:

continue

try:

wl = float(parts[0])

refl_pct = float(parts[1])

rows.append((wl, refl_pct))

except ValueError:

continue

rows.sort()

wl = np.array([r[0] for r in rows])

refl = np.array([r[1] for r in rows]) / 100.0

return wl, refl

def upsample_spatial_2d(profile_2d, target_len):

"""Upsample a (bands, spatial_bins) map to (bands, target_len) via nearest-neighbor repeat."""

bins = profile_2d.shape[1]

factor = target_len / bins

if float(factor).is_integer():

return np.repeat(profile_2d, int(factor), axis=1)

x_old = np.linspace(0, 1, bins)

x_new = np.linspace(0, 1, target_len)

out = np.empty((profile_2d.shape[0], target_len), dtype=np.float32)

for b in range(profile_2d.shape[0]):

out[b] = np.interp(x_new, x_old, profile_2d[b])

return out

def write_envi_hdr(out_hdr_path, samples, lines, bands, wavelengths, wl_units,

interleave="bil", description="Calibrated reflectance cube"):

with open(out_hdr_path, "w") as f:

f.write("ENVI\n")

f.write(f"description = {{{description}}}\n")

f.write(f"samples = {samples}\n")

f.write(f"lines = {lines}\n")

f.write(f"bands = {bands}\n")

f.write("header offset = 0\n")

f.write("file type = ENVI Standard\n")

f.write("data type = 4\n") # float32

f.write(f"interleave = {interleave}\n")

f.write("byte order = 0\n")

f.write(f"wavelength units = {wl_units}\n")

wl_str = ", ".join(f"{w:.4f}" for w in wavelengths)

f.write(f"wavelength = {{{wl_str}}}\n")

f.write("reflectance scale factor = 1.0\n")

if __name__ == "__main__":

main()

 

Stima Keystone in FigSpec F60-C

In accoppiata al precedente post qui si effettua una stima dell'errore Keystone sul sensore del FigSpec F60-CL (viene stimato l'errore spaziale al posto di quello spettrale indicato da Smile Effect)

Il grafico sottostante e' il risultato dell'elaborazione  

In sintesi sembra che da 390 a 750 nm l'errore sia accettabile e stabile intorno a 3 px.Da 750 nm in poi l'errore esplode ma siamo in zona di assorbimento atmosferico per cui il segnale e' molto basso e quindi diciamo che e' ragionevole trovare questo tipo di errore ed non e' dipendente dallo strumento ma dalla fisica dell'atmosfera

Rimuovendo la coda si ha questo grafico 

l'interpretazione potrebbe essere che il sensore (ma non ho documentazione della ditta in merito) sia internamente diviso in due sotto sensori (uno visibile ed uno NIR) con una zona di giunzione intorno a 550 nm  

 

 

#!/usr/bin/env python3

"""

Estimates keystone effect: for a fixed physical spatial edge in the scene,

tracks how its detected column (sample) position shifts across bands

(wavelengths). No keystone -> position stays constant. Keystone present ->

position drifts smoothly with wavelength.

Method:

1. Average many lines together per band to build a stable (bands, samples)

spatial profile (reduces along-track noise, keeps the across-track

structure).

2. Pick a reference band (highest-contrast by default) and find its

strongest spatial edge (max |gradient|), sub-pixel refined.

3. For every other band, search a small window around that same position

for the local edge, sub-pixel refined, so we always track the SAME

physical edge rather than jumping to a different one.

4. Plot/report edge column position vs wavelength.

Usage:

python3 detect_keystone.py figspec.hdr figspec.spe --out keystone_diagnostic.png

Works on raw DN or radiance data (recommended) - reflectance is fine too

since spatial edges are a scene/optics property, not atmosphere-dependent

like the smile-detection absorption feature.

"""

import os

import re

import argparse

import numpy as np

def auto_offset_and_lines(path, bands, samples, itemsize, forced_lines=None,

forced_offset=None):

per_line_bytes = bands * samples * itemsize

file_size = os.path.getsize(path)

if forced_offset is not None and forced_lines is not None:

return forced_offset, forced_lines

expected_lines = forced_lines if forced_lines is not None else file_size // per_line_bytes

expected = per_line_bytes * expected_lines

if file_size == expected:

return 0, expected_lines

remainder = file_size % per_line_bytes

lines = file_size // per_line_bytes

print(f" auto-detected offset={remainder} bytes, lines={lines} "

f"(file size {file_size} doesn't match header dims exactly)")

return remainder, lines

def subpixel_extreme(y, i0, find_max=True):

"""Parabolic sub-pixel refinement around discrete extreme at index i0."""

if i0 <= 0 or i0 >= len(y) - 1:

return float(i0)

y0, y1, y2 = y[i0 - 1], y[i0], y[i0 + 1]

denom = (y0 - 2 * y1 + y2)

if abs(denom) < 1e-9:

return float(i0)

delta = 0.5 * (y0 - y2) / denom

delta = np.clip(delta, -1, 1)

return i0 + delta

def main():

ap = argparse.ArgumentParser()

ap.add_argument("hdr_path")

ap.add_argument("data_path")

ap.add_argument("--out", default="keystone_diagnostic.png")

ap.add_argument("--line-stride", type=int, default=5,

help="Use every Nth line when averaging (speed/memory tradeoff)")

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

help="Reference band index to define the target edge "

"(default: auto-pick band with strongest overall edge contrast)")

ap.add_argument("--search-window-px", type=int, default=8,

help="Half-width (in samples) of the search window around the "

"reference edge position, used for every other band")

ap.add_argument("--wl-min", type=float, default=None,

help="Ignore bands below this wavelength (nm) in the fit/report")

ap.add_argument("--wl-max", type=float, default=None,

help="Ignore bands above this wavelength (nm) in the fit/report")

ap.add_argument("--split-nm", type=float, default=None,

help="If the data looks like a step (not a smooth curve), compare mean "

"positions before/after this wavelength instead of trusting the "

"quadratic fit (e.g. --split-nm 570 to test a detector splice)")

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

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

args = ap.parse_args()

meta = parse_hdr(args.hdr_path)

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

dtype = meta["dtype"]

interleave = meta["interleave"]

wavelengths = meta["wavelengths"]

pixel_size = meta["pixel_size"]

if interleave != "bil":

raise ValueError("This script currently assumes BIL interleave.")

itemsize = dtype.itemsize

offset, lines = auto_offset_and_lines(

args.data_path, bands, samples, itemsize,

forced_lines=args.lines, forced_offset=args.offset)

print(f"Using offset={offset}, lines={lines}")

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

shape=(lines, bands, samples))

print(f"Averaging every {args.line_stride}th line "

f"({lines // args.line_stride} lines used) to build a per-band spatial profile...")

idx_lines = np.arange(0, lines, args.line_stride)

accum = np.zeros((bands, samples), dtype=np.float64)

for i, li in enumerate(idx_lines):

accum += cube[li, :, :]

if i % 200 == 0:

print(f" {i}/{len(idx_lines)}", end="\r")

mean_spatial = (accum / len(idx_lines)).astype(np.float32) # (bands, samples)

print(f"\nDone averaging. Shape: {mean_spatial.shape}")

# spatial gradient per band (along samples axis)

grad = np.gradient(mean_spatial, axis=1) # (bands, samples)

if args.ref_band is None:

contrast = np.sum(np.abs(grad), axis=1)

ref_band = int(np.argmax(contrast))

print(f"Auto-selected reference band: {ref_band} "

f"({wavelengths[ref_band]:.1f} nm), contrast={contrast[ref_band]:.1f}")

else:

ref_band = args.ref_band

ref_profile_grad = np.abs(grad[ref_band])

ref_i0 = int(np.argmax(ref_profile_grad))

ref_pos = subpixel_extreme(ref_profile_grad, ref_i0, find_max=True)

print(f"Reference edge position (band {ref_band}): column {ref_pos:.2f}")

win = args.search_window_px

edge_positions = np.full(bands, np.nan)

saturated = np.zeros(bands, dtype=bool)

for b in range(bands):

lo = max(0, int(round(ref_pos)) - win)

hi = min(samples, int(round(ref_pos)) + win + 1)

local_grad = np.abs(grad[b, lo:hi])

if len(local_grad) < 3 or local_grad.max() < 1e-6:

continue

i0_local = int(np.argmax(local_grad))

if i0_local == 0 or i0_local == len(local_grad) - 1:

saturated[b] = True # hit the window boundary - unreliable

pos_local = subpixel_extreme(local_grad, i0_local, find_max=True)

edge_positions[b] = lo + pos_local

valid = ~np.isnan(edge_positions) & ~saturated

if args.wl_min is not None:

valid &= wavelengths >= args.wl_min

if args.wl_max is not None:

valid &= wavelengths <= args.wl_max

n_saturated = int(saturated.sum())

if n_saturated > 0:

print(f"\nWARNING: {n_saturated}/{bands} bands hit the search-window edge "

f"(unreliable detections, likely low SNR) - excluded from the fit.")

sat_wls = wavelengths[saturated]

if len(sat_wls) > 0:

print(f" saturated band wavelength range: {sat_wls.min():.1f} - {sat_wls.max():.1f} nm")

wl_valid = wavelengths[valid]

pos_valid = edge_positions[valid]

if len(pos_valid) < 10:

print("Too few valid bands to estimate keystone reliably.")

return

coeffs = np.polyfit(wl_valid, pos_valid, 2)

fit_curve = np.polyval(coeffs, wl_valid)

ptp_measured = pos_valid.max() - pos_valid.min()

ptp_fit = fit_curve.max() - fit_curve.min()

print(f"\nEdge column position across {bands} bands:")

print(f" raw measured range: {pos_valid.min():.2f} - {pos_valid.max():.2f} px "

f"(peak-to-valley = {ptp_measured:.2f} px)")

print(f" smooth quadratic fit range: {fit_curve.min():.2f} - {fit_curve.max():.2f} px "

f"(peak-to-valley = {ptp_fit:.2f} px)")

print(f" quadratic fit coefficients (a*wl^2+b*wl+c): {coeffs}")

if args.split_nm is not None:

seg1 = pos_valid[wl_valid < args.split_nm]

seg2 = pos_valid[wl_valid >= args.split_nm]

if len(seg1) > 3 and len(seg2) > 3:

print(f"\nTwo-segment comparison (split at {args.split_nm} nm) - use this if the "

f"data looks like a STEP rather than a smooth curve (e.g. a detector splice):")

print(f" segment 1 (< {args.split_nm}nm, n={len(seg1)}): "

f"mean={seg1.mean():.2f} px, std={seg1.std():.2f} px")

print(f" segment 2 (>= {args.split_nm}nm, n={len(seg2)}): "

f"mean={seg2.mean():.2f} px, std={seg2.std():.2f} px")

print(f" step offset between segments: {seg1.mean() - seg2.mean():.2f} px")

if pixel_size:

print(f" pixel size from header: {pixel_size} -> "

f"keystone p2v ~ {ptp_fit * pixel_size:.3f} (same units as pixel size)")

try:

import matplotlib

matplotlib.use("Agg")

import matplotlib.pyplot as plt

plt.figure(figsize=(11, 5))

excluded = ~valid & ~np.isnan(edge_positions)

if excluded.sum() > 0:

plt.scatter(wavelengths[excluded], edge_positions[excluded], s=6, alpha=0.3,

color="lightgray", label="excluded (window-edge / out of range)")

plt.scatter(wl_valid, pos_valid, s=6, alpha=0.4, label="measured (sub-pixel)", color="steelblue")

plt.plot(wl_valid, fit_curve, color="darkorange", linewidth=2,

label=f"quadratic fit (p2v={ptp_fit:.2f} px)")

plt.axhline(ref_pos, color="gray", linestyle="--", alpha=0.5,

label=f"reference position (band {ref_band})")

plt.xlabel("Wavelength (nm)")

plt.ylabel("Detected edge column (samples)")

plt.title("Keystone estimate: spatial edge position vs wavelength")

plt.legend()

plt.tight_layout()

plt.savefig(args.out, dpi=120)

print(f"\nSaved plot to {args.out}")

except ImportError:

print("\n(matplotlib not installed - numeric results above still usable)")

if __name__ == "__main__":

main()