Add anomaly mask: automatic threshold detection across all viz layers

lidar_pipeline/pipeline.py

@@ -61,9 +61,10 @@ from .visualizations import (
     generate_openness,
     generate_mslrm, generate_sailore,
     generate_roughness, generate_wavelet,
-    generate_svf, generate_aniso_open,
+       generate_svf, generate_aniso_open,
-    generate_flow_accumulation,
+        generate_flow_accumulation,
-)
+        generate_anomaly_mask,
+    )
 from .gpu import gpu_cleanup, num_gpus, restrict_gpus, safe_gpu_call
 from .ign import generate_ign_overlay
 from .rendering import tif_to_png
@@ -84,6 +85,7 @@ VIZ_STEPS = [
     ('roughness',  generate_roughness),
     ('wavelet',    generate_wavelet),
     ('flow_acc',   generate_flow_accumulation),
+    ('anomaly',    generate_anomaly_mask),
     ('ortho',      lambda d, b, v, r: generate_ign_overlay(
                     d, b, v, r,
                     layer='ORTHOIMAGERY.ORTHOPHOTOS',

lidar_pipeline/rendering.py

@@ -168,6 +168,14 @@ COLORMAPS = {
         'vmin_mode': 'fixed', 'vmin_val': 0,
         'vmax_mode': 'percentile', 'vmax_pct': 98,
     },
+    'anomaly': {
+        'cmap': 'YlOrRd',
+        'title': 'Carte d\'Anomalies (détection automatique)',
+        'legend': 'Score d\'anomalie composite (0–1)\nRouge = Haute anomalie (structures suspectes)\nJaune = Anomalie modérée\nBlanc = Aucun signal (terrain naturel)\n\nSeuil auto: pixels > 2σ de la moyenne locale\nCombiné: MSRM, SVF, Ondelette, Openness, Rugosité',
+        'description': 'Détection automatique — cible à vérifier sur terrain',
+        'vmin_mode': 'fixed', 'vmin_val': 0,
+        'vmax_mode': 'fixed', 'vmax_val': 1,
+    },
 }
 
 # RGB entries (ortho/topo) are handled specially

lidar_pipeline/visualizations.py

@@ -1112,3 +1112,146 @@ def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None):
     except Exception as e:
         logger.error(f"  ✗ Erreur openness anisotropique: {e}", exc_info=True)
         return None
+
+
+# ============================================================
+# Anomaly Mask — automatic threshold detection
+# ============================================================
+
+def generate_anomaly_mask(dem_file, basename, vis_dir, resolution, shared=None, n_sigma=2.0):
+    """Composite anomaly mask — automatic threshold detection (GPU if available).
+
+    Reads pre-computed visualization layers (MSRM, SVF, Wavelet, Openness Neg,
+    Roughness), normalizes each to z-scores, and combines them into a composite
+    anomaly score. Pixels beyond `n_sigma` standard deviations of the local mean
+    are flagged as suspicious.
+
+    The output is a continuous score (0–1) where:
+      - 0 = no anomaly (flat/natural terrain)
+      - 1 = high anomaly (potential archaeological structure)
+
+    This mask is directly usable in GIS for polygon extraction and field survey
+    planning.
+
+    Args:
+        n_sigma: Number of standard deviations for the anomaly threshold.
+                 Lower = more sensitive (more false positives).
+                 Default 2.0 (good balance for archaeological detection).
+    """
+    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
+    logger.info(f"  → Détection automatique d'anomalies (seuil {n_sigma}σ){gpu_tag}...")
+    t0 = time.time()
+    output = vis_dir / f"{basename}_anomaly.tif"
+
+    try:
+        if shared:
+            transform = shared.transform
+            crs = shared.crs
+            nan_mask = shared.nan_mask
+        else:
+            dem_np, transform, crs = _read_dem(dem_file)
+            nan_mask = np.isnan(dem_np)
+
+        rows, cols = nan_mask.shape
+
+        # Collect available visualization layers from disk
+        # Each is loaded, normalized to z-score, and contributes to the composite
+        layer_configs = [
+            # (filename_pattern, weight)
+            ("mslrm", 2.5),       # Multi-scale relief — strongest signal
+            ("negative_openness", 1.8),  # Fossés, dolines
+            ("roughness", 1.5),   # Surface irregularity
+            ("wavelet", 1.5),     # Circular + linear structures
+            ("svf", 1.3),         # Sky-view depressions
+            ("flow_acc", 1.2),    # Drainage channels / ditches
+            ("aniso_open", 1.0),  # Anisotropic structures
+            ("positive_openness", 0.8),  # Surélevations
+        ]
+
+        layers = []
+        for pattern, weight in layer_configs:
+            layer_path = vis_dir / f"{basename}_{pattern}.tif"
+            if not layer_path.exists():
+                layer_path = vis_dir / f"{basename}_negative_openness.tif" if "neg" in pattern else None
+                if layer_path is None or not layer_path.exists():
+                    continue
+            try:
+                with rasterio.open(layer_path) as src:
+                    data = src.read(1).astype(np.float64)
+                # Z-score normalization
+                valid = data[~nan_mask]
+                if len(valid) == 0:
+                    continue
+                mean_val = np.nanmean(valid)
+                std_val = max(np.nanstd(valid), 0.01)
+                zscore = (data - mean_val) / std_val
+                zscore[nan_mask] = 0.0
+                layers.append((zscore, weight))
+            except Exception as e:
+                logger.debug(f"    Couche {pattern} non disponible: {e}")
+                continue
+
+        if not layers:
+            # Fallback: use MSRM + roughness computed on the fly
+            logger.info("    Aucune couche trouvée — calcul MSRM + rugosité en direct...")
+            dem_np_safe = shared.dem_np if shared else dem_np
+
+            # Quick MSRM (single scale 10m for speed)
+            sigma_px = max(5, 10.0 / resolution)
+            local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px) if shared else \
+                _filter_nanaware(dem_np_safe, xp_gaussian_filter, sigma=sigma_px)
+            quick_mslrm = dem_np_safe - local_mean
+            quick_mslrm[nan_mask] = np.nan
+            valid = quick_mslrm[~nan_mask]
+            std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
+            quick_mslrm = quick_mslrm / std_val
+            layers.append((quick_mslrm, 2.5))
+
+            # Quick roughness
+            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)
+            else:
+                fine_mean = _filter_nanaware(dem_np_safe.astype(np.float64), xp_uniform_filter, size=fine_size)
+                fine_mean_sq = _filter_nanaware(dem_np_safe.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
+            roughness = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0))
+            roughness[nan_mask] = np.nan
+            valid_r = roughness[~nan_mask]
+            std_val_r = max(np.nanstd(valid_r), 0.01) if len(valid_r) > 0 else 0.01
+            roughness = roughness / std_val_r
+            layers.append((roughness, 1.5))
+
+        # Weighted RMS combination (unsigned — all deviations are suspicious)
+        combined = np.zeros((rows, cols), dtype=np.float64)
+        total_weight = 0.0
+        for layer_data, weight in layers:
+            combined += (layer_data ** 2) * weight
+            total_weight += weight
+
+        if total_weight > 0:
+            combined = np.sqrt(combined / total_weight)
+
+        # Apply n_sigma threshold: pixels below n_sigma are suppressed
+        combined = np.maximum(combined - n_sigma, 0.0)
+
+        # Rescale to 0–1 for visualization (percentile-based stretch)
+        valid_combined = combined[~nan_mask]
+        if len(valid_combined) > 0 and np.nanmax(valid_combined) > 0:
+            p99 = np.percentile(valid_combined, 99)
+            if p99 > 0:
+                combined = np.clip(combined / p99, 0, 1)
+
+        combined[nan_mask] = np.nan
+
+        _save_tif(output, combined.astype(np.float32), transform, crs)
+        n_anomaly = int(np.sum(combined > 0.1)) if np.any(combined > 0) else 0
+        pct_anomaly = n_anomaly / max(np.sum(~nan_mask), 1) * 100
+        logger.info(f"  ✓ Détection anomalies terminée ({time.time()-t0:.1f}s) — "
+                     f"{pct_anomaly:.1f}% de la zone ({n_anomaly} px) au-delà de {n_sigma}σ")
+        return output
+    except Exception as e:
+        logger.error(f"  ✗ Erreur détection anomalies: {e}", exc_info=True)
+        return None