Recaler les lignes de balayage des passes et accélérer le rendu
Troisième passe du calage vertical : chaque ligne de balayage (décalage et inclinaison due au roulis) est recalée contre le consensus des autres faisceaux, à toutes les échelles, avec un profil d'étalonnage par faisceau et par degré d'angle qui retire les écarts non linéaires en travers de la fauchée. Les lignes sans recouvrement sont corrigées contre leur propre faisceau. Efface les lignes en creux et la marche au bord de fauchée mesurées sur LHD_FXX_0999_6882 (validé sur des blocs jamais vus). Calcul vectorisé, CuPy si GPU ; la gigue par fenêtres de temps devient inutile quand scan_angle existe. Rendu plus rapide : encodage AVIF speed 9 (0,6 s au lieu de 4 s par dalle), classification IGN par extraction directe laspy au lieu de PDAL (4,9 s au lieu de 13,5 s), comblement des trous et gradients sur GPU, cache numba persistant dans l'image. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Antoine Jacquin
@ -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`.
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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:
|
||||
|
||||
@ -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())
|
||||
|
||||
@ -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()
|
||||
|
||||
@ -125,6 +125,16 @@ class SharedDEM:
|
||||
"""Compute gradient components lazily on first access."""
|
||||
if self._gradient is None:
|
||||
logger.debug(" → Calcul gradient...")
|
||||
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))
|
||||
@ -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()
|
||||
|
||||
Reference in New Issue
Block a user