diff --git a/docker-compose.worker.yml b/docker-compose.worker.yml index 2821e3c..6be3c12 100644 --- a/docker-compose.worker.yml +++ b/docker-compose.worker.yml @@ -34,6 +34,12 @@ services: # Les générations lancées depuis une webapp distante utilisent le GPU - LIDAR_GPU=1 - LIDAR_WORKERS=auto + # Une ligne par worker : BLAS/OpenMP mono-thread, sinon 12 workers × N + # threads écrasent les 14 cœurs (load 68+ observé pendant les runs) + - OMP_NUM_THREADS=1 + - OPENBLAS_NUM_THREADS=1 + - MKL_NUM_THREADS=1 + - NUMEXPR_NUM_THREADS=1 # Protéger l'API si le réseau n'est pas de confiance : même valeur que # LIDAR_REMOTE_TOKEN sur chaque webapp distante (sinon, laisser commenté) # - LIDAR_API_TOKEN=change-moi diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index 275c852..8f32f28 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -2,7 +2,8 @@ Handles ground classification via PDAL (IGN supplier pre-classification, SMRF or CSF) and DTM rasterisation -using scipy binned_statistic_2d. Gaps without LiDAR data (common in +using scipy binned_statistic_2d (GPU-first via gpu.bin_mean_2d, scipy as +fallback). Gaps without LiDAR data (common in complex/rocky terrain) are filled with a terrain-aware interpolation so the DTM stays continuous. """ @@ -10,6 +11,7 @@ DTM stays continuous. import json import logging import subprocess +import time from pathlib import Path import numpy as np @@ -17,6 +19,8 @@ import rasterio from rasterio.transform import from_bounds from scipy.stats import binned_statistic_2d +from .gpu import bin_mean_2d + logger = logging.getLogger("lidar") # Classes LAS exploitables de la pré-classification LiDAR HD (noms → codes) @@ -1221,7 +1225,10 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, logger.info(" → Génération DTM...") try: + t_read = time.perf_counter() las = laspy.read(str(las_file)) + logger.info(f" Lecture {len(las.points):,} points " + f"({time.perf_counter() - t_read:.1f}s)") except Exception as e: # laspy can't read COPC v1.1 — try PDAL conversion logger.warning(f" laspy: {e}") @@ -1240,6 +1247,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, strip_offsets = {} strip_jitter = {} gps_time = None + t_align = time.perf_counter() if strip_align: try: strip_offsets = _strip_offsets_for_file(las_file, las) @@ -1270,6 +1278,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, f"{len(jt)} fenêtres de {STRIP_JITTER_BIN:g} s)") else: logger.debug(" Gigue intra-faisceau : rien à corriger") + logger.info(f" Calage faisceaux : {time.perf_counter() - t_align:.1f}s") try: @@ -1324,23 +1333,36 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, # non calés par faisceau : bande de contexte, l'image finale est # recadrée sur la dalle avant livraison). if used_edge_buffer > 0 and source_laz is not None: + t_neigh = time.perf_counter() nx, ny, nz = _neighbor_ground_points( source_laz, (min_x, min_y, max_x, max_y), neighbor_classes if neighbor_classes is not None else [2]) + logger.info(f" Voisines (raccord) : {len(nx):,} points " + f"({time.perf_counter() - t_neigh:.1f}s)") if len(nx): xs = np.concatenate([xs, nx]) ys = np.concatenate([ys, ny]) zs = np.concatenate([zs, nz]) - stat = binned_statistic_2d( - xs, ys, zs, - statistic='mean', - bins=[width, height], - range=[[min_x, max_x], [min_y, max_y]] - ) - - dtm = stat.statistic.T - dtm = dtm[::-1, :] # Flip Y so north is at top + t_raster = time.perf_counter() + dtm = bin_mean_2d(xs, ys, zs, width, height, + (min_x, max_x), (min_y, max_y)) + if dtm is not None: + # Sortie déjà en (height, width) : seul le flip Y reste à faire + dtm = dtm[::-1, :] # nord en haut + logger.info(f" ✓ Rasterisation {width}x{height} GPU " + f"({time.perf_counter() - t_raster:.1f}s)") + else: + stat = binned_statistic_2d( + xs, ys, zs, + statistic='mean', + bins=[width, height], + range=[[min_x, max_x], [min_y, max_y]] + ) + dtm = stat.statistic.T + dtm = dtm[::-1, :] # Flip Y so north is at top + logger.info(f" ✓ Rasterisation {width}x{height} CPU " + f"({time.perf_counter() - t_raster:.1f}s)") # Comblement « historique » (fonctionnement d'origine, rétabli) : # seuls les petits trous proches des données sont remplis ; les grands @@ -1357,6 +1379,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, logger.info(f" {nan_count:,} pixels sans données ({nan_pct:.1f}%)") max_gap_pixels = max(1, int(1.0 / resolution)) + t_fill = time.perf_counter() from rasterio.fill import fillnodata valid_mask = ~np.isnan(dtm) dtm_filled = fillnodata(dtm, mask=valid_mask, max_search_distance=max_gap_pixels) @@ -1364,12 +1387,14 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, filled_count = np.count_nonzero(small_gap_mask) if filled_count > 0: dtm = np.where(small_gap_mask, dtm_filled, dtm) - logger.info(f" {filled_count:,} petits trous comblés (< {max_gap_pixels}px)") + logger.info(f" {filled_count:,} petits trous comblés " + f"(< {max_gap_pixels}px, {time.perf_counter() - t_fill:.1f}s)") # Save as GeoTIFF output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif" transform = from_bounds(min_x, min_y, max_x, max_y, width, height) + t_write = time.perf_counter() with rasterio.open( output_tif, 'w', driver='GTiff', height=height, width=width, @@ -1383,6 +1408,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, # Tampon de raccord inscrit dans le fichier : changement de # --edge-buffer ⇒ invalidation automatique du cache DTM. dst.update_tags(**{EDGE_BUFFER_TAG: f"{used_edge_buffer:g}"}) + logger.info(f" Écriture GeoTIFF : {time.perf_counter() - t_write:.1f}s") if strip_align: _write_strip_align_sidecar(dtm_dir, basename, output_suffix, diff --git a/lidar_pipeline/gpu.py b/lidar_pipeline/gpu.py index fe8d0a2..85e3a7d 100644 --- a/lidar_pipeline/gpu.py +++ b/lidar_pipeline/gpu.py @@ -419,6 +419,64 @@ def xp_maximum_filter(arr, footprint=None, size=None): return ndimage.maximum_filter(arr, footprint=footprint, size=size) +# --------------------------------------------------------------------------- +# Rasterisation MNT — moyenne z par cellule (bincount 2D) +# --------------------------------------------------------------------------- + +def _bin_mean_core(lib, xs, ys, zs, width, height, x_range, y_range): + """Moyenne z par cellule d'une grille régulière, backend-agnostique. + + `lib` = numpy ou cupy (mêmes primitives). Sémantique calquée sur + scipy.stats.binned_statistic_2d(statistic='mean', bins=[width, height], + range=[[xmin, xmax], [ymin, ymax]]) : points hors emprise ignorés, valeur + exactement sur le bord droit/haut rangée dans la dernière cellule. + Retourne (height, width) float64, NaN sur les cellules vides. + + float64 obligatoire pour x/y : à des coordonnées Lambert 93 (~1e6 m) la + résolution float32 est ~6 cm — grossière devant un pixel de 0,2 m. + """ + (xmin, xmax), (ymin, ymax) = x_range, y_range + x = lib.asarray(xs, dtype=lib.float64) + y = lib.asarray(ys, dtype=lib.float64) + z = lib.asarray(zs, dtype=lib.float64) + inside = (x >= xmin) & (x <= xmax) & (y >= ymin) & (y <= ymax) + x, y, z = x[inside], y[inside], z[inside] + # floor((v - min) / pas) ; les valeurs intérieures sont >= 0 donc la + # troncature astype == floor. Le clip range le bord droit (et l'arrondi + # float adjacent) dans la dernière cellule, comme scipy. + ix = ((x - xmin) * (float(width) / (xmax - xmin))).astype(lib.int64) + iy = ((y - ymin) * (float(height) / (ymax - ymin))).astype(lib.int64) + ix = lib.clip(ix, 0, width - 1) + iy = lib.clip(iy, 0, height - 1) + n = int(width) * int(height) + idx = iy * int(width) + ix + sums = lib.bincount(idx, weights=z, minlength=n) + counts = lib.bincount(idx, minlength=n) + mean = sums / lib.where(counts == 0, 1, counts) + return lib.where(counts == 0, float("nan"), mean).reshape(int(height), int(width)) + + +def bin_mean_2d(xs, ys, zs, width, height, x_range, y_range): + """Rasterisation « moyenne par cellule » sur GPU (appelant : dtm.create_dtm_fast). + + Retourne un tableau numpy (height, width) float64 (NaN = cellule vide), + ou None si le GPU est indisponible ou échoue (OOM le plus souvent) — + l'appelant retombe alors sur scipy. Un échec ici n'appelle PAS + disable_gpu() : la rastérisation est ponctuelle, les visualisations qui + suivent doivent garder leur accélérateur. + """ + if not _gpu_available(): + return None + try: + result = _bin_mean_core(_cp, xs, ys, zs, width, height, + x_range, y_range) + return to_cpu(result) + except Exception as e: + logger.warning(f"Rasterisation GPU échouée ({e}) — repli scipy") + gpu_cleanup() + return None + + # --------------------------------------------------------------------------- # Misc # --------------------------------------------------------------------------- @@ -458,9 +516,14 @@ def safe_gpu_call(func, *args, **kwargs): try: return func(*args, **kwargs) except Exception as e: - err_msg = str(e) - if _cp is not None and ('CUDA' in err_msg or 'cuda' in err_msg or 'GPU' in err_msg or 'Out of memory' in err_msg): - logger.warning(f"Erreur GPU ({e.__class__.__name__}), retry en CPU...") + # GPU actif : on TOUTE erreur (OOM, types numpy/cupy mêlés après un + # échec de transfert en cours de run, ...) le GPU est désactivé et le + # calcul est retranché en CPU — une panne partielle ne doit pas + # faire échouer la visualisation entière. En mode CPU, on relance + # l'erreur d'origine (déjà en CPU, rien à retrancher). + if _cp is not None: + logger.warning(f"Erreur GPU ({e.__class__.__name__}: {e}), " + f"retry en CPU...") disable_gpu() cpu_args = tuple(to_cpu(a) for a in args) cpu_kwargs = {k: to_cpu(v) for k, v in kwargs.items()} diff --git a/lidar_pipeline/tests/test_gpu.py b/lidar_pipeline/tests/test_gpu.py index b4f22bf..0d65ed5 100644 --- a/lidar_pipeline/tests/test_gpu.py +++ b/lidar_pipeline/tests/test_gpu.py @@ -90,4 +90,113 @@ def test_to_gpu_roundtrip(): assert isinstance(gpu_arr, cp.ndarray) result = to_cpu(gpu_arr) assert isinstance(result, np.ndarray) - np.testing.assert_array_almost_equal(result, [1.0, 2.0, 3.0]) \ No newline at end of file + np.testing.assert_array_almost_equal(result, [1.0, 2.0, 3.0]) + + +# --------------------------------------------------------------------------- +# bin_mean_2d — rasterisation MNT (moyenne z par cellule) +# --------------------------------------------------------------------------- + +def _scipy_reference(xs, ys, zs, width, height, x_range, y_range): + from scipy.stats import binned_statistic_2d + stat = binned_statistic_2d( + xs, ys, zs, statistic='mean', bins=[width, height], + range=[list(x_range), list(y_range)] + ) + return stat.statistic.T # (height, width) + + +def test_bin_mean_core_parity_scipy(): + """Le cœur numpy de bin_mean_2d reproduit binned_statistic_2d (mean).""" + from lidar_pipeline.gpu import _bin_mean_core + rng = np.random.default_rng(42) + x_range, y_range = (1000.0, 1010.0), (6800.0, 6808.0) + w, h = 40, 32 + n = 5000 + xs = rng.uniform(*x_range, n) + ys = rng.uniform(*y_range, n) + zs = rng.uniform(50.0, 150.0, n) + # Cas limites : bord gauche/bas (inclus), bord droit/haut (inclus dans la + # dernière cellule), hors emprise (ignorés), points exactement sur une + # arête intérieure + xs = np.concatenate([xs, [1000.0, 1010.0, 999.9, 1010.1, 1002.5]]) + ys = np.concatenate([ys, [6800.0, 6808.0, 6808.1, 6807.9, 6804.0]]) + zs = np.concatenate([zs, [99.0, 101.0, 777.0, 777.0, 103.0]]) + got = _bin_mean_core(np, xs, ys, zs, w, h, x_range, y_range) + ref = _scipy_reference(xs, ys, zs, w, h, x_range, y_range) + assert got.shape == (h, w) + np.testing.assert_allclose(got, ref, equal_nan=True) + + +def test_bin_mean_core_empty_and_single(): + """Cellules vides → NaN ; un seul point → sa valeur partout où présent.""" + from lidar_pipeline.gpu import _bin_mean_core + x_range, y_range = (0.0, 10.0), (0.0, 10.0) + # Aucun point dans l'emprise + out = _bin_mean_core(np, [-5.0], [-5.0], [1.0], 5, 5, x_range, y_range) + assert out.shape == (5, 5) + assert np.isnan(out).all() + # Un point au centre exact : cellule (2, 2) + out = _bin_mean_core(np, [5.0], [5.0], [7.5], 5, 5, x_range, y_range) + assert out[2, 2] == 7.5 + assert np.isnan(out).sum() == 24 + # Cellule avec plusieurs points : moyenne exacte + out = _bin_mean_core(np, [5.0, 5.1, 5.2], [5.0, 5.0, 5.0], + [10.0, 20.0, 30.0], 5, 5, x_range, y_range) + assert out[2, 2] == 20.0 + + +def test_bin_mean_2d_gpu_or_none(): + """bin_mean_2d : None sans GPU, sinon sortie identique au cœur numpy.""" + from lidar_pipeline.gpu import bin_mean_2d, _bin_mean_core + rng = np.random.default_rng(7) + x_range, y_range = (1000.0, 1200.0), (6800.0, 7000.0) + w, h = 100, 100 + n = 20000 + xs = rng.uniform(*x_range, n) + ys = rng.uniform(*y_range, n) + zs = rng.uniform(50.0, 150.0, n) + got = bin_mean_2d(xs, ys, zs, w, h, x_range, y_range) + ref = _bin_mean_core(np, xs, ys, zs, w, h, x_range, y_range) + if got is None: + assert isinstance(ref, np.ndarray) # repli scipy assuré par l'appelant + else: + assert isinstance(got, np.ndarray) + np.testing.assert_allclose(got, ref, equal_nan=True, rtol=1e-9) + + +def test_safe_gpu_call_retries_on_cpu_for_any_error(monkeypatch): + """GPU actif : TOUTE erreur (pas seulement les messages CUDA) désactive + le GPU et retranche le calcul en CPU. + + Cas réel : après un échec de transfert GPU, des types numpy/cupy mêlés + (« Unsupported type ») ne contenaient aucun + mot-clé CUDA et propageaient l'erreur — la visualisation entière + échouait alors pour rien. + """ + from lidar_pipeline import gpu + + class _FakeCP: # GPU « actif » sans CuPy + pass + + monkeypatch.setattr(gpu, "_cp", _FakeCP()) + monkeypatch.setattr(gpu, "HAS_GPU", True) + calls = [] + + def f(x): + calls.append(x) + if len(calls) == 1: + raise TypeError("Unsupported type ") + return x + 1 + + assert gpu.safe_gpu_call(f, 21) == 22 + assert gpu.HAS_GPU is False # GPU désactivé après l'erreur + assert gpu._cp is None + # Mode CPU : l'erreur se relance telle quelle (rien à retrancher) + def g(x): + raise ValueError("boom") + try: + gpu.safe_gpu_call(g, 1) + raise AssertionError("devait relancer l'erreur") + except ValueError: + pass \ No newline at end of file