diff --git a/AGENTS.md b/AGENTS.md index 1935cd9..d1ba12e 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -29,15 +29,20 @@ - **`generate_*` signature is strict**: `(dem_file, basename, vis_dir, resolution, shared=None)` returning `Path` on success, `None` on failure. IGN overlays (`ortho`, `topo`) omit `shared`. - **Return `None` on failure, never raise**: `dtm.py`, `visualizations.py`, and `ign.py` all return `None` to let the pipeline continue. Raising aborts the entire file. - **Logger is always `logging.getLogger("lidar")`**, never `__name__`. All modules route through this single logger so worker processes can configure it. +- **Ajustement conjoint des lignes de balayage (3ᵉ passe du calage, `STRIP_ALIGN_VERSION` 3)** : chaque ligne de balayage (~150 lignes/s, découpée par `_scan_line_ids` : retour de la dent de scie de `scan_angle`) peut être décalée ET inclinée (roulis) de 1 à 15 cm — lignes en creux isolées, passe entière basculée visible en bande (mesuré sur LHD_FXX_0999_6882). `_joint_line_corrections` (`dtm.py`) recale chaque ligne (a + b·u) contre le consensus des AUTRES faisceaux de sa maille 1 m (plan local par la pente), moindres carrés tronqués à 5 MAD par `bincount`, pas amorti 0,5, arrêt sous 1 mm, recentrage (pas de dérive d'ensemble), estimation sur 1 point sur 3 ; CuPy si GPU actif, repli numpy. S'y ajoute un **profil d'étalonnage par faisceau et par degré d'angle** (`STRIP_ANGLE_BIN`, commun à toutes les lignes du faisceau, moyenne retirée, pente conservée car l'inclinaison par ligne est indéterminable quand la fauchée ne traverse la dalle qu'en partie ; classes pauvres = valeur voisine) : il retire les écarts non linéaires en travers de la fauchée (±1,5 cm mesurés au bord), sources de marches parallèles au vol au bord de fauchée. Lignes sans recouvrement : `_scan_line_corrections_beam` (contre la surface de leur propre faisceau, composante ligne à ligne seule, σ 3 lignes). Validé sur blocs jamais vus (damier 10 m). La gigue par fenêtres de temps (2ᵉ passe) n'est plus calculée quand `scan_angle` existe (redondante). Coût CPU ~26 s par dalle pour tout le calage (73 s avant). Sidecar : `lines`, `line_window`, `line_cell`, `line_model`. - **Openness à échelle fixe** : `generate_openness` normalise par des références figées `OPENNESS_POS_REF` / `OPENNESS_NEG_REF` (degrés, médianes mesurées sur 15 dalles réparties sur le territoire) et plus par z-score de dalle — même ouverture = même couleur, mosaïque jointive. - **Filename special-cases** in `_expected_output_path()`: `pos_open` → `positive_openness`, `neg_open` → `negative_openness`, `hillshade` → `hillshade_multi`. +- **Génération imposée** : la carte ne propose plus de réglage de classification ni de raccord — `_build_command` (`mapserve.py`) passe toujours `--ground-classification ign --ign-classes sol --edge-buffer 100` ; les champs `ground_class`/`ign_classes`/`reclassify`/`edge_buffer` éventuellement envoyés sont ignorés. Défauts CLI alignés (`ign`, `100`). +- **Classification IGN par extraction directe** : `_extract_ign_ground` (`dtm.py`) filtre les classes avec laspy (retours ≥ 1) et écrit le LAS sol sans passer par PDAL (4,9 s au lieu de 13,5 s par dalle) ; la lecture de la détection automatique est réutilisée (`_LAST_READ`). PDAL reste le repli. +- **Encodage AVIF rapide** : `AVIF_SPEED = 9` (`rendering.py`, `tiles.py`, `_SUBTILE_AVIF_SPEED` dans `index.py`) — dalle 5000 × 5000 px encodée en 0,6 s au lieu de 4 s (+3 % de taille, −0,3 dB) ; l'encodage était l'étape la plus longue du rendu d'une couche. +- **Préparation sur GPU** : `_fill_nans` (transformée de distance `cupyx`) et les gradients de `SharedDEM` passent sur GPU quand il est actif (repli CPU). `NUMBA_CACHE_DIR=/tmp/numba-cache` (Dockerfile) évite la recompilation des noyaux numba à chaque worker. - **Default output is AVIF**, not WebP. Use `--format webp` for WebP. Quality default is 60 (visually lossless on smooth color ramps, ~÷3 vs q98). - **Calage vertical des faisceaux de vol** : chaque tuile mélange plusieurs passes (1-2 `PointSourceId` par passe) parfois biaisées verticalement de quelques cm (±2,5 cm mesurés sur 1000_6882). `create_dtm_fast` mesure l'offset robuste de chaque faisceau (points sol, maille 1 m, surface médiane itérée 3×) et retranche les offsets ≥ 0,5 cm (`STRIP_ALIGN_THRESHOLD` dans `dtm.py`) avant rastérisation. Offsets calculés **par tuile** (ils dérivent le long d'une ligne de vol : jamais de table globale), mémoïsés par LAS sol, consignés dans `DTM/*_dtm*_stripalign.json` (sidecar de cache : absent, ou version/seuil/paramètres différents ⇒ régénération du DTM). Désactivable : `--no-strip-align`. - **Gigue intra-faisceau (2ᵉ passe du calage)** : les lignes de balayage successives d'une MÊME passe peuvent être décalées verticalement de façon aléatoire (vibration capteur / bruit haute fréquence de trajectoire) — un offset constant par faisceau n'y suffit pas. `_strip_jitter_offsets` découpe chaque faisceau en fenêtres de temps GPS (`STRIP_JITTER_BIN` = 0,1 s, origine de temps propre à chaque faisceau), mesure l'offset robuste de chaque fenêtre contre la surface médiane des AUTRES faisceaux (maille 1 m partagée, ≥ `STRIP_JITTER_MIN_CELLS` = 40 cellules), lisse la série (médiane glissante `STRIP_JITTER_SMOOTH` = 5 fenêtres), la borne à ± `STRIP_JITTER_MAX` (10 cm) puis l'interpole au temps GPS de chaque point (`_apply_strip_jitter`) ; en recouvrement à deux faisceaux, chacun reçoit une série (chacun absorbe sa part). Requiert la dimension `gps_time` (silencieusement ignorée sinon). Sidecar version 2 (séries dans `jitter`), couverte par `--no-strip-align`. - **Couches produites et servies = `PANEL_VIZ` (`index.py`, aujourd'hui `('relief_oriente',)`)** : pipeline sans `--only`/`--skip` (`panel_steps()`), génération lancée depuis la carte (`_panel_viz_steps` dans `mapserve.py`) et couches servies (`tiles.available_layers` filtre : panneau, `/tiles/…`, TileJSON, WMTS, JOSM). Les autres visualisations restent calculables avec `--only`. Les tests de la carte lèvent la restriction via `_setup(..., panel=None)`. - **Relief orienté (`relief_oriente`, seule couche de la carte)** : image RGB unique (GeoTIFF uint8 3 bandes, rendue telle quelle comme ortho/topo via `RGB_KEYWORDS` dans `rendering.py`) — clarté CIELAB = openness positive **locale** (MNT − gaussienne `RELIEF_DETREND_M` = 10 m, rayons `RELIEF_RADII_M` = 5/10/20 m, 16 directions) 65 % + ombrage 35 % ; teinte = aspect, chroma fixe (`RELIEF_CHROMA`). Échelle log fixe `RELIEF_OPEN_RANGE` (pas de statistique par dalle) et support total 40 m < bande de raccord 100 m : dalles jointives. Rapide : détendance + rayons sur grille décimée ~0,8 m (`RELIEF_GRID_M`), noyau dédié qui n'accumule que la moyenne des angles (`_mean_horizon_*` : CuPy `RawKernel` sur GPU, numba parallèle sur CPU, numpy en repli), colorisation fusionnée (numba) ou vectorisée sans trigonométrie (CuPy) via une table L* × teinte (`_relief_lut`). ~5 s de calcul hors préparation sur CPU 12 cœurs. Tout changement de constante change le rendu : régénérer les dalles (`--only relief_oriente --force`). - **Openness sous-échantillonnée** : `generate_openness` calcule le lancé de rayons (l'étape la plus coûteuse : 532 s/tuile à 0,2 m sur CPU) sur une grille décimée par blocs (`OPENNESS_DOWNSAMPLE = 2` : max par bloc en positive, min en négative — préserve les reliefs qui bornent l'horizon) puis rééchantillonne en bilinéaire. Coût ÷ facteur³ : 532 s → 40 s (×13). Signal archéologique préservé (corr. 0,93 après lissage) ; la texture de bruit sub-métrique disparaît. `--openness-downsample 1` = pleine résolution. SVF et openness anisotrope ne sont PAS concernées. -- **Raccord des bords entre tuiles** : les rendus à grand noyau (openness/SVF : rayons 100 m ; LRM : 15 m) tronquent leur fenêtre au bord de dalle — bandes d'artefacts à chaque changement de tuile. `--edge-buffer N` (défaut 0 = off ; case « Raccord des bords » de la carte, `EDGE_BUFFER_METERS` = 100 m dans `mapserve.py`) fait rastériser le DTM sur la **dalle nominale 1 km alignée sur la grille** plus une bande de N m remplie avec les points sol des 8 LAZ voisines (`_neighbor_ground_points` dans `dtm.py` : PDAL en flux, découpe + filtre de classes IGN ; voisine absente = téléchargement automatique depuis le catalogue IGN avant le run, **isolée dans `input/edge_neighbors/`** pour ne pas gonfler le corpus des passes globales (dédupliqué sur tout le lot, `_fetch_edge_neighbors` dans `pipeline.py`) ; introuvable ou échec = bande vide). Les visualisations calculent sur l'emprise étendue puis `rendering.py` (`_core_tile_window`, via `tif_to_crop`/`tif_to_png`) recadre les sorties sur la dalle 1 km exacte lue dans le nom LHD — les AVIF restent des carrés 1 km alignés dans la mosaïque. Tampon consigné dans le tag GeoTIFF `LIDAR_EDGE_BUFFER` du DTM : changer `--edge-buffer` invalide le cache DTM automatiquement (tag absent = 0). Bandes voisines non calées par faisceaux (contexte seul, recadrée hors image finale). Coût : ~7 s de lecture par voisine + ~44 % de pixels en plus à 100 m/0,2 m. Nom hors pattern LHD : option ignorée (bornes d'en-tête, pas de recadrage). +- **Raccord des bords entre tuiles** : les rendus à grand noyau (openness/SVF : rayons 100 m ; LRM : 15 m) tronquent leur fenêtre au bord de dalle — bandes d'artefacts à chaque changement de tuile. `--edge-buffer N` (défaut 100 ; toujours appliqué par la génération depuis la carte, `EDGE_BUFFER_METERS` = 100 m dans `mapserve.py`, plus de case à cocher ; 0 = off en ligne de commande) fait rastériser le DTM sur la **dalle nominale 1 km alignée sur la grille** plus une bande de N m remplie avec les points sol des 8 LAZ voisines (`_neighbor_ground_points` dans `dtm.py` : PDAL en flux, découpe + filtre de classes IGN ; voisine absente = téléchargement automatique depuis le catalogue IGN avant le run, **isolée dans `input/edge_neighbors/`** pour ne pas gonfler le corpus des passes globales (dédupliqué sur tout le lot, `_fetch_edge_neighbors` dans `pipeline.py`) ; introuvable ou échec = bande vide). Les visualisations calculent sur l'emprise étendue puis `rendering.py` (`_core_tile_window`, via `tif_to_crop`/`tif_to_png`) recadre les sorties sur la dalle 1 km exacte lue dans le nom LHD — les AVIF restent des carrés 1 km alignés dans la mosaïque. Tampon consigné dans le tag GeoTIFF `LIDAR_EDGE_BUFFER` du DTM : changer `--edge-buffer` invalide le cache DTM automatiquement (tag absent = 0). Bandes voisines non calées par faisceaux (contexte seul, recadrée hors image finale). Coût : ~7 s de lecture par voisine + ~44 % de pixels en plus à 100 m/0,2 m. Nom hors pattern LHD : option ignorée (bornes d'en-tête, pas de recadrage). - **Tests use lazy imports inside each test function**, never at module top, to avoid importing CuPy/GDAL at import time. - **`_`-prefixed names are critical private**: `_create_ground_pipeline`, `_fallback_to_smrf`, `_fill_nans`, `_init_gpu`, `_process_file_standalone` — do not call from outside their module. - **`build_index()` écrit l'inventaire + les paliers sources** : `output/index_tiles.json` (dalles, couches, URLs versionnées — servi par `/api/tiles`), vignettes `index_thumbs/` (≈3,9 m/px + `_mid` 1,56 m/px) et quadrants `index_subtiles/` — paliers de la pyramide XYZ (`tiles.py`). Chaque tuile du run en cours porte ses coins WGS84 pour les cadres de progression. Plus d'HTML : l'interface vit dans `mapui.py`. diff --git a/Dockerfile b/Dockerfile index 42ccd28..d541bc3 100644 --- a/Dockerfile +++ b/Dockerfile @@ -30,6 +30,9 @@ RUN pip3 install --no-cache-dir -r requirements.txt # the pre-built wheel (sm_89 = RTX 4060 Ti, sm_120 = RTX 5060 Ti). # nvcc must be in PATH at runtime for JIT compilation. ENV CUPY_CUDA_COMPILE_WITH_CACHE=1 +# Cache des noyaux numba (relief orienté, accumulation) : sans lui, chaque +# worker (processus spawn par dalle) recompile ses noyaux (~1 s par dalle). +ENV NUMBA_CACHE_DIR=/tmp/numba-cache ENV PATH=/usr/local/cuda/bin:${PATH} RUN pip3 install --no-cache-dir cupy-cuda12x diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index 8f32f28..7407cb5 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -55,7 +55,30 @@ IGN_CLASS_NAMES = { # la surface médiane des autres faisceaux (même maille 1 m), la série est lissée # (médiane glissante) puis interpolée au temps GPS de chaque point. Requiert la # dimension gps_time (ignorée silencieusement sinon). -STRIP_ALIGN_VERSION = 2 +# +# 3ᵉ passe — décalage ligne à ligne : les fenêtres de 0,1 s (lissées sur +# 0,5 s) regroupent ~15 lignes de balayage (~150 lignes/s) et ne travaillent +# qu'en recouvrement. Or deux lignes SUCCESSIVES d'une même passe peuvent +# différer de 1 à 2 cm (motif alterné, mesuré sur LHD_FXX_0999_6882) : stries +# fines perpendiculaires au vol sur tout le MNT. Chaque faisceau est découpé +# en lignes (sauts de scan_angle en dents de scie), le décalage robuste de +# chaque ligne (décalage ET inclinaison le long de la ligne : le roulis +# bascule les lignes) est mesuré contre la surface de SON faisceau (maille 0,5 m, +# boîte 1,5 m, plan local ajusté à la position réelle du point, 3 itérations +# pour que la ligne ne fausse pas sa propre référence), et seule la composante +# ligne à ligne est retirée (série − lissage gaussien σ 3 lignes ; une médiane +# glissante suivrait un motif alterné au lieu de l'effacer) : les variations +# lentes restent aux passes précédentes. Fonctionne sans recouvrement ; requiert +# gps_time et scan_angle. +# +# Ajustement conjoint (prioritaire) : la surface d'un faisceau absorbe toute +# erreur plus large que sa boîte de référence ; là où plusieurs faisceaux se +# recouvrent, chaque ligne (décalage + inclinaison) est donc recalée contre le +# consensus des AUTRES faisceaux, à toutes les échelles (roulis lent d'une +# passe entière, lignes isolées très décalées), par pas amortis et itérés. +# La correction sur son propre faisceau ne sert plus qu'aux lignes sans +# recouvrement. +STRIP_ALIGN_VERSION = 3 STRIP_ALIGN_THRESHOLD = 0.005 # m : écart mini pour corriger un faisceau (0,5 cm) STRIP_ALIGN_CELL = 1.0 # m : maille de comparaison des faisceaux STRIP_ALIGN_MIN_SHARED = 500 # cellules sol communes mini pour valider un offset @@ -63,11 +86,31 @@ STRIP_JITTER_BIN = 0.1 # s : durée d'une fenêtre de temps GPS (gigue) STRIP_JITTER_SMOOTH = 5 # fenêtres : largeur de la médiane glissante STRIP_JITTER_MIN_CELLS = 40 # cellules sol partagées mini pour valider une fenêtre STRIP_JITTER_MAX = 0.10 # m : amplitude maxi d'une correction de gigue (garde-fou) +STRIP_LINE_CELL = 0.5 # m : maille de la surface de référence d'un faisceau +STRIP_LINE_BOX = 3 # mailles : lissage de la référence (boîte 1,5 m) +STRIP_LINE_WINDOW = 3 # lignes : σ du lissage gaussien retiré (garde le ligne-à-ligne) +STRIP_LINE_ITERS = 3 # itérations (atténuation par la ligne elle-même < 1 %) +STRIP_LINE_MAX = 0.05 # m : correction maxi d'une ligne (garde-fou) +STRIP_LINE_MIN_POINTS = 100 # points sol mini pour mesurer une ligne +STRIP_LINE_GAP = 0.05 # s : trou de temps qui coupe une ligne (fin de passe) +STRIP_LINE_MIN_RMS = 0.002 # m : faisceau laissé tel quel sous ce niveau +STRIP_LINE_MODEL = "conjoint+decalage+inclinaison+profil-angle" # modèle (consigné, invalide le cache) +_SCAN_ANGLE_UNIT = 0.006 # ° par unité de scan_angle (LAS 1.4) +STRIP_ANGLE_BIN = 1.0 # ° : pas du profil de correction par faisceau selon l'angle +STRIP_ANGLE_MIN_POINTS = 500 # points en recouvrement mini par classe d'angle (sinon valeur voisine) +STRIP_LINE_SUBSAMPLE = 3 # estimation sur 1 point sur 3 (correction appliquée à tous) +STRIP_JOINT_CELL = 1.0 # m : maille du consensus des autres faisceaux +STRIP_JOINT_ITERS = 8 # itérations maxi de l'ajustement conjoint +STRIP_JOINT_DAMPING = 0.5 # pas amorti : deux faisceaux se rapprochent sans se croiser +STRIP_JOINT_TOL = 0.001 # m : arrêt quand le pas moyen passe sous 1 mm +STRIP_JOINT_MIN_POINTS = 30 # points en recouvrement mini pour recaler une ligne +STRIP_JOINT_MAX = 0.15 # m : correction maxi d'un point (garde-fou) # Mémo des offsets par fichier : la classification est partagée entre # résolutions, le même LAS sol est rasterisé à 0,5 m puis 0,2 m. _STRIP_OFFSETS_CACHE = {} _STRIP_JITTER_CACHE = {} +_STRIP_LINES_CACHE = {} def _strip_surface_grid(x, y, z, inv, n_sources, cell=STRIP_ALIGN_CELL): @@ -278,6 +321,381 @@ def _strip_jitter_offsets(x, y, z, psid, t, cell=STRIP_ALIGN_CELL, return result +def _module_of(arr): + """Module (numpy ou cupy) d'un tableau.""" + mod = type(arr).__module__ + if mod.startswith("cupy"): + import cupy + return cupy + return np + + +def _array_module(use_gpu): + """cupy si le GPU est actif et demandé, sinon numpy.""" + if use_gpu: + from . import gpu as _gpu + if _gpu.is_gpu_active() and _gpu._cp is not None: + return _gpu._cp + return np + + +def _scan_line_ids(t, angle, gap=STRIP_LINE_GAP): + """Identifiant de ligne de balayage de points d'UN faisceau triés par temps. + + Nouvelle ligne à chaque saut de scan_angle de plus de la moitié de + l'amplitude dans le sens opposé au balayage (retour de la dent de scie) + ou à chaque trou de temps > gap (fin de passe). Les trous laissés par la + végétation retirée dans une ligne ne la coupent pas. + """ + m = _module_of(angle) + angle = m.asarray(angle, dtype=m.float64) + if len(angle) == 0: + return m.zeros(0, dtype=m.int64) + amp = float(m.percentile(angle, 99) - m.percentile(angle, 1)) + d = m.diff(angle) + moving = d[d != 0] + sweep = float(m.sign(m.median(moving))) if len(moving) else 1.0 + # Retour de la dent de scie : grand saut de sens OPPOSÉ au balayage. Un + # trou de végétation fait aussi sauter l'angle, mais dans le sens du + # balayage : il ne coupe pas la ligne. + new = m.concatenate([m.ones(1, dtype=bool), (d * sweep < -0.5 * amp) | (m.diff(t) > gap)]) + return m.cumsum(new) - 1 + + +def _group_median(values, groups, n_groups, min_count): + """Médiane des valeurs finies par groupe (vectorisée), NaN sous min_count.""" + finite = np.isfinite(values) + order = np.lexsort((np.where(finite, values, np.inf), groups)) + g = groups[order] + v = values[order] + nf = np.bincount(groups[finite], minlength=n_groups) + start = np.searchsorted(g, np.arange(n_groups)) + med = np.full(n_groups, np.nan) + ok = nf >= max(1, min_count) + lo = start[ok] + (nf[ok] - 1) // 2 + hi = start[ok] + nf[ok] // 2 + med[ok] = 0.5 * (v[lo] + v[hi]) + return med + + +def _line_fit(r, w, line, u, n_lines, min_points): + """Moindres carrés pondérés par ligne : r ≈ a + b·u (sommes par bincount). + + Returns: + (a, b, ok) ; b vaut 0 là où l'étendue de u ne permet pas d'estimer + une inclinaison (a seul, moyenne pondérée). + """ + m = _module_of(r) + rw = m.where(w, r, 0.0) + wf = w.astype(m.float64) + s0 = m.bincount(line, weights=wf, minlength=n_lines) + s1 = m.bincount(line, weights=wf * u, minlength=n_lines) + s2 = m.bincount(line, weights=wf * u * u, minlength=n_lines) + sr = m.bincount(line, weights=rw, minlength=n_lines) + sur = m.bincount(line, weights=rw * u, minlength=n_lines) + det = s0 * s2 - s1 * s1 + ok = s0 >= min_points + tilt = ok & (det > 1e-2 * m.maximum(s0, 1) ** 2) + a = m.where(ok, sr / m.maximum(s0, 1), 0.0) + b = m.zeros(n_lines) + safe = m.where(tilt, det, 1.0) + a = m.where(tilt, (s2 * sr - s1 * sur) / safe, a) + b = m.where(tilt, (s0 * sur - s1 * sr) / safe, b) + return a, b, ok + + +def _robust_mask(r, floor=0.03, k=5.0): + """Résidus finis sous k MAD (plancher floor) : écarte végétation basse, + points mal classés et bords de trous sans trier chaque ligne.""" + m = _module_of(r) + f = m.isfinite(r) + if not bool(f.any()): + return f + mad = 1.4826 * float(m.median(m.abs(r[f]))) + return f & (m.abs(r) < max(floor, k * mad)) + + +def _scan_line_corrections_beam(x, y, z, t, angle, cell=STRIP_LINE_CELL, + box=STRIP_LINE_BOX, window=STRIP_LINE_WINDOW, + iters=STRIP_LINE_ITERS, max_corr=STRIP_LINE_MAX, + min_points=STRIP_LINE_MIN_POINTS, + subsample=STRIP_LINE_SUBSAMPLE): + """Corrections ligne à ligne d'un faisceau contre SA surface (à SOUSTRAIRE). + + Chaque ligne est modélisée par un décalage ET une inclinaison le long de + la ligne (a + b·u, u = scan_angle normalisé) : une erreur de roulis + bascule la ligne, un bout plus haut que l'autre. Seule la composante + ligne à ligne est retirée (série − gaussienne σ window lignes). Estimation + sur 1 point sur subsample, moindres carrés tronqués (5 MAD) vectorisés. + + Returns: + (corr par point dans l'ordre d'entrée, décalages par ligne, + inclinaisons par ligne au bord de fauchée). + """ + from scipy.ndimage import gaussian_filter1d, uniform_filter + n = len(z) + o = np.argsort(t, kind="stable") + line_all = np.empty(n, dtype=np.int64) + line_all[o] = _scan_line_ids(t[o], angle[o]) + nl = int(line_all.max()) + 1 if n else 0 + tot_a, tot_b = np.zeros(nl), np.zeros(nl) + if nl < 10 * window: + return np.zeros(n), tot_a, tot_b + u_all = np.asarray(angle, dtype=np.float64) / max(np.percentile(np.abs(angle), 99), 1e-9) + sel = np.arange(n) % max(1, int(subsample)) == 0 + xs, ys, zs = x[sel], y[sel], np.asarray(z, dtype=np.float64)[sel] + line, u = line_all[sel], u_all[sel] + # Référence = plan local : moyennes (x, y, z) par boîte de mailles, pente + # de la surface lissée ; évalué à la position réelle du point (la moyenne + # d'une maille n'est pas en son centre : sur une pente, une interpolation + # au centre crée un biais qui dépend de la position de la ligne). + x0, y0 = xs.min(), ys.min() + ix = np.floor((xs - x0) / cell).astype(np.int64) + iy = np.floor((ys - y0) / cell).astype(np.int64) + W, H = int(ix.max()) + 1, int(iy.max()) + 1 + flat = iy * W + ix + count_b = uniform_filter(np.bincount(flat, minlength=W * H).reshape(H, W).astype(np.float64), box) + valid = count_b > 0 + inv_c = np.where(valid, 1.0 / np.maximum(count_b, 1e-12), np.nan) + mean_x = uniform_filter(np.bincount(flat, weights=xs - x0, minlength=W * H).reshape(H, W), box) * inv_c + mean_y = uniform_filter(np.bincount(flat, weights=ys - y0, minlength=W * H).reshape(H, W), box) * inv_c + dxp = (xs - x0) - mean_x.ravel()[flat] + dyp = (ys - y0) - mean_y.ravel()[flat] + min_pts = max(10, min_points // max(1, int(subsample))) + z_work = zs.copy() + for _ in range(iters): + mean_z = uniform_filter(np.bincount(flat, weights=z_work, minlength=W * H).reshape(H, W), box) * inv_c + gz_y, gz_x = np.gradient(np.where(valid, mean_z, np.nanmean(mean_z)), cell) + r = z_work - (mean_z.ravel()[flat] + gz_x.ravel()[flat] * dxp + gz_y.ravel()[flat] * dyp) + a, b, ok = _line_fit(r, _robust_mask(r), line, u, nl, min_pts) + if ok.sum() < 10 * window: + break + idx = np.flatnonzero(ok) + a = np.interp(np.arange(nl), idx, a[ok]) + b = np.interp(np.arange(nl), idx, b[ok]) + # Seule la composante ligne à ligne est retirée (une médiane glissante + # suivrait un motif alterné au lieu de l'effacer) + da = a - gaussian_filter1d(a, window, mode="nearest") + db = b - gaussian_filter1d(b, window, mode="nearest") + da[~ok] = 0.0 + db[~ok] = 0.0 + tot_a += da + tot_b += db + z_work = zs - np.clip(tot_a[line] + tot_b[line] * u, -max_corr, max_corr) + corr = np.clip(tot_a[line_all] + tot_b[line_all] * u_all, -max_corr, max_corr) + return corr, tot_a, tot_b + + +def _joint_line_corrections(x, y, z, psid, t, angle, cell=STRIP_JOINT_CELL, + iters=STRIP_JOINT_ITERS, damping=STRIP_JOINT_DAMPING, + tol=STRIP_JOINT_TOL, min_points=STRIP_JOINT_MIN_POINTS, + subsample=STRIP_LINE_SUBSAMPLE, max_corr=STRIP_JOINT_MAX, + use_gpu=True): + """Ajustement conjoint des lignes de tous les faisceaux contre le + consensus des AUTRES faisceaux (décalage + inclinaison par ligne). + + À chaque itération, le résidu d'un point est mesuré contre la moyenne des + autres faisceaux de sa maille (ramenée à sa position par la pente de la + surface), chaque ligne est ajustée par moindres carrés tronqués, et une + fraction damping du pas est appliquée à toutes les lignes à la fois ; + l'altitude moyenne est recentrée (pas de dérive d'ensemble). Sur GPU + (CuPy) si disponible, repli numpy sur toute erreur. + + Returns: + (corr à SOUSTRAIRE par point, id de ligne global par point, + lignes recalées (bool par ligne), nombre d'itérations) — numpy. + """ + m = _array_module(use_gpu) + if m is not np: + try: + return _joint_line_corrections_impl(m, x, y, z, psid, t, angle, cell, iters, + damping, tol, min_points, subsample, max_corr) + except Exception as e: + logger.warning(f" Ajustement conjoint GPU impossible ({e}) — repli CPU") + return _joint_line_corrections_impl(np, x, y, z, psid, t, angle, cell, iters, + damping, tol, min_points, subsample, max_corr) + + +def _joint_line_corrections_impl(m, x, y, z, psid, t, angle, cell, iters, damping, + tol, min_points, subsample, max_corr): + def host(a): + return a.get() if m is not np else a + n = len(z) + psid_h = np.asarray(psid) + beams, bidx_h = np.unique(psid_h, return_inverse=True) + # Classe d'angle par faisceau : profil d'étalonnage commun à toutes les + # lignes d'un faisceau (écart non linéaire en travers de la fauchée). + abin_h = np.round(np.asarray(angle, dtype=np.float64) * _SCAN_ANGLE_UNIT / STRIP_ANGLE_BIN).astype(np.int64) + amin = int(abin_h.min()) if n else 0 + n_ab = int(abin_h.max()) - amin + 1 if n else 1 + abin_h = bidx_h * n_ab + (abin_h - amin) + x, y = m.asarray(x, dtype=m.float64), m.asarray(y, dtype=m.float64) + z, t = m.asarray(z, dtype=m.float64), m.asarray(t, dtype=m.float64) + angle = m.asarray(angle, dtype=m.float64) + bidx = m.asarray(bidx_h) + gl = m.zeros(n, dtype=m.int64) + u = m.zeros(n) + base = 0 + for i in range(len(beams)): + idx = m.flatnonzero(bidx == i) + o = m.argsort(t[idx]) + gl[idx[o]] = _scan_line_ids(t[idx][o], angle[idx][o]) + base + base = int(gl[idx].max()) + 1 + u[idx] = angle[idx] / max(float(m.percentile(m.abs(angle[idx]), 99)), 1e-9) + n_lines = base + touched = m.zeros(n_lines, dtype=bool) + if len(beams) < 2 or n == 0: + return np.zeros(n), host(gl), host(touched), 0 + abin = m.asarray(abin_h) + n_bins = len(beams) * n_ab + sel = m.arange(0, n, max(1, int(subsample))) + xs, ys, zs = x[sel], y[sel], z[sel] + gs, us, bs, cs = gl[sel], u[sel], bidx[sel], abin[sel] + x0, y0 = float(xs.min()), float(ys.min()) + ix = m.floor((xs - x0) / cell).astype(m.int64) + iy = m.floor((ys - y0) / cell).astype(m.int64) + W, H = int(ix.max()) + 1, int(iy.max()) + 1 + nk = W * H + key = iy * W + ix + bkey = bs * nk + key + dxc = (xs - x0) - (ix + 0.5) * cell + dyc = (ys - y0) - (iy + 0.5) * cell + n_all = m.bincount(key, minlength=nk) + n_own = m.bincount(bkey, minlength=len(beams) * nk) + n_other = n_all[key] - n_own[bkey] + overlap = n_other >= 2 + min_pts = max(10, min_points // max(1, int(subsample))) + A = m.zeros(n_lines) + B = m.zeros(n_lines) + C = m.zeros(n_bins) + min_bin = max(30, STRIP_ANGLE_MIN_POINTS // max(1, int(subsample))) + it = 0 + for it in range(1, iters + 1): + zc = zs - m.clip(A[gs] + B[gs] * us + C[cs], -max_corr, max_corr) + s_all = m.bincount(key, weights=zc, minlength=nk) + s_own = m.bincount(bkey, weights=zc, minlength=len(beams) * nk) + mean = m.where(n_all > 0, s_all / m.maximum(n_all, 1), m.nan).reshape(H, W) + gy, gx = m.gradient(m.where(m.isfinite(mean), mean, m.nanmean(mean)), cell) + ref = ((s_all[key] - s_own[bkey]) / m.maximum(n_other, 1) + + gx.ravel()[key] * dxc + gy.ravel()[key] * dyc) + r = m.where(overlap, zc - ref, m.nan) + w = m.zeros(len(r), dtype=bool) + for i in range(len(beams)): + mb = bs == i + w[mb] = _robust_mask(r[mb]) + a, b, ok = _line_fit(r, w, gs, us, n_lines, min_pts) + touched |= ok + # Profil par faisceau et classe d'angle : moyenne tronquée du résidu + # restant après le pas de ligne ; sa moyenne (portée par les lignes) + # est retirée faisceau par faisceau. Sa pente est conservée : quand la + # fauchée ne traverse la dalle qu'en partie, l'inclinaison des lignes + # est indéterminée et seul le profil peut la porter. + r2 = m.where(w, r - (a[gs] + b[gs] * us), 0.0) + wf = w.astype(m.float64) + cnt = m.bincount(cs, weights=wf, minlength=n_bins) + prof = m.where(cnt >= min_bin, m.bincount(cs, weights=r2, minlength=n_bins) / m.maximum(cnt, 1), 0.0) + su = m.bincount(cs, weights=wf * us, minlength=n_bins) + ub = m.where(cnt > 0, su / m.maximum(cnt, 1), 0.0) + for i in range(len(beams)): + sl = slice(i * n_ab, (i + 1) * n_ab) + k = cnt[sl] >= min_bin + if int(k.sum()) >= 3: + wk = cnt[sl][k] + uk, pk = ub[sl][k], prof[sl][k] + um, pm = float((wk * uk).sum() / wk.sum()), float((wk * pk).sum() / wk.sum()) + fitted = prof[sl] - pm + # Classe trop pauvre (bord de fauchée, classe incomplète) : + # valeur de la classe valide la plus proche, pas zéro — c'est + # précisément au bord que l'écart est le plus fort. + idx = m.arange(n_ab, dtype=m.float64) + prof[sl] = m.interp(idx, idx[k], fitted[k]) + else: + prof[sl] = 0.0 + A += damping * a + B += damping * b + C += damping * prof + A -= float(m.mean(A[gs] + B[gs] * us + C[cs])) # pas de dérive d'ensemble + step_lines = float(m.mean(m.abs(a[ok]))) * damping if bool(ok.any()) else 0.0 + step_prof = float(m.max(m.abs(prof))) * damping + if max(step_lines, step_prof) < tol: + break + corr = m.clip(A[gl] + B[gl] * u + C[abin], -max_corr, max_corr) + return host(corr), host(gl), host(touched), it + + +def _rolling_median_fast(values, window): + """Médiane glissante centrée (bords répétés), vectorisée.""" + from numpy.lib.stride_tricks import sliding_window_view + half = window // 2 + return np.median(sliding_window_view(np.pad(values, half, mode="edge"), window), axis=1) + + +def _scan_line_corrections(x, y, z, psid, t, angle, min_rms=STRIP_LINE_MIN_RMS): + """Corrections ligne à ligne de tous les faisceaux d'une tuile. + + 1. Ajustement conjoint contre les autres faisceaux (lignes en recouvrement). + 2. Lignes sans recouvrement : composante ligne à ligne contre la surface + de leur propre faisceau. + z doit être déjà calé (offsets constants et gigue temporelle). + + Returns: + (corrections à SOUSTRAIRE par point, {psid: (lignes, rms, max)} des + faisceaux corrigés). + """ + psid = np.asarray(psid) + z = np.asarray(z, dtype=np.float64) + corr, gl, touched, _ = _joint_line_corrections(x, y, z, psid, t, angle) + stats = {} + for p in np.unique(psid): + m = psid == p + if int(m.sum()) < 50 * STRIP_LINE_MIN_POINTS: + continue + alone = ~touched[gl[m]] + if alone.mean() > 0.05: + c_self, per_line, _ = _scan_line_corrections_beam( + x[m], y[m], z[m] - corr[m], t[m], angle[m]) + c_beam = corr[m] + np.where(alone, c_self, 0.0) + else: + c_beam = corr[m] + rms = float(np.sqrt(np.mean(c_beam ** 2))) if m.any() else 0.0 + if rms < min_rms: + corr[m] = 0.0 + continue + corr[m] = c_beam + n_lines = int(len(np.unique(gl[m]))) + stats[int(p)] = (n_lines, rms, float(np.max(np.abs(c_beam)))) + return corr, stats + + +def _scan_angle(las): + """scan_angle (LAS 1.4) ou scan_angle_rank (LAS ≤ 1.3), None si absent.""" + for name in ("scan_angle", "scan_angle_rank"): + try: + return np.asarray(getattr(las, name), dtype=np.float64) + except AttributeError: + continue + return None + + +def _scan_lines_for_file(las_file, las, z_aligned, t): + """Corrections ligne à ligne d'un LAS sol, mémoïsées par (chemin, mtime).""" + try: + p = Path(las_file) + cache_key = (str(p), p.stat().st_mtime_ns) + except OSError: + cache_key = (str(las_file), 0) + if cache_key in _STRIP_LINES_CACHE: + return _STRIP_LINES_CACHE[cache_key] + angle = _scan_angle(las) + result = (np.zeros(len(z_aligned)), {}) + if angle is not None and len(angle) == len(z_aligned): + result = _scan_line_corrections( + np.asarray(las.x, dtype=np.float64), np.asarray(las.y, dtype=np.float64), + z_aligned, np.asarray(las.point_source_id), t, angle) + _STRIP_LINES_CACHE[cache_key] = result + return result + + def _apply_strip_jitter(psid, t, jitter): """Corrections de gigue interpolées au temps GPS de chaque point. @@ -357,7 +775,7 @@ def _strip_offsets_for_file(las_file, las): def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets, - jitter=None): + jitter=None, lines=None): """Consigne les calages appliqués (version, seuil, offsets et gigue). Le sidecar sert de suivi de cache : un DTM sans sidecar, ou produit avec @@ -381,6 +799,11 @@ def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets, "jitter_bin": STRIP_JITTER_BIN, "jitter_smooth": STRIP_JITTER_SMOOTH, "jitter": jitter_payload, + "line_window": STRIP_LINE_WINDOW, + "line_cell": STRIP_LINE_CELL, + "line_model": STRIP_LINE_MODEL, + "lines": {str(p): {"lines": n, "rms_m": round(r, 4), "max_m": round(m, 4)} + for p, (n, r, m) in (lines or {}).items()}, } try: sidecar = Path(dtm_dir) / f"{basename}_dtm{output_suffix}_stripalign.json" @@ -742,6 +1165,8 @@ def detect_ground_method(laz_file): if las is None: logger.info(f" → Méthode: SMRF (défaut — lecture impossible)") return 'smrf' + _LAST_READ.clear() + _LAST_READ[str(laz_file)] = las # réutilisé par l'extraction IGN total_points = len(las.points) if total_points == 0: @@ -801,6 +1226,39 @@ def detect_ground_method(laz_file): return method +_LAST_READ = {} + + +def _extract_ign_ground(laz_file, output_las, codes): + """Extraction directe (laspy) des classes IGN choisies vers un LAS. + + Même résultat que le pipeline PDAL (filtres ReturnNumber ≥ 1, + NumberOfReturns ≥ 1, classes) mais sans relecture ni conversion : + ~4 s au lieu de ~9 s par dalle, et la lecture de la détection automatique + est réutilisée. Format de points et échelles du fichier source conservés. + + Returns: + True si le fichier a été écrit avec au moins un point. + """ + import laspy + las = _LAST_READ.pop(str(laz_file), None) + if las is None: + las = laspy.read(str(laz_file)) + keep = ((np.asarray(las.return_number) >= 1) + & (np.asarray(las.number_of_returns) >= 1) + & np.isin(np.asarray(las.classification), np.asarray(sorted(codes)))) + if not keep.any(): + return False + header = laspy.LasHeader(point_format=las.header.point_format, + version=las.header.version) + header.scales = las.header.scales + header.offsets = las.header.offsets + out = laspy.LasData(header) + out.points = las.points[keep] + out.write(str(output_las)) + return True + + def classify_ground(laz_file, temp_dir, method='auto', force=False, ign_classes="sol"): """Classify ground points using PDAL ground classification filter. @@ -843,6 +1301,20 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False, ign_classes= logger.info(f" Reclassification forcée — suppression de {output_las.name}") output_las.unlink() + # Pré-classification IGN : extraction directe, PDAL en secours. + if method == 'ign': + try: + if _extract_ign_ground(laz_file, output_las, ign_codes or [2]): + logger.info(f" ✓ Classification sol IGN terminée (extraction directe)") + return output_las + logger.warning(" Aucun point des classes IGN demandées — repli SMRF") + output_las.unlink(missing_ok=True) + return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label) + except Exception as e: + logger.warning(f" Extraction IGN directe impossible ({e}) — pipeline PDAL") + output_las.unlink(missing_ok=True) + _LAST_READ.clear() + pipeline_json = _create_ground_pipeline(laz_file, output_las, method, ign_codes=ign_codes) pipeline_file = temp_dir / f"pipeline_{method_label}.json" @@ -1265,7 +1737,11 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, gps_time = np.asarray(las.gps_time, dtype=np.float64) except AttributeError: gps_time = None - if gps_time is not None and len(gps_time) == len(las.points): + # L'ajustement conjoint ligne à ligne (3ᵉ passe) recale déjà chaque + # ligne contre les autres faisceaux, à toutes les échelles : la gigue + # par fenêtres de temps n'est calculée que si scan_angle manque. + lines_possible = _scan_angle(las) is not None + if gps_time is not None and len(gps_time) == len(las.points) and not lines_possible: try: strip_jitter = _strip_jitter_for_file(las_file, las, strip_offsets) except Exception as e: @@ -1328,6 +1804,20 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, zs = zs - lut[np.asarray(las.point_source_id, dtype=np.int64)] if strip_jitter and gps_time is not None: zs = zs - _apply_strip_jitter(las.point_source_id, gps_time, strip_jitter) + # 3ᵉ passe : décalage ligne à ligne (après les deux premières). + strip_lines = {} + if strip_align and gps_time is not None and len(gps_time) == len(zs): + t_lines = time.perf_counter() + try: + line_corr, strip_lines = _scan_lines_for_file(las_file, las, zs, gps_time) + zs = zs - line_corr + except Exception as e: + logger.warning(f" Mesure du décalage ligne à ligne impossible ({e}) — non corrigé") + strip_lines = {} + for p_, (nl_, rms_, max_) in sorted(strip_lines.items()): + logger.info(f" Lignes PSID {p_} : rms {rms_ * 100:.1f} cm, max {max_ * 100:.1f} cm " + f"({nl_} lignes)") + logger.info(f" Décalage ligne à ligne : {time.perf_counter() - t_lines:.1f}s") # Points sol des tuiles voisines dans la bande de raccord (best-effort, # non calés par faisceau : bande de contexte, l'image finale est @@ -1412,7 +1902,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, if strip_align: _write_strip_align_sidecar(dtm_dir, basename, output_suffix, - strip_offsets, strip_jitter) + strip_offsets, strip_jitter, strip_lines) logger.info(f" ✓ DTM créé: {output_tif.name}") return output_tif diff --git a/lidar_pipeline/index.py b/lidar_pipeline/index.py index cd5c253..d638637 100644 --- a/lidar_pipeline/index.py +++ b/lidar_pipeline/index.py @@ -294,7 +294,7 @@ _CARTO_SUBTILED_VIZ = () # speed 0 (lent/meilleur) → 10 (rapide) ; 6+8 quasi identique en taille, # on garde 8 (2× plus rapide) pour 2500×2500 px. _SUBTILE_AVIF_QUALITY = 55 -_SUBTILE_AVIF_SPEED = 8 +_SUBTILE_AVIF_SPEED = 9 # Vignette de sous-tuile (px) : à la vue d'ensemble chaque sous-tuile active # décode sa vignette en navigateur — 160 px couvre l'affichage jusqu'à ~220 px diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index ea5bda0..22dc894 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -76,7 +76,8 @@ _file_filter = FilePrefixFilter() from .progress import report_event from .dtm import (classify_ground, create_dtm_fast, STRIP_ALIGN_VERSION, - STRIP_ALIGN_THRESHOLD, STRIP_JITTER_BIN, STRIP_JITTER_SMOOTH) + STRIP_ALIGN_THRESHOLD, STRIP_JITTER_BIN, STRIP_JITTER_SMOOTH, + STRIP_LINE_WINDOW, STRIP_LINE_CELL, STRIP_LINE_MODEL) from .visualizations import ( SharedDEM, generate_hillshade, generate_slope, generate_aspect, @@ -435,7 +436,10 @@ class LidarArchaeoPipeline: return (data.get("version") == STRIP_ALIGN_VERSION and abs(float(data.get("threshold", -1)) - STRIP_ALIGN_THRESHOLD) < 1e-9 and abs(float(data.get("jitter_bin", -1)) - STRIP_JITTER_BIN) < 1e-9 - and int(data.get("jitter_smooth", -1)) == STRIP_JITTER_SMOOTH) + and int(data.get("jitter_smooth", -1)) == STRIP_JITTER_SMOOTH + and int(data.get("line_window", -1)) == STRIP_LINE_WINDOW + and abs(float(data.get("line_cell", -1)) - STRIP_LINE_CELL) < 1e-9 + and data.get("line_model") == STRIP_LINE_MODEL) except Exception: return False diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index daaaa81..6beb02d 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -170,6 +170,12 @@ RGB_LEGENDS = { } RGB_KEYWORDS = tuple(RGB_LEGENDS) +# Vitesse d'encodage AVIF (libavif, 0 = lent/compact … 10 = rapide). Mesuré +# sur une dalle 5000 × 5000 px (q60) : défaut 4,1 s ; speed 9 0,6 s pour +3 % +# de taille et −0,3 dB de PSNR, invisible. L'encodage était l'étape la plus +# longue du rendu d'une couche. +AVIF_SPEED = 9 + # Fusion des textes de légende (titre / lecture du rendu / méthode de calcul) # depuis la source unique VIZ_LEGENDS (index.py, sans dépendance lourde). from .index import VIZ_LEGENDS, parse_basename_coords @@ -804,7 +810,8 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, if quality >= 100: img.save(str(output_file), format=pil_format, lossless=True) else: - img.save(str(output_file), format=pil_format, quality=quality) + img.save(str(output_file), format=pil_format, quality=quality, + **({'speed': AVIF_SPEED} if pil_format == 'AVIF' else {})) # Delete source TIFF (unless --keep-tif) if not keep_tif: @@ -888,7 +895,8 @@ def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, outpu if quality >= 100: img.save(str(output_file), format=pil_format, lossless=True) else: - img.save(str(output_file), format=pil_format, quality=quality) + img.save(str(output_file), format=pil_format, quality=quality, + **({'speed': AVIF_SPEED} if pil_format == 'AVIF' else {})) # Delete source TIFF (unless --keep-tif) if not keep_tif: diff --git a/lidar_pipeline/tests/test_dtm.py b/lidar_pipeline/tests/test_dtm.py index 9becae8..8d45fd2 100644 --- a/lidar_pipeline/tests/test_dtm.py +++ b/lidar_pipeline/tests/test_dtm.py @@ -661,6 +661,13 @@ class TestStripAlignSidecar: assert not p._strip_align_matches("TILE", "_r0p2") # Calage désactivé + DTM calé : à régénérer assert not P(tmp_path, strip_align=False)._strip_align_matches("TILE", "_r0p2") + # Paramètres ligne à ligne différents : à régénérer + _write_strip_align_sidecar(p.dtm_dir, "TILE", "_r0p2", {}) + data = json.loads(bad.read_text()) + data["line_window"] = 99 + bad.write_text(json.dumps(data)) + assert not p._strip_align_matches("TILE", "_r0p2") + class TestEdgeBuffer: @@ -784,3 +791,159 @@ class TestEdgeBuffer: with rasterio.open(str(dtm)) as src: assert abs(src.bounds.left - 0.05) < 1e-6 # bornes de l'en-tête assert read_dtm_edge_buffer(dtm) == 0.0 + + +def _synthetic_beam(n_lines=240, spacing=0.4, seed=0, offsets=None, tilts=None): + """Faisceau synthétique : lignes de balayage (dents de scie de scan_angle) + sur un terrain en pente traversé par un fossé ; offsets verticaux par ligne.""" + rng = np.random.default_rng(seed) + xs, ys, zs, ts, angs, ids = [], [], [], [], [], [] + x_line = np.arange(0, 100, 0.12) + for k in range(n_lines): + y = k * spacing + rng.normal(0, 0.02, x_line.size) + z = 100 + 0.03 * x_line + 0.05 * y + z = z - 0.5 * (np.abs(x_line - 41) < 1.0) # fossé perpendiculaire aux lignes + u = np.linspace(-1, 1, x_line.size) + z = z + offsets[k] + (0 if tilts is None else tilts[k]) * u + rng.normal(0, 0.01, x_line.size) + xs.append(x_line); ys.append(y); zs.append(z); ids.append(np.full(x_line.size, k)) + ts.append(k * 0.0067 + np.linspace(0, 0.006, x_line.size)) + angs.append(np.linspace(-3300, 3300, x_line.size)) + return (np.concatenate(xs), np.concatenate(ys), np.concatenate(zs), + np.concatenate(ts), np.concatenate(angs), np.concatenate(ids)) + + +class TestScanLineAlignment: + """3ᵉ passe du calage : décalage vertical entre lignes de balayage successives.""" + + def test_line_ids_follow_sawtooth_not_vegetation_gaps(self): + from lidar_pipeline.dtm import _scan_line_ids + t = np.arange(50) * 0.0001 + ang = np.tile(np.linspace(-3000, 3000, 10), 5) + keep = np.ones(50, bool); keep[13:17] = False # trou de végétation dans la ligne 2 + ids = _scan_line_ids(t[keep], ang[keep]) + assert ids.max() == 4 + t2 = t.copy(); t2[30:] += 1.0 # fin de passe : nouvelle ligne + assert _scan_line_ids(t2, np.zeros(50)).max() == 1 + + def test_group_median_matches_numpy(self): + from lidar_pipeline.dtm import _group_median + rng = np.random.default_rng(3) + g = rng.integers(0, 20, 2000) + v = rng.normal(size=2000); v[::17] = np.nan + med = _group_median(v, g, 20, 1) + for k in range(20): + vals = v[(g == k) & np.isfinite(v)] + assert med[k] == pytest.approx(np.median(vals)) + + def test_removes_alternating_line_offsets_and_keeps_ditch(self): + from lidar_pipeline.dtm import _scan_line_corrections_beam + n = 240 + rng = np.random.default_rng(1) + offsets = 0.015 * (-1.0) ** np.arange(n) + rng.normal(0, 0.006, n) + x, y, z, t, ang, ids = _synthetic_beam(n, offsets=offsets) + corr, per_line, _ = _scan_line_corrections_beam(x, y, z, t, ang) + assert len(per_line) == n + from scipy.ndimage import gaussian_filter1d + hp = lambda v: v - gaussian_filter1d(v, 3, mode="nearest") + before, after = hp(offsets)[10:-10].std(), hp(offsets - per_line)[10:-10].std() + assert after < 0.2 * before, f"{before*1000:.1f} → {after*1000:.1f} mm" + ditch = np.abs(x - 41) < 0.8 + depth = lambda zz: np.median(zz[~ditch & (np.abs(x - 41) < 4)]) - np.median(zz[ditch]) + assert depth(z - corr) == pytest.approx(depth(z), abs=0.01) + + def test_removes_alternating_line_tilts(self): + """Roulis : lignes basculées alternativement (un bout haut, l'autre bas).""" + from scipy.ndimage import gaussian_filter1d + from lidar_pipeline.dtm import _scan_line_corrections_beam + n = 240 + rng = np.random.default_rng(4) + tilts = 0.02 * (-1.0) ** np.arange(n) + rng.normal(0, 0.008, n) + x, y, z, t, ang, ids = _synthetic_beam(n, offsets=np.zeros(n), tilts=tilts) + corr, _, per_tilt = _scan_line_corrections_beam(x, y, z, t, ang) + hp = lambda v: v - gaussian_filter1d(v, 3, mode="nearest") + before, after = hp(tilts)[10:-10].std(), hp(tilts - per_tilt)[10:-10].std() + assert after < 0.2 * before, f"{before*1000:.1f} → {after*1000:.1f} mm" + + def test_joint_adjustment_fixes_all_scales_against_other_beam(self): + """Deux faisceaux superposés : le faisceau fautif (dérive lente, roulis, + ligne isolée à −8 cm) est recalé sur l'autre à toutes les échelles, + sans dérive de l'altitude d'ensemble.""" + from lidar_pipeline.dtm import _joint_line_corrections + n = 240 + k = np.arange(n) + bad_off = 0.03 * np.sin(2 * np.pi * k / 120) # dérive lente (non vue sur sa propre surface) + bad_off[100] -= 0.08 # ligne isolée très décalée + bad_tilt = 0.02 * np.cos(2 * np.pi * k / 60) # roulis lent + xa, ya, za, ta, aa, _ = _synthetic_beam(n, seed=1, offsets=bad_off, tilts=bad_tilt) + xb, yb, zb, tb, ab, _ = _synthetic_beam(n, seed=2, offsets=np.zeros(n)) + tb = tb + 1000.0 + x = np.r_[xa, xb]; y = np.r_[ya, yb]; z = np.r_[za, zb]; t = np.r_[ta, tb] + ang = np.r_[aa, ab]; psid = np.r_[np.full(len(za), 1), np.full(len(zb), 2)] + corr, gl, touched, iters = _joint_line_corrections(x, y, z, psid, t, ang) + truth = np.r_[bad_off[np.repeat(k, len(za) // n)] + bad_tilt[np.repeat(k, len(za) // n)] + * np.tile(np.linspace(-1, 1, len(za) // n), n), np.zeros(len(zb))] + # Sans vérité terrain, l'écart est partagé entre les faisceaux : c'est + # l'écart ENTRE faisceaux (points homologues, même géométrie) qui doit + # disparaître. + na = len(za) + before = truth[:na] - truth[:na].mean() + rel = (truth[:na] - corr[:na]) - (0.0 - corr[na:]) + after = rel - rel.mean() + assert np.std(after) < 0.25 * np.std(before), \ + f"{np.std(before)*1000:.1f} → {np.std(after)*1000:.1f} mm" + assert abs(np.mean(corr)) < 0.002 # pas de dérive d'ensemble + assert touched.mean() > 0.9 and iters <= 8 + + def test_joint_adjustment_removes_static_angle_profile(self): + """Étalonnage en arc selon l'angle (même pour toutes les lignes) : + non linéaire, invisible pour décalage + inclinaison, retiré par le + profil par faisceau et classe d'angle.""" + from lidar_pipeline.dtm import _joint_line_corrections + n = 240 + xa, ya, za, ta, aa, _ = _synthetic_beam(n, seed=5, offsets=np.zeros(n)) + uu = aa / 3300.0 + prof = 0.02 * (uu ** 2 - np.mean(uu ** 2)) + za = za + prof + xb, yb, zb, tb, ab, _ = _synthetic_beam(n, seed=6, offsets=np.zeros(n)) + x = np.r_[xa, xb]; y = np.r_[ya, yb]; z = np.r_[za, zb]; t = np.r_[ta, tb + 1000.0] + ang = np.r_[aa, ab]; psid = np.r_[np.full(len(za), 1), np.full(len(zb), 2)] + corr, _, _, _ = _joint_line_corrections(x, y, z, psid, t, ang) + na = len(za) + rel = (prof - corr[:na]) - (0.0 - corr[na:]) + assert np.std(rel - rel.mean()) < 0.3 * np.std(prof), \ + f"{np.std(prof)*1000:.1f} → {np.std(rel - rel.mean())*1000:.1f} mm" + + def test_clean_beam_left_untouched(self): + from lidar_pipeline.dtm import _scan_line_corrections + x, y, z, t, ang, ids = _synthetic_beam(240, offsets=np.zeros(240)) + corr, stats = _scan_line_corrections(x, y, z, np.full(len(z), 7), t, ang) + assert stats == {} and not corr.any() + + +class TestIgnDirectExtraction: + """Classification IGN : extraction directe par laspy (PDAL en secours).""" + + def test_keeps_requested_classes_and_valid_returns(self, tmp_path): + import laspy + from lidar_pipeline.dtm import _extract_ign_ground + n = 1000 + rng = np.random.default_rng(0) + hdr = laspy.LasHeader(point_format=6, version="1.4") + hdr.scales = np.array([0.01, 0.01, 0.001]); hdr.offsets = np.array([1000.0, 6800000.0, 0.0]) + las = laspy.LasData(hdr) + las.x = 1000 + rng.uniform(0, 100, n); las.y = 6800000 + rng.uniform(0, 100, n) + las.z = rng.uniform(100, 110, n) + cls = rng.choice([1, 2, 3, 6, 9], n); las.classification = cls + rn = np.ones(n, dtype=np.uint8); rn[:10] = 0; las.return_number = rn + las.number_of_returns = np.ones(n, dtype=np.uint8) + las.point_source_id = np.full(n, 42, dtype=np.uint16) + src = tmp_path / "t.las"; las.write(str(src)) + out = tmp_path / "g.las" + assert _extract_ign_ground(src, out, [2]) + g = laspy.read(str(out)) + expected = (cls == 2) & (rn >= 1) + assert len(g.points) == int(expected.sum()) + assert set(np.unique(np.asarray(g.classification))) == {2} + assert g.header.point_format.id == 6 and np.all(np.asarray(g.point_source_id) == 42) + assert _extract_ign_ground(src, tmp_path / "e.las", [1, 2]) + assert len(laspy.read(str(tmp_path / "e.las")).points) == int(((np.isin(cls, [1, 2])) & (rn >= 1)).sum()) diff --git a/lidar_pipeline/tiles.py b/lidar_pipeline/tiles.py index 727d021..2946a0f 100644 --- a/lidar_pipeline/tiles.py +++ b/lidar_pipeline/tiles.py @@ -30,6 +30,7 @@ TILE_MAX_NATIVE_Z = 19 # 0,2 m/px ≈ résolution du z19 à la latitude TILE_DIRNAME = "index_xyz" # cache disque, sous le dossier de sortie WEBP_QUALITY = 78 AVIF_QUALITY = 60 +AVIF_SPEED = 9 # encodage rapide (cf. rendering.AVIF_SPEED : ×7, +3 % de taille) # PNG palettisé (PNG8 + alpha) : ~5× plus léger (170 → 32 Ko sur une dalle # réelle) pour un écart moyen de ~4 niveaux sur une rampe de couleur. Laissé # DÉSACTIVÉ par défaut : le PNG canonique reste sans perte, la fidélité prime @@ -727,7 +728,7 @@ def _encode(img, fmt): elif fmt == "webp": img.save(buf, format="WEBP", quality=WEBP_QUALITY) elif fmt == "avif": - img.save(buf, format="AVIF", quality=AVIF_QUALITY) + img.save(buf, format="AVIF", quality=AVIF_QUALITY, speed=AVIF_SPEED) else: raise ValueError(f"format de tuile inconnu : {fmt}") return buf.getvalue() diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 218c807..5dfb823 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -125,8 +125,18 @@ class SharedDEM: """Compute gradient components lazily on first access.""" if self._gradient is None: logger.debug(" → Calcul gradient...") - dy = np.gradient(self.filled, self.resolution, axis=0) - dx = np.gradient(self.filled, self.resolution, axis=1) + dy = dx = None + if _gpu_mod.is_gpu_active() and self.filled_gpu is not None: + try: + g = self.filled_gpu + dy = to_cpu(xp.gradient(g, self.resolution, axis=0)) + dx = to_cpu(xp.gradient(g, self.resolution, axis=1)) + except Exception as e: + logger.warning(f" Gradient GPU impossible ({e}) — repli CPU") + dy = dx = None + if dy is None: + dy = np.gradient(self.filled, self.resolution, axis=0) + dx = np.gradient(self.filled, self.resolution, axis=1) slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) slope_deg = np.degrees(slope_rad) self._gradient = (dy, dx, slope_rad, slope_deg) @@ -223,6 +233,21 @@ def _fill_nans(arr): nan_mask = np.isnan(arr) if not np.any(nan_mask): return arr, nan_mask + # GPU (cupyx) : la transformée de distance sur 25 M de pixels coûte ~3 s + # sur CPU, l'étape la plus longue de la préparation une fois le reste sur GPU. + if _gpu_mod.is_gpu_active(): + try: + from cupyx.scipy.ndimage import distance_transform_edt as edt_gpu + cp = _gpu_mod._cp + _, idx = edt_gpu(cp.asarray(nan_mask), return_distances=False, return_indices=True) + iy, ix = cp.asnumpy(idx[0]), cp.asnumpy(idx[1]) + del idx + gpu_cleanup() + filled = arr.copy() + filled[nan_mask] = arr[iy[nan_mask], ix[nan_mask]] + return filled, nan_mask + except Exception as e: + logger.warning(f" Comblement GPU impossible ({e}) — repli CPU") from scipy.ndimage import distance_transform_edt _, (iy, ix) = distance_transform_edt(nan_mask, return_indices=True) filled = arr.copy()