Accélérer la rugosité 2.5x via écarts-types par sommes intégrales

lidar_pipeline/visualizations.py

@@ -843,11 +843,55 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
 # Roughness
 # ============================================================
 
+def _integral_sums(x):
+    """Sommes intégrales 2D : S[i,j] = somme de x[0:i, 0:j] (float64).
+
+    Ligne/colonne 0 remplies de zéros — permet la somme d'une fenêtre
+    quelconque par 4 coins, y compris contre le bord (indices 0).
+    """
+    rows, cols = x.shape
+    S = xp.zeros((rows + 1, cols + 1), dtype=np.float64)
+    S[1:, 1:] = xp.cumsum(xp.cumsum(x.astype(np.float64), axis=0), axis=1)
+    return S
+
+
+def _box_std_from_integral(Sx, Sx2, size):
+    """Écart-type local sur fenêtre size×size via sommes intégrales.
+
+    Coût indépendant de la taille de fenêtre (4 accès par pixel) — contre un
+    uniform_filter dont le coût croît avec la fenêtre (75 px à 0,2 m pour
+    l'échelle large). Aux bords, la fenêtre est tronquée et normalisée par le
+    nombre réel d'éléments (les dalles se recouvrent, le bord est sans effet
+    visuel).
+    """
+    rows = Sx.shape[0] - 1
+    cols = Sx.shape[1] - 1
+    r = size // 2
+    iy0 = xp.maximum(xp.arange(rows) - r, 0)
+    iy1 = xp.minimum(xp.arange(rows) + r + 1, rows)
+    ix0 = xp.maximum(xp.arange(cols) - r, 0)
+    ix1 = xp.minimum(xp.arange(cols) + r + 1, cols)
+    # Nombre d'éléments réels de la fenêtre (tronquée aux bords)
+    counts = ((iy1 - iy0)[:, None] * (ix1 - ix0)[None, :]).astype(np.float64)
+
+    def box_sum(S):
+        return (S[iy1][:, ix1] - S[iy0][:, ix1]
+                - S[iy1][:, ix0] + S[iy0][:, ix0])
+
+    mean = box_sum(Sx) / counts
+    mean_sq = box_sum(Sx2) / counts
+    return xp.sqrt(xp.maximum(mean_sq - mean * mean, 0))
+
+
 def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
     """Surface roughness - multi-scale standard deviation (GPU-accelerated).
 
     Combines fine (3m) and broad (15m) roughness for better detection
-    of archaeological features at multiple scales.
+    of archaeological features at multiple scales. Les écarts-types locaux
+    sont calculés par sommes intégrales : deux cumsum partagés entre les
+    deux échelles, extraction par 4 coins — coût constant quelle que soit
+    la fenêtre (un uniform_filter coûte proportionnellement à sa taille,
+    75 px à 0,2 m pour l'échelle large).
     """
     gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
     logger.info(f"  → Rugosité de surface{gpu_tag}...")
@@ -860,38 +904,33 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
             crs = shared.crs
             dem_np = shared.dem_np
             nan_mask = shared.nan_mask
+            if _gpu_mod.HAS_GPU:
+                filled = shared.filled_gpu
+            else:
+                filled = shared.filled
         else:
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
+            filled, _ = _fill_nans(dem_np)
+            if _gpu_mod.HAS_GPU:
+                filled = to_gpu(filled)
+
+        # Sommes intégrales partagées par les deux échelles (X et X²)
+        Sx = _integral_sums(filled)
+        Sx2 = _integral_sums(filled.astype(np.float64) ** 2)
 
-        # Fine roughness (3m window)
         fine_size = max(3, int(3 / resolution))
         if fine_size % 2 == 0:
             fine_size += 1
-
-        if shared:
-            fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size)
-            fine_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
-            fine_mean_sq[shared.nan_mask] = np.nan
-        else:
-            fine_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=fine_size)
-            fine_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
-        roughness_fine = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0))
-        roughness_fine[nan_mask] = np.nan
-
-        # Broad roughness (15m window)
         broad_size = max(3, int(15 / resolution))
         if broad_size % 2 == 0:
             broad_size += 1
 
-        if shared:
-            broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size)
-            broad_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
-            broad_mean_sq[shared.nan_mask] = np.nan
-        else:
-            broad_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=broad_size)
-            broad_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
-        roughness_broad = np.sqrt(np.maximum(broad_mean_sq - broad_mean * broad_mean, 0))
+        roughness_fine = to_cpu(_box_std_from_integral(Sx, Sx2, fine_size))
+        roughness_broad = to_cpu(_box_std_from_integral(Sx, Sx2, broad_size))
+        del Sx, Sx2
+        gpu_cleanup()
+        roughness_fine[nan_mask] = np.nan
         roughness_broad[nan_mask] = np.nan
 
         # Std normalization per scale then weighted combination