Continua l'analisi dello strumento FigSpec FS60-CL
Come si vede dal truecolor sottostante le bande laterali dell'immagine risultano decisamente meno luminose della parte centrale
Si tratta di un effetto di "vignettatura"
La cosa piu' fastidiosa e' pur soltante on RGB si ha che la distribuzione della vignettatura e' funzione della lunghezza d'onda ed e' asimmetrica anche tra il lato destro ed il lato sinistro
La vignettatura potrebbe essere corretta tramite il pannello di bianco ma questo dovrebbe essere abbastanza grande da includere tutta la larghezza del sensore
File analizzato: true_color.png
Dimensioni immagine: 480 colonne x 2845 righe
Soglia nodata (luma): 3.0
--- Profilo across-track (luma) ---
Fit lineare: luma = -0.04739 * colonna + 98.642 (R^2=0.0451)
Variazione stimata dal bordo sinistro al bordo destro: -22.70 (in unita' DN 0-255)
--- Confronto meta' sinistra vs meta' destra (split al 50%) ---
Sinistra: media=93.36 mediana=79.00 std=56.71 n_pixel=676032
Destra: media=82.23 mediana=70.67 std=55.43 n_pixel=663230
Differenza (destra - sinistra): -11.13 DN (-11.9% relativo)
Canale R: slope=-0.03882 DN/colonna, R^2=0.0387, variazione totale=-18.59 DN
Canale G: slope=-0.06856 DN/colonna, R^2=0.0748, variazione totale=-32.84 DN
Canale B: slope=-0.03480 DN/colonna, R^2=0.0244, variazione totale=-16.67 DN
Dimensioni immagine: 480 colonne x 2845 righe
Soglia nodata (luma): 3.0
--- Profilo across-track (luma) ---
Fit lineare: luma = -0.04739 * colonna + 98.642 (R^2=0.0451)
Variazione stimata dal bordo sinistro al bordo destro: -22.70 (in unita' DN 0-255)
--- Confronto meta' sinistra vs meta' destra (split al 50%) ---
Sinistra: media=93.36 mediana=79.00 std=56.71 n_pixel=676032
Destra: media=82.23 mediana=70.67 std=55.43 n_pixel=663230
Differenza (destra - sinistra): -11.13 DN (-11.9% relativo)
Canale R: slope=-0.03882 DN/colonna, R^2=0.0387, variazione totale=-18.59 DN
Canale G: slope=-0.06856 DN/colonna, R^2=0.0748, variazione totale=-32.84 DN
Canale B: slope=-0.03480 DN/colonna, R^2=0.0244, variazione totale=-16.67 DN
in modo abbastanza prevedibile il problema e' su tutte le bande con un netto gradino a 680 nm circa
script per il solo dato truecolor rgb
#!/usr/bin/env python3
"""
analyze_crosstrack_illumination.py
Analizza un'immagine true-color generata da un sensore pushbroom per
evidenziare eventuali differenze di illuminazione lungo la direzione
across-track (cioè lungo le colonne, che nel dato pushbroom corrispondono
alla direzione perpendicolare al volo).
Cosa fa:
1. Carica l'immagine (PNG/JPEG/TIFF a 3 canali).
2. Maschera i pixel di nodata (neri, dovuti a rotazione in georeferenziazione).
3. Calcola il profilo di luminosità medio/mediano per ogni colonna (RGB e luma).
4. Confronta statisticamente la metà sinistra vs la metà destra dell'immagine.
5. Salva:
- un grafico del profilo across-track (con eventuale fit lineare/polinomiale
per quantificare il gradiente)
- una mappa "flat-field corrected" (normalizzata per colonna) per vedere se
l'effetto sparisce
- un report testuale con i numeri chiave
Uso:
python3 analyze_crosstrack_illumination.py input.png [--out out_dir]
"""
import argparse
import sys
from pathlib import Path
import numpy as np
from PIL import Image
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
def load_image(path):
img = Image.open(path).convert("RGB")
arr = np.asarray(img).astype(np.float64) # (rows, cols, 3)
return arr
def nodata_mask(arr, thresh=3.0):
"""
Ritorna una maschera booleana (rows, cols) True dove il pixel e' dato valido.
I bordi neri dovuti alla rotazione in georeferenziazione hanno tipicamente
tutti i canali vicino a 0.
"""
luma = arr.mean(axis=2)
return luma > thresh
def column_profile(arr, mask):
"""
Calcola media e mediana per colonna, per ogni canale RGB e per la luma,
considerando solo i pixel validi (mask True). Colonne senza dati validi
vengono restituite come NaN.
"""
rows, cols, _ = arr.shape
mean_rgb = np.full((cols, 3), np.nan)
median_rgb = np.full((cols, 3), np.nan)
mean_luma = np.full(cols, np.nan)
valid_count = mask.sum(axis=0)
luma = arr.mean(axis=2)
for c in range(cols):
col_mask = mask[:, c]
if col_mask.sum() == 0:
continue
for ch in range(3):
vals = arr[col_mask, c, ch]
mean_rgb[c, ch] = vals.mean()
median_rgb[c, ch] = np.median(vals)
mean_luma[c] = luma[col_mask, c].mean()
return mean_rgb, median_rgb, mean_luma, valid_count
def robust_linear_fit(x, y):
"""Fit lineare semplice ignorando i NaN, ritorna (slope, intercept, r2)."""
good = ~np.isnan(y)
if good.sum() < 2:
return np.nan, np.nan, np.nan
xs, ys = x[good], y[good]
A = np.vstack([xs, np.ones_like(xs)]).T
slope, intercept = np.linalg.lstsq(A, ys, rcond=None)[0]
pred = slope * xs + intercept
ss_res = np.sum((ys - pred) ** 2)
ss_tot = np.sum((ys - ys.mean()) ** 2)
r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
return slope, intercept, r2
def analyze(path, out_dir, nodata_thresh=3.0, split_fraction=0.5):
arr = load_image(path)
rows, cols, _ = arr.shape
mask = nodata_mask(arr, nodata_thresh)
mean_rgb, median_rgb, mean_luma, valid_count = column_profile(arr, mask)
x = np.arange(cols)
# Fit lineare sulla luma per quantificare il gradiente sinistra->destra
slope, intercept, r2 = robust_linear_fit(x, mean_luma)
# Confronto meta' sinistra vs meta' destra (solo pixel validi)
split_col = int(cols * split_fraction)
left_mask = mask.copy()
left_mask[:, split_col:] = False
right_mask = mask.copy()
right_mask[:, :split_col] = False
luma = arr.mean(axis=2)
left_vals = luma[left_mask]
right_vals = luma[right_mask]
report_lines = []
report_lines.append(f"File analizzato: {path}")
report_lines.append(f"Dimensioni immagine: {cols} colonne x {rows} righe")
report_lines.append(f"Soglia nodata (luma): {nodata_thresh}")
report_lines.append("")
report_lines.append("--- Profilo across-track (luma) ---")
report_lines.append(f"Fit lineare: luma = {slope:.5f} * colonna + {intercept:.3f} (R^2={r2:.4f})")
variazione_totale = slope * (cols - 1)
report_lines.append(f"Variazione stimata dal bordo sinistro al bordo destro: {variazione_totale:+.2f} (in unita' DN 0-255)")
report_lines.append("")
report_lines.append("--- Confronto meta' sinistra vs meta' destra (split al 50%) ---")
report_lines.append(f"Sinistra: media={left_vals.mean():.2f} mediana={np.median(left_vals):.2f} std={left_vals.std():.2f} n_pixel={left_vals.size}")
report_lines.append(f"Destra: media={right_vals.mean():.2f} mediana={np.median(right_vals):.2f} std={right_vals.std():.2f} n_pixel={right_vals.size}")
diff_media = right_vals.mean() - left_vals.mean()
report_lines.append(f"Differenza (destra - sinistra): {diff_media:+.2f} DN ({100*diff_media/left_vals.mean():+.1f}% relativo)")
report_lines.append("")
for ch, name in enumerate(["R", "G", "B"]):
s, i, r2c = robust_linear_fit(x, mean_rgb[:, ch])
report_lines.append(f"Canale {name}: slope={s:.5f} DN/colonna, R^2={r2c:.4f}, variazione totale={s*(cols-1):+.2f} DN")
report_text = "\n".join(report_lines)
print(report_text)
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
(out_dir / "report.txt").write_text(report_text, encoding="utf-8")
# ---- Figura 1: profilo across-track ----
fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
ax = axes[0]
colors = {"R": "red", "G": "green", "B": "blue"}
for ch, name in enumerate(["R", "G", "B"]):
ax.plot(x, mean_rgb[:, ch], color=colors[name], alpha=0.7, label=f"media {name}")
ax.plot(x, mean_luma, color="black", linewidth=2, label="luma media")
fit_line = slope * x + intercept
ax.plot(x, fit_line, color="orange", linestyle="--", linewidth=2,
label=f"fit lineare luma (slope={slope:.4f})")
ax.axvline(split_col, color="gray", linestyle=":", label="split 50%")
ax.set_ylabel("DN medio (0-255)")
ax.set_title("Profilo di illuminazione across-track (per colonna)")
ax.legend(loc="best", fontsize=8)
ax.grid(alpha=0.3)
ax2 = axes[1]
ax2.plot(x, valid_count, color="purple")
ax2.set_ylabel("Pixel validi per colonna")
ax2.set_xlabel("Colonna (direzione across-track)")
ax2.set_title("Numero di pixel non-nodata usati per colonna")
ax2.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "crosstrack_profile.png", dpi=150)
plt.close(fig)
# ---- Figura 2: istogrammi sinistra vs destra ----
fig, ax = plt.subplots(figsize=(8, 5))
bins = np.linspace(0, 255, 60)
ax.hist(left_vals, bins=bins, alpha=0.6, label=f"Sinistra (media={left_vals.mean():.1f})", color="steelblue", density=True)
ax.hist(right_vals, bins=bins, alpha=0.6, label=f"Destra (media={right_vals.mean():.1f})", color="darkorange", density=True)
ax.axvline(left_vals.mean(), color="steelblue", linestyle="--")
ax.axvline(right_vals.mean(), color="darkorange", linestyle="--")
ax.set_xlabel("Luma (DN)")
ax.set_ylabel("Densita'")
ax.set_title("Distribuzione luminosita': meta' sinistra vs meta' destra")
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "left_vs_right_histogram.png", dpi=150)
plt.close(fig)
# ---- Figura 3: mappa originale vs flat-field corrected ----
# Normalizza ogni colonna dividendo per il proprio valore medio di luma,
# poi riporta alla media globale. Questo "appiattisce" un eventuale
# gradiente sistematico across-track (vignettatura), lasciando invece
# intatte le vere variazioni di scena (fiume, vegetazione, ecc.) che
# non dipendono dalla colonna.
global_mean = np.nanmean(mean_luma)
correction = np.where(mean_luma > 1e-6, global_mean / mean_luma, 1.0)
correction = np.clip(correction, 0.5, 2.0) # evita correzioni estreme sui bordi rumorosi
corrected = arr.copy()
for ch in range(3):
corrected[:, :, ch] = arr[:, :, ch] * correction[np.newaxis, :]
corrected = np.clip(corrected, 0, 255).astype(np.uint8)
corrected[~mask] = 0
original_uint8 = np.clip(arr, 0, 255).astype(np.uint8)
original_uint8[~mask] = 0
fig, axes = plt.subplots(1, 2, figsize=(10, 12))
axes[0].imshow(original_uint8)
axes[0].set_title("Originale")
axes[0].axis("off")
axes[1].imshow(corrected)
axes[1].set_title("Flat-field corrected\n(normalizzato per colonna)")
axes[1].axis("off")
fig.tight_layout()
fig.savefig(out_dir / "original_vs_corrected.png", dpi=150)
plt.close(fig)
# Salva anche l'immagine corretta a se' stante, utile come confronto diretto
Image.fromarray(corrected).save(out_dir / "true_color_flatfield_corrected.png")
print(f"\nOutput salvati in: {out_dir.resolve()}")
print(" - report.txt")
print(" - crosstrack_profile.png")
print(" - left_vs_right_histogram.png")
print(" - original_vs_corrected.png")
print(" - true_color_flatfield_corrected.png")
def main():
parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("input", help="Percorso dell'immagine true-color (PNG/JPEG/TIFF)")
parser.add_argument("--out", default="crosstrack_analysis", help="Cartella di output")
parser.add_argument("--nodata-thresh", type=float, default=3.0,
help="Soglia luma sotto la quale un pixel e' considerato nodata (default 3.0)")
parser.add_argument("--split-fraction", type=float, default=0.5,
help="Frazione di colonne che definisce il confine sinistra/destra (default 0.5)")
args = parser.parse_args()
analyze(args.input, args.out, args.nodata_thresh, args.split_fraction)
if __name__ == "__main__":
main()
script per il cubo iperspettrale
#!/usr/bin/env python3
"""
analyze_crosstrack_cube.py
Analizza un cubo iperspettrale pushbroom (formato ENVI: .hdr + .bil/.bsq/.bip,
oppure un array numpy .npy/.npz) per individuare e quantificare differenze di
illuminazione lungo la direzione across-track (le colonne del sensore),
banda per banda.
Perche' questo script serve
---------------------------
Nel true-color renderizzato avevi notato una fascia destra piu' scura della
sinistra. Su un solo composito RGB non si puo' distinguere se il problema
nasce:
(a) nel sensore/ottica (vignettatura radiometrica, gia' presente nella
riflettanza calibrata), oppure
(b) nel rendering RGB (stretch, white balance, ecc.), oppure
(c) da un vero effetto fisico di scena (BRDF/anisotropia, ombre, angolo
di illuminazione solare rispetto alla direzione di volo).
Lavorando sul cubo (idealmente gia' calibrato in riflettanza con dark/white/
pannello) e guardando il profilo across-track banda per banda, si puo' capire
se il gradiente e' costante su tutte le lunghezze d'onda (probabile
vignettatura/ottica del sensore, o dark/white non uniformi lungo lo swath)
oppure cambia con la banda (piu' probabile un effetto fisico/spettrale di
scena o un residuo di calibrazione spettrale).
Cosa fa lo script
-----------------
1. Carica il cubo (rows=along-track, cols=across-track/swath, bands=spettrale).
Auto-rileva l'interleave ENVI (BIL tipico per pushbroom) tramite 'spectral'.
2. Maschera i pixel nodata (es. tutti zero o NaN, o sotto una soglia).
3. Per ogni banda calcola il profilo medio/mediano per colonna (across-track).
4. Quantifica il gradiente sinistra/destra per banda con un fit lineare e
con il confronto diretto delle due meta'.
5. Produce:
- una heatmap (banda x colonna) del profilo across-track normalizzato,
per vedere a colpo d'occhio se il gradiente e' presente su tutte le bande
o solo in alcune (utile per capire la causa fisica),
- un grafico dello slope (pendenza) del gradiente in funzione della
lunghezza d'onda/banda,
- un grafico di alcune bande rappresentative (es. blu/verde/rosso/NIR se
riconoscibili, altrimenti prima/centro/ultima banda),
- un report testuale banda per banda,
- opzionalmente un cubo corretto (flat-field per colonna, per banda) in
formato .npy, utile per verificare se il problema sparisce.
Uso
---
# Cubo ENVI (serve il file .hdr, il dato binario associato viene trovato
# automaticamente in base al campo 'file' nell'header o stesso nome)
python3 analyze_crosstrack_cube.py cubo.hdr --out out_cube
# Cubo numpy gia' in memoria come array (righe, colonne, bande)
python3 analyze_crosstrack_cube.py cubo.npy --out out_cube --shape rows,cols,bands
# Salva anche il cubo corretto
python3 analyze_crosstrack_cube.py cubo.hdr --out out_cube --save-corrected
"""
import argparse
import sys
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# --------------------------------------------------------------------------
# Caricamento cubo
# --------------------------------------------------------------------------
def _parse_envi_header_datafile(hdr_path):
"""
Legge il campo 'file =' / 'data file =' dall'header ENVI, se presente,
per capire come si chiama il file binario associato. Ritorna None se non
trovato o non parsabile.
"""
try:
text = Path(hdr_path).read_text(encoding="utf-8", errors="ignore")
except Exception:
return None
for line in text.splitlines():
low = line.strip().lower()
if low.startswith("file =") or low.startswith("data file ="):
value = line.split("=", 1)[1].strip()
value = value.strip("{}").strip()
if value:
return value
return None
def _find_envi_data_file(hdr_path):
"""
Cerca il file binario associato a un header ENVI quando 'spectral' non
riesce a trovarlo automaticamente (lo fa solo per estensioni .img, .dat,
.sli o nessuna estensione). Strategie, in ordine:
1. Campo 'file =' dentro l'header stesso.
2. Qualunque file nella stessa cartella con lo stesso nome base
(stesso nome del .hdr senza estensione) e estensione diversa da
.hdr, .png, .jpg, .txt, .csv, .py (per evitare falsi positivi).
Ritorna il Path del file trovato, o None.
"""
hdr_path = Path(hdr_path)
base = hdr_path.with_suffix("") # rimuove .hdr
declared = _parse_envi_header_datafile(hdr_path)
if declared:
candidate = (hdr_path.parent / declared)
if candidate.exists():
return candidate
candidate2 = hdr_path.parent / Path(declared).name
if candidate2.exists():
return candidate2
exclude_suffixes = {".hdr", ".png", ".jpg", ".jpeg", ".txt", ".csv", ".py", ".md"}
candidates = [p for p in hdr_path.parent.glob(base.name + "*")
if p.is_file() and p.suffix.lower() not in exclude_suffixes and p != hdr_path]
if candidates:
candidates.sort(key=lambda p: p.stat().st_size, reverse=True)
return candidates[0]
return None
def load_cube(path, data_file_override=None):
"""
Carica un cubo iperspettrale da:
- file ENVI header (.hdr) -> usa la libreria 'spectral'
- file binario ENVI (qualunque estensione, es. .spe, .bil, .bsq, .bip,
.raw, .dat...) che ha un .hdr associato nella stessa cartella
- file numpy (.npy / .npz)
Ritorna array numpy (rows, cols, bands) in float64, e un dict di metadati
(wavelengths se disponibili, nomi banda, ecc.)
"""
path = Path(path)
meta = {"wavelengths": None, "band_names": None, "source": str(path)}
if path.suffix.lower() == ".npy":
cube = np.load(str(path))
return cube.astype(np.float64), meta
if path.suffix.lower() == ".npz":
data = np.load(str(path))
for key in ("cube", "data", "arr_0"):
if key in data:
return data[key].astype(np.float64), meta
first_key = list(data.keys())[0]
return data[first_key].astype(np.float64), meta
hdr_path = None
data_path = Path(data_file_override) if data_file_override else None
if path.suffix.lower() == ".hdr":
hdr_path = path
else:
candidate_hdrs = [
path.with_suffix(".hdr"),
Path(str(path) + ".hdr"),
]
for c in candidate_hdrs:
if c.exists():
hdr_path = c
if data_path is None:
data_path = path
break
if hdr_path is None or not hdr_path.exists():
raise ValueError(
f"Formato non riconosciuto o header ENVI non trovato per {path}. "
f"Usa un file .hdr (ENVI), il file binario con .hdr associato, oppure .npy/.npz"
)
import spectral
if data_path is None:
try:
img = spectral.envi.open(str(hdr_path))
except spectral.io.envi.EnviDataFileNotFoundError:
found = _find_envi_data_file(hdr_path)
if found is None:
raise ValueError(
f"Impossibile trovare il file dati binario associato a {hdr_path}. "
f"Passa esplicitamente il file dati con --data-file, oppure rinomina/"
f"copia il binario accanto all'header con estensione .img/.dat/.bin."
)
print(f"[info] File dati non auto-rilevato da 'spectral'; uso: {found}")
img = spectral.envi.open(str(hdr_path), image=str(found))
else:
img = spectral.envi.open(str(hdr_path), image=str(data_path))
cube = img.load().astype(np.float64)
wl = getattr(img, "bands", None)
if wl is not None and getattr(wl, "centers", None) is not None:
meta["wavelengths"] = np.array(wl.centers, dtype=float)
# conserviamo l'header originale (dict) e l'interleave, utili per
# riesportare il cubo corretto in formato ENVI mantenendo wavelength,
# unita', sensor type, ecc.
meta["envi_header"] = dict(getattr(img, "metadata", {}) or {})
meta["interleave"] = getattr(img, "interleave", None)
return np.asarray(cube), meta
def ensure_rows_cols_bands(cube, shape_hint=None):
"""
Se il cubo non e' gia' (rows, cols, bands), e l'utente ha fornito
--shape rows,cols,bands, lo reshape. Altrimenti assume che sia gia'
nell'ordine corretto (che e' il default per 'spectral' e per la
convenzione usata nello script true-color della sessione precedente).
"""
if shape_hint is not None:
rows, cols, bands = shape_hint
cube = cube.reshape(rows, cols, bands)
return cube
# --------------------------------------------------------------------------
# Maschera nodata
# --------------------------------------------------------------------------
def nodata_mask(cube, thresh=1e-6):
"""
True dove il pixel e' valido. Un pixel e' nodata se e' NaN su tutte le
bande oppure se il valore assoluto medio su tutte le bande e' sotto
soglia (tipico bordo nero da rotazione/georeferenziazione).
"""
finite = np.isfinite(cube).all(axis=2)
mean_abs = np.nanmean(np.abs(cube), axis=2)
valid = finite & (mean_abs > thresh)
return valid
# --------------------------------------------------------------------------
# Profilo across-track per banda
# --------------------------------------------------------------------------
def crosstrack_profile_per_band(cube, mask):
"""
Ritorna:
mean_profile: (cols, bands) media per colonna e banda (NaN se nessun
pixel valido in quella colonna)
valid_count: (cols,) numero di pixel validi per colonna
"""
rows, cols, bands = cube.shape
mean_profile = np.full((cols, bands), np.nan)
valid_count = mask.sum(axis=0)
for c in range(cols):
col_mask = mask[:, c]
if col_mask.sum() == 0:
continue
mean_profile[c, :] = cube[col_mask, c, :].mean(axis=0)
return mean_profile, valid_count
def linear_fit(x, y):
good = np.isfinite(y)
if good.sum() < 2:
return np.nan, np.nan, np.nan
xs, ys = x[good], y[good]
A = np.vstack([xs, np.ones_like(xs)]).T
slope, intercept = np.linalg.lstsq(A, ys, rcond=None)[0]
pred = slope * xs + intercept
ss_res = np.sum((ys - pred) ** 2)
ss_tot = np.sum((ys - ys.mean()) ** 2)
r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
return slope, intercept, r2
def export_envi_cube(out_path, cube, meta, interleave="bil"):
"""
Salva un cubo (rows, cols, bands) in formato ENVI (coppia .hdr + file dati),
riusando i metadati originali (wavelength, unita', sensor type, ecc.) se
disponibili in meta['envi_header'], cosi' il cubo esportato resta
apribile in ENVI/QGIS/altri tool con le stesse informazioni spettrali
dell'originale.
`out_path` puo' essere passato con o senza estensione: verra' sempre
scritto come coppia <out_path>.hdr + <out_path>.dat (interleave BIL di
default, tipico per dati pushbroom).
"""
import spectral
out_path = Path(out_path)
if out_path.suffix.lower() == ".hdr":
out_path = out_path.with_suffix("")
orig_header = dict(meta.get("envi_header") or {})
# Ripulisce campi che 'spectral' ricalcola da solo o che diventerebbero
# inconsistenti (dimensioni, nome file, offset) per evitare header corrotti
for key in ("samples", "lines", "bands", "header offset", "file type",
"byte order", "data type", "interleave"):
orig_header.pop(key, None)
metadata = orig_header # contiene ad es. wavelength, wavelength units, sensor type, map info...
spectral.envi.save_image(
str(out_path.with_suffix(".hdr")),
cube.astype(np.float32),
dtype=np.float32,
force=True,
interleave=interleave,
metadata=metadata,
)
return out_path.with_suffix(".hdr")
# --------------------------------------------------------------------------
# Analisi principale
# --------------------------------------------------------------------------
def analyze(cube_path, out_dir, shape_hint=None, nodata_thresh=1e-6,
split_fraction=0.5, save_corrected=False, data_file_override=None,
output_format="envi"):
cube, meta = load_cube(cube_path, data_file_override=data_file_override)
cube = ensure_rows_cols_bands(cube, shape_hint)
rows, cols, bands = cube.shape
print(f"Cubo caricato: {rows} righe x {cols} colonne x {bands} bande")
mask = nodata_mask(cube, nodata_thresh)
mean_profile, valid_count = crosstrack_profile_per_band(cube, mask) # (cols, bands)
x = np.arange(cols)
wavelengths = meta.get("wavelengths")
if wavelengths is not None and len(wavelengths) == bands:
band_axis = wavelengths
band_label = "Lunghezza d'onda (nm)"
else:
band_axis = np.arange(bands)
band_label = "Indice banda"
# --- Fit lineare per banda: slope, intercept, r2 ---
slopes = np.full(bands, np.nan)
r2s = np.full(bands, np.nan)
for b in range(bands):
s, i, r2 = linear_fit(x, mean_profile[:, b])
slopes[b] = s
r2s[b] = r2
# variazione percentuale bordo-bordo per banda, rispetto al valore medio globale della banda
global_mean_per_band = np.nanmean(mean_profile, axis=0)
variation_total = slopes * (cols - 1)
variation_pct = 100 * variation_total / np.where(global_mean_per_band != 0, global_mean_per_band, np.nan)
# --- Confronto meta' sinistra vs destra per banda ---
split_col = int(cols * split_fraction)
left_mask = mask.copy()
left_mask[:, split_col:] = False
right_mask = mask.copy()
right_mask[:, :split_col] = False
left_mean = np.array([cube[left_mask, b].mean() if left_mask.any() else np.nan for b in range(bands)])
right_mean = np.array([cube[right_mask, b].mean() if right_mask.any() else np.nan for b in range(bands)])
diff_pct = 100 * (right_mean - left_mean) / np.where(left_mean != 0, left_mean, np.nan)
# --- Report testuale ---
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
lines = []
lines.append(f"Cubo: {cube_path}")
lines.append(f"Dimensioni: {rows} righe x {cols} colonne x {bands} bande")
lines.append(f"Soglia nodata: {nodata_thresh}")
lines.append("")
lines.append("Riepilogo globale (mediato su tutte le bande):")
lines.append(f" Slope medio: {np.nanmean(slopes):.6f} per colonna")
lines.append(f" Variazione media bordo-bordo: {np.nanmean(variation_total):+.4f} ({np.nanmean(variation_pct):+.2f}% rispetto alla media di banda)")
lines.append(f" Differenza media destra-sinistra: {np.nanmean(diff_pct):+.2f}%")
lines.append("")
n_bands_show = min(bands, 30)
lines.append(f"Dettaglio per banda (prime {n_bands_show} mostrate nel report; tutte le bande nel file .csv):")
lines.append(f"{'banda':>6} {'asse':>10} {'slope':>12} {'R^2':>8} {'var_tot':>10} {'var_%':>8} {'diff_dx-sx_%':>14}")
for b in range(n_bands_show):
lines.append(f"{b:6d} {band_axis[b]:10.2f} {slopes[b]:12.6f} {r2s[b]:8.3f} {variation_total[b]:10.3f} {variation_pct[b]:8.2f} {diff_pct[b]:14.2f}")
if bands > n_bands_show:
lines.append(f"... (+{bands - n_bands_show} bande, vedi crosstrack_per_band_stats.csv)")
report_text = "\n".join(lines)
print(report_text)
(out_dir / "report.txt").write_text(report_text, encoding="utf-8")
# CSV completo per tutte le bande
csv_path = out_dir / "crosstrack_per_band_stats.csv"
with open(csv_path, "w", encoding="utf-8") as f:
f.write("band_index,band_axis,slope,r2,variation_total,variation_pct,diff_right_minus_left_pct,global_mean\n")
for b in range(bands):
f.write(f"{b},{band_axis[b]},{slopes[b]},{r2s[b]},{variation_total[b]},{variation_pct[b]},{diff_pct[b]},{global_mean_per_band[b]}\n")
print(f"CSV per-banda salvato in {csv_path}")
# ---------------------------------------------------------------
# Figura 1: heatmap banda x colonna, normalizzata per banda
# (ogni riga della heatmap divisa per la propria media, cosi' il
# confronto tra bande con intensita' molto diverse resta leggibile)
# ---------------------------------------------------------------
norm_profile = mean_profile / np.where(global_mean_per_band != 0, global_mean_per_band, np.nan)[np.newaxis, :]
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(
norm_profile.T,
aspect="auto",
origin="lower",
extent=[0, cols, band_axis[0], band_axis[-1]],
cmap="RdBu_r",
vmin=0.7, vmax=1.3,
)
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel(band_label)
ax.set_title("Profilo across-track normalizzato per banda\n(rosso=piu' luminoso della media di banda, blu=piu' scuro)")
fig.colorbar(im, ax=ax, label="valore / media di banda")
fig.tight_layout()
fig.savefig(out_dir / "heatmap_band_vs_column.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Figura 2: slope (gradiente) in funzione della banda
# ---------------------------------------------------------------
fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
axes[0].plot(band_axis, variation_pct, color="darkred")
axes[0].axhline(0, color="gray", linestyle=":")
axes[0].set_ylabel("Variazione bordo-bordo (%)")
axes[0].set_title("Gradiente across-track in funzione della banda")
axes[0].grid(alpha=0.3)
axes[1].plot(band_axis, diff_pct, color="darkblue")
axes[1].axhline(0, color="gray", linestyle=":")
axes[1].set_ylabel("Diff. destra-sinistra (%)")
axes[1].set_xlabel(band_label)
axes[1].grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "gradient_vs_wavelength.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Figura 3: profili across-track per alcune bande rappresentative
# ---------------------------------------------------------------
sample_idx = sorted(set([0, bands // 4, bands // 2, (3 * bands) // 4, bands - 1]))
fig, ax = plt.subplots(figsize=(9, 5))
cmap = plt.get_cmap("viridis")
for k, b in enumerate(sample_idx):
ax.plot(x, mean_profile[:, b], color=cmap(k / max(1, len(sample_idx) - 1)),
label=f"banda {b} ({band_axis[b]:.0f}{'nm' if band_label.startswith('Lunghezza') else ''})")
ax.axvline(split_col, color="gray", linestyle=":", label="split 50%")
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel("Valore medio (DN o riflettanza)")
ax.set_title("Profilo across-track per bande rappresentative")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "sample_bands_profile.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Correzione flat-field per colonna, per banda (opzionale)
# ---------------------------------------------------------------
if save_corrected:
print("Applico correzione flat-field per colonna (per banda)...")
correction = np.where(mean_profile > 1e-9, global_mean_per_band[np.newaxis, :] / mean_profile, 1.0)
correction = np.clip(correction, 0.3, 3.0) # evita amplificazioni estreme ai bordi rumorosi
corrected = cube * correction[np.newaxis, :, :]
corrected[~mask] = 0
if output_format == "envi":
corrected_path = export_envi_cube(out_dir / "cube_flatfield_corrected", corrected, meta,
interleave="bil")
print(f"Cubo corretto salvato in formato ENVI: {corrected_path} (+ file dati associato)")
else:
corrected_path = out_dir / "cube_flatfield_corrected.npy"
np.save(corrected_path, corrected)
print(f"Cubo corretto salvato in formato numpy: {corrected_path}")
# Verifica: ricalcola il profilo sul cubo corretto e stampa il nuovo slope medio
mean_profile_corr, _ = crosstrack_profile_per_band(corrected, mask)
slopes_corr = np.full(bands, np.nan)
for b in range(bands):
s, i, r2 = linear_fit(x, mean_profile_corr[:, b])
slopes_corr[b] = s
print(f"Slope medio PRIMA della correzione: {np.nanmean(slopes):.6f}")
print(f"Slope medio DOPO la correzione: {np.nanmean(slopes_corr):.6f}")
print(f"\nOutput salvati in: {out_dir.resolve()}")
print(" - report.txt")
print(" - crosstrack_per_band_stats.csv")
print(" - heatmap_band_vs_column.png")
print(" - gradient_vs_wavelength.png")
print(" - sample_bands_profile.png")
if save_corrected:
if output_format == "envi":
print(" - cube_flatfield_corrected.hdr (+ file dati binario associato)")
else:
print(" - cube_flatfield_corrected.npy")
def main():
parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("input", help="Cubo iperspettrale: file .hdr (ENVI), .npy o .npz")
parser.add_argument("--out", default="crosstrack_cube_analysis", help="Cartella di output")
parser.add_argument("--shape", default=None,
help="Solo se il .npy non e' gia' (rows,cols,bands): 'rows,cols,bands'")
parser.add_argument("--nodata-thresh", type=float, default=1e-6,
help="Soglia sotto la quale un pixel e' nodata (default 1e-6)")
parser.add_argument("--split-fraction", type=float, default=0.5,
help="Frazione di colonne che definisce il confine sinistra/destra")
parser.add_argument("--save-corrected", action="store_true",
help="Salva anche il cubo corretto (flat-field per colonna, per banda) come .npy")
parser.add_argument("--data-file", default=None,
help="Percorso esplicito del file binario ENVI associato all'header, "
"da usare se l'auto-rilevamento fallisce (es. estensioni non standard come .spe)")
parser.add_argument("--output-format", choices=["envi", "npy"], default="envi",
help="Formato del cubo corretto salvato con --save-corrected (default: envi)")
args = parser.parse_args()
shape_hint = None
if args.shape:
shape_hint = tuple(int(v) for v in args.shape.split(","))
if len(shape_hint) != 3:
print("ERRORE: --shape deve essere 'rows,cols,bands'", file=sys.stderr)
sys.exit(1)
analyze(args.input, args.out, shape_hint, args.nodata_thresh,
args.split_fraction, args.save_corrected, data_file_override=args.data_file,
output_format=args.output_format)
if __name__ == "__main__":
main()
