"""Terrain visualization functions for LiDAR archaeological analysis.

Each function takes (dem_file, basename, vis_dir, resolution) as explicit
parameters and returns the path to the output GeoTIFF file, or None on error.

When a SharedDEM object is provided via the `shared` parameter, pre-computed
data (gradient, NaN mask, LRM) is reused across visualizations to avoid
redundant I/O and computation.
"""

import logging
import math
import time
import warnings
from pathlib import Path

import numpy as np
import rasterio
from .gpu import HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, xp_minimum_filter, xp_maximum_filter, gpu_cleanup
from . import gpu as _gpu_mod

logger = logging.getLogger("lidar")

# CuPy module reference — lazily imported on first GPU use.
# If disable_gpu() is called at runtime, HAS_GPU becomes False
# and xp delegates to numpy instead.
_cp = None


class _XPProxy:
    """Proxy that delegates array operations to cupy or numpy.

    Checks HAS_GPU on every attribute access so that disable_gpu()
    (called on CUDA errors) takes effect immediately, without needing
    to change every call site in visualizations.py.
    """

    def __getattr__(self, name):
        global _cp
        from . import gpu as _gpu_mod
        if _gpu_mod.HAS_GPU:
            if _cp is None:
                try:
                    import cupy
                    _cp = cupy
                except ImportError:
                    pass
            if _cp is not None:
                return getattr(_cp, name)
        return getattr(np, name)


xp = _XPProxy()


class SharedDEM:
    """Pre-computed DEM data shared across all visualizations.

    Reads the DEM once and lazily computes on first access:
    - NaN mask and filled DEM (avoids 20+ calls to _fill_nans)
    - Gradient components (shared by hillshade, slope)
    - LRM at 15m kernel (shared by mslrm + sailore)

    Attributes are computed lazily on first access to avoid computing
    data that is never used (e.g. LRM when only hillshade needs generation).
    """

    def __init__(self, dem_file, resolution):
        dem_np, transform, crs = _read_dem(dem_file)
        self.dem_file = dem_file
        self.resolution = resolution
        self.transform = transform
        self.crs = crs
        self.nan_mask = np.isnan(dem_np)
        self.dem_np = dem_np.astype(np.float32)

        # Lazy caches — computed on first access
        self._filled = None
        self._gradient = None  # (dy, dx, slope_rad, slope_deg)
        self._lrm_15 = None

        # GPU lazy caches
        self._filled_gpu = None
        self._dem_gpu = None

    @property
    def filled(self):
        """Filled DEM (NaN interpolated) — computed lazily."""
        if self._filled is None:
            logger.debug("    → Calcul filled DEM (interpolation NaN)...")
            self._filled, _ = _fill_nans(self.dem_np)
        return self._filled

    @property
    def dy(self):
        self._ensure_gradient()
        return self._gradient[0]

    @property
    def dx(self):
        self._ensure_gradient()
        return self._gradient[1]

    @property
    def slope_rad(self):
        self._ensure_gradient()
        return self._gradient[2]

    @property
    def slope_deg(self):
        self._ensure_gradient()
        return self._gradient[3]

    @property
    def lrm_15(self):
        """LRM at 15m kernel — computed lazily."""
        if self._lrm_15 is None:
            logger.debug("    → Calcul LRM 15m...")
            sigma_15 = 15.0 / self.resolution
            local_mean_15 = _filter_nanaware_from_filled(self, xp_gaussian_filter, sigma=sigma_15)
            self._lrm_15 = self.dem_np - local_mean_15
            self._lrm_15[self.nan_mask] = np.nan
        return self._lrm_15

    def _ensure_gradient(self):
        """Compute gradient components lazily on first access."""
        if self._gradient is None:
            logger.debug("    → Calcul gradient...")
            dy = np.gradient(self.filled, self.resolution, axis=0)
            dx = np.gradient(self.filled, self.resolution, axis=1)
            slope_rad = np.arctan(np.sqrt(dx**2 + dy**2))
            slope_deg = np.degrees(slope_rad)
            self._gradient = (dy, dx, slope_rad, slope_deg)

    @property
    def filled_gpu(self):
        """Lazy GPU copy of the filled DEM."""
        if self._filled_gpu is None and _gpu_mod.HAS_GPU:
            self._filled_gpu = to_gpu(self.filled)
        return self._filled_gpu

    @property
    def dem_gpu(self):
        """Lazy GPU copy of the DEM."""
        if self._dem_gpu is None and _gpu_mod.HAS_GPU:
            self._dem_gpu = to_gpu(self.dem_np)
        return self._dem_gpu


def _filter_nanaware_from_filled(shared, filter_func, *args, **kwargs):
    """Apply filter on pre-filled DEM data (skips expensive _fill_nans).

    Uses the SharedDEM.filled array directly, then restores NaN mask.
    If GPU is available, reuses the lazy GPU copy to avoid redundant transfers.
    """
    if _gpu_mod.HAS_GPU:
        filled_gpu = shared.filled_gpu
    else:
        filled_gpu = None

    if filled_gpu is not None:
        result_gpu = filter_func(filled_gpu, *args, **kwargs)
        result = to_cpu(result_gpu)
        gpu_cleanup()
    else:
        result = filter_func(shared.filled, *args, **kwargs)
    result[shared.nan_mask] = np.nan
    return result


def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None, nan_mask=None):
    """Helper to save a 2D or 3D array as GeoTIFF.

    Args:
        nan_mask: Optional boolean mask (True=NaN) to apply before saving.
                  Restores NaN zones in gradient-derived products that were
                  computed on the filled DEM.
    """
    if nan_mask is not None:
        data = np.array(data, dtype=dtype, copy=True)
        data[nan_mask] = np.nan

    # Auto-detect nodata for float types with NaN
    if nodata is None and dtype.startswith('float') and np.any(np.isnan(data)):
        nodata = float('nan')

    if data.ndim == 2:
        height, width = data.shape
        with rasterio.open(
            output_path, 'w', driver='GTiff',
            height=height, width=width, count=count,
            dtype=dtype, crs=crs, transform=transform,
            compress='deflate', nodata=nodata
        ) as dst:
            dst.write(data.astype(dtype), 1)
    elif data.ndim == 3:
        bands, height, width = data.shape
        with rasterio.open(
            output_path, 'w', driver='GTiff',
            height=height, width=width, count=bands,
            dtype=dtype, crs=crs, transform=transform,
            compress='deflate', nodata=nodata
        ) as dst:
            for i in range(bands):
                dst.write(data[i].astype(dtype), i + 1)


def _read_dem(dem_file):
    """Read DEM file and return (data, transform, crs)."""
    with rasterio.open(dem_file) as src:
        return src.read(1), src.transform, src.crs


def _fill_nans(arr):
    """Fill NaN values using nearest-neighbor interpolation.

    Returns (filled_array, nan_mask) so the caller can restore NaN after filtering.
    """
    from scipy.interpolate import NearestNDInterpolator
    nan_mask = np.isnan(arr)
    if not np.any(nan_mask):
        return arr, nan_mask
    valid = ~nan_mask
    y_coords, x_coords = np.where(valid)
    if len(y_coords) == 0:
        return np.zeros_like(arr), nan_mask
    z_values = arr[valid]
    interp = NearestNDInterpolator(
        np.column_stack((y_coords, x_coords)), z_values
    )
    y_missing, x_missing = np.where(nan_mask)
    filled = arr.copy()
    filled[y_missing, x_missing] = interp(y_missing, x_missing)
    return filled, nan_mask


def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
    """Apply a filter to an array while preserving NaN zones.

    1. Fill NaN with nearest-neighbor interpolation
    2. Apply the filter
    3. Restore original NaN mask on the result

    Args:
        arr: Input array (numpy or cupy).
        filter_func: Function that takes (array, *args, **kwargs) and returns filtered array.
        use_gpu: If True, apply filter on GPU (send filled array to GPU first).
    Returns:
        Filtered array with original NaN positions preserved.
    """
    is_gpu_arr = _gpu_mod.HAS_GPU and _cp is not None and isinstance(arr, _cp.ndarray)
    arr_np = to_cpu(arr) if is_gpu_arr else arr

    filled, nan_mask = _fill_nans(arr_np)

    if use_gpu and _gpu_mod.HAS_GPU:
        filled_gpu = to_gpu(filled)
        result_gpu = filter_func(filled_gpu, *args, **kwargs)
        result = to_cpu(result_gpu)
        gpu_cleanup()
    else:
        result = filter_func(filled, *args, **kwargs)

    result[nan_mask] = np.nan
    return result


# ============================================================
# Shared ray-tracing core
# ============================================================

def _prepare_dem_for_raycast(dem_file, shared, resolution):
    """Load DEM and prepare padded array for ray-tracing.

    Returns (dem_gpu_or_cpu, dem_np, rows, cols, res, nan_mask,
             transform, crs, padded) ready for ray-tracing.
    """
    if shared:
        dem_np = shared.dem_np
        nan_mask = shared.nan_mask
        transform = shared.transform
        crs = shared.crs
        dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled
    else:
        dem_np, transform, crs = _read_dem(dem_file)
        nan_mask = np.isnan(dem_np)
        filled, _ = _fill_nans(dem_np)
        dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled
    res = resolution
    rows, cols = dem_np.shape
    return dem, dem_np, rows, cols, res, nan_mask, transform, crs


def _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m=None):
    """Core ray-tracing: compute max zenith/nadir angles per direction and radius.

    For each pixel, in each direction, traces rays outward up to max_dist steps,
    recording the max upward angle (positive openness) and max downward angle
    (negative openness) reached at each radius checkpoint.

    Padding is done on CPU (numpy) to avoid GPU memory pressure and
    pre-compiled kernel mismatches (CUDA_ERROR_NO_BINARY_FOR_GPU on sm_89).
    The padded array is transferred to GPU once, then each direction is
    processed and results are streamed back to CPU.

    Args:
        dem: GPU or CPU filled DEM array (rows, cols).
        rows, cols: dimensions.
        res: resolution in m/px.
        n_dirs: number of directions.
        max_dist: max ray steps.
        radii_m: list of radii in meters to record checkpoints.
                 If None, records only at max_dist.

    Returns:
        pos_angles: array of shape (n_dirs, n_radii, rows, cols) — max zenith angles
        neg_angles: array of shape (n_dirs, n_radii, rows, cols) — max nadir angles
    """
    angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
    dx_dir = np.cos(angles)
    dy_dir = np.sin(angles)

    if radii_m is not None:
        radii_steps = [min(int(r / res), max_dist) for r in radii_m]
        n_radii = len(radii_m)
    else:
        radii_steps = [max_dist]
        n_radii = 1

    # Pad on CPU (numpy) — avoids GPU memory pressure and
    # pre-compiled kernel issues (NO_BINARY_FOR_GPU on sm_89).
    dem_np = to_cpu(dem)
    padded_np = np.pad(dem_np, max_dist, mode='constant', constant_values=np.nan)

    # Transfer padded DEM to GPU for computation
    padded = to_gpu(padded_np)
    # Free the CPU copy — we don't need it anymore
    del dem_np, padded_np

    # Process one direction at a time to limit GPU memory.
    # Store results as flat CPU arrays — transfer back to GPU at the end.
    pos_results = [None] * n_dirs
    neg_results = [None] * n_dirs

    for d_idx in range(n_dirs):
        ddx, ddy = dx_dir[d_idx], dy_dir[d_idx]

        # Pre-compute valid steps for this direction
        valid_steps = []
        for step in range(1, max_dist + 1):
            px = int(round(ddx * step))
            py = int(round(ddy * step))
            dist_m = math.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
            if dist_m < res * 0.5:
                continue
            valid_steps.append((step, px, py, dist_m))

        # Running max angles per radius
        running_pos = xp.zeros((n_radii, rows, cols))
        running_neg = xp.zeros((n_radii, rows, cols))
        # Track which radius checkpoints have been passed
        radii_remaining = set(range(n_radii))

        for step, px, py, dist_m in valid_steps:
            # Slice from padded array, subtract original dem
            view = padded[max_dist + py:max_dist + py + rows,
                         max_dist + px:max_dist + px + cols]
            elev_diff = view - dem
            del view  # free slice reference

            # Positive: angle to terrain above viewer
            pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m)
            # Negative: angle to terrain below viewer
            neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m)
            del elev_diff  # free intermediate

            # Update running max for all radius checkpoints still active
            for r_idx in radii_remaining:
                pos_angle_safe = xp.nan_to_num(pos_angle, nan=0)
                neg_angle_safe = xp.nan_to_num(neg_angle, nan=0)
                running_pos[r_idx] = xp.where(xp.isnan(pos_angle), running_pos[r_idx],
                                              xp.maximum(running_pos[r_idx], pos_angle_safe))
                running_neg[r_idx] = xp.where(xp.isnan(neg_angle), running_neg[r_idx],
                                              xp.maximum(running_neg[r_idx], neg_angle_safe))

            # Check which radii have been passed
            new_remaining = set()
            for r_idx in radii_remaining:
                if step < radii_steps[r_idx]:
                    new_remaining.add(r_idx)
            radii_remaining = new_remaining
            if not radii_remaining:
                break

        # Store results on CPU, free GPU memory before next direction
        pos_results[d_idx] = to_cpu(running_pos)
        neg_results[d_idx] = to_cpu(running_neg)
        del running_pos, running_neg
        gpu_cleanup()

    # Free the large padded array
    del padded
    gpu_cleanup()

    # Reassemble into final arrays (on CPU to avoid GPU memory pressure)
    pos_angles = np.array(pos_results)
    neg_angles = np.array(neg_results)

    return pos_angles, neg_angles


# ============================================================
# Core terrain visualizations
# ============================================================

def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
    """Generate multi-directional hillshade with contrast enhancement — GPU if available.

    Combines 8-direction hillshade with slope shading for balanced illumination.
    Applies percentile normalization and gamma correction to restore
    contrast lost by averaging multiple azimuths.
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Hillshade multidirectionnel{gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_hillshade_multi.tif"
    used_gpu = _gpu_mod.HAS_GPU

    try:
        if shared:
            transform = shared.transform
            crs = shared.crs
            dem = to_gpu(shared.dem_np)
            dy = to_gpu(shared.dy) if _gpu_mod.HAS_GPU else shared.dy
            dx = to_gpu(shared.dx) if _gpu_mod.HAS_GPU else shared.dx
            slope = to_gpu(shared.slope_rad) if _gpu_mod.HAS_GPU else shared.slope_rad
            aspect = xp.arctan2(dy, dx)
            sin_slope = xp.sin(slope)
            cos_slope = xp.cos(slope)
        else:
            dem_np, transform, crs = _read_dem(dem_file)
            dem = to_gpu(dem_np)
            dy, dx = xp.gradient(dem)
            slope = xp.arctan(xp.sqrt(dx**2 + dy**2))
            aspect = xp.arctan2(dy, dx)
            sin_slope = xp.sin(slope)
            cos_slope = xp.cos(slope)

        # 8 azimuths for balanced illumination (eliminates directional bias)
        azimuts = [0, 45, 90, 135, 180, 225, 270, 315]
        altitude = 35  # Higher altitude for better micro-relief detection
        hillshades = []

        alt_rad = xp.radians(xp.array(altitude))
        sin_alt = xp.sin(alt_rad)
        cos_alt = xp.cos(alt_rad)

        for az in azimuts:
            az_rad = xp.radians(xp.array(az))
            hs = sin_alt * sin_slope + cos_alt * cos_slope * xp.cos(az_rad - aspect)
            hillshades.append(xp.clip(hs, 0, 1))

        combined_hillshade = xp.mean(xp.array(hillshades), axis=0)
        slope_shaded = cos_slope
        combined = 0.7 * combined_hillshade + 0.3 * slope_shaded

        # Contrast enhancement: percentile stretch + gamma
        combined_np = to_cpu(combined)
        nan_mask = shared.nan_mask if shared else np.isnan(dem_np)
        valid = combined_np[~nan_mask]
        if len(valid) > 0:
            p2, p98 = np.percentile(valid, 2), np.percentile(valid, 98)
            if p98 - p2 > 0.01:
                combined_np = np.clip((combined_np - p2) / (p98 - p2), 0, 1)
        # Gamma correction to enhance shadows
        gamma = 0.8
        combined_np = np.power(combined_np, gamma)

        _save_tif(output, combined_np.astype(np.float32), transform, crs, nan_mask=nan_mask)
        logger.info(f"  ✓ Hillshade terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur hillshade: {e}", exc_info=True)
        return None


def generate_slope(dem_file, basename, vis_dir, resolution, shared=None):
    """Generate slope map (degrees) — GPU if available."""
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Pente (Slope){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_slope.tif"

    try:
        if shared:
            transform = shared.transform
            crs = shared.crs
            slope = shared.slope_deg
            nan_mask = shared.nan_mask
            if _gpu_mod.HAS_GPU:
                slope = to_gpu(slope)
        else:
            dem_np, transform, crs = _read_dem(dem_file)
            dem = to_gpu(dem_np)
            dy, dx = xp.gradient(dem)
            slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) * 180 / xp.pi
            nan_mask = np.isnan(dem_np)
        _save_tif(output, to_cpu(slope) if _gpu_mod.HAS_GPU else slope, transform, crs, nan_mask=nan_mask)
        logger.info(f"  ✓ Pente terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur slope: {e}", exc_info=True)
        return None


# ============================================================
# GPU-accelerated visualizations
# ============================================================

def generate_svf(dem_file, basename, vis_dir, resolution, shared=None):
    """Sky-View Factor - ray-tracing on 16 azimuths, multi-radius (GPU if available).

    Traces rays in 16 directions at 3 radii (25, 50, 100m) and combines
    with weights favoring medium range for archaeological feature detection.
    SVF = (1/N) * sum(cos²(horizon_angle)). Valleys/crevices have low SVF,
    ridges/peaks have high SVF.
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Sky-View Factor (ray-tracing multi-rayon){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_svf.tif"

    try:
        dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
            _prepare_dem_for_raycast(dem_file, shared, resolution)

        radii_m = [25, 50, 100]
        radius_weights = [0.3, 0.4, 0.3]  # Medium range weighted more
        max_dist = min(int(100 / res), 300)
        n_dirs = 16

        pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)

        # pos/neg are now numpy arrays (CPU) — combine on CPU
        svf_combined = np.zeros((rows, cols), dtype=np.float32)
        for r_idx in range(len(radii_m)):
            horizon = np.maximum(pos_angles[:, r_idx], neg_angles[:, r_idx])
            svf_r = np.mean(np.cos(horizon) ** 2, axis=0)
            svf_combined += svf_r * radius_weights[r_idx]

        svf_np = svf_combined
        svf_np[nan_mask] = np.nan
        _save_tif(output, svf_np, transform, crs)
        logger.info(f"  ✓ SVF terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur SVF: {e}", exc_info=True)
        return None


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.
    """
    name = "positive_openness" if positive else "negative_openness"
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → {name.replace('_', ' ').title()} (ray-tracing multi-rayon){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_{name}.tif"

    try:
        dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
            _prepare_dem_for_raycast(dem_file, shared, resolution)

        radii_m = [25, 50, 100]
        max_dist = min(int(100 / res), 300)
        n_dirs = 8

        pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)

        # Select positive or negative
        if positive:
            angles = pos_angles
        else:
            angles = neg_angles

        # 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)
        openness_result[nan_mask] = np.nan

        # Std normalization for cross-tile comparability
        valid = openness_result[~nan_mask]
        if len(valid) > 0:
            std_val = max(np.nanstd(valid), 0.01)
            openness_result = openness_result / std_val

        _save_tif(output, openness_result, transform, crs)
        logger.info(f"  ✓ {name} terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur openness: {e}", exc_info=True)
        return None


def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
    """Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available).

    Scales adapt to resolution. Std normalization per scale.
    Weighted combination favoring archaeologically relevant scales (5-25m).
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Multi-Scale Relief Model (MSRM){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_mslrm.tif"

    try:
        if shared:
            transform = shared.transform
            crs = shared.crs
            dem_np = shared.dem_np
            nan_mask = shared.nan_mask
        else:
            dem_np, transform, crs = _read_dem(dem_file)
            nan_mask = np.isnan(dem_np)

        # Adaptive scales: finer at higher resolution
        min_scale = max(2.0, resolution * 4)
        candidate_scales = [2, 5, 10, 20, 50, 100, 200]
        sigmas = [s for s in candidate_scales if s >= min_scale]

        # Archaeological weights: favor 5-25m range (ditches, enclosures, tumulus)
        scale_weights = {
            2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6, 200: 0.4,
        }
        weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])

        logger.info(f"    MSRM échelles: {sigmas}m")
        lrm_stack = []

        for sigma in sigmas:
            sigma_px = sigma / resolution
            if shared:
                local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px)
            else:
                local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
            lrm = dem_np - local_mean
            lrm[nan_mask] = np.nan
            # Std normalization: x / std — preserves sign and contrast better than z-score
            valid_lrm = lrm[~nan_mask]
            lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
            lrm = lrm / lrm_std
            lrm_stack.append(lrm.astype(np.float32))

        # Weighted combination — preserve sign for RdBu_r colormap
        # Positive = elevated (red), Negative = depression (blue)
        lrm_array = np.array(lrm_stack)
        weights_3d = weights[:, np.newaxis, np.newaxis]
        with np.errstate(invalid='ignore', divide='ignore'):
            with warnings.catch_warnings():
                warnings.filterwarnings('ignore', message='Mean of empty slice')
                # Signed RMS: magnitude from RMS, sign from weighted mean
                signed_mean = np.nansum(lrm_array * weights_3d, axis=0) / np.sum(weights)
                rms_magnitude = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
                mslrm = np.sign(signed_mean) * rms_magnitude
        mslrm[nan_mask] = np.nan
        _save_tif(output, mslrm.astype(np.float32), transform, crs)
        logger.info(f"  ✓ MSRM terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur MSRM: {e}", exc_info=True)
        return None


# ============================================================
# SAILORE
# ============================================================

def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
    """SAILORE - Self-Adaptive Improved Local Relief Model (GPU if available).

    Kernel size adapts to local slope: flat areas get larger kernels,
    steep areas get smaller kernels. Scales adapt to resolution.
    Reuses shared.lrm_15 when available to avoid recomputation.
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → SAILORE (LRM adaptatif){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_sailore.tif"

    try:
        if shared:
            transform = shared.transform
            crs = shared.crs
            dem_np = shared.dem_np
            nan_mask = shared.nan_mask
            slope_deg = shared.slope_deg
        else:
            dem_np, transform, crs = _read_dem(dem_file)
            nan_mask = np.isnan(dem_np)
            gy, gx = np.gradient(dem_np, resolution)
            slope = np.arctan(np.sqrt(gx**2 + gy**2))
            slope_deg = np.degrees(slope)
            slope_deg[nan_mask] = np.nan

        # Fixed physical scales (independent of resolution)
        sigma_min_m = 2.0   # 2m — fine detail
        sigma_max_m = 25.0  # 25m — broad relief
        sigma_min = sigma_min_m / resolution
        sigma_max = sigma_max_m / resolution
        slope_norm = np.clip(slope_deg / 30.0, 0, 1)

        # LRM fine (2m) — always compute
        if shared:
            lrm_fine = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_min)
        else:
            lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min)
        lrm_fine[nan_mask] = np.nan

        # LRM medium (13.5m) — reuse shared.lrm_15 (σ=15m) when available
        sigma_mid = (sigma_min + sigma_max) / 2
        if shared and abs(15.0 / resolution - sigma_mid) < 2.0 / resolution:
            # shared.lrm_15 is close enough to medium scale
            lrm_medium = shared.lrm_15.copy()
        else:
            if shared:
                lrm_medium = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_mid)
            else:
                lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_mid)
        lrm_medium[nan_mask] = np.nan

        # LRM coarse (25m) — always compute
        if shared:
            lrm_coarse = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_max)
        else:
            lrm_coarse = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_max)
        lrm_coarse[nan_mask] = np.nan

        w_fine = slope_norm
        w_medium = 1 - 2 * np.abs(slope_norm - 0.5)
        w_coarse = 1 - slope_norm
        w_total = w_fine + w_medium + w_coarse
        w_total[w_total == 0] = 1

        sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total
        sailore[nan_mask] = np.nan

        _save_tif(output, sailore.astype(np.float32), transform, crs)
        logger.info(f"  ✓ SAILORE terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur SAILORE: {e}", exc_info=True)
        return None


# ============================================================
# Roughness
# ============================================================

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.
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Rugosité de surface{gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_roughness.tif"

    try:
        if shared:
            transform = shared.transform
            crs = shared.crs
            dem_np = shared.dem_np
            nan_mask = shared.nan_mask
        else:
            dem_np, transform, crs = _read_dem(dem_file)
            nan_mask = np.isnan(dem_np)

        # 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_broad[nan_mask] = np.nan

        # Std normalization per scale then weighted combination
        fine_valid = roughness_fine[~nan_mask]
        broad_valid = roughness_broad[~nan_mask]
        fine_std = max(np.nanstd(fine_valid), 0.01) if len(fine_valid) > 0 else 0.01
        broad_std = max(np.nanstd(broad_valid), 0.01) if len(broad_valid) > 0 else 0.01

        roughness = 0.7 * roughness_fine / fine_std + 0.3 * roughness_broad / broad_std
        roughness[nan_mask] = np.nan

        _save_tif(output, roughness, transform, crs)
        logger.info(f"  ✓ Rugosité terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur rugosité: {e}", exc_info=True)
        return None


# ============================================================
# Wavelet (Mexican Hat + Directional Gabor)
# ============================================================

def _gabor_kernel_2d(size, sigma, wavelength, theta):
    """Create a 2D Gabor kernel.

    Args:
        size: kernel size (odd integer)
        sigma: standard deviation
        wavelength: wavelength of sinusoid
        theta: orientation angle in radians (0 = horizontal)
    """
    center = size // 2
    y, x = np.ogrid[-center:center+1, -center:center+1]
    # Rotate coordinates
    x_theta = x * np.cos(theta) + y * np.sin(theta)
    y_theta = -x * np.sin(theta) + y * np.cos(theta)
    sigma_sq = 2 * sigma * sigma
    kernel = np.exp(-(x_theta**2 + y_theta**2) / sigma_sq)
    kernel *= np.cos(2 * np.pi * x_theta / wavelength)
    return kernel


def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
    """Multi-scale wavelet analysis: Mexican Hat + Directional Gabor (GPU if available).

    Mexican Hat (radial): detects circular features (tumulus, enclos ronds).
    Gabor (directional): detects linear features (chemins, murs, fossés).
    4 Gabor orientations (0°, 45°, 90°, 135°) at key archaeological scales.
    Both combined with RMS for orientation-invariant detection.
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Ondelette Mexican Hat + Gabor directionnelle{gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_wavelet.tif"

    try:
        if shared:
            transform = shared.transform
            crs = shared.crs
            dem_np = shared.dem_np
            nan_mask = shared.nan_mask
            filled = shared.filled.astype(np.float64)
        else:
            dem_np, transform, crs = _read_dem(dem_file)
            nan_mask = np.isnan(dem_np)
            filled, _ = _fill_nans(dem_np.astype(np.float64))

        # --- Mexican Hat scales ---
        min_scale = max(resolution * 2, 1.0)
        candidate_scales = [0.5, 1, 2, 5, 10, 20, 50, 100]
        mex_scales = [s for s in candidate_scales if s >= min_scale]

        mex_weights_map = {
            0.5: 0.6, 1.0: 0.8, 2.0: 1.5, 5.0: 2.0,
            10.0: 1.8, 20.0: 1.5, 50.0: 1.0, 100.0: 0.6,
        }
        mex_weights = np.array([mex_weights_map.get(s, 1.0) for s in mex_scales])

        logger.info(f"    Échelles CWT: {mex_scales}m (résolution {resolution}m/px)")

        from scipy.ndimage import gaussian_laplace, convolve

        wavelet_stack = []

        # Mexican Hat (radial) — multi-scale
        for scale_m in mex_scales:
            sigma_px = scale_m / resolution
            if _gpu_mod.HAS_GPU:
                try:
                    from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace
                    response = -gpu_gaussian_laplace(to_gpu(filled), sigma=sigma_px)
                    response = to_cpu(response)
                except Exception:
                    response = -gaussian_laplace(filled, sigma=sigma_px)
            else:
                response = -gaussian_laplace(filled, sigma=sigma_px)
            response[nan_mask] = np.nan
            valid = response[~nan_mask]
            std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
            response = response / std_val
            wavelet_stack.append(response)

        # Gabor (directional) — 4 orientations at 3 key scales
        gabor_scales_m = [5, 10, 20]  # Key archaeological scales for linear features
        gabor_orientations = [0, np.pi/4, np.pi/2, 3*np.pi/4]  # 0°, 45°, 90°, 135°
        gabor_scale_weights = {5: 2.0, 10: 1.8, 20: 1.5}

        for scale_m in gabor_scales_m:
            sigma_px = max(2, scale_m / resolution / 3)  # Sigma relative to wavelength
            wavelength_px = max(3, scale_m / resolution)
            kernel_size = max(5, int(wavelength_px * 2.5))
            if kernel_size % 2 == 0:
                kernel_size += 1

            for theta in gabor_orientations:
                kernel = _gabor_kernel_2d(kernel_size, sigma_px, wavelength_px, theta)
                # Normalize kernel
                kernel = kernel / max(np.abs(kernel).max(), 1e-10)
                response = convolve(filled, kernel, mode='constant', cval=0)
                response[nan_mask] = np.nan
                valid = response[~nan_mask]
                std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
                response = response / std_val
                wavelet_stack.append(response)

        # Build weights: Mexican Hat weights + Gabor weights
        gabor_weights_list = []
        for scale_m in gabor_scales_m:
            w = gabor_scale_weights.get(scale_m, 1.0)
            gabor_weights_list.extend([w] * len(gabor_orientations))
        gabor_weights = np.array(gabor_weights_list)
        all_weights = np.concatenate([mex_weights, gabor_weights])

        # Weighted RMS combination
        stack = np.array(wavelet_stack)
        weights_3d = all_weights[:, np.newaxis, np.newaxis]
        with np.errstate(invalid='ignore', divide='ignore'):
            with warnings.catch_warnings():
                warnings.filterwarnings('ignore', message='Mean of empty slice')
                combined = np.sqrt(np.nansum((stack ** 2) * weights_3d, axis=0) / np.sum(all_weights))
        combined[nan_mask] = np.nan

        _save_tif(output, combined.astype(np.float32), transform, crs)
        logger.info(f"  ✓ Ondelette terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur ondelette: {e}", exc_info=True)
        return None


# ============================================================
# Flow Accumulation
# ============================================================

def generate_flow_accumulation(dem_file, basename, vis_dir, resolution, shared=None):
    """Flow Accumulation — priority-flood sink filling + D8 accumulation (GPU if available).

    Detects channels, ditches, and drainage paths by computing how many
    upstream cells flow through each cell. Archaeological ditches and
    natural drainage features both accumulate high flow values.

    Uses log10 transformation for visualization.
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Accumulation d'écoulement (flow accumulation){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_flow_acc.tif"

    try:
        if shared:
            transform = shared.transform
            crs = shared.crs
            dem_np = shared.dem_np
            nan_mask = shared.nan_mask
            filled = shared.filled
        else:
            dem_np, transform, crs = _read_dem(dem_file)
            nan_mask = np.isnan(dem_np)
            filled, _ = _fill_nans(dem_np)

        # Priority-flood sink filling (Wang & Liu 2006, O(n log n))
        from heapq import heappush, heappop
        rows, cols = filled.shape
        dem_filled = filled.copy()

        # Use heap for priority-flood
        visited = np.zeros((rows, cols), dtype=bool)
        heap = []
        # Seed with all border cells
        for x in range(cols):
            heappush(heap, (dem_filled[0, x], 0, x))
            heappush(heap, (dem_filled[rows-1, x], rows-1, x))
            visited[0, x] = True
            visited[rows-1, x] = True
        for y in range(1, rows-1):
            heappush(heap, (dem_filled[y, 0], y, 0))
            heappush(heap, (dem_filled[y, cols-1], y, cols-1))
            visited[y, 0] = True
            visited[y, cols-1] = True

        while heap:
            elev, cy, cx = heappop(heap)
            for dy, dx in [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]:
                ny, nx = cy + dy, cx + dx
                if 0 <= ny < rows and 0 <= nx < cols and not visited[ny, nx]:
                    if dem_filled[ny, nx] > elev:
                        dem_filled[ny, nx] = elev
                    heappush(heap, (dem_filled[ny, nx], ny, nx))
                    visited[ny, nx] = True

        logger.info(f"    ✓ Sink filling terminé ({time.time()-t0:.1f}s)")

        # D8 flow direction and accumulation
        flow_acc = np.zeros((rows, cols), dtype=np.int64)
        flow_dir = np.zeros((rows, cols), dtype=np.int8) - 1

        # D8 neighbors (ordered by angle)
        neighbors = [
            (-1,  0), (-1,  1), ( 0,  1), ( 1,  1),
            ( 1,  0), ( 1, -1), ( 0, -1), (-1, -1)
        ]

        # Process cells in ascending elevation order for correct accumulation
        flat_idx = np.argsort(dem_filled.ravel())
        flow_acc_flat = flow_acc.ravel()

        for idx in flat_idx:
            y = idx // cols
            x = idx % cols

            # Find steepest downslope neighbor
            max_slope = -np.inf
            best_dir = -1
            for di, (dy, dx) in enumerate(neighbors):
                ny, nx = y + dy, x + dx
                if 0 <= ny < rows and 0 <= nx < cols:
                    slope = (dem_filled[y, x] - dem_filled[ny, nx])
                    dist = math.sqrt(dy**2 + dx**2)
                    slope_per_m = slope / dist
                    if slope_per_m > max_slope:
                        max_slope = slope_per_m
                        best_dir = di

            flow_dir[y, x] = best_dir
            flow_acc[y, x] = 1  # Count self

        # Accumulate flow (upstream to downstream)
        # Process in reverse elevation order (highest first)
        for idx in reversed(flat_idx):
            y = idx // cols
            x = idx % cols
            d = flow_dir[y, x]
            if d >= 0:
                dy, dx = neighbors[d]
                ny, nx = y + dy, x + dx
                if 0 <= ny < rows and 0 <= nx < cols:
                    flow_acc[ny, nx] += flow_acc[y, x]

        logger.info(f"    ✓ D8 accumulation terminé ({time.time()-t0:.1f}s)")

        # Log transform for visualization
        flow_result = np.log10(np.maximum(flow_acc.astype(np.float32), 1.0))
        flow_result[nan_mask] = np.nan

        _save_tif(output, flow_result, transform, crs)
        logger.info(f"  ✓ Flow accumulation terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur flow accumulation: {e}", exc_info=True)
        return None


# ============================================================
# Anisotropic Openness
# ============================================================

def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None):
    """Anisotropic Openness - weighted directional openness, multi-radius (GPU if available).

    Computes positive and negative openness with anisotropic weighting:
    NW/SE directions weighted more heavily to enhance detection of structures
    aligned NE-SW (common in French archaeological sites: villas, enclosures).

    Multi-radius (25, 50, 100m) with equal weight, results std-normalized
    for cross-tile comparability.
    """
    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
    logger.info(f"  → Openness Anisotropique (multi-rayon){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_aniso_open.tif"

    try:
        dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
            _prepare_dem_for_raycast(dem_file, shared, resolution)

        n_dirs = 8
        # Anisotropic weights: emphasize NW-SE and NE-SW directions
        weights = np.array([1.0, 1.5, 1.0, 1.5, 1.0, 1.5, 1.0, 1.5])

        radii_m = [25, 50, 100]
        max_dist = min(int(100 / res), 300)

        pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)

        # Weighted combination across directions and radii
        weight_total = np.sum(weights)
        n_radii = len(radii_m)

        pos_combined = np.zeros((rows, cols), dtype=np.float64)
        neg_combined = np.zeros((rows, cols), dtype=np.float64)

        for r_idx in range(n_radii):
            for d_idx in range(n_dirs):
                w = weights[d_idx]
                pos_combined += pos_angles[d_idx, r_idx] * w / (n_radii * weight_total)
                neg_combined += neg_angles[d_idx, r_idx] * w / (n_radii * weight_total)

        aniso_result = np.degrees(pos_combined - neg_combined).astype(np.float32)
        aniso_result[nan_mask] = np.nan

        # Std normalization for cross-tile comparability
        valid = aniso_result[~nan_mask]
        if len(valid) > 0:
            std_val = max(np.nanstd(valid), 0.01)
            aniso_result = aniso_result / std_val

        _save_tif(output, aniso_result, transform, crs)
        logger.info(f"  ✓ Openness anisotropique terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
        return output
    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