From 41206f65de7b7e9bbe7c590b6360a679ecbdcbd5 Mon Sep 17 00:00:00 2001 From: Antoine Jacquin Date: Sat, 19 Sep 2026 00:53:22 +0200 Subject: [PATCH] =?UTF-8?q?Caler=20les=20faisceaux=20de=20vol,=20d=C3=A9ci?= =?UTF-8?q?mer=20l'openness=20et=20cibler=20les=20couches=20affich=C3=A9es?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Trois chantiers liés à la qualité et au coût des rendus : - Calage vertical des faisceaux : les passes d'une tuile peuvent être biaisées de quelques cm (±2,5 cm mesurés sur 1000_6882), créant des marches aux recouvrements. Les offsets par PointSourceId sont mesurés sur les points sol de la tuile (réf. médiane itérée) et retranchés ≥ 0,5 cm avant rastérisation, avec sidecar de cache et application au plancher bare-earth. - Openness décimée ×2 : lancé de rayons sur grille par blocs (max/min) puis rééchantillonnage bilinéaire — 532 s → 40 s par tuile à 0,2 m sur CPU, signal archéologique préservé. Réglable --openness-downsample. - Génération webapp concentrée sur les couches affichées : les défauts /api/generate, /api/preview et le sélecteur génèrent le panneau complet (slope, aspect, pos_open) au lieu d'aspect seul. --- .gitignore | 1 + AGENTS.md | 2 + lidar_pipeline/cli.py | 19 ++ lidar_pipeline/dtm.py | 194 +++++++++++++++++++- lidar_pipeline/index.py | 4 +- lidar_pipeline/pipeline.py | 47 ++++- lidar_pipeline/tests/test_dtm.py | 86 +++++++++ lidar_pipeline/tests/test_visualizations.py | 46 +++++ lidar_pipeline/tests/test_webapp.py | 11 +- lidar_pipeline/visualizations.py | 45 ++++- lidar_pipeline/webapp.py | 25 ++- 11 files changed, 455 insertions(+), 25 deletions(-) diff --git a/.gitignore b/.gitignore index cb49126..24831de 100644 --- a/.gitignore +++ b/.gitignore @@ -60,3 +60,4 @@ notebooks/ # Éventuels fichiers de cache matplotlib matplotlibrc docker-compose.webapp.override.yml +output-test/ diff --git a/AGENTS.md b/AGENTS.md index 2cc345f..7304721 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -28,6 +28,8 @@ - **Logger is always `logging.getLogger("lidar")`**, never `__name__`. All modules route through this single logger so worker processes can configure it. - **Filename special-cases** in `_expected_output_path()`: `pos_open` → `positive_openness`, `neg_open` → `negative_openness`, `hillshade` → `hillshade_multi`. - **Default output is AVIF**, not WebP. Use `--format webp` for WebP. Quality default is 98. +- **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 différents ⇒ régénération du DTM), appliqués aussi au plancher `--bare-earth`. Désactivable : `--no-strip-align`. +- **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és. - **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()` writes 3 files**: `output/index.html` (data shell, `const TILES` embedded), `output/assets/app.css` and `output/assets/app.js` (source: `_APP_CSS`/`_APP_JS` constants in `index.py`). `webapp.py` serves `/assets` with no-cache headers. Each tile carries `meta` — ground method read from `DTM/*_dtm{_rXpY}_method.txt` (falls back to the primary-resolution sidecar) + per-viz dates/sizes. diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py index f74f671..0c4f90e 100644 --- a/lidar_pipeline/cli.py +++ b/lidar_pipeline/cli.py @@ -150,6 +150,23 @@ def main(): "Requalifie le point le plus bas de chaque colonne en terrain — utile sous " "végétation dense ou en relief raide où la classification du sol sous-couvre le terrain." ) + parser.add_argument( + "--openness-downsample", + type=int, + default=None, + metavar="FACTEUR", + help="Facteur de sous-échantillonnage du calcul d'openness (lancé de rayons) : " + "grille décimée par blocs puis rééchantillonnage. Défaut : 2 (~8× plus " + "rapide, rendu quasi identique) ; 1 = pleine résolution" + ) + parser.add_argument( + "--no-strip-align", + action="store_true", + help="Désactiver le calage vertical des faisceaux de vol. Par défaut, les écarts " + "verticaux ≥ 0,5 cm entre lignes de vol (PointSourceId) d'une tuile sont mesurés " + "sur les points sol et corrigés avant rastérisation du MNT (offsets consignés " + "dans DTM/*_stripalign.json)" + ) parser.add_argument( "--keep-tif", action="store_true", @@ -344,6 +361,8 @@ def main(): force_classify=args.force_classification, keep_tif=args.keep_tif, bare_earth=args.bare_earth, + strip_align=not args.no_strip_align, + openness_downsample=args.openness_downsample, quality=quality, only_viz=only_viz, skip_viz=skip_viz, diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index c65aee1..0da25a9 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -31,6 +31,144 @@ IGN_CLASS_NAMES = { } +# Calage vertical des faisceaux de vol (strip alignment). Une tuile LiDAR HD +# est couverte par plusieurs passes d'acquisition, chacune portée par un ou +# deux PointSourceId (faisceaux). Les passes sont bien alignées horizontalement +# mais certains faisceaux portent un biais vertical de quelques cm (mesuré +# jusqu'à ~5 cm sur LHD_FXX_1000_6882 : les deux faisceaux d'une passe écartés +# de ±2,5 cm). À 0,2 m/px ces écarts créent des marches et du bruit aux +# coutures des zones de recouvrement. On mesure l'offset robuste de chaque +# faisceau sur les points sol de la tuile elle-même (les PointSourceId changent +# selon la campagne, rien n'est codé en dur) et on le retranche avant +# rastérisation. Les offsets dérivent le long d'une ligne de vol (signe inversé +# entre tuiles voisines mesuré) : le calcul est donc par tuile, jamais global. +STRIP_ALIGN_VERSION = 1 +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 + +# 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 = {} + + +def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL, + threshold=STRIP_ALIGN_THRESHOLD, + min_shared=STRIP_ALIGN_MIN_SHARED): + """Mesure les offsets verticaux relatifs entre faisceaux d'une tuile. + + Méthode : surface sol par faisceau (moyenne des points à moins de 0,5 m du + minimum de chaque maille de 1 m, robuste à la végétation résiduelle), puis + offset de chaque faisceau = médiane de l'écart à la surface médiane de + référence (itéré 3 fois), sur les seules cellules couvertes par au moins + deux faisceaux. Ne corrige rien d'autre que le vertical : le calage + horizontal mesuré sur les données est excellent (≤ 1 cm). + + Args: + x, y, z, psid: coordonnées et PointSourceId des points sol. + cell: taille de maille de comparaison (m). + threshold: seuil (m) en dessous duquel un offset est ignoré. + min_shared: cellules communes minimales pour valider un faisceau. + + Returns: + dict {psid: offset} des offsets à SOUSTRAIRE (z - offset), ne + contenant que les |offset| >= threshold ; vide si rien à corriger + (faisceau unique, pas de recouvrement, tuile déjà alignée). + """ + us, inv = np.unique(psid, return_inverse=True) + if len(us) < 2: + return {} + x0 = np.floor(np.min(x) / cell) * cell + y0 = np.floor(np.min(y) / cell) * cell + xi = ((x - x0) / cell).astype(np.int64) + yi = ((y - y0) / cell).astype(np.int64) + ny = int(yi.max()) + 1 + key = xi * ny + yi + + def _surface(k): + m = inv == k + kk, zz = key[m], z[m] + order = np.argsort(kk, kind='stable') + k_s, z_s = kk[order], zz[order] + starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]]) + mins = np.minimum.reduceat(z_s, starts) + ukey = k_s[starts] + thr = mins[np.searchsorted(ukey, kk)] + sel = np.flatnonzero(zz <= thr + 0.5) + k2, z2 = kk[sel], zz[sel] + cnt = np.bincount(k2) + sums = np.bincount(k2, weights=z2) + v = np.flatnonzero(cnt >= 2) + return v, sums[v] / cnt[v] + + surfaces = [_surface(k) for k in range(len(us))] + populated = [c for c, _ in surfaces if len(c)] + if not populated: + return {} + allcells = np.unique(np.concatenate(populated)) + grid = np.full((len(us), len(allcells)), np.nan) + for k, (c, zs) in enumerate(surfaces): + grid[k, np.searchsorted(allcells, c)] = zs + comparable = np.sum(~np.isnan(grid), axis=0) >= 2 + + offsets = np.zeros(len(us)) + for _ in range(3): + ref = np.nanmedian(grid + offsets[:, None], axis=0) + for k in range(len(us)): + m = ~np.isnan(grid[k]) & comparable + if int(m.sum()) >= min_shared: + offsets[k] = np.median(grid[k][m] - ref[m]) + + return {int(p): round(float(offsets[k]), 3) + for k, p in enumerate(us) if abs(offsets[k]) >= threshold} + + +def _strip_offsets_for_file(las_file, las): + """Offsets de calage d'un LAS sol, mémoïsés 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_OFFSETS_CACHE: + return _STRIP_OFFSETS_CACHE[cache_key] + try: + psid = np.asarray(las.point_source_id) + except AttributeError: + psid = None # dimension absente (producteur tiers) : pas de calage + if psid is not None and len(psid) == len(las.points): + offsets = _strip_vertical_offsets( + np.asarray(las.x, dtype=np.float64), + np.asarray(las.y, dtype=np.float64), + np.asarray(las.z, dtype=np.float64), + psid) + else: + offsets = {} + _STRIP_OFFSETS_CACHE[cache_key] = offsets + return offsets + + +def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets): + """Consigne les offsets de calage appliqués (version et seuil inclus). + + Le sidecar sert de suivi de cache : un DTM sans sidecar, ou produit avec + une version/un seuil différents, est régénéré. Il est écrit même quand + aucun offset n'a été appliqué, pour ne pas re-mesurer une tuile déjà + connue comme bien alignée. + """ + payload = { + "version": STRIP_ALIGN_VERSION, + "threshold": STRIP_ALIGN_THRESHOLD, + "offsets": offsets, + } + try: + sidecar = Path(dtm_dir) / f"{basename}_dtm{output_suffix}_stripalign.json" + sidecar.write_text(json.dumps(payload, ensure_ascii=False), + encoding="utf-8") + except Exception as e: + logger.warning(f" Écriture sidecar calage faisceaux impossible: {e}") + + def _strip_lidar_ext(path): """Extract base name from a LAZ/LAS file (mirrors pipeline._file_basename).""" name = Path(path).name @@ -675,7 +813,8 @@ def _interpolate_holes(dtm, downsample=8): return filled, int(holes.sum()) -def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000): +def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000, + strip_offsets=None): """Rasterize the per-cell minimum z (lowest return) of the full point cloud. In complex/forested terrain the ground is under-classified, leaving DTM @@ -689,6 +828,9 @@ def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000): width, height: Output grid dimensions (pixels). bounds: (min_x, min_y, max_x, max_y) the grid covers. chunk_size: Points per streaming chunk. + strip_offsets: Optionnel : dict {point_source_id: offset} issu du + calage des faisceaux, retranché aux Z du nuage complet pour + rester cohérent avec le MNT calé. Returns: (height, width) float32 array of per-cell min z (NaN where no point). @@ -698,12 +840,23 @@ def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000): grid = np.full((height, width), np.nan, dtype=np.float32) rng = [[min_x, max_x], [min_y, max_y]] + lut = None + if strip_offsets: + lut = np.zeros(65536) + for p, off in strip_offsets.items(): + lut[int(p) & 0xFFFF] = off + def process(points): if len(points) == 0: return x = np.asarray(points.x, dtype=np.float64) y = np.asarray(points.y, dtype=np.float64) z = np.asarray(points.z, dtype=np.float64) + if lut is not None: + try: + z = z - lut[np.asarray(points.point_source_id, dtype=np.int64)] + except AttributeError: + pass st = binned_statistic_2d(x, y, z, statistic='min', bins=[width, height], range=rng) # Match the DTM convention: .T then flip Y (north at top). @@ -726,7 +879,7 @@ def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000): def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output_suffix="", source_laz=None, bare_earth=False, - pure=False): + pure=False, strip_align=True): """Create DTM using fast binning method with gap filling. Args: @@ -745,6 +898,10 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, pure: Sans effet (conservé pour compatibilité). Fonctionnement historique rétabli : petits trous comblés par fillnodata, grands trous laissés en nodata (rendus en noir dans les rendus). + strip_align: Si True (défaut), mesure et corrige les écarts verticaux + entre faisceaux de vol (PointSourceId) avant rastérisation ; les + offsets ≥ STRIP_ALIGN_THRESHOLD (0,5 cm) sont consignés dans un + sidecar *_dtm*_stripalign.json. Returns: Path to output DTM GeoTIFF, or None on failure. @@ -774,6 +931,22 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, logger.error(f" ✗ Fichier vide (0 points): {las_file.name}") return None + # Calage vertical des faisceaux de vol avant rastérisation : best-effort, + # en cas d'échec de la mesure on continue non calé (jamais d'abort). + strip_offsets = {} + if strip_align: + try: + strip_offsets = _strip_offsets_for_file(las_file, las) + except Exception as e: + logger.warning(f" Mesure du calage faisceaux impossible ({e}) — MNT non calé") + strip_offsets = {} + if strip_offsets: + logger.info(" Calage faisceaux : " + ", ".join( + f"PSID {p} {off:+.3f} m" for p, off in sorted(strip_offsets.items()))) + else: + logger.debug(" Calage faisceaux : aucun écart >= " + f"{STRIP_ALIGN_THRESHOLD * 100:.1f} cm, rien à corriger") + try: min_x, max_x = float(las.header.min[0]), float(las.header.max[0]) @@ -786,8 +959,17 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, logger.debug(f" Grid: {width}x{height} pixels ({len(las.points):,} points)") logger.info(f" Rasterisation {width}x{height} ({len(las.points):,} points)...") + xs = np.asarray(las.x, dtype=np.float64) + ys = np.asarray(las.y, dtype=np.float64) + zs = np.asarray(las.z, dtype=np.float64) + if strip_offsets: + lut = np.zeros(65536) + for p, off in strip_offsets.items(): + lut[int(p) & 0xFFFF] = off + zs = zs - lut[np.asarray(las.point_source_id, dtype=np.int64)] + stat = binned_statistic_2d( - las.x, las.y, las.z, + xs, ys, zs, statistic='mean', bins=[width, height], range=[[min_x, max_x], [min_y, max_y]] @@ -803,7 +985,8 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, # explicite (--bare-earth). if bare_earth and source_laz is not None: min_grid = _min_return_grid(source_laz, width, height, - (min_x, min_y, max_x, max_y)) + (min_x, min_y, max_x, max_y), + strip_offsets=strip_offsets) # Cellules sans sol mesuré (NaN) : le retour le plus bas devient # la mesure — sinon la comparaison NaN est fausse et la cellule # retombe sur l'interpolation en fin de passe. @@ -844,6 +1027,9 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, ) as dst: dst.write(dtm.astype('float32'), 1) + if strip_align: + _write_strip_align_sidecar(dtm_dir, basename, output_suffix, + strip_offsets) logger.info(f" ✓ DTM créé: {output_tif.name}") return output_tif diff --git a/lidar_pipeline/index.py b/lidar_pipeline/index.py index 7f2f3c6..61bb79d 100644 --- a/lidar_pipeline/index.py +++ b/lidar_pipeline/index.py @@ -3063,9 +3063,11 @@ function clearGhosts() { ghostCells.forEach(g => g.remove()); ghostCells.length // Visualisations choisies dans la barre de génération (noms d'étapes --only). // Sert au preview (détection des tuiles incomplètes) et au lancement : une // tuile existante mais privée d'une de ces visualisations est incluse. +// Rien de sélectionné → toutes les couches du sélecteur (= panneau affiché). function selectedGenViz() { + const all = genViz ? Array.from(genViz.options).map(o => o.value) : []; const sel = genViz ? Array.from(genViz.selectedOptions).map(o => o.value) : []; - return sel.length ? sel : ['aspect']; + return sel.length ? sel : (all.length ? all : ['aspect']); } function exitGenMode() { diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 444521d..6c7e127 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -55,7 +55,7 @@ class FilePrefixFilter(logging.Filter): _file_filter = FilePrefixFilter() from .progress import report_event -from .dtm import classify_ground, create_dtm_fast +from .dtm import classify_ground, create_dtm_fast, STRIP_ALIGN_VERSION, STRIP_ALIGN_THRESHOLD from .visualizations import ( SharedDEM, generate_hillshade, generate_slope, generate_aspect, @@ -109,7 +109,7 @@ VIZ_STEPS = [ class LidarArchaeoPipeline: """Orchestrates the LiDAR archaeological analysis pipeline.""" - def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', ign_classes="sol", force_classify=False, keep_tif=False, bare_earth=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_ids=None, no_index=False, incremental_index=False): + def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', ign_classes="sol", force_classify=False, keep_tif=False, bare_earth=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_ids=None, no_index=False, incremental_index=False, strip_align=True, openness_downsample=None): self.input_dir = Path(input_dir) self.output_dir = Path(output_dir) # Accept single float or comma-separated string for multi-resolution @@ -134,6 +134,8 @@ class LidarArchaeoPipeline: self.gpu_ids = gpu_ids self.no_index = no_index self.incremental_index = incremental_index + self.strip_align = strip_align + self.openness_downsample = openness_downsample self._last_index_rebuild = 0.0 self.temp_dir = self.output_dir / "temp" @@ -368,6 +370,27 @@ class LidarArchaeoPipeline: recorded = self._dtm_method_name(basename, res_suffix) return recorded is None or recorded == self._effective_ground_method() + def _strip_align_matches(self, basename, res_suffix): + """True si le sidecar de calage des faisceaux correspond à la config. + + Un DTM sans sidecar (antérieur au calage) est régénéré pour mesurer + et consigner ses offsets ; un sidecar de version ou de seuil + différents aussi. Calage désactivé : tout DTM porteur d'un sidecar + (donc calé) est régénéré non calé. + """ + sidecar = self.dtm_dir / f"{basename}_dtm{res_suffix}_stripalign.json" + if not self.strip_align: + return not sidecar.exists() + if not sidecar.exists(): + return False + try: + import json + data = json.loads(sidecar.read_text(encoding="utf-8")) + return (data.get("version") == STRIP_ALIGN_VERSION + and abs(float(data.get("threshold", -1)) - STRIP_ALIGN_THRESHOLD) < 1e-9) + except Exception: + return False + def _write_dtm_method(self, basename, res_suffix): """Record the ground classification method used to build a DTM.""" try: @@ -416,7 +439,10 @@ class LidarArchaeoPipeline: res_suffix = self._res_suffix(res) dtm_path = self.dtm_dir / f"{basename}_dtm{res_suffix}.tif" if dtm_path.exists() and not self.force_classify: - if method_matches: + if not self._strip_align_matches(basename, res_suffix): + logger.info(f" DTM{res_suffix} sans calage de faisceaux conforme — régénération (offsets verticaux mesurés et appliqués)") + dtm_path.unlink() + elif method_matches: import rasterio try: with rasterio.open(dtm_path) as src: @@ -468,7 +494,8 @@ class LidarArchaeoPipeline: output_suffix=res_suffix, source_laz=laz_file, bare_earth=self.bare_earth, - pure=pure_ign) + pure=pure_ign, + strip_align=self.strip_align) t_dtm = time.time() - t2 if not dtm_file: logger.error(f" ✗ Échec DTM {res}m/px ({t_dtm:.1f}s)") @@ -483,6 +510,12 @@ class LidarArchaeoPipeline: self._write_dtm_method(basename, res_suffix) # Process each resolution: visualizations + PDF + # Option de calcul (surcharge le défaut du module) appliquée ICI car les + # workers (spawn) réimportent les modules à froid : c'est le seul endroit + # qui s'exécute dans le processus qui fait le calcul. + if self.openness_downsample is not None: + from . import visualizations as _viz_mod + _viz_mod.OPENNESS_DOWNSAMPLE = max(1, int(self.openness_downsample)) all_vis_results = {} for res in self.resolutions: res_suffix = self._res_suffix(res) @@ -580,7 +613,7 @@ class LidarArchaeoPipeline: active_ids = self.gpu_ids if self.gpu_ids else available_gpu_ids() resolutions_str = ','.join(str(r) for r in self.resolutions) future_to_file = { - executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), resolutions_str, self.force, self.ground_method, self.ign_classes, self.force_classify, self.keep_tif, self.bare_earth, self.quality, self.only_viz, self.skip_viz, self.output_format, active_ids[file_idx % len(active_ids)] if active_ids else None): laz_file + executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), resolutions_str, self.force, self.ground_method, self.ign_classes, self.force_classify, self.keep_tif, self.bare_earth, self.quality, self.only_viz, self.skip_viz, self.output_format, active_ids[file_idx % len(active_ids)] if active_ids else None, self.openness_downsample): laz_file for file_idx, laz_file in enumerate(files) } done = 0 @@ -682,7 +715,7 @@ class LidarArchaeoPipeline: logger.warning(f" Note: Impossible de supprimer les fichiers temporaires: {e}") -def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', ign_classes="sol", force_classify=False, keep_tif=False, bare_earth=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_id=None): +def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', ign_classes="sol", force_classify=False, keep_tif=False, bare_earth=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_id=None, openness_downsample=None): """Standalone function for multiprocessing — creates its own pipeline instance. Each worker gets its own temp directory to avoid file conflicts. @@ -709,7 +742,7 @@ def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, fo worker_logger.addHandler(handler) worker_logger.addFilter(_file_filter) - pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, ign_classes=ign_classes, force_classify=force_classify, keep_tif=keep_tif, bare_earth=bare_earth, quality=quality, only_viz=only_viz, skip_viz=skip_viz, output_format=output_format) + pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, ign_classes=ign_classes, force_classify=force_classify, keep_tif=keep_tif, bare_earth=bare_earth, quality=quality, only_viz=only_viz, skip_viz=skip_viz, output_format=output_format, openness_downsample=openness_downsample) basename = _file_basename(laz_file_str) pipeline.temp_dir = pipeline.output_dir / "temp" / basename pipeline.temp_dir.mkdir(exist_ok=True) diff --git a/lidar_pipeline/tests/test_dtm.py b/lidar_pipeline/tests/test_dtm.py index 86db990..549c637 100644 --- a/lidar_pipeline/tests/test_dtm.py +++ b/lidar_pipeline/tests/test_dtm.py @@ -567,3 +567,89 @@ class TestStripLidarExt: from lidar_pipeline.dtm import _strip_lidar_ext from pathlib import Path assert _strip_lidar_ext(Path("/data/input/file.copc.laz")) == "file" + + +class TestStripVerticalOffsets: + """Vertical de-bias of flight-line point sources (strip alignment).""" + + def _synthetic(self, offsets, n=1_500_000, extent=300.0, seed=0): + """Interleaved sources over one tile; each carries a known Z bias.""" + rng = np.random.default_rng(seed) + x = rng.uniform(0, extent, n) + y = rng.uniform(0, extent, n) + z = 100 + 0.02 * x - 0.01 * y + rng.normal(0, 0.01, n) + psid = rng.integers(0, len(offsets), n) + return x, y, z + np.asarray(offsets)[psid], psid.astype(np.uint16) + + def test_recovers_known_offsets(self): + """±5 cm biases are recovered; the aligned source stays untouched.""" + from lidar_pipeline.dtm import _strip_vertical_offsets + x, y, z, psid = self._synthetic((0.0, 0.05, -0.05)) + offs = _strip_vertical_offsets(x, y, z, psid) + assert abs(offs.get(1, 0.0) - 0.05) < 0.01 + assert abs(offs.get(2, 0.0) + 0.05) < 0.01 + assert 0 not in offs # biais ~0 < seuil : pas de correction + + def test_small_bias_below_threshold_ignored(self): + """Biases under the 0.5 cm threshold trigger no correction.""" + from lidar_pipeline.dtm import _strip_vertical_offsets + x, y, z, psid = self._synthetic((0.0, 0.003, -0.003)) + assert _strip_vertical_offsets(x, y, z, psid) == {} + + def test_single_source_returns_empty(self): + """A single point source cannot be compared: no offsets.""" + from lidar_pipeline.dtm import _strip_vertical_offsets + x, y, z, psid = self._synthetic((0.05,)) + assert _strip_vertical_offsets(x, y, z, psid) == {} + + def test_no_shared_cells_returns_empty(self): + """Sources covering disjoint areas (no overlap) are not corrected.""" + from lidar_pipeline.dtm import _strip_vertical_offsets + rng = np.random.default_rng(1) + n = 200_000 + x = np.concatenate([rng.uniform(0, 100, n), rng.uniform(200, 300, n)]) + y = rng.uniform(0, 300, 2 * n) + z = 100 + rng.normal(0, 0.01, 2 * n) + psid = np.concatenate([np.zeros(n, np.uint16), np.ones(n, np.uint16)]) + z[psid == 1] += 0.10 + assert _strip_vertical_offsets(x, y, z, psid) == {} + + +class TestStripAlignSidecar: + def test_sidecar_roundtrip_and_threshold(self, tmp_path): + """Sidecar records version/threshold/offsets and matches config.""" + from lidar_pipeline.dtm import ( + _write_strip_align_sidecar, STRIP_ALIGN_VERSION, STRIP_ALIGN_THRESHOLD) + import json + offsets = {1049: 0.026, 1147: -0.026} + _write_strip_align_sidecar(tmp_path, "TILE", "_r0p2", offsets) + data = json.loads((tmp_path / "TILE_dtm_r0p2_stripalign.json").read_text()) + assert data["version"] == STRIP_ALIGN_VERSION + assert data["threshold"] == STRIP_ALIGN_THRESHOLD + assert data["offsets"] == {"1049": 0.026, "1147": -0.026} # clés JSON en chaînes + + def test_pipeline_match_logic(self, tmp_path): + """_strip_align_matches invalidates legacy DTMs and config changes.""" + from lidar_pipeline.pipeline import LidarArchaeoPipeline + from lidar_pipeline.dtm import _write_strip_align_sidecar + import json + + class P(LidarArchaeoPipeline): + def __init__(self, out, strip_align): + self.output_dir = out + self.dtm_dir = out / "DTM" + self.dtm_dir.mkdir(exist_ok=True) + self.strip_align = strip_align + + p = P(tmp_path, strip_align=True) + # DTM hérité sans sidecar : à régénérer + assert not p._strip_align_matches("TILE", "_r0p2") + # Sidecar conforme : valide + _write_strip_align_sidecar(p.dtm_dir, "TILE", "_r0p2", {}) + assert p._strip_align_matches("TILE", "_r0p2") + # Seuil différent : à régénérer + bad = p.dtm_dir / "TILE_dtm_r0p2_stripalign.json" + bad.write_text(json.dumps({"version": 1, "threshold": 0.02, "offsets": {}})) + 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") diff --git a/lidar_pipeline/tests/test_visualizations.py b/lidar_pipeline/tests/test_visualizations.py index 7e0319d..f7fecf8 100644 --- a/lidar_pipeline/tests/test_visualizations.py +++ b/lidar_pipeline/tests/test_visualizations.py @@ -80,6 +80,52 @@ class TestOpenness: assert result is not None assert result.exists() + def test_downsampled_matches_full_resolution(self, tmp_path, tmp_output_dir): + """Factor 2: same output shape, highly correlated with factor 1. + + MNT dédié à 1 m/px (structures de plusieurs dizaines de pixels, régime + de production 0,2 m/px) : la fixture synthétique partagée à 5 m/px a un + mur large de 2 px, hors régime pour valider une décimation ×2. + """ + import rasterio + from rasterio.transform import from_bounds + from lidar_pipeline import visualizations + from lidar_pipeline.visualizations import generate_openness + + size = 240 + x = np.linspace(0, size, size) + X, Y = np.meshgrid(x, x) + dem = 100.0 + 0.01 * X + 0.005 * Y + dem += 5.0 * np.exp(-((X - 120)**2 + (Y - 120)**2) / (2 * 40**2)) + dist_wall = np.abs((X - 40) * 0.707 + (Y - 60) * 0.707) / np.sqrt(2) + dem += 1.5 * np.exp(-dist_wall**2 / (2 * 8**2)) + dem -= 2.0 * np.exp(-np.abs(X - 200)**2 / (2 * 12**2)) + # Sans bruit blanc : à 1 m/px un bruit σ=5 cm domine l'angle + # d'horizon à courte distance (max le long du rayon) et masque la + # géométrie que ce test cherche à valider (décimation + zoom). + dem_file = tmp_path / "dem_1m.tif" + with rasterio.open(dem_file, 'w', driver='GTiff', height=size, width=size, + count=1, dtype='float32', crs='EPSG:2154', + transform=from_bounds(660000, 6700000, 660240, 6700240, size, size)) as dst: + dst.write(dem.astype('float32'), 1) + + saved = visualizations.OPENNESS_DOWNSAMPLE + try: + visualizations.OPENNESS_DOWNSAMPLE = 1 + r1 = generate_openness(dem_file, "full", tmp_output_dir, 1.0, positive=True) + visualizations.OPENNESS_DOWNSAMPLE = 2 + r2 = generate_openness(dem_file, "dec", tmp_output_dir, 1.0, positive=True) + finally: + visualizations.OPENNESS_DOWNSAMPLE = saved + + with rasterio.open(r1) as s1, rasterio.open(r2) as s2: + a, b = s1.read(1), s2.read(1) + assert a.shape == b.shape + m = ~np.isnan(a) & ~np.isnan(b) + assert m.sum() > 0 + corr = np.corrcoef(a[m], b[m])[0, 1] + assert corr > 0.97, f"corrélation openness décimée/native trop faible : {corr:.3f}" + class TestMSLRM: def test_generates_tif(self, synthetic_dem, tmp_output_dir): diff --git a/lidar_pipeline/tests/test_webapp.py b/lidar_pipeline/tests/test_webapp.py index e4a0c19..ba299c6 100644 --- a/lidar_pipeline/tests/test_webapp.py +++ b/lidar_pipeline/tests/test_webapp.py @@ -558,12 +558,15 @@ def test_start_next_queued_launches_after_run(tmp_path, monkeypatch): webapp._queue[:] = saved_queue -def test_build_command_default_viz_aspect(): - """Sans choix de visualisation, la commande génère uniquement aspect.""" - from lidar_pipeline.webapp import _build_command +def test_build_command_default_viz_panel(): + """Sans choix de visualisation, la commande génère les couches affichées.""" + from lidar_pipeline.webapp import _build_command, _panel_viz_steps cmd = _build_command([(1054, 6882)]) i = cmd.index("--only") - assert cmd[i + 1] == "aspect" + expected = _panel_viz_steps() + n = len(expected) + assert cmd[i + 1:i + 1 + n] == expected + assert "pos_open" in expected # panneau : slope, aspect, positive_openness def test_build_command_custom_viz(): diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 7470d5c..276cf5b 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -628,12 +628,24 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): return None +# Sous-échantillonnage du calcul d'openness : le lancé de rayons est l'étape +# la plus coûteuse du pipeline (rayons de 100 m = 500 pas à 0,2 m/px). L'openness +# étant un champ angulaire lisse (moyenne d'horizons jusqu'à 100 m), on la +# calcule sur une grille décimée par blocs — max local pour l'openness +# positive, min pour la négative, afin de préserver les reliefs qui bornent +# l'angle d'horizon (talus, murs) — puis on la rééchantillonne à la +# résolution demandée. Coût ÷ facteur³ (cellules ÷ facteur², rayons ÷ facteur) +# ; facteur 2 ≈ ×8 plus rapide, rendu quasi identique. 1 = pleine résolution. +OPENNESS_DOWNSAMPLE = 2 + + def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None): """Positive/Negative Openness - multi-radius ray-tracing with std normalization. - Traces rays in 8 directions at 3 radii (25, 50, 100m). - Results are combined with equal weight across radii, then normalized - by standard deviation for cross-tile comparability. + Traces rays in 8 directions at 3 radii (25, 50, 100m) on a block-decimated + grid (cf. OPENNESS_DOWNSAMPLE), then bilinearly resamples the result back + to the requested resolution. Results are combined with equal weight across + radii, then normalized by standard deviation for cross-tile comparability. """ name = "positive_openness" if positive else "negative_openness" gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" @@ -646,8 +658,21 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh _prepare_dem_for_raycast(dem_file, shared, resolution) radii_m = [25, 50, 100] - max_dist = int(max(radii_m) / res) # rayon réel en pixels, non tronqué n_dirs = 8 + full_rows, full_cols = rows, cols + + # Grille décimée par blocs (cf. OPENNESS_DOWNSAMPLE) + factor = max(1, int(OPENNESS_DOWNSAMPLE)) + if factor > 1 and rows >= 2 * factor and cols >= 2 * factor: + r2 = (rows // factor) * factor + c2 = (cols // factor) * factor + blocks = dem[:r2, :c2].reshape(r2 // factor, factor, c2 // factor, factor) + dem = blocks.max(axis=(1, 3)) if positive else blocks.min(axis=(1, 3)) + rows, cols = dem.shape + res = res * factor + logger.info(f" Grille décimée ×{factor} ({rows}×{cols}) — rééchantillonnage final") + + max_dist = int(max(radii_m) / res) # rayon réel en pixels, non tronqué pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m) @@ -660,6 +685,18 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh # Mean across directions and radii (equal weight) — on CPU now openness = np.mean(angles, axis=(0, 1)) openness_result = np.degrees(openness).astype(np.float32) + + # Retour à la grille pleine résolution (bilinéaire ; bord ajusté au plus + # proche voisin si les dimensions n'étaient pas divisibles par le facteur) + if (rows, cols) != (full_rows, full_cols): + from scipy.ndimage import zoom + zoomed = zoom(openness_result, factor, order=1) + if zoomed.shape != (full_rows, full_cols): + zoomed = np.pad(zoomed, + ((0, full_rows - zoomed.shape[0]), + (0, full_cols - zoomed.shape[1])), + mode='edge') + openness_result = zoomed.astype(np.float32) openness_result[nan_mask] = np.nan # Z-score (écarts locaux en sigmas) : unités comparables entre tuiles, diff --git a/lidar_pipeline/webapp.py b/lidar_pipeline/webapp.py index bc1eed6..1cd0c32 100644 --- a/lidar_pipeline/webapp.py +++ b/lidar_pipeline/webapp.py @@ -463,6 +463,20 @@ def _fallback_viz_step_names(): return [KEYWORD_TO_STEP.get(k, k) for k in VIZ_LABELS] +def _panel_viz_steps(): + """Couches réellement affichées par la webapp, en noms d'étapes --only. + + Source unique : PANEL_VIZ (index.py) — le panneau de couches et le + sélecteur de génération sont pilotés par ce registre, les valeurs par + défaut de /api/preview et /api/generate doivent donc l'être aussi, sans + quoi une tuile générée « par défaut » serait incomplète pour le panneau. + """ + from .index import PANEL_VIZ, KEYWORD_TO_STEP + if not PANEL_VIZ: + return ["aspect"] + return [KEYWORD_TO_STEP.get(v, v) for v in PANEL_VIZ] + + def _viz_step_labels(): """Libellés français des étapes de visualisation (clé = nom d'étape). @@ -480,7 +494,8 @@ class PreviewRequest(BaseModel): regenerate: bool = Field(False, description="Inclure les tuiles déjà générées") viz: Optional[list] = Field(None, description="Visualisations demandées, noms d'étapes du " - "pipeline (ex: aspect, wavelet) ; défaut : aspect. " + "pipeline (ex: aspect, wavelet) ; défaut : les " + "couches affichées dans le panneau. " "Une tuile existante est considérée faite " "seulement si elle les possède toutes") all_missing: bool = Field(False, @@ -517,7 +532,7 @@ class GenerateRequest(BaseModel): viz: list = Field(None, description="Visualisations à générer, noms d'étapes du " "pipeline (ex: aspect, wavelet, slope) ; " - "défaut : aspect") + "défaut : les couches affichées dans le panneau") all_missing: bool = Field(False, description="Ignorer la liste de tuiles : traiter toutes " "les dalles LHD présentes dans input/ qui " @@ -1152,7 +1167,7 @@ def preview(req: PreviewRequest): # Webapp légère : le distant connaît input/ et output/ complets return _proxy_api("POST", "/api/preview", json.loads(req.json())) names = _viz_step_names() - viz = [v for v in (req.viz or []) if v] or ["aspect"] + viz = [v for v in (req.viz or []) if v] or _panel_viz_steps() invalid = [v for v in viz if v not in names] if invalid: raise HTTPException( @@ -1217,7 +1232,7 @@ def _build_command(tiles, regenerate=False, ground_class="ign", bare_earth=False cmd = [sys.executable, "-u", "-m", "lidar_pipeline", str(INPUT_DIR), "-o", str(OUTPUT_DIR), "-r", ",".join(str(r) for r in GENERATE_RESOLUTIONS), - "--only", *(viz or ["aspect"]), + "--only", *(viz or _panel_viz_steps()), "--ground-classification", ground_class, "--ign-classes", ign_classes] if bare_earth: @@ -1246,7 +1261,7 @@ def _resolve_request(req): input/ peuvent avoir changé entre la mise en file et le départ du run. """ names = _viz_step_names() - viz = [v for v in (req.viz or []) if v] or ["aspect"] + viz = [v for v in (req.viz or []) if v] or _panel_viz_steps() invalid = [v for v in viz if v not in names] if invalid: raise HTTPException(