lidar_rendu/lidar_pipeline/tests/test_visualizations.py

"""Tests for visualization functions.

Each test creates a small synthetic DEM and runs a visualization function,
checking that it produces a valid output file.
"""

import numpy as np
import pytest
from pathlib import Path


# --- Core terrain visualizations (no GPU required) ---

class TestHillshade:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_hillshade
        result = generate_hillshade(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()
        assert result.suffix == ".tif"

    def test_output_values_valid(self, synthetic_dem, tmp_output_dir):
        import rasterio
        from lidar_pipeline.visualizations import generate_hillshade
        result = generate_hillshade(synthetic_dem, "test", tmp_output_dir, 5.0)
        with rasterio.open(result) as src:
            data = src.read(1)
            assert data.shape[0] > 0
            assert np.nanmin(data) >= 0
            assert np.nanmax(data) <= 1


class TestSlope:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_slope
        result = generate_slope(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()

    def test_slope_values_degrees(self, synthetic_dem, tmp_output_dir):
        import rasterio
        from lidar_pipeline.visualizations import generate_slope
        result = generate_slope(synthetic_dem, "test", tmp_output_dir, 5.0)
        with rasterio.open(result) as src:
            data = src.read(1)
            assert np.nanmin(data) >= 0
            assert np.nanmax(data) <= 90


# --- GPU-accelerated visualizations ---

class TestSVF:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_svf
        result = generate_svf(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()

    def test_svf_values_0_1(self, synthetic_dem, tmp_output_dir):
        import rasterio
        from lidar_pipeline.visualizations import generate_svf
        result = generate_svf(synthetic_dem, "test", tmp_output_dir, 5.0)
        with rasterio.open(result) as src:
            data = src.read(1)
            valid = data[~np.isnan(data)]
            assert np.nanmin(valid) >= 0
            assert np.nanmax(valid) <= 1


class TestOpenness:
    def test_positive_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_openness
        result = generate_openness(synthetic_dem, "test", tmp_output_dir, 5.0, positive=True)
        assert result is not None
        assert result.exists()

    def test_negative_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_openness
        result = generate_openness(synthetic_dem, "test", tmp_output_dir, 5.0, positive=False)
        assert result is not None
        assert result.exists()

    def test_downsampled_matches_full_resolution(self, tmp_path, tmp_output_dir):
        """Factor 2: same output shape, highly correlated with factor 1.

        Dedicated 1 m/px DTM (structures spanning tens of pixels, as in the
        0.2 m/px production regime): the shared 5 m/px synthetic fixture has
        a wall only 2 px wide, outside the regime needed to validate a ×2
        decimation.
        """
        import rasterio
        from rasterio.transform import from_bounds
        from lidar_pipeline import visualizations
        from lidar_pipeline.visualizations import generate_openness

        size = 240
        x = np.linspace(0, size, size)
        X, Y = np.meshgrid(x, x)
        dem = 100.0 + 0.01 * X + 0.005 * Y
        dem += 5.0 * np.exp(-((X - 120)**2 + (Y - 120)**2) / (2 * 40**2))
        dist_wall = np.abs((X - 40) * 0.707 + (Y - 60) * 0.707) / np.sqrt(2)
        dem += 1.5 * np.exp(-dist_wall**2 / (2 * 8**2))
        dem -= 2.0 * np.exp(-np.abs(X - 200)**2 / (2 * 12**2))
        # No white noise: at 1 m/px a σ=5 cm noise dominates the short-range
        # horizon angle (max along the ray) and hides the geometry this test
        # validates (decimation + zoom).
        dem_file = tmp_path / "dem_1m.tif"
        with rasterio.open(dem_file, 'w', driver='GTiff', height=size, width=size,
                           count=1, dtype='float32', crs='EPSG:2154',
                           transform=from_bounds(660000, 6700000, 660240, 6700240, size, size)) as dst:
            dst.write(dem.astype('float32'), 1)

        saved = visualizations.OPENNESS_DOWNSAMPLE
        try:
            visualizations.OPENNESS_DOWNSAMPLE = 1
            r1 = generate_openness(dem_file, "full", tmp_output_dir, 1.0, positive=True)
            visualizations.OPENNESS_DOWNSAMPLE = 2
            r2 = generate_openness(dem_file, "dec", tmp_output_dir, 1.0, positive=True)
        finally:
            visualizations.OPENNESS_DOWNSAMPLE = saved

        with rasterio.open(r1) as s1, rasterio.open(r2) as s2:
            a, b = s1.read(1), s2.read(1)
        assert a.shape == b.shape
        m = ~np.isnan(a) & ~np.isnan(b)
        assert m.sum() > 0
        corr = np.corrcoef(a[m], b[m])[0, 1]
        assert corr > 0.97, f"decimated/native openness correlation too low: {corr:.3f}"

    @pytest.mark.parametrize("positive", [True, False])
    def test_scale_independent_of_rest_of_tile(self, tmp_path, tmp_output_dir, positive):
        """Same local relief = same value, whatever the rest of the tile.

        Two DTMs identical on their western half; one also carries a hill to
        the east, more than 100 m (max radius) from the compared area. A
        per-tile normalization (z-score) shifted the whole scale and made
        mosaics non-seamless; with the frozen references, the western half
        must come out identical.
        """
        import rasterio
        from rasterio.transform import from_bounds
        from lidar_pipeline.visualizations import generate_openness

        size = 300
        x = np.arange(size, dtype=float)
        X, Y = np.meshgrid(x, x)
        flat = 100.0 + 2.0 * np.exp(-((X - 60)**2 + (Y - 150)**2) / (2 * 15**2))
        hill = flat + 40.0 * np.exp(-((X - 260)**2 + (Y - 150)**2) / (2 * 20**2))
        results = []
        for name, dem in (("flat", flat), ("hill", hill)):
            f = tmp_path / f"{name}.tif"
            with rasterio.open(f, 'w', driver='GTiff', height=size, width=size,
                               count=1, dtype='float32', crs='EPSG:2154',
                               transform=from_bounds(660000, 6700000, 660300, 6700300,
                                                     size, size)) as dst:
                dst.write(dem.astype('float32'), 1)
            out = generate_openness(f, name, tmp_output_dir, 1.0, positive=positive)
            with rasterio.open(out) as src:
                results.append(src.read(1))
        # Tolerance: interpolation residual of the decimated grid (≈0.01);
        # a per-tile z-score would shift the whole area by several tenths.
        west = np.s_[:, :100]
        np.testing.assert_allclose(results[0][west], results[1][west], atol=0.05)


def _write_dem(path, dem, res=1.0):
    import rasterio
    from rasterio.transform import from_origin
    with rasterio.open(path, 'w', driver='GTiff', height=dem.shape[0], width=dem.shape[1],
                       count=1, dtype='float32', crs='EPSG:2154',
                       transform=from_origin(660000, 6700300, res, res)) as dst:
        dst.write(dem.astype('float32'), 1)
    return path


class TestReliefOriente:
    def test_generates_rgb_uint8(self, synthetic_dem, tmp_output_dir):
        import rasterio
        from lidar_pipeline.visualizations import generate_relief_oriente
        result = generate_relief_oriente(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None and result.name == "test_relief_oriente.tif"
        with rasterio.open(result) as src, rasterio.open(synthetic_dem) as dem:
            assert src.count == 3 and src.dtypes[0] == 'uint8'
            assert (src.height, src.width) == (dem.height, dem.width)
            rgb = src.read()
        assert rgb.std() > 5, "uniform image: neither relief nor orientation rendered"

    def test_nodata_gets_fixed_color(self, tmp_path, tmp_output_dir):
        import rasterio
        from lidar_pipeline.visualizations import generate_relief_oriente, RELIEF_NODATA_RGB
        x = np.arange(120, dtype=float)
        dem = 100 + 0.05 * x[None, :] + 0.02 * x[:, None]
        dem[10:20, 10:20] = np.nan
        out = generate_relief_oriente(_write_dem(tmp_path / "d.tif", dem), "t", tmp_output_dir, 1.0)
        with rasterio.open(out) as src:
            rgb = src.read()
        assert tuple(rgb[:, 15, 15]) == RELIEF_NODATA_RGB

    def test_horizon_kernels_agree(self):
        """numba and numpy (and the CUDA kernel, same code) give the same angle."""
        pytest.importorskip("numba")
        from lidar_pipeline.visualizations import (
            _horizon_rays, _mean_horizon_numba, _mean_horizon_numpy)
        rng = np.random.default_rng(0)
        dem = rng.normal(0, 0.3, (60, 70)).astype(np.float32)
        dem[20:30, 30:40] += 2.0
        offs, dist, cps = _horizon_rays(0.8, 16, (5, 10, 20))
        a = _mean_horizon_numba(dem, offs, dist, cps)
        b = _mean_horizon_numpy(dem, offs, dist, cps)
        np.testing.assert_allclose(a, b, atol=1e-5)

    def test_matches_legacy_ray_trace_geometry(self):
        """Same ray geometry as _ray_trace_horizons (8 directions)."""
        from lidar_pipeline.visualizations import (
            _horizon_rays, _mean_horizon_numpy, _ray_trace_horizons)
        rng = np.random.default_rng(1)
        dem = rng.normal(0, 0.5, (50, 50)).astype(np.float32)
        radii = (5, 10, 20)
        pos, _ = _ray_trace_horizons(dem, 50, 50, 1.0, 8, 20, list(radii))
        legacy = np.mean(pos, axis=(0, 1))
        offs, dist, cps = _horizon_rays(1.0, 8, radii)
        np.testing.assert_allclose(_mean_horizon_numpy(dem, offs, dist, cps), legacy, atol=1e-5)

    def test_seamless_independent_of_rest_of_tile(self, tmp_path, tmp_output_dir):
        """Same local relief = same colors, whatever the rest of the tile
        (support: 20 m radius + 4σ = 40 m detrend window, ~60 m in all, well
        under the 100 m edge buffer)."""
        import rasterio
        from lidar_pipeline.visualizations import generate_relief_oriente
        size = 300
        X, Y = np.meshgrid(np.arange(size, dtype=float), np.arange(size, dtype=float))
        flat = 100.0 + 1.5 * np.exp(-((X - 60)**2 + (Y - 150)**2) / (2 * 6**2))
        hill = flat + 40.0 * np.exp(-((X - 260)**2 + (Y - 150)**2) / (2 * 10**2))
        imgs = []
        for name, dem in (("flat", flat), ("hill", hill)):
            out = generate_relief_oriente(_write_dem(tmp_path / f"{name}.tif", dem),
                                          name, tmp_output_dir, 1.0)
            with rasterio.open(out) as src:
                imgs.append(src.read().astype(int))
        diff = np.abs(imgs[0][:, :, :150] - imgs[1][:, :, :150])
        assert diff.max() <= 1, f"color difference far from the hill: {diff.max()}"

    def test_numba_colorize_matches_vectorized(self):
        """Fused CPU kernel = vectorized path (CuPy/numpy) up to rounding."""
        pytest.importorskip("numba")
        from lidar_pipeline.gpu import xp_zoom
        from lidar_pipeline.visualizations import (
            _relief_colorize_numba, _relief_colorize_xp, _relief_lut, _pad_to)
        rng = np.random.default_rng(2)
        open_c = rng.uniform(0.3, 15, (30, 25)).astype(np.float32)
        dx = rng.normal(0, 0.3, (120, 100)).astype(np.float32)
        dy = rng.normal(0, 0.3, (120, 100)).astype(np.float32)
        lut = _relief_lut()
        a = _relief_colorize_numba(open_c, 4, dx, dy, lut).astype(int)
        b = _relief_colorize_xp(np, _pad_to(np, xp_zoom(open_c, 4), 120, 100), dx, dy, lut).astype(int)
        diff = np.abs(a - b).max(axis=2)
        assert (diff > 3).mean() < 0.01, f"{(diff > 3).mean():.3%} pixels diverge"

    def test_lut_lightness_is_monotonic_and_hue_neutral(self):
        """Lightness increases with L*; at fixed L*, all hues have the same
        perceived luminance (no false relief caused by color)."""
        from lidar_pipeline.visualizations import _relief_lut
        lut = _relief_lut().astype(float) / 255
        lin = np.where(lut <= 0.04045, lut / 12.92, ((lut + 0.055) / 1.055) ** 2.4)
        Y = lin @ np.array([0.2126, 0.7152, 0.0722])          # (256, 360)
        assert np.all(np.diff(Y.mean(axis=1)) >= -1e-6)
        Lstar = 116 * np.cbrt(Y[100:180]) - 16                # L* ≈ 40–70
        # Residual spread: sRGB gamut clipping of a few hues at C* 60
        assert (Lstar.max(axis=1) - Lstar.min(axis=1)).max() < 6


class TestMSLRM:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_mslrm
        result = generate_mslrm(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()


class TestSAILORE:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_sailore
        result = generate_sailore(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()


class TestRoughness:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_roughness
        result = generate_roughness(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()

    def test_roughness_non_negative(self, synthetic_dem, tmp_output_dir):
        import rasterio
        from lidar_pipeline.visualizations import generate_roughness
        result = generate_roughness(synthetic_dem, "test", tmp_output_dir, 5.0)
        with rasterio.open(result) as src:
            data = src.read(1)
            # Standard deviation is always >= 0
            assert np.nanmin(data) >= 0


class TestWavelet:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_wavelet
        result = generate_wavelet(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()

    def test_output_median_centered(self, synthetic_dem, tmp_output_dir):
        """Robust rescaling: median ≈ 1 whatever the terrain.

        This is what lets a single fixed color stretch give homogeneous
        colors across tiles (see knots in rendering.py).
        """
        import rasterio
        from lidar_pipeline.visualizations import generate_wavelet
        result = generate_wavelet(synthetic_dem, "test", tmp_output_dir, 5.0)
        with rasterio.open(result) as src:
            data = src.read(1)
        valid = data[np.isfinite(data)]
        assert abs(np.median(valid) - 1.0) < 0.05

    def test_ditch_on_hilltop_not_amplified(self, tmp_path, tmp_output_dir):
        """A ditch on a hilltop does not stand out more than on flat ground.

        Detrending (Gaussian local mean ~35 m) must neutralize the
        topographic position:
        - the index peak on the hilltop ditch stays comparable to the same
          ditch on flat ground;
        - the hilltop background (no structure) stays comparable to the flat
          background.

        Noisy synthetic DTM (σ=8 cm): without detrending, the hilltop
        background comes out ~1.7× the flat background (yellow hilltops on
        the map).
        """
        import rasterio
        from rasterio.transform import from_bounds
        from lidar_pipeline.visualizations import generate_wavelet

        size = 600
        res = 1.0
        x = np.arange(size) * res
        y = np.arange(size) * res
        X, Y = np.meshgrid(x, y)
        dem = 100.0 + 0.01 * X

        # Realistic hill (sigma 80 m, 25 m high) on the left + sensor noise
        dem += 25.0 * np.exp(-((X - 150)**2 + (Y - 300)**2) / (2 * 80**2))
        rng = np.random.default_rng(42)
        dem += rng.normal(0, 0.08, dem.shape)

        # Two identical ditches: one on the hilltop, the other on flat ground
        for xc in (150, 480):
            dem -= 1.5 * np.exp(-((X - xc)**2) / (2 * 1.2**2))

        dem_file = tmp_path / "ditch_dem.tif"
        transform = from_bounds(660000, 6700000, 660600, 6700600, size, size)
        with rasterio.open(
            dem_file, 'w', driver='GTiff', height=size, width=size,
            count=1, dtype='float32', crs='EPSG:2154', transform=transform,
        ) as dst:
            dst.write(dem.astype('float32'), 1)

        result = generate_wavelet(dem_file, "ditch", tmp_output_dir, res)
        assert result is not None and result.exists()
        with rasterio.open(result) as src:
            data = src.read(1)

        # Peaks over a vertical band around each ditch
        peak_hilltop = np.nanmax(data[:, 145:156])
        peak_flat = np.nanmax(data[:, 475:486])
        assert peak_flat > 5.0  # the ditch clearly stands out above the background
        # No amplification of the ditch by the topographic position
        assert peak_hilltop < 1.5 * peak_flat

        # Hilltop background (35 m east of the hilltop ditch) vs flat background
        bg_hilltop = np.nanmedian(data[250:350, 185:226])
        bg_flat = np.nanmedian(data[250:350, 500:561])
        assert bg_hilltop < 1.4 * bg_flat  # without detrending: ~1.7×


class TestFlowAccumulation:
    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
        from lidar_pipeline.visualizations import generate_flow_accumulation
        result = generate_flow_accumulation(synthetic_dem, "test", tmp_output_dir, 5.0)
        assert result is not None
        assert result.exists()

    def test_flow_log_values(self, synthetic_dem, tmp_output_dir):
        import rasterio
        from lidar_pipeline.visualizations import generate_flow_accumulation
        result = generate_flow_accumulation(synthetic_dem, "test", tmp_output_dir, 5.0)
        with rasterio.open(result) as src:
            data = src.read(1)
            # log1p(x) >= 0 for x >= 0
            valid = data[~np.isnan(data)]
            assert np.nanmin(valid) >= 0


class TestRayTrace:
    def test_rays_are_traced(self, synthetic_dem, tmp_output_dir):
        """Verify _ray_trace_horizons returns expected shapes."""
        from lidar_pipeline.visualizations import _ray_trace_horizons, _prepare_dem_for_raycast
        import rasterio
        with rasterio.open(synthetic_dem) as src:
            dem_np = src.read(1)
        rows, cols = dem_np.shape
        # Create a simple filled DEM for testing
        import numpy as np
        filled = np.nan_to_num(dem_np, nan=0)
        # Test with numpy (no GPU)
        pos, neg = _ray_trace_horizons(
            filled, rows, cols, 5.0, n_dirs=4, max_dist=10, radii_m=[25, 50]
        )
        assert pos.shape == (4, 2, rows, cols)
        assert neg.shape == (4, 2, rows, cols)


class TestNodataPreserved:
    """Nodata is preserved in the renderings (legacy behaviour).

    DTM holes are filled upstream (create_dtm_fast, all modes); if some
    nodata remains anyway, hillshade/slope/aspect keep it (rendered black on
    the map) instead of inventing interpolated values.
    """

    @staticmethod
    def _dem_with_hole(synthetic_dem, tmp_path):
        import rasterio
        with rasterio.open(synthetic_dem) as src:
            arr = src.read(1).copy()
            profile = src.profile.copy()
        arr[80:120, 80:120] = np.nan
        dem_hole = tmp_path / "dem_hole.tif"
        profile.update(dtype='float32', nodata=float('nan'))
        with rasterio.open(dem_hole, 'w', **profile) as dst:
            dst.write(arr.astype('float32'), 1)
        return dem_hole

    def test_aspect_solo_preserves_nodata(self, synthetic_dem, tmp_path):
        from lidar_pipeline.visualizations import generate_aspect
        dem_hole = self._dem_with_hole(synthetic_dem, tmp_path)
        out = generate_aspect(dem_hole, "solo", tmp_path, 5.0)
        assert out is not None and out.exists()
        import rasterio
        with rasterio.open(out) as src:
            data = src.read(1)
        assert np.isnan(data[80:120, 80:120]).all(), "the hole must stay nodata"
        # The gradient at the hole's edge spreads NaN over a 1 px ring:
        # check an area far from the hole
        assert not np.isnan(data[0:40, 0:40]).any(), "NaN far from the hole"

    def test_aspect_shared_preserves_nodata(self, synthetic_dem, tmp_path):
        from lidar_pipeline.visualizations import SharedDEM, generate_aspect
        dem_hole = self._dem_with_hole(synthetic_dem, tmp_path)
        shared = SharedDEM(dem_hole, 5.0)
        out = generate_aspect(dem_hole, "shared", tmp_path, 5.0, shared=shared)
        assert out is not None and out.exists()
        import rasterio
        with rasterio.open(out) as src:
            data = src.read(1)
        assert np.isnan(data[80:120, 80:120]).all(), "the hole must stay nodata"
        assert not np.isnan(data[0:40, 0:40]).any(), "NaN far from the hole"

    def test_slope_and_hillshade_preserve_nodata(self, synthetic_dem, tmp_path):
        from lidar_pipeline.visualizations import generate_slope, generate_hillshade
        dem_hole = self._dem_with_hole(synthetic_dem, tmp_path)
        import rasterio
        for gen, name in ((generate_slope, "p"), (generate_hillshade, "h")):
            out = gen(dem_hole, name, tmp_path, 5.0)
            assert out is not None and out.exists()
            with rasterio.open(out) as src:
                data = src.read(1)
            assert np.isnan(data[80:120, 80:120]).any(), f"{out.name}: hole disappeared"


def test_ray_trace_horizons_cpu_fallback_on_oom(monkeypatch):
    """On GPU OOM, ray-tracing disables the GPU then retries on CPU."""
    import lidar_pipeline.visualizations as viz
    import lidar_pipeline.gpu as gpu_mod

    calls = {"n": 0}
    disabled = []

    def fake_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=None):
        calls["n"] += 1
        if calls["n"] == 1:
            raise RuntimeError("Out of memory allocating 600,000,000 bytes")
        return ("pos", "neg")

    monkeypatch.setattr(viz, "_ray_trace_horizons_core", fake_core)
    monkeypatch.setattr(gpu_mod, "is_gpu_active", lambda: True)
    monkeypatch.setattr(gpu_mod, "disable_gpu", lambda: disabled.append(True))
    result = viz._ray_trace_horizons(None, 4, 4, 0.5, 8, 10)
    assert result == ("pos", "neg")
    assert calls["n"] == 2
    assert disabled == [True]


def test_ray_trace_horizons_reraises_non_oom(monkeypatch):
    """A non-OOM error is not masked by the CPU fallback."""
    import lidar_pipeline.visualizations as viz
    import lidar_pipeline.gpu as gpu_mod

    def fake_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=None):
        raise ValueError("other error")

    monkeypatch.setattr(viz, "_ray_trace_horizons_core", fake_core)
    monkeypatch.setattr(gpu_mod, "is_gpu_active", lambda: True)
    try:
        viz._ray_trace_horizons(None, 4, 4, 0.5, 8, 10)
        assert False, "ValueError expected"
    except ValueError:
        pass


class TestPriorityFlood:
    def test_numba_matches_python(self):
        """numba and Python results are identical on a DEM with a pit."""
        from lidar_pipeline.visualizations import _priority_flood_numba, _priority_flood_python

        dem = np.zeros((20, 20), dtype=np.float64)
        dem[10, 10] = -5.0
        dem[9:12, 9:12] = -3.0
        nodata = np.zeros((20, 20), dtype=bool)

        result_numba = _priority_flood_numba(dem.copy(), nodata)
        result_python = _priority_flood_python(dem.copy(), nodata)

        if result_numba is not None:
            assert np.allclose(result_numba, result_python)

    def test_pit_is_filled(self):
        """An isolated pit is raised to the level of its rim."""
        from lidar_pipeline.visualizations import _priority_flood

        dem = np.full((10, 10), 5.0, dtype=np.float64)
        dem[5, 5] = 1.0
        nodata = np.zeros((10, 10), dtype=bool)

        result = _priority_flood(dem, nodata)
        assert result[5, 5] == 5.0

    def test_nodata_cells_untouched(self):
        """Nodata cells are never modified."""
        from lidar_pipeline.visualizations import _priority_flood

        dem = np.full((10, 10), 5.0, dtype=np.float64)
        dem[5, 5] = 1.0
        nodata = np.zeros((10, 10), dtype=bool)
        nodata[2, 2] = True
        dem[2, 2] = 999.0

        result = _priority_flood(dem, nodata)
        assert result[2, 2] == 999.0


class TestDensiteSol:
    def test_density_levels_log_scale(self):
        from lidar_pipeline.visualizations import density_levels
        d = np.array([0.0, np.nan, 0.25, 0.36, 0.5, 1.0, 45.0, 45.3, 1000.0])
        assert density_levels(d).tolist() == [0, 0, 0, 1, 2, 4, 14, 15, 15]

    def test_generate_reads_density_sidecar(self, tmp_path):
        import rasterio
        from rasterio.transform import from_bounds
        from lidar_pipeline.dtm import density_path
        from lidar_pipeline.visualizations import generate_densite_sol
        dem = tmp_path / "T_dtm_r0p2.tif"
        dem.touch()
        with rasterio.open(density_path(dem), 'w', driver='GTiff', width=4, height=1,
                           count=1, dtype='float32', crs='EPSG:2154',
                           transform=from_bounds(0, 0, 4, 1, 4, 1)) as dst:
            dst.write(np.array([[0.0, 1.0, 4.0, 64.0]], dtype='float32'), 1)
        out = generate_densite_sol(dem, "T", tmp_path, 0.2)
        with rasterio.open(out) as src:
            assert src.read(1).tolist() == [[0, 4, 8, 15]]
            assert src.width == 4  # 1 m grid kept

    def test_missing_sidecar_returns_none(self, tmp_path):
        from lidar_pipeline.visualizations import generate_densite_sol
        assert generate_densite_sol(tmp_path / "X_dtm.tif", "X", tmp_path, 0.2) is None