Rasterisation GPU du DTM (bin_mean_2d) + GPU fallback renforcé et workers mono-thread

Rastérisation des points sol via gpu.bin_mean_2d (bincount 2D CuPy, sémantique
identique au repli scipy binned_statistic_2d) : ~x7 plus rapide sur cette
étape, logs de durée par phase dans create_dtm_fast. safe_gpu_call retente en
CPU sur toute erreur GPU (pas seulement les messages CUDA) : un transfert
échoué en cours de run mélangeait types numpy/cupy et faisait échouer la
visualisation entière. Workers en mono-thread BLAS/OpenMP sur la machine de
traitement (plus de saturation des cœurs pendant les runs).

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

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])
 
+        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,

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)
+        # GPU actif : on TOUTE erreur (OOM, types numpy/cupy mêlés après un
-        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):
+        # échec de transfert en cours de run, ...) le GPU est désactivé et le
-            logger.warning(f"Erreur GPU ({e.__class__.__name__}), retry en CPU...")
+        # 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()}

lidar_pipeline/tests/test_gpu.py

@@ -91,3 +91,112 @@ def test_to_gpu_roundtrip():
     result = to_cpu(gpu_arr)
     assert isinstance(result, np.ndarray)
     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 <class 'numpy.ndarray'> ») 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 <class 'numpy.ndarray'>")
+        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