Vignettatura FigSpec FS60-C (2)

TLDR: volare sempre orientati con l'azimuth del Sole 

l'azimuth del Sole cambia con il giorno e l'ora.All'alba ed al tramonto e' molto vicino a Est ed Ovest ma a mezzogiorno e' orientato verso Sud per la latitudini settentrionali 

Frugando sulla chiavetta USB del software del sensore ho trovato un file denominato FS6261009镜头系数 che tradotto suona come "coefficiente dell'obbiettivo"....aprendolo si vede che si tratta di un file json in cui sono inseriti due array (uno monodimensionale di 480 valori ed uno bidimensionale 72x480)

 

il grafico assomiglia ad una correzione di vignettatura ma e' solo per una banda (sembra la 121) 

Ho provato ad usare questa informazione per correggere la vignettatura ma nonostante la correzione l'effetto al bordo non e' stato eliminato ..l'ipotesi successiva e' che non sia qualcosa di relativo al sensore ma alla geometria di ripresa

Dai dati GPS si ha 

Punti usati: 166
Lunghezza traccia: 134.1 m
Azimuth linea di volo: 251.86 deg (reciproco: 71.86 deg)
Direzione across-track (perpendicolare al volo): 341.86 deg / 161.86 deg
 

Al giorno e l'ora di ripresa il Sole aveva una elevazione di 66 gradi ed un azimuth, quindi la traiettoria di volo era praticamente parallela al direzione del Sole

Il sensore della camera ha una FOV di 24.75 gradi (IFOV 0.9 mrad × 480 canali spaziali) 

  

Il fenomeno ai bordi quindi puo' essere non solo la vignettatura ma anche il forward scattering ed il backscattering che generano il fenomeno di hotspot BRDF

In pratica la direzione di ripresa e' tale che in una direzione l'ombra dell'oggetto e' coperta dall'oggetto stesso (back scattering, massima riflettanza hotspot propriamente detto, guardando dalla parte opposta dal Sole) mentre nella direzione opposta si osserva l'ombra piuttosto che l'oggetto stesso (forward scattering, guardando verso il Sole)

Il fenomeno affligge tutti i sensori iperspettrali sia da satellite che airborne ed e' particolarmente importante quanto piu' grande e' l'angolo di field of view. Un satellite che ha tipicamente questo problema e' Modis con una FOV di 55 gradi con una correzione basata sul modello a kernel RossThick-LiSparse. Lo stesso fenomeno si puo' verificare in analisi di immagini multitemporali anche di satelliti a bassa FOV (tipo Landsat 15 gradi) pur essendo eliosincroni

Altrimenti per una correzione piu' semplice viene utilizzato il modello RPV — Rahman-Pinty-Verstraete  

 

 #!/usr/bin/env python3
"""
brdf_hotspot_validation.py

Conferma quantitativa (non solo qualitativa) dell'ipotesi che l'asimmetria
sinistra/destra osservata nel profilo across-track del true-color sia
spiegabile con l'effetto hotspot/anti-hotspot BRDF, usando il modello
semi-empirico RPV (Rahman-Pinty-Verstraete) — lo stesso tipo di modello
usato in letteratura (es. MISR, POLDER) per descrivere la riflettanza
bidirezionale di superfici naturali includendo l'hotspot.

Idea generale
-------------
1. Dalla traccia GPS (o da un azimuth di volo fornito manualmente) si
   ricava la direzione across-track (perpendicolare al volo).
2. Nota la posizione del sole (azimuth, elevazione) e il FOV del sensore,
   si calcola per ogni colonna dello swath la geometria di vista (zenith e
   azimuth relativo al sole) e si predice con RPV il fattore BRDF atteso
   (normalizzato, cioe' solo la FORMA della variazione, non il valore
   assoluto di riflettanza).
3. Dal true-color si estrae il profilo across-track empirico (luma media
   per colonna), e se disponibile un file di calibrazione vignettatura
   (vendor, vedi analyze_crosstrack_illumination.py --vignette-json) lo si
   usa per "scorporare" la componente ottica (vignettatura) da quella di
   scena (BRDF), isolando il residuo asimmetrico da confrontare col modello.
4. Si calcola la correlazione tra la forma prevista da RPV e il residuo
   osservato: un'alta correlazione conferma quantitativamente che
   l'asimmetria e' spiegabile con l'hotspot/anti-hotspot, non con un
   artefatto strumentale.

Il modello RPV (riassunto)
---------------------------
R(theta_i, theta_r, phi) = rho0 * M(theta_i,theta_r) * F(g,Theta) * H(G)

  M(theta_i,theta_r) = (cos(theta_i) cos(theta_r) (cos(theta_i)+cos(theta_r)))^(k-1)
      -> termine di forma "bowl/dome"; k=1 lo rende neutro (default)

  F(g,Theta) = (1-Theta^2) / (1 + 2*Theta*cos(g) + Theta^2)^1.5
      -> funzione di fase di Henyey-Greenstein; g = angolo di scattering;
         Theta<0 favorisce il backscattering (tipico canopy vegetale)

  H(G) = 1 + (1-h) / (1 + G/h)
      -> funzione hotspot; G = sqrt(tan^2(theta_i)+tan^2(theta_r)
         -2*tan(theta_i)*tan(theta_r)*cos(phi)) e' la "distanza" in spazio
         tangente tra le direzioni sole e sensore (G=0 esattamente
         all'hotspot); h controlla quanto e' "stretto" il picco (tipico
         0.05-0.3 per canopy vegetali)

rho0 viene posto a 1 perche' qui interessa solo la FORMA relativa, non il
valore assoluto di riflettanza.

Uso
---
    python3 brdf_hotspot_validation.py true_color.png \\
        --gps-csv figspec.gps \\
        --sun-azimuth 154 --sun-elevation 66 \\
        --vignette-json test.json \\
        --out out_brdf

    # oppure con azimuth di volo manuale invece del GPS:
    python3 brdf_hotspot_validation.py true_color.png \\
        --flight-azimuth 251.86 --sun-azimuth 154 --sun-elevation 66 --out out_brdf
"""

import argparse
import json
import sys
from pathlib import Path

import numpy as np
from PIL import Image
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt


# ---------------------------------------------------------------------
# Geometria di volo (azimuth) dai punti GPS
# ---------------------------------------------------------------------

def flight_azimuth_from_gps(csv_path, lat_col="Latitude", lon_col="Longitude", delimiter=","):
    import pandas as pd
    df = pd.read_csv(csv_path, delimiter=delimiter)
    lat = df[lat_col].values.astype(float)
    lon = df[lon_col].values.astype(float)

    lat0 = np.radians(lat.mean())
    R = 6371000.0
    x = np.radians(lon - lon.mean()) * R * np.cos(lat0)
    y = np.radians(lat - lat.mean()) * R
    pts = np.column_stack([x, y])
    pts_centered = pts - pts.mean(axis=0)

    cov = np.cov(pts_centered.T)
    eigvals, eigvecs = np.linalg.eigh(cov)
    main_dir = eigvecs[:, np.argmax(eigvals)]
    az = np.degrees(np.arctan2(main_dir[0], main_dir[1])) % 360
    return az


# ---------------------------------------------------------------------
# Modello RPV
# ---------------------------------------------------------------------

def rpv_reflectance(theta_i_deg, theta_r_deg, phi_deg, k=1.0, theta_hg=-0.15, h=0.1, rho0=1.0):
    """
    Calcola il fattore di riflettanza bidirezionale RPV (non normalizzato
    in valore assoluto, rho0 e' un fattore di scala arbitrario).

    theta_i_deg : zenith solare (scalare, gradi)
    theta_r_deg : zenith di vista (array, gradi)
    phi_deg     : azimuth relativo vista-sole (array, gradi) = az_vista - az_sole
    k           : parametro di forma bowl/dome (1.0 = neutro)
    theta_hg    : parametro di asimmetria della fase di Henyey-Greenstein (-1..1)
    h           : parametro di nitidezza dell'hotspot (piu' piccolo = picco piu' stretto)
    rho0        : fattore di scala (irrilevante se si normalizza il risultato)
    """
    ti = np.radians(theta_i_deg)
    tr = np.radians(theta_r_deg)
    phi = np.radians(phi_deg)

    cos_ti, cos_tr = np.cos(ti), np.cos(tr)
    sin_ti, sin_tr = np.sin(ti), np.sin(tr)

    # angolo di scattering (fase)
    cos_g = cos_ti * cos_tr + sin_ti * sin_tr * np.cos(phi)
    cos_g = np.clip(cos_g, -1.0, 1.0)

    # funzione di fase Henyey-Greenstein
    F = (1 - theta_hg**2) / (1 + 2 * theta_hg * cos_g + theta_hg**2) ** 1.5

    # funzione hotspot (spazio tangente)
    tan_ti, tan_tr = np.tan(ti), np.tan(tr)
    G = np.sqrt(np.clip(tan_ti**2 + tan_tr**2 - 2 * tan_ti * tan_tr * np.cos(phi), 0, None))
    H = 1 + (1 - h) / (1 + G / h)

    # termine di forma bowl/dome
    M = (cos_ti * cos_tr * (cos_ti + cos_tr)) ** (k - 1)

    R = rho0 * M * F * H
    return R, np.degrees(np.arccos(cos_g)), G


# ---------------------------------------------------------------------
# Profilo across-track empirico dal true-color
# ---------------------------------------------------------------------

def empirical_column_profile(image_path, nodata_thresh=3.0):
    img = Image.open(image_path).convert("RGB")
    arr = np.asarray(img).astype(np.float64)
    luma = arr.mean(axis=2)
    mask = luma > nodata_thresh
    cols = arr.shape[1]

    profile = np.full(cols, np.nan)
    for c in range(cols):
        if mask[:, c].any():
            profile[c] = luma[mask[:, c], c].mean()
    return profile


def load_vignette_gain(vignette_json, band_key, cols):
    with open(vignette_json, "r", encoding="utf-8") as f:
        data = json.load(f)
    coeffs = {e["Key"]: np.array(e["Value"], dtype=float) for e in data.get("Coefficients", [])}
    if not coeffs:
        raise ValueError(f"Nessun coefficiente trovato in {vignette_json}")
    key = band_key if band_key is not None and band_key in coeffs else list(coeffs.keys())[0]
    gain = coeffs[key]
    if len(gain) != cols:
        x_orig = np.linspace(0, 1, len(gain))
        x_new = np.linspace(0, 1, cols)
        gain = np.interp(x_new, x_orig, gain)
    return gain, key


def fit_symmetric_vignetting(profile):
    """
    Fallback se non e' disponibile un file di calibrazione vendor: stima
    la componente di vignettatura come un fit quadratico del profilo (bella
    approssimazione per una caduta di illuminazione a "campana" simmetrica),
    e ritorna il "gain" equivalente (1/vignetting_shape_normalizzata) cosi'
    l'interfaccia e' la stessa di load_vignette_gain.
    """
    cols = len(profile)
    x = np.arange(cols)
    valid = np.isfinite(profile)
    coeffs = np.polyfit(x[valid], profile[valid], deg=2)
    fitted = np.polyval(coeffs, x)
    fitted_norm = fitted / np.nanmean(fitted)
    gain = 1.0 / fitted_norm
    return gain, None


# ---------------------------------------------------------------------
# Analisi principale
# ---------------------------------------------------------------------

def analyze(image_path, out_dir, sun_azimuth, sun_elevation, fov_deg,
            flight_azimuth=None, gps_csv=None, vignette_json=None, band_key=None,
            rpv_k=1.0, rpv_theta=-0.15, rpv_h=0.1, nodata_thresh=3.0):

    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)

    # --- 1. azimuth di volo ---
    if flight_azimuth is None:
        if gps_csv is None:
            raise ValueError("Serve --flight-azimuth oppure --gps-csv per determinare la geometria di volo")
        flight_azimuth = flight_azimuth_from_gps(gps_csv)
        print(f"Azimuth di volo stimato dai punti GPS: {flight_azimuth:.2f} deg")
    else:
        print(f"Azimuth di volo (fornito manualmente): {flight_azimuth:.2f} deg")

    across_1 = (flight_azimuth + 90) % 360
    across_2 = (flight_azimuth - 90) % 360
    print(f"Direzioni across-track: {across_1:.2f} deg / {across_2:.2f} deg")

    theta_i = 90 - sun_elevation
    print(f"Sole: azimuth={sun_azimuth:.1f} deg, elevazione={sun_elevation:.1f} deg (zenith={theta_i:.1f} deg)")

    def fold180(a):
        a = a % 360
        return min(a, 360 - a)

    d1 = fold180(across_1 - sun_azimuth)
    d2 = fold180(across_2 - sun_azimuth)
    print(f"Scarto tra direzioni across-track e azimuth solare: {d1:.1f} deg / {d2:.1f} deg")

    # --- 2. profilo empirico dal true-color ---
    profile = empirical_column_profile(image_path, nodata_thresh)
    cols = len(profile)
    x = np.arange(cols)
    center = (cols - 1) / 2.0

    # --- 3. geometria di vista per colonna ---
    half_fov = fov_deg / 2.0
    view_zenith = np.abs((x - center) / center) * half_fov  # 0 al centro, half_fov ai bordi

    # due possibili assegnazioni lato<->azimuth (non sappiamo a priori quale
    # lato del sensore corrisponde a quale bearing fisico, dipende dal
    # montaggio/scansione del sensore: testiamo entrambe e scegliamo quella
    # che correla meglio con l'osservato, riportando comunque i risultati
    # di entrambe per trasparenza)
    offset = x - center
    view_az_A = np.where(offset < 0, across_1, across_2)
    view_az_B = np.where(offset < 0, across_2, across_1)

    results = {}
    for label, view_az in [("A", view_az_A), ("B", view_az_B)]:
        phi = view_az - sun_azimuth
        R, phase_angle, G = rpv_reflectance(theta_i, view_zenith, phi, k=rpv_k, theta_hg=rpv_theta, h=rpv_h)
        results[label] = {
            "view_az": view_az,
            "R": R,
            "phase_angle": phase_angle,
            "G": G,
        }

    # --- 4. vignettatura (per isolare il residuo BRDF dall'empirico) ---
    if vignette_json is not None:
        gain, used_key = load_vignette_gain(vignette_json, band_key, cols)
        vign_source = f"file vendor ({vignette_json}, Key={used_key})"
    else:
        gain, used_key = fit_symmetric_vignetting(profile)
        vign_source = "fit quadratico simmetrico (fallback, nessun file vendor fornito)"
    vignetting_shape = 1.0 / gain
    vignetting_shape_norm = vignetting_shape / np.nanmean(vignetting_shape)
    print(f"Componente di vignettatura stimata da: {vign_source}")

    empirical_norm = profile / np.nanmean(profile)
    residual_empirical = empirical_norm / vignetting_shape_norm
    residual_empirical_norm = residual_empirical / np.nanmean(residual_empirical)

    # --- 5. confronto modello vs osservato, per entrambe le assegnazioni ---
    def compare(model_R):
        model_norm = model_R / np.nanmean(model_R)
        valid = np.isfinite(residual_empirical_norm) & np.isfinite(model_norm)
        r = np.corrcoef(residual_empirical_norm[valid], model_norm[valid])[0, 1]
        rmse = np.sqrt(np.nanmean((residual_empirical_norm[valid] - model_norm[valid]) ** 2))
        # combinato: vignettatura * BRDF, confrontato con il profilo grezzo (non scorporato)
        combined = model_norm * vignetting_shape_norm
        combined_norm = combined / np.nanmean(combined)
        r_combined = np.corrcoef(empirical_norm[valid], combined_norm[valid])[0, 1]
        return model_norm, r, rmse, combined_norm, r_combined

    lines = []
    lines.append(f"Immagine: {image_path}")
    lines.append(f"Colonne (samples): {cols}")
    lines.append(f"FOV assunto: {fov_deg:.2f} deg (default: da datasheet FS-60, IFOV 0.9 mrad x 480 samples)")
    lines.append(f"Azimuth di volo: {flight_azimuth:.2f} deg")
    lines.append(f"Direzioni across-track: {across_1:.2f} / {across_2:.2f} deg")
    lines.append(f"Sole: azimuth={sun_azimuth:.1f} deg, elevazione={sun_elevation:.1f} deg (zenith={theta_i:.1f} deg)")
    lines.append(f"Scarto tra direzioni across-track e azimuth solare: {d1:.1f} deg / {d2:.1f} deg")
    lines.append(f"Parametri RPV: k={rpv_k}, Theta={rpv_theta}, h={rpv_h}")
    lines.append(f"Vignettatura stimata da: {vign_source}")
    lines.append("")

    best_label, best_r = None, -2
    for label in ("A", "B"):
        view_az = results[label]["view_az"]
        model_norm, r, rmse, combined_norm, r_combined = compare(results[label]["R"])
        results[label]["model_norm"] = model_norm
        results[label]["r_residual"] = r
        results[label]["rmse_residual"] = rmse
        results[label]["combined_norm"] = combined_norm
        results[label]["r_combined"] = r_combined

        side_desc = f"colonne basse -> {view_az[0]:.1f} deg, colonne alte -> {view_az[-1]:.1f} deg"
        lines.append(f"--- Assegnazione {label} ({side_desc}) ---")
        lines.append(f"  Correlazione modello BRDF vs residuo empirico (vignettatura scorporata): r={r:.3f} (RMSE={rmse:.4f})")
        lines.append(f"  Correlazione modello combinato (vignettatura x BRDF) vs profilo grezzo: r={r_combined:.3f}")
        lines.append(f"  Range angolo di fase coperto dal FOV: {results[label]['phase_angle'].min():.1f} - {results[label]['phase_angle'].max():.1f} deg")
        lines.append(f"  Range fattore hotspot G (0=hotspot esatto): {results[label]['G'].min():.3f} - {results[label]['G'].max():.3f}")
        lines.append("")

        if r > best_r:
            best_r = r
            best_label = label

    lines.append(f"Assegnazione con correlazione migliore: {best_label} (r={best_r:.3f})")
    hotspot_zenith_in_fov = theta_i <= half_fov
    lines.append(f"Zenith solare ({theta_i:.1f} deg) {'RIENTRA' if hotspot_zenith_in_fov else 'NON rientra'} "
                 f"nel range di zenith di vista coperto dal FOV (0-{half_fov:.1f} deg): "
                 f"{'ci si aspetta un picco hotspot pieno vicino al bordo' if hotspot_zenith_in_fov else 'ci si aspetta solo un aumento progressivo verso il bordo, senza raggiungere il picco pieno'}")

    report_text = "\n".join(lines)
    print("\n" + report_text)
    (out_dir / "report.txt").write_text(report_text, encoding="utf-8")

    # --- Figura 1: confronto residuo empirico vs modello BRDF (entrambe le assegnazioni) ---
    fig, ax = plt.subplots(figsize=(10, 5.5))
    ax.plot(x, residual_empirical_norm, color="black", linewidth=1.8, label="residuo empirico (vignettatura scorporata)")
    ax.plot(x, results["A"]["model_norm"], color="crimson", linestyle="--",
            label=f"modello RPV - assegnazione A (r={results['A']['r_residual']:.3f})")
    ax.plot(x, results["B"]["model_norm"], color="steelblue", linestyle="--",
            label=f"modello RPV - assegnazione B (r={results['B']['r_residual']:.3f})")
    ax.set_xlabel("Colonna (direzione across-track)")
    ax.set_ylabel("Valore normalizzato (media=1)")
    ax.set_title("Residuo empirico (dopo rimozione vignettatura) vs modello BRDF (RPV)")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig(out_dir / "residual_vs_brdf_model.png", dpi=150)
    plt.close(fig)

    # --- Figura 2: profilo grezzo vs modello combinato (vignettatura x BRDF) ---
    fig, ax = plt.subplots(figsize=(10, 5.5))
    ax.plot(x, empirical_norm, color="black", linewidth=1.8, label="profilo empirico grezzo (normalizzato)")
    ax.plot(x, vignetting_shape_norm, color="gray", linestyle=":", label="sola componente vignettatura")
    ax.plot(x, results[best_label]["combined_norm"], color="darkorange", linewidth=2,
            label=f"modello combinato vignettatura x BRDF (assegn. {best_label}, r={results[best_label]['r_combined']:.3f})")
    ax.set_xlabel("Colonna (direzione across-track)")
    ax.set_ylabel("Valore normalizzato (media=1)")
    ax.set_title("Profilo grezzo vs modello combinato (vignettatura ottica x effetto BRDF hotspot)")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig(out_dir / "raw_vs_combined_model.png", dpi=150)
    plt.close(fig)

    print(f"\nOutput salvati in: {out_dir.resolve()}")
    print(" - report.txt")
    print(" - residual_vs_brdf_model.png")
    print(" - raw_vs_combined_model.png")


def main():
    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument("image", help="Immagine true-color (PNG/JPEG/TIFF)")
    parser.add_argument("--out", default="out_brdf", help="Cartella di output")
    parser.add_argument("--gps-csv", default=None, help="File CSV GPS per stimare l'azimuth di volo automaticamente")
    parser.add_argument("--flight-azimuth", type=float, default=None, help="Azimuth di volo manuale (gradi), alternativa a --gps-csv")
    parser.add_argument("--sun-azimuth", type=float, required=True, help="Azimuth solare (gradi da nord)")
    parser.add_argument("--sun-elevation", type=float, required=True, help="Elevazione solare sopra l'orizzonte (gradi)")
    parser.add_argument("--fov", type=float, default=24.75,
                         help="FOV across-track totale del sensore (gradi). Default 24.75 deg, "
                              "stimato da datasheet FS-60 (IFOV 0.9 mrad x 480 samples)")
    parser.add_argument("--vignette-json", default=None, help="File JSON di calibrazione vignettatura vendor (opzionale)")
    parser.add_argument("--band-key", type=int, default=None, help="Key della banda da usare dal file vignette-json")
    parser.add_argument("--rpv-k", type=float, default=1.0, help="Parametro RPV k (forma bowl/dome, default 1.0=neutro)")
    parser.add_argument("--rpv-theta", type=float, default=-0.15, help="Parametro RPV Theta (asimmetria fase, default -0.15)")
    parser.add_argument("--rpv-h", type=float, default=0.1, help="Parametro RPV h (nitidezza hotspot, default 0.1)")
    parser.add_argument("--nodata-thresh", type=float, default=3.0, help="Soglia luma per nodata (default 3.0)")
    args = parser.parse_args()

    analyze(args.image, args.out, args.sun_azimuth, args.sun_elevation, args.fov,
            flight_azimuth=args.flight_azimuth, gps_csv=args.gps_csv,
            vignette_json=args.vignette_json, band_key=args.band_key,
            rpv_k=args.rpv_k, rpv_theta=args.rpv_theta, rpv_h=args.rpv_h,
            nodata_thresh=args.nodata_thresh)


if __name__ == "__main__":
    main()

 

 

 

Vignettatura FigSpec FS60-C

Aggiornamento

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 

 

 

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()