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
Rimuovendo la coda si ha questo grafico
#!/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 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
m = re.search(r"pixel size\s*=\s*([\d.]+)", txt, re.IGNORECASE)
pixel_size = float(m.group(1)) if m else None
dtype_map = {1: np.uint8, 2: np.int16, 3: np.int32, 4: np.float32,
5: np.float64, 12: np.uint16, 13: np.uint32,
14: np.int64, 15: np.uint64}
np_dtype = np.dtype(dtype_map[dtype_code])
np_dtype = np_dtype.newbyteorder(">" if byte_order == 1 else "<")
return dict(samples=samples, lines=lines, bands=bands, dtype=np_dtype,
interleave=interleave, wavelengths=np.array(wavelengths),
pixel_size=pixel_size)
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()
Nessun commento:
Posta un commento