martedì 6 ottobre 2026

Icosaedro Aruco Tags

Questo test prende un icosaedro sul quale sono stati apposti degli aruco tags (in modo da avere almeno 3 facce sempre visibili alla camera) e tramite le normali calcola la posizione del centro dell'icosaedro tramite opencv

Un limite di questo metodo nel costruirselo da soli e' che e' difficile centrare il centro dell'aruco tag sul centro della faccia...e questo limita molto la precisione del posizionamento del centro dell'icosaedro dato che e' calcolata dall'intersezione delle normali agli aruco tags 

 


 

 

"""
Calibrazione della camera RGB della RealSense D415 con scacchiera.

Requisiti:
pip install pyrealsense2 opencv-contrib-python numpy

Uso:
1. Stampa checkerboard_A4.pdf A SCALA 100% (non "adatta alla pagina").
2. Misura col righello il quadratino di verifica stampato sul foglio:
se non è esattamente SQUARE_SIZE_M metri di lato, correggi la costante
qui sotto con la misura reale.
3. Esegui lo script, inquadra la scacchiera da tante angolazioni/distanze/
inclinazioni diverse (coprendo tutta l'inquadratura: centro, angoli,
vicino, lontano, inclinata). Premi SPAZIO per catturare un frame quando
gli angoli vengono rilevati (disegnati a colori), ESC/'q' per terminare
e calibrare.
4. Servono almeno 15-20 catture buone, distribuite in tutta l'immagine,
per una calibrazione affidabile.

Il risultato (camera_matrix, dist_coeffs) viene salvato in
d415_calibration.npz, pronto per essere caricato nello script ArUco.
"""

import cv2
import numpy as np
import pyrealsense2 as rs

# ------------------------- CONFIGURAZIONE -------------------------

# Angoli interni della scacchiera (colonne, righe) = (num_quadrati_x - 1, num_quadrati_y - 1)
PATTERN_SIZE = (9, 6)

# Lato di un quadrato in METRI: usa il valore misurato col righello sul foglio stampato
SQUARE_SIZE_M = 0.020

FRAME_W, FRAME_H, FPS = 1280, 720, 30

OUT_NPZ = "d415_calibration.npz"

MIN_CAPTURES = 15


def build_object_points(pattern_size, square_size):
cols, rows = pattern_size
objp = np.zeros((rows * cols, 3), np.float32)
objp[:, :2] = np.mgrid[0:cols, 0:rows].T.reshape(-1, 2) * square_size
return objp


def main():
objp_template = build_object_points(PATTERN_SIZE, SQUARE_SIZE_M)

pipeline = rs.pipeline()
config = rs.config()
config.enable_stream(rs.stream.color, FRAME_W, FRAME_H, rs.format.bgr8, FPS)
profile = pipeline.start(config)

criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

obj_points = [] # punti 3D reali (nel piano della scacchiera)
img_points = [] # punti 2D corrispondenti nell'immagine
img_size = None

print("SPAZIO = cattura frame corrente, 'q'/ESC = termina e calibra")

try:
while True:
frames = pipeline.wait_for_frames()
color_frame = frames.get_color_frame()
if not color_frame:
continue

frame = np.asanyarray(color_frame.get_data())
gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
img_size = gray.shape[::-1] # (w, h)

found, corners = cv2.findChessboardCorners(
gray, PATTERN_SIZE,
flags=cv2.CALIB_CB_ADAPTIVE_THRESH + cv2.CALIB_CB_NORMALIZE_IMAGE
)

display = frame.copy()
if found:
corners_refined = cv2.cornerSubPix(
gray, corners, (11, 11), (-1, -1), criteria
)
cv2.drawChessboardCorners(display, PATTERN_SIZE, corners_refined, found)
else:
corners_refined = None

cv2.putText(display, f"Catture: {len(obj_points)}/{MIN_CAPTURES}+",
(10, 30), cv2.FONT_HERSHEY_SIMPLEX, 0.7,
(0, 255, 0) if found else (0, 0, 255), 2)

cv2.imshow("Calibrazione D415 - SPAZIO=cattura, q=fine", display)
key = cv2.waitKey(1) & 0xFF

if key == ord(' ') and found:
obj_points.append(objp_template.copy())
img_points.append(corners_refined)
print(f"Cattura #{len(obj_points)} registrata")

elif key in (ord('q'), 27):
break

finally:
pipeline.stop()
cv2.destroyAllWindows()

if len(obj_points) < 4:
print("Troppo poche catture per calibrare (minimo consigliato 15-20).")
return

print(f"Calibrazione in corso con {len(obj_points)} catture...")
ret, camera_matrix, dist_coeffs, rvecs, tvecs = cv2.calibrateCamera(
obj_points, img_points, img_size, None, None
)

# Errore di riproiezione medio, indicatore della qualità della calibrazione.
# Uso numpy invece di cv2.norm: su alcune build di OpenCV (es. 5.0) cv2.norm
# può lamentare un mismatch di tipo (CV_32FC1 vs CV_32FC2) anche se le shape
# sono corrette, quindi si aggira il problema calcolando la norma a mano.
total_error = 0
for i in range(len(obj_points)):
proj, _ = cv2.projectPoints(obj_points[i], rvecs[i], tvecs[i], camera_matrix, dist_coeffs)
proj = np.asarray(proj, dtype=np.float64).reshape(-1, 2)
observed = np.asarray(img_points[i], dtype=np.float64).reshape(-1, 2)
error = np.linalg.norm(observed - proj, axis=1).mean()
total_error += error
mean_error = total_error / len(obj_points)

print("RMS calibrazione:", ret)
print("Errore di riproiezione medio (px, minore è meglio, <0.5 è ottimo):", mean_error)
print("Camera matrix:\n", camera_matrix)
print("Distortion coeffs:\n", dist_coeffs.ravel())

np.savez(OUT_NPZ, camera_matrix=camera_matrix, dist_coeffs=dist_coeffs,
image_size=img_size, mean_reproj_error=mean_error)
print(f"Salvato in {OUT_NPZ}")


if __name__ == "__main__":
main()


 

"""
Stima della posizione 3D del centro di un icosaedro con marker ArUco sulle facce,
tramite intersezione ai minimi quadrati delle normali dei marker rilevati
in uno stream video dalla RealSense D415 (canale RGB), via pyrealsense2 -
stesso approccio usato in calibrate_d415.py, che negozia la risoluzione con
la camera in modo affidabile (a differenza di cv2.VideoCapture generico via
UVC/V4L2, che spesso resta bloccato a risoluzioni basse). Il punto stimato
viene mediato sugli ultimi 10 frame per ridurre il jitter.

Richiede: pyrealsense2, opencv-contrib-python (per il modulo cv2.aruco), numpy
pip install pyrealsense2 opencv-contrib-python numpy
"""

from collections import deque

import cv2
import numpy as np
import pyrealsense2 as rs

# ------------------------- CONFIGURAZIONE -------------------------

# Calibrazione camera: SOSTITUISCI con i valori reali della tua camera
# (ottenuti con calibrate_d415.py). Senza calibrazione corretta la stima
# del centro sarà sbagliata.
CAMERA_MATRIX = np.array([
[931.62977526, 0.0, 648.9939992],
[ 0.0, 935.13433783, 366.09135937],
[ 0.0, 0.0, 1.0]
], dtype=np.float64)

DIST_COEFFS = np.zeros((5, 1), dtype=np.float64)

# Risoluzione/fps del canale RGB richiesti alla D415. Deve combaciare con la
# risoluzione usata in calibrate_d415.py (FRAME_W/FRAME_H), altrimenti
# fx/fy/cx/cy non sono più validi.
FRAME_W, FRAME_H, FPS = 1280, 720, 30

MARKER_LENGTH = 0.026 # lato del marker in metri: misuralo con un calibro
ARUCO_DICT = cv2.aruco.getPredefinedDictionary(cv2.aruco.DICT_4X4_50)
DETECTOR_PARAMS = cv2.aruco.DetectorParameters()
DETECTOR = cv2.aruco.ArucoDetector(ARUCO_DICT, DETECTOR_PARAMS)

# Se true, inverte il segno della normale (dipende da come è orientato
# il marker rispetto alla faccia: prova entrambi i valori e guarda quale
# fa convergere il punto rosso dentro il solido)
FLIP_NORMAL = True

# Numero di stime consecutive su cui calcolare la media mobile del centro
SMOOTHING_WINDOW = 10

WINDOW_NAME = "Icosaedro - stima centro 3D"


def marker_object_points(length):
"""Punti 3D degli angoli del marker nel suo sistema locale, piano z=0,
ordine: alto-sx, alto-dx, basso-dx, basso-sx (ordine standard ArUco)."""
h = length / 2.0
return np.array([
[-h, h, 0],
[ h, h, 0],
[ h, -h, 0],
[-h, -h, 0]
], dtype=np.float64)


OBJ_POINTS = marker_object_points(MARKER_LENGTH)


def estimate_marker_pose(corners, camera_matrix):
"""Pose di un singolo marker via solvePnP (sostituisce l'estimatePoseSingleMarkers
ormai deprecato). Ritorna rvec, tvec."""
ok, rvec, tvec = cv2.solvePnP(
OBJ_POINTS, corners.reshape(-1, 2), camera_matrix, DIST_COEFFS,
flags=cv2.SOLVEPNP_IPPE_SQUARE
)
return rvec, tvec


def closest_point_to_lines(points, directions):
"""
Trova il punto 3D che minimizza la somma delle distanze quadratiche da un
insieme di rette skew, ciascuna definita da (punto, direzione unitaria).
Risolve A x = b con A = sum(I - d d^T), b = sum((I - d d^T) p).
"""
A = np.zeros((3, 3))
b = np.zeros(3)
I = np.eye(3)
for p, d in zip(points, directions):
d = d / np.linalg.norm(d)
M = I - np.outer(d, d)
A += M
b += M @ p
x, *_ = np.linalg.lstsq(A, b, rcond=None)
return x


def open_realsense_color_stream():
"""Apre il pipeline RealSense sul solo canale colore, come in
calibrate_d415.py. Ritorna il pipeline avviato."""
pipeline = rs.pipeline()
config = rs.config()
config.enable_stream(rs.stream.color, FRAME_W, FRAME_H, rs.format.bgr8, FPS)
profile = pipeline.start(config)

color_profile = rs.video_stream_profile(profile.get_stream(rs.stream.color))
intr = color_profile.get_intrinsics()
print(f"Stream RGB D415 avviato: {intr.width}x{intr.height}@{FPS}fps")
if (intr.width, intr.height) != (FRAME_W, FRAME_H):
print("ATTENZIONE: la risoluzione effettiva riportata dalla camera "
f"({intr.width}x{intr.height}) non combacia con quella richiesta "
f"({FRAME_W}x{FRAME_H}). Verifica che FRAME_W/FRAME_H siano supportati "
"da questo sensore (realsense-viewer per l'elenco modalità).")

return pipeline


def main():
pipeline = open_realsense_color_stream()

# Finestra ridimensionabile dall'utente (trascinando i bordi)
cv2.namedWindow(WINDOW_NAME, cv2.WINDOW_NORMAL)

sign = -1.0 if FLIP_NORMAL else 1.0

# Buffer circolare con gli ultimi N centri stimati, per la media mobile
center_history = deque(maxlen=SMOOTHING_WINDOW)

try:
while True:
frames = pipeline.wait_for_frames()
color_frame = frames.get_color_frame()
if not color_frame:
continue

frame = np.asanyarray(color_frame.get_data())
gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
corners, ids, _ = DETECTOR.detectMarkers(gray)

line_points, line_dirs = [], []

if ids is not None:
cv2.aruco.drawDetectedMarkers(frame, corners, ids)

for c in corners:
rvec, tvec = estimate_marker_pose(c, CAMERA_MATRIX)
R, _ = cv2.Rodrigues(rvec)

# Normale della faccia in coordinate camera (asse Z locale del marker)
normal = R @ np.array([0.0, 0.0, 1.0])
center = tvec.reshape(3)

# Il centro dell'icosaedro sta lungo la normale, dalla parte
# opposta a dove punta il marker (verso "dentro" al solido)
line_points.append(center)
line_dirs.append(sign * normal)

cv2.drawFrameAxes(frame, CAMERA_MATRIX, DIST_COEFFS, rvec, tvec, MARKER_LENGTH * 0.5)

if len(line_points) >= 2:
center3d = closest_point_to_lines(line_points, line_dirs)
center_history.append(center3d)

# Media mobile sugli ultimi SMOOTHING_WINDOW valori disponibili
center_avg = np.mean(center_history, axis=0)

img_pt, _ = cv2.projectPoints(
center_avg.reshape(1, 3), np.zeros(3), np.zeros(3),
CAMERA_MATRIX, DIST_COEFFS
)
x, y = img_pt.ravel().astype(int)
cv2.circle(frame, (x, y), 6, (0, 0, 255), -1)
cv2.putText(
frame,
f"centro (media {len(center_history)}/{SMOOTHING_WINDOW}): "
f"{center_avg[0]:.3f}, {center_avg[1]:.3f}, {center_avg[2]:.3f} m",
(10, 30), cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 0, 255), 2
)
else:
cv2.putText(frame, "servono >=2 marker visibili", (10, 30),
cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 165, 255), 2)

cv2.imshow(WINDOW_NAME, frame)
if cv2.waitKey(1) & 0xFF == ord('q'):
break

finally:
pipeline.stop()
cv2.destroyAllWindows()


if __name__ == "__main__":
main()

 

 

 

martedì 29 settembre 2026

Nvidia Dev Env con Docker e Visual Code

Mi sono comprato un vecchio portatile MSI con una NVidia GeForce GTX 860M che supporta i driver 580.178.04 Cuda 13

La scheda video e' molto vecchia per cui la soluzione migliore per fare lo sviluppo e' utilizzare un docker container con Cuda SDK 12 gia' pronto 


 

Si installa Visual Code con estensione Dev Containers e C/C++ 

 Da notare che il container non ha installato gdb... per questo motivo in postCreateCommand lo installa tramite apt

 in ./devcontainer/devcontainer.json

{
"name": "CUDA Development Environment",
"image": "nvidia/cuda:12.0.0-devel-ubuntu22.04",
"customizations": {
"vscode": {
"extensions": [
"ms-vscode.cpptools",
"ms-vscode.cmake-tools"
]
}
},
"containerEnv": {
"NVCC_FLAGS": "-arch=sm_50"
},
"runArgs": [
"--gpus=all"
],
"postCreateCommand": "apt-get update && apt-get install -y gdb && nvcc --version"
}


 in .vscode/launch.json

{
"version": "0.2.0",
"configurations": [
{
"name": "Debug CUDA (GDB)",
"type": "cppdbg",
"request": "launch",
"program": "${fileDirname}/${fileBasenameNoExtension}",
"args": [],
"stopAtEntry": false,
"cwd": "${fileDirname}",
"environment": [],
"externalConsole": false,
"MIMode": "gdb",
"miDebuggerPath": "/usr/bin/gdb",
"setupCommands": [
{
"description": "Abilita pretty-printing per gdb",
"text": "enable pretty-printing",
"ignoreFailures": true
}
],
"preLaunchTask": "Compila CUDA con nvcc"
}
]
}

 

in .vscode/tasks.json

si deve impostare il target sm_50 perche' la scheda e' vecchia ed ha cuda capabilities 5 

{
"version": "2.0.0",
"tasks": [
{
"type": "shell",
"label": "Compila CUDA con nvcc",
"command": "nvcc",
"args": [
"-g",
"-arch=sm_50",
"${file}",
"-o",
"${fileDirname}/${fileBasenameNoExtension}"
],
"group": {
"kind": "build",
"isDefault": true
},
"problemMatcher": [
"$gcc"
],
"detail": "Compilatore CUDA NVCC per GTX 860M"
}
]
}

 

test.cu

#include <iostream>

__global__ void helloFromGPU() {
printf("Hello from GPU! Thread index: %d\n", threadIdx.x);
}

int main() {
std::cout << "Hello from CPU!" << std::endl;
helloFromGPU<<<1, 5>>>();
cudaDeviceSynchronize();
return 0;
}

 per iniziare lo sviluppo si deve aprire CTRL+SHIFT+P Dev Containers: Reopen in container

A questo punto in basso a sinistra si ha un box azzurro

ed il terminale punta alla shell del container root@02cb6c6a1ae0:/workspaces/progetto_cuda#

 

Per fare il debug del kernel Cuda si usa compute-sanitizer

compute-sanitizer ./test_cuda
========= COMPUTE-SANITIZER
Hello from CPU!
Hello from GPU! Thread index: 0
Hello from GPU! Thread index: 1
Hello from GPU! Thread index: 2
Hello from GPU! Thread index: 3
Hello from GPU! Thread index: 4
========= ERROR SUMMARY: 0 errors 

 

 

 

 

 

Cronaca di un fallimento: crisotilo su Enmap (2)

Prosegue la serie degli schiaffi

Mappa random forest contro verita' a terra


Ho provato a cambiare approccio passando random forest.

Per fare questa cosa pero' per avere  un classificatore random forest oltre ad avere un dataset di positivi devo avere anche un dataset negativo

Per fare questo e' stato creato un tema puntuale di punti random non coincidenti con quelli del censimento amianto

 

import numpy as np
import geopandas as gpd
import shapely
from shapely.geometry import box
from scipy.spatial import cKDTree

# ---------------- parametri ----------------
SRC = "amianto_PRA_pubblico.shp"   # punti amianto
MASK_RASTER = None                 # es. "mask_built.tif" (1 = edificato valido, stessa griglia EnMAP)
AREA = None                        # es. "area_campionamento.shp" (poligoni); ignorato se MASK_RASTER è impostato
MIN_DIST = 50                     # metri minimi dai punti amianto
N_RATIO = 4                        # negativi per ogni positivo
SEED = 0
OUT = "background_random.shp"
CRS_DEFAULT = 32632
# --------------------------------------------

rng = np.random.default_rng(SEED)

# punti amianto
pts = gpd.read_file(SRC)
pts = pts[pts.geometry.notna() & ~pts.geometry.is_empty]
if pts.crs is None:
    pts = pts.set_crs(CRS_DEFAULT)
elif pts.crs.is_geographic:
    pts = pts.to_crs(CRS_DEFAULT)

xy = np.column_stack([pts.geometry.x, pts.geometry.y])
n = len(xy)
n_bg = N_RATIO * n
tree = cKDTree(xy)

if MASK_RASTER:
    # ---- campionamento su pixel della maschera (centri pixel, senza duplicati) ----
    import rasterio
    with rasterio.open(MASK_RASTER) as src:
        if src.crs != pts.crs:
            raise ValueError(f"CRS diverso: raster {src.crs} vs punti {pts.crs}. Riproietta uno dei due.")
        mask = src.read(1) == 1
        T = src.transform
    rows, cols = np.nonzero(mask)
    pick = rng.permutation(len(rows))
    cx, cy = rasterio.transform.xy(T, rows[pick], cols[pick], offset="center")
    cand = np.column_stack([cx, cy])
    d, _ = tree.query(cand, k=1)
    cand = cand[d > MIN_DIST]
    if len(cand) < n_bg:
        print(f"ATTENZIONE: disponibili solo {len(cand)} punti su {n_bg} richiesti. "
              f"Riduci MIN_DIST o N_RATIO, oppure allarga la maschera.")
    out = cand[:n_bg]

else:
    # ---- campionamento casuale continuo dentro un poligono / bounding box ----
    if AREA:
        poly = gpd.read_file(AREA).to_crs(pts.crs).union_all()
    else:
        poly = box(*pts.total_bounds)
    minx, miny, maxx, maxy = poly.bounds

    out = np.empty((0, 2))
    for _ in range(200):                      # limite di sicurezza contro cicli infiniti
        if len(out) >= n_bg:
            break
        m = max(2 * (n_bg - len(out)), 10_000)
        cand = np.column_stack([rng.uniform(minx, maxx, m),
                                rng.uniform(miny, maxy, m)])
        cand = cand[shapely.contains_xy(poly, cand[:, 0], cand[:, 1])]
        d, _ = tree.query(cand, k=1)
        out = np.vstack([out, cand[d > MIN_DIST]])
    if len(out) < n_bg:
        print(f"ATTENZIONE: generati solo {len(out)} punti su {n_bg} richiesti. "
              f"Riduci MIN_DIST o allarga l'area.")
    out = out[:n_bg]

bg = gpd.GeoDataFrame(
    {"id": np.arange(len(out)), "classe": 0},
    geometry=gpd.points_from_xy(out[:, 0], out[:, 1]),
    crs=pts.crs,
)
bg.to_file(OUT)
print(f"Positivi: {n} | Negativi salvati: {len(bg)} | file: {OUT}")

fatto questo sono stati estratti gli spettri dei rispettivi punti e creato un campo boolean false in caso di assenza di amianto. La libreria spettrale e' stata fusa con quella dei tetti in amianto ...quindi ho ottenuto una libreria con circa 14000 spettro con circa 50% presenza di amianto, 50% assenza di amianto

Questa libreria e' stata usata per addestrare una rete random forest con

import json, pickle
import numpy as np
import geopandas as gpd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedGroupKFold
from sklearn.metrics import (average_precision_score, roc_auc_score,
                             precision_recall_curve)
import joblib

# ---------------- parametri ----------------
GPKG = "C:/Users/l.innocenti/Documents/MF/random_forest.gpkg"
LAYER = "random_forest"                      # None = prima tabella; altrimenti nome tabella
FIELD_PROFILES = "profile"
FIELD_BG = "background"
CELL = 3000                       # metri, cella per i gruppi della CV spaziale
# intervalli (nm) da escludere: vapore acqueo e bordi rumorosi
DROP_RANGES = [(0, 550), (1300, 1500), (1750, 2150), (2400, 3000)]
N_TREES = 500
SEED = 0
OUT_MODEL = "C:/Users/l.innocenti/Documents/MF/rf_amianto.joblib"
# --------------------------------------------


def parse_profile(v):
    """Decodifica il campo profile di EnMAP-Box/QPS (dict, JSON testo o blob)."""
    if v is None:
        return None
    if isinstance(v, dict):                  # già decodificato (campo JSON del gpkg)
        return v
    if isinstance(v, (bytes, bytearray, memoryview)):
        b = bytes(v)
        try:
            return json.loads(b.decode("utf-8"))
        except Exception:
            return pickle.loads(b)           # solo su file tuoi, di cui ti fidi
    if isinstance(v, str):
        return json.loads(v)
    return None


df = gpd.read_file(GPKG, layer=LAYER)
print("Record letti:", len(df), "| campi:", list(df.columns))

# --- spettri ---
prof = [parse_profile(v) for v in df[FIELD_PROFILES]]
lengths = [len(p["y"]) if p and p.get("y") is not None else 0 for p in prof]
n_bands = max(set(lengths), key=lengths.count)          # lunghezza più frequente
ok = np.array([l == n_bands for l in lengths])
print(f"Bande: {n_bands} | profili scartati per lunghezza/vuoti: {(~ok).sum()}")

first = next(p for p, o in zip(prof, ok) if o)
wl = np.array(first["x"], float) if first.get("x") is not None else None
bbl = np.array(first["bbl"], bool) if first.get("bbl") is not None else np.ones(n_bands, bool)

band_mask = bbl.copy()
if wl is not None:
    for lo, hi in DROP_RANGES:
        band_mask &= ~((wl >= lo) & (wl <= hi))
else:
    print("ATTENZIONE: lunghezze d'onda assenti, uso solo bbl (nessun filtro per intervalli)")

X = np.array([np.asarray(p["y"], float) for p, o in zip(prof, ok) if o])
keep = np.where(ok)[0]

# --- etichette: 1 = amianto, 0 = background ---
bg = df[FIELD_BG].iloc[keep]
valid = bg.notna().to_numpy()
y = (~bg.fillna(False).astype(bool)).to_numpy().astype(int)

# --- pulizia spettri ---
X = X[:, band_mask]
good = valid & np.isfinite(X).all(axis=1) & (np.abs(X).sum(axis=1) > 0)
X, y, keep = X[good], y[good], keep[good]
print(f"Campioni: {len(y)} | amianto: {y.sum()} | background: {(y == 0).sum()} | bande usate: {X.shape[1]}")

# --- gruppi spaziali per la CV ---
groups = np.arange(len(y))                                # fallback: ogni campione un gruppo
if hasattr(df, "geometry") and df.geometry.notna().any():
    g = df.iloc[keep]
    if g.crs is not None and g.crs.is_geographic:
        g = g.to_crs(32632)
    has = g.geometry.notna().to_numpy()
    cx = (g.geometry.x.to_numpy()[has] // CELL).astype(np.int64)
    cy = (g.geometry.y.to_numpy()[has] // CELL).astype(np.int64)
    groups[has] = -(cx * 100000 + cy) - 1                 # id negativi, distinti dal fallback
else:
    print("ATTENZIONE: nessuna geometria, la CV NON è spaziale e le metriche saranno ottimistiche")

# --- modello ---
def make_rf():
    return RandomForestClassifier(
        n_estimators=N_TREES, max_features="sqrt", min_samples_leaf=3,
        class_weight="balanced_subsample", n_jobs=-1, random_state=SEED)

# --- cross-validation spaziale ---
cv = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=SEED)
oof = np.zeros(len(y))
for k, (tr, te) in enumerate(cv.split(X, y, groups), 1):
    m = make_rf().fit(X[tr], y[tr])
    oof[te] = m.predict_proba(X[te])[:, 1]
    print(f"fold {k}: PR-AUC = {average_precision_score(y[te], oof[te]):.3f}")

print(f"\nPR-AUC globale: {average_precision_score(y, oof):.3f} | ROC-AUC: {roc_auc_score(y, oof):.3f}")

prec, rec, thr = precision_recall_curve(y, oof)
f1 = 2 * prec[:-1] * rec[:-1] / np.clip(prec[:-1] + rec[:-1], 1e-9, None)
i = f1.argmax()
print(f"Soglia F1 ottimale: {thr[i]:.3f} (precision {prec[i]:.3f}, recall {rec[i]:.3f})")
for t in (0.5, 0.7, 0.9):
    j = np.searchsorted(thr, t)
    print(f"soglia {t}: precision {prec[j]:.3f}, recall {rec[j]:.3f}")

# --- modello finale su tutti i dati ---
rf = make_rf().fit(X, y)
joblib.dump({"model": rf, "band_mask": band_mask, "wavelengths": wl,
             "n_bands_total": n_bands, "threshold_f1": float(thr[i])}, OUT_MODEL)
print("Modello salvato in", OUT_MODEL)

# --- bande più importanti ---
imp = rf.feature_importances_
labels = wl[band_mask] if wl is not None else np.arange(X.shape[1])
for j in np.argsort(imp)[::-1][:10]:
    print(f"banda {labels[j]:.0f}: {imp[j]:.4f}")

 a questo punto ho preso una immagine sulla quale la rete non e' stata addestrata e fatto inferenza tramite

import numpy as np
import joblib
import rasterio
import geopandas as gpd
import xml.etree.ElementTree as ET
from rasterio.windows import Window
from scipy.ndimage import maximum_filter
from scipy.spatial import cKDTree

# ---------------- parametri ----------------
IMG = "C:/Users/l.innocenti/Documents/MF/ENMAP_SPECTRAL_IMAGE.TIF"    # stesso prodotto da cui vengono gli spettri
METADATA_XML = "C:/Users/l.innocenti/Documents/MF/ENMAP_METADATA.XML"  # METADATA.XML dello stesso prodotto
MODEL = "C:/Users/l.innocenti/Documents/MF/rf_amianto.joblib"
OUT_PROB = "C:/Users/l.innocenti/Documents/MF/prob_amianto.tif"

CENSUS = "amianto_PRA_pubblico.shp"                                    # punti del censimento
LIBRARY = "C:/Users/l.innocenti/Documents/MF/random_forest.gpkg"       # per escludere i punti usati nel training
LIBRARY_LAYER = "random_forest"

BUILT_MASK = None      # raster 0/1 sulla stessa griglia dell'immagine (1 = edificato): rende il "caso" onesto
MTMF_SCORE = None      # raster punteggio MTMF sulla stessa griglia, per il confronto a parità di pixel
NODATA_DEFAULT = -32768
BLOCK = 256
TOL_NM = 1.0           # scarto massimo ammesso tra lunghezza d'onda del modello e della banda abbinata
NEIGH = 2              # un punto è "hit" se c'è un pixel segnalato entro 2 pixel
EXCL_DIST = 150        # m: esclude i punti del censimento vicini a spettri usati nel training
TOP_FRACTIONS = [0.001, 0.005, 0.01, 0.02, 0.05]   # quota di pixel validi segnalati
# --------------------------------------------


def wavelengths_from_xml(xml_path, n):
    """Legge un solo blocco da n lunghezze d'onda (nm) dal METADATA.XML di EnMAP."""
    root = ET.parse(xml_path).getroot()
    for tag in ("wavelengthCenterOfBand", "waveLength"):
        arr = np.array([float(e.text) for e in root.iter()
                        if e.tag.split("}")[-1] == tag and e.text])
        if len(arr) >= n and len(arr) % n == 0:
            arr = arr[:n]
            if 400 < arr[0] < 440 and 2400 < arr[-1] < 2500:
                print(f"Lunghezze d'onda lette da <{tag}>: {arr[0]:.1f} ... {arr[-1]:.1f} nm")
                return arr
    raise SystemExit("Blocco di lunghezze d'onda non riconosciuto nell'XML")


# ===== 1) predizione a blocchi =====
d = joblib.load(MODEL)
rf, band_mask, n_model = d["model"], d["band_mask"], d["n_bands_total"]
wl_model = d["wavelengths"]
if wl_model is None:
    raise SystemExit("Il modello non ha lunghezze d'onda salvate: impossibile abbinare le bande.")

with rasterio.open(IMG) as src:
    wl_img = wavelengths_from_xml(METADATA_XML, src.count)

    # abbinamento banda-modello -> banda-immagine, con indice crescente
    # (il modello è un sottoinsieme ordinato delle bande dell'immagine)
    img_idx = np.empty(len(wl_model), dtype=int)
    prev = -1
    for k, w in enumerate(wl_model):
        cand = np.arange(prev + 1, len(wl_img))
        if len(cand) == 0:
            raise SystemExit("Abbinamento impossibile: finite le bande dell'immagine.")
        j = cand[np.abs(wl_img[cand] - w).argmin()]
        img_idx[k] = j
        prev = j
    err = np.abs(wl_img[img_idx] - wl_model)
    print(f"Bande modello: {n_model} | bande immagine: {src.count} | "
          f"scarto max: {err.max():.4f} nm | medio: {err.mean():.4f} nm")
    if err.max() > TOL_NM:
        raise SystemExit(f"Scarto oltre {TOL_NM} nm: le lunghezze d'onda non corrispondono, "
                         "libreria e immagine vengono da sensori/versioni diverse.")
    if len(np.unique(img_idx)) != len(img_idx):
        raise SystemExit("Due bande del modello cadono sulla stessa banda dell'immagine.")

    missing = np.setdiff1d(np.arange(len(wl_img)), img_idx)
    print(f"Bande dell'immagine non usate ({len(missing)}):", np.round(wl_img[missing], 1))

    read_idx = (img_idx + 1).tolist()          # rasterio: indici da 1

    nodata = src.nodata if src.nodata is not None else NODATA_DEFAULT
    H, W = src.height, src.width
    T, CRS = src.transform, src.crs

    prof = src.profile.copy()
    prof.update(count=1, dtype="float32", nodata=-1.0, compress="lzw",
                tiled=True, blockxsize=256, blockysize=256)
    prob = np.full((H, W), -1.0, dtype=np.float32)

    with rasterio.open(OUT_PROB, "w", **prof) as dst:
        for r0 in range(0, H, BLOCK):
            for c0 in range(0, W, BLOCK):
                h, w = min(BLOCK, H - r0), min(BLOCK, W - c0)
                win = Window(c0, r0, w, h)
                arr = src.read(read_idx, window=win)                  # (n_model, h, w)
                X = arr.reshape(len(read_idx), -1).T[:, band_mask].astype(np.float32)

                good = (np.isfinite(X).all(axis=1)
                        & (X != nodata).all(axis=1)
                        & (np.abs(X).sum(axis=1) > 0))
                p = np.full(X.shape[0], -1.0, dtype=np.float32)
                if good.any():
                    p[good] = rf.predict_proba(X[good])[:, 1]

                blk = p.reshape(h, w)
                prob[r0:r0 + h, c0:c0 + w] = blk
                dst.write(blk, 1, window=win)
            print(f"righe {min(r0 + BLOCK, H)}/{H}")
print("Probabilità salvata in", OUT_PROB)

# ===== 2) confronto con il censimento =====
valid = prob >= 0
eval_mask = valid.copy()
if BUILT_MASK:
    with rasterio.open(BUILT_MASK) as m:
        built = m.read(1) == 1
    if built.shape != valid.shape:
        raise SystemExit("BUILT_MASK non ha la stessa griglia dell'immagine")
    eval_mask &= built
else:
    print("ATTENZIONE: senza BUILT_MASK il 'caso' è calcolato su tutti i pixel validi "
          "(campagna, boschi...), quindi gli arricchimenti assoluti sono gonfiati")

cen = gpd.read_file(CENSUS)
cen = cen[cen.geometry.notna()].to_crs(CRS)
xs, ys = cen.geometry.x.to_numpy(), cen.geometry.y.to_numpy()

# esclude i punti vicini a spettri della libreria (positivi e background usati nel training)
lib = gpd.read_file(LIBRARY, layer=LIBRARY_LAYER)
lib = lib[lib.geometry.notna()].to_crs(CRS)
tree = cKDTree(np.column_stack([lib.geometry.x, lib.geometry.y]))
dist, _ = tree.query(np.column_stack([xs, ys]), k=1)
far = dist > EXCL_DIST

cols_f, rows_f = ~T * (xs, ys)
rows, cols = np.floor(rows_f).astype(int), np.floor(cols_f).astype(int)
inside = (rows >= 0) & (rows < H) & (cols >= 0) & (cols < W)
sel = inside & far
sel[sel] = eval_mask[rows[sel], cols[sel]]           # solo pixel valutabili
rows, cols = rows[sel], cols[sel]
n_pts = len(rows)
print(f"\nPunti del censimento: {len(cen)} | nell'immagine: {inside.sum()} "
      f"| lontani dal training e su pixel validi: {n_pts}")


def hit_stats(score, frac):
    thr = np.quantile(score[eval_mask], 1 - frac)
    flag = (score >= thr) & eval_mask
    dil = maximum_filter(flag.astype(np.uint8), size=2 * NEIGH + 1) > 0
    obs = dil[rows, cols].mean()
    chance = dil[eval_mask].mean()
    ci = 1.96 * np.sqrt(obs * (1 - obs) / n_pts)
    return thr, int(flag.sum()), obs, ci, chance


scores = {"RF": prob}
if MTMF_SCORE:
    with rasterio.open(MTMF_SCORE) as m:
        s = m.read(1).astype(np.float32)
    if s.shape != valid.shape:
        raise SystemExit("MTMF_SCORE non ha la stessa griglia dell'immagine")
    s[~np.isfinite(s)] = -np.inf
    scores["MTMF"] = s

print(f"\n{'metodo':6} {'% pixel':>8} {'n pixel':>9} {'soglia':>9} {'hit %':>14} {'caso %':>8} {'arricch.':>9}")
for name, sc in scores.items():
    for f in TOP_FRACTIONS:
        thr, npx, obs, ci, chance = hit_stats(sc, f)
        print(f"{name:6} {100*f:8.2f} {npx:9d} {thr:9.3f} "
              f"{100*obs:6.1f}±{100*ci:4.1f} {100*chance:8.1f} {obs/max(chance,1e-9):9.2f}")

 Punti del censimento: 56865 | nell'immagine: 4189 | lontani dal training e su pixel validi: 2925

metodo  % pixel   n pixel    soglia          hit %   caso %  arricch.
RF         0.10      1048     0.994    3.8± 0.7      1.2      3.23
RF         0.50      5238     0.970   14.9± 1.3      4.1      3.59
RF         1.00     10476     0.933   27.1± 1.6      6.7      4.03
RF         2.00     20952     0.876   47.4± 1.8     11.3      4.20
RF         5.00     52380     0.738   68.4± 1.7     21.3      3.22

Nella riga 2,00 segnali 20952 pixel, e il 47,4% dei punti del censimento ha un pixel segnalato a meno di 60 m. Se li avessi segnalati a caso, ti aspetteresti l'11,3%. Il rapporto è 47,4 / 11,3 ≈ 4,2.

Prendendo il file prob_amianto.tif e guardando quanti punti di verita' a terra cadono nell'intorno dei pixel 

import numpy as np
import rasterio
import geopandas as gpd
from scipy.spatial import cKDTree

# ---------------- parametri ----------------
PROB = "C:/Users/l.innocenti/Documents/MF/prob_amianto.tif"      # probabilità 0-1, nodata = -1
CENSUS = "C:/Users/l.innocenti/Documents/MF/amianto_PRA_pubblico.shp"
OUT = "C:/Users/l.innocenti/Documents/MF/punti_distanza.gpkg"    # None per non salvare

DIST = 60.0                                   # metri
THRESHOLDS = [0.7, 0.8, 0.9, 0.95, 0.99]      # un pixel è "segnalato" se prob >= soglia
N_CHANCE = 200_000                            # pixel validi campionati per stimare il "caso"
SEED = 0

# opzionale: escludi i punti entro EXCL_DIST m da spettri usati nel training
LIBRARY = None                                # es. "C:/.../random_forest.gpkg"
LIBRARY_LAYER = "random_forest"
EXCL_DIST = 150.0
# --------------------------------------------

with rasterio.open(PROB) as src:
    prob = src.read(1)
    T, CRS = src.transform, src.crs
    H, W = src.height, src.width
    px = abs(T.a)
    nodata = src.nodata if src.nodata is not None else -1.0

valid = np.isfinite(prob) & (prob != nodata) & (prob >= 0)
print(f"Immagine {W}x{H} | pixel {px:.0f} m | pixel validi: {valid.sum()}")

# ---------- punti ----------
cen = gpd.read_file(CENSUS)
cen = cen[cen.geometry.notna()].to_crs(CRS).reset_index(drop=True)
xs, ys = cen.geometry.x.to_numpy(), cen.geometry.y.to_numpy()
n_tot = len(cen)

cf, rf_ = ~T * (xs, ys)
cols, rows = np.floor(cf).astype(int), np.floor(rf_).astype(int)
inside = (rows >= 0) & (rows < H) & (cols >= 0) & (cols < W)
use = inside.copy()
use[inside] = valid[rows[inside], cols[inside]]          # solo punti su pixel validi
print(f"Punti totali: {n_tot} | dentro l'immagine: {inside.sum()} | su pixel validi: {use.sum()}")

if LIBRARY:
    lib = gpd.read_file(LIBRARY, layer=LIBRARY_LAYER)
    lib = lib[lib.geometry.notna()].to_crs(CRS)
    d_lib, _ = cKDTree(np.column_stack([lib.geometry.x, lib.geometry.y])).query(
        np.column_stack([xs, ys]), k=1)
    use &= d_lib > EXCL_DIST
    print(f"Dopo esclusione punti vicini al training (< {EXCL_DIST:.0f} m): {use.sum()}")

pts = np.column_stack([xs[use], ys[use]])
n = len(pts)
if n == 0:
    raise SystemExit("Nessun punto utilizzabile.")

# ---------- "caso": pixel validi casuali ----------
rng = np.random.default_rng(SEED)
vr, vc = np.nonzero(valid)
pick = rng.choice(len(vr), size=min(N_CHANCE, len(vr)), replace=False)
cx, cy = T * (vc[pick] + 0.5, vr[pick] + 0.5)
rand_pts = np.column_stack([cx, cy])

# ---------- distanza dal pixel segnalato più vicino (centro pixel) ----------
res = cen.loc[use, ["geometry"]].copy()
print(f"\nDistanza dal centro del pixel segnalato più vicino, soglia {DIST:.0f} m "
      f"(n punti = {n})\n")
print(f"{'soglia':>7} {'pixel segn.':>12} {'< %dm' % DIST:>8} {'>= %dm' % DIST:>8} "
      f"{'% entro':>8} {'caso %':>8} {'arricch.':>9}")

for t in THRESHOLDS:
    fr, fc = np.nonzero(valid & (prob >= t))
    if len(fr) == 0:
        print(f"{t:7.2f} {0:12d}   (nessun pixel sopra soglia)")
        continue
    fx, fy = T * (fc + 0.5, fr + 0.5)
    tree = cKDTree(np.column_stack([fx, fy]))

    d, _ = tree.query(pts, k=1)
    d_rand, _ = tree.query(rand_pts, k=1)
    near = d < DIST
    obs = near.mean()
    chance = (d_rand < DIST).mean()
    print(f"{t:7.2f} {len(fr):12d} {near.sum():8d} {(~near).sum():8d} "
          f"{100 * obs:8.1f} {100 * chance:8.1f} {obs / max(chance, 1e-9):9.2f}")
    res[f"d_{t:.2f}"] = d

# ---------- salvataggio ----------
if OUT:
    res.to_file(OUT, driver="GPKG")
    print("\nSalvato:", OUT, "(distanza in metri per ogni soglia)")


si ha che con una soglia di 0.7 si ha una percentuale di corretta detection entro 60 m del 65% dei casi (tirando a caso la percentuale sarebbe stata del 14.5% quindi il segnale e' ben presente) 

Rendendo la soglia piu' stringente la percentuale crolla miseramente 

soglia  pixel segn.    < 60m   >= 60m  % entro   caso %  arricch.
   0.70        60332     1912     1013     65.4     14.5      4.51
   0.80        38932     1590     1335     54.4     10.3      5.29
   0.90        16122      959     1966     32.8      5.2      6.28
   0.95         8183      514     2411     17.6      3.0      5.82
   0.99         1747      109     2816      3.7      0.8      4.43 

 

CEDev e CEmu TI-84 Plus CE

Programmazione di TI-84 Plus CE (identica alla TI-83 Premium CE, versione francese) con un emulatore    Per prima cosa si inizia con il comp...