lidar_rendu/lidar_pipeline/tests/test_dtm.py

"""Tests for DTM module."""

import json
import numpy as np
import pytest
from pathlib import Path
from unittest.mock import patch, MagicMock


class TestSMRFPipeline:
    def test_pipeline_json_valid(self):
        """create_smrf_pipeline produces valid JSON with expected stages."""
        from lidar_pipeline.dtm import create_smrf_pipeline
        result = create_smrf_pipeline("/data/input/test.laz", "/data/output/test_ground.las")
        pipeline = json.loads(result)

        assert "pipeline" in pipeline
        stages = pipeline["pipeline"]

        # Should have: reader, range filter (ReturnNumber), assign, ELM, outlier, SMRF, range filter (Classification), writer
        stage_types = [s.get("type") if isinstance(s, dict) else None for s in stages]

        # First stage is the filename string (reader)
        assert isinstance(stages[0], str)
        assert "test.laz" in stages[0]

        # Must contain preprocessing steps
        assert "filters.assign" in stage_types
        assert "filters.elm" in stage_types
        assert "filters.outlier" in stage_types

        # Must contain SMRF filter
        assert "filters.smrf" in stage_types

        # Must contain ReturnNumber filter
        range_stages = [s for s in stages if isinstance(s, dict) and s.get("type") == "filters.range"]
        assert len(range_stages) >= 1
        # At least one should filter ReturnNumber
        assert any("ReturnNumber" in str(s.get("limits", "")) for s in range_stages)

    def test_pipeline_elm_parameters(self):
        """ELM filter has terrain-adapted parameters."""
        from lidar_pipeline.dtm import create_smrf_pipeline
        result = create_smrf_pipeline("/input/a.laz", "/output/a_ground.las")
        pipeline = json.loads(result)

        elm_stage = [s for s in pipeline["pipeline"] if isinstance(s, dict) and s.get("type") == "filters.elm"][0]
        assert elm_stage["cell"] == 5.0
        assert elm_stage["threshold"] == 2.0

    def test_pipeline_outlier_parameters(self):
        """Outlier filter uses statistical method."""
        from lidar_pipeline.dtm import create_smrf_pipeline
        result = create_smrf_pipeline("/input/a.laz", "/output/a_ground.las")
        pipeline = json.loads(result)

        outlier_stage = [s for s in pipeline["pipeline"] if isinstance(s, dict) and s.get("type") == "filters.outlier"][0]
        assert outlier_stage["method"] == "statistical"
        assert outlier_stage["mean_k"] == 8
        assert outlier_stage["multiplier"] == 3.0

    def test_pipeline_output_path(self):
        """Pipeline output path is set correctly."""
        from lidar_pipeline.dtm import create_smrf_pipeline
        result = create_smrf_pipeline("/input/a.laz", "/output/a_ground.las")
        pipeline = json.loads(result)
        # Last stage should be writer with correct output path
        writer = [s for s in pipeline["pipeline"] if isinstance(s, dict) and s.get("type") == "writers.las"][0]
        assert writer["filename"] == "/output/a_ground.las"


class TestCSFPipeline:
    def test_pipeline_json_valid(self):
        """create_csf_pipeline produces valid JSON with CSF filter."""
        from lidar_pipeline.dtm import create_csf_pipeline
        result = create_csf_pipeline("/data/input/test.laz", "/data/output/test_ground.las")
        pipeline = json.loads(result)

        assert "pipeline" in pipeline
        stages = pipeline["pipeline"]
        stage_types = [s.get("type") if isinstance(s, dict) else None for s in stages]

        # Must contain CSF filter
        assert "filters.csf" in stage_types

        # Must contain ReturnNumber filter
        range_stages = [s for s in stages if isinstance(s, dict) and s.get("type") == "filters.range"]
        assert any("ReturnNumber" in str(s.get("limits", "")) for s in range_stages)

    def test_csf_parameters(self):
        """CSF pipeline has expected parameters."""
        from lidar_pipeline.dtm import create_csf_pipeline
        result = create_csf_pipeline("/input/a.laz", "/output/a_ground.las")
        pipeline = json.loads(result)

        csf_stage = [s for s in pipeline["pipeline"] if isinstance(s, dict) and s.get("type") == "filters.csf"][0]
        assert csf_stage["resolution"] == 1.0  # cloth 1 m : ~4× plus rapide, MNT inchangé
        assert csf_stage["rigidness"] == 3
        assert csf_stage["smooth"] is True
        assert "hdiff" not in csf_stage  # hdiff is not a valid PDAL CSF parameter


class TestInterpolateHoles:
    def test_fills_interior_hole_with_surface(self):
        """Large interior NaN hole is filled (no NaN left, value is plausible)."""
        from lidar_pipeline.dtm import _interpolate_holes
        # Linear-in-column surface z = 0.02 * x, with a large square hole in the middle.
        x = np.arange(40, dtype=float) * 0.02
        dtm = np.tile(x, (40, 1))
        dtm[16:24, 16:24] = np.nan
        filled, count = _interpolate_holes(dtm)
        assert count == 64
        assert not np.isnan(filled).any()
        # Filled values stay within the surrounding z range (no wild extrapolation).
        zmin, zmax = np.nanmin(dtm), np.nanmax(dtm)
        hole_vals = filled[16:24, 16:24]
        assert np.all(hole_vals >= zmin - 1e-6)
        assert np.all(hole_vals <= zmax + 1e-6)
        # A linear surface is interpolated near-exactly in the interior.
        expected = np.tile(x[16:24], (8, 1))
        assert np.allclose(hole_vals, expected, atol=0.02)
        # Original valid cells are untouched.
        valid = ~np.isnan(dtm)
        assert np.allclose(filled[valid], dtm[valid])

    def test_no_holes_returns_unchanged(self):
        """No NaN → returns same array and zero count."""
        from lidar_pipeline.dtm import _interpolate_holes
        dtm = np.arange(64, dtype=float).reshape(8, 8)
        filled, count = _interpolate_holes(dtm)
        assert count == 0
        assert np.shares_memory(filled, dtm)

    def test_all_nan_returns_unchanged(self):
        """No valid data → cannot interpolate, returns zeros-free NaN array."""
        from lidar_pipeline.dtm import _interpolate_holes
        dtm = np.full((8, 8), np.nan)
        filled, count = _interpolate_holes(dtm)
        assert count == 0
        assert np.isnan(filled).all()


class TestDetectGroundMethod:
    def _make_mock_las(self, num_returns, z_values):
        """Create a mock laspy object with specified NumberOfReturns and z."""
        mock_las = MagicMock()
        mock_las.NumberOfReturns = np.array(num_returns)
        mock_las.z = np.array(z_values)
        mock_las.points = MagicMock()
        mock_las.points.__len__ = lambda self: len(num_returns)
        return mock_las

    @patch('lidar_pipeline.dtm._read_with_pdal')
    @patch('laspy.read')
    def test_urban_terrain_returns_csf(self, mock_read, mock_pdal):
        """High single-return ratio (>0.6) should select CSF."""
        from lidar_pipeline.dtm import detect_ground_method

        # 70% single returns = urban
        n = 10000
        num_returns = np.ones(n, dtype=int)
        num_returns[:int(n * 0.3)] = 2  # 30% multi-return
        z_values = np.random.normal(100, 5, n)  # Low variance = flat terrain

        mock_read.return_value = self._make_mock_las(num_returns, z_values)

        result = detect_ground_method(Path("/data/input/test.laz"))
        assert result == 'csf'

    @patch('lidar_pipeline.dtm._read_with_pdal')
    @patch('laspy.read')
    def test_natural_terrain_returns_smrf(self, mock_read, mock_pdal):
        """Low single-return ratio and moderate variance should select SMRF."""
        from lidar_pipeline.dtm import detect_ground_method

        # 40% single returns, moderate variance
        n = 10000
        num_returns = np.ones(n, dtype=int)
        num_returns[:int(n * 0.6)] = 2  # 60% multi-return (forest)
        z_values = np.random.normal(100, 15, n)  # Moderate variance

        mock_read.return_value = self._make_mock_las(num_returns, z_values)

        result = detect_ground_method(Path("/data/input/test.laz"))
        assert result == 'smrf'

    @patch('lidar_pipeline.dtm._read_with_pdal')
    @patch('laspy.read')
    def test_mountainous_terrain_returns_csf(self, mock_read, mock_pdal):
        """High variance terrain (>30m std) selects CSF for complex terrain."""
        from lidar_pipeline.dtm import detect_ground_method

        # Moderate single-return ratio but very high height variance
        n = 10000
        num_returns = np.ones(n, dtype=int)
        num_returns[:int(n * 0.5)] = 2
        z_values = np.random.normal(100, 50, n)  # Very high variance = mountainous

        mock_read.return_value = self._make_mock_las(num_returns, z_values)

        result = detect_ground_method(Path("/data/input/test.laz"))
        assert result == 'csf'


class TestIGNPipeline:
    def test_pipeline_keeps_supplier_classification(self):
        """create_ign_pipeline réutilise la pré-classification (classe 2) sans refiltrer."""
        from lidar_pipeline.dtm import create_ign_pipeline
        result = create_ign_pipeline("/input/a.laz", "/output/a_ground.las")
        pipeline = json.loads(result)

        stages = pipeline["pipeline"]
        stage_types = [s.get("type") if isinstance(s, dict) else None for s in stages]

        # Aucun algorithme de classification, pas de remise à zéro, pas de filtres de bruit
        assert "filters.smrf" not in stage_types
        assert "filters.csf" not in stage_types
        assert "filters.assign" not in stage_types
        assert "filters.elm" not in stage_types
        assert "filters.outlier" not in stage_types

        # Filtre ReturnNumber conservé + extraction des points classe 2
        range_stages = [s for s in stages if isinstance(s, dict) and s.get("type") == "filters.range"]
        assert any("ReturnNumber" in str(s.get("limits", "")) for s in range_stages)
        assert any(s.get("limits") == "Classification[2:2]" for s in range_stages)

        writer = [s for s in stages if isinstance(s, dict) and s.get("type") == "writers.las"][0]
        assert writer["filename"] == "/output/a_ground.las"


class TestDetectIGN:
    def _make_mock_las(self, classification, num_returns, z_values):
        mock_las = MagicMock()
        mock_las.classification = classification
        mock_las.NumberOfReturns = np.array(num_returns)
        mock_las.z = np.array(z_values)
        mock_las.points = MagicMock()
        mock_las.points.__len__ = lambda self: len(num_returns)
        return mock_las

    @patch('lidar_pipeline.dtm._read_with_pdal')
    @patch('laspy.read')
    def test_preclassified_returns_ign(self, mock_read, mock_pdal):
        """Fichier pré-classifié (majorité classe 2) → méthode IGN."""
        from lidar_pipeline.dtm import detect_ground_method

        n = 10000
        num_returns = np.ones(n, dtype=int)
        cls = np.zeros(n, dtype=np.uint8)
        cls[int(n * 0.15):] = 2  # 85 % de points classe 2
        z_values = np.random.normal(100, 5, n)

        mock_read.return_value = self._make_mock_las(cls, num_returns, z_values)
        assert detect_ground_method(Path("/data/input/test.laz")) == 'ign'

    @patch('lidar_pipeline.dtm._read_with_pdal')
    @patch('laspy.read')
    def test_unclassified_falls_back_to_smrf_or_csf(self, mock_read, mock_pdal):
        """Sans classification exploitable → détection SMRF/CSF classique."""
        from lidar_pipeline.dtm import detect_ground_method

        n = 10000
        num_returns = np.ones(n, dtype=int)
        num_returns[:int(n * 0.6)] = 2  # 60 % multi-retours (forêt) → non urbain
        cls = np.zeros(n, dtype=np.uint8)  # aucun point classe 2
        z_values = np.random.normal(100, 5, n)

        mock_read.return_value = self._make_mock_las(cls, num_returns, z_values)
        assert detect_ground_method(Path("/data/input/test.laz")) == 'smrf'


class TestClassifyGroundMethod:
    @patch('lidar_pipeline.dtm.subprocess')
    def test_classify_ground_auto_calls_detect(self, mock_subprocess):
        """classify_ground with method='auto' should call detect_ground_method."""
        from lidar_pipeline.dtm import classify_ground

        # Mock detect_ground_method to return 'csf'
        with patch('lidar_pipeline.dtm.detect_ground_method', return_value='csf') as mock_detect:
            mock_subprocess.run.return_value = MagicMock(returncode=0)
            result = classify_ground(Path("/data/input/test.laz"), Path("/tmp"), method='auto')

            mock_detect.assert_called_once()

    @patch('lidar_pipeline.dtm.subprocess')
    def test_classify_ground_smrf_uses_smrf_pipeline(self, mock_subprocess):
        """classify_ground with method='smrf' should create SMRF pipeline."""
        from lidar_pipeline.dtm import classify_ground, _create_ground_pipeline

        mock_subprocess.run.return_value = MagicMock(returncode=0)

        with patch('lidar_pipeline.dtm.detect_ground_method'):
            # Create temp dir
            import tempfile
            with tempfile.TemporaryDirectory() as tmpdir:
                result = classify_ground(Path("/data/input/test.laz"), Path(tmpdir), method='smrf')

                # Check the pipeline JSON was written with SMRF
                pipeline_file = Path(tmpdir) / "pipeline_smrf.json"
                if pipeline_file.exists():
                    pipeline = json.loads(pipeline_file.read_text())
                    stage_types = [s.get("type") if isinstance(s, dict) else None for s in pipeline["pipeline"]]
                    assert "filters.smrf" in stage_types

    @patch('lidar_pipeline.dtm.subprocess')
    def test_classify_ground_csf_uses_csf_pipeline(self, mock_subprocess):
        """classify_ground with method='csf' should create CSF pipeline."""
        from lidar_pipeline.dtm import classify_ground

        mock_subprocess.run.return_value = MagicMock(returncode=0)

        with patch('lidar_pipeline.dtm.detect_ground_method'):
            import tempfile
            with tempfile.TemporaryDirectory() as tmpdir:
                result = classify_ground(Path("/data/input/test.laz"), Path(tmpdir), method='csf')

                pipeline_file = Path(tmpdir) / "pipeline_csf.json"
                if pipeline_file.exists():
                    pipeline = json.loads(pipeline_file.read_text())
                    stage_types = [s.get("type") if isinstance(s, dict) else None for s in pipeline["pipeline"]]
                    assert "filters.csf" in stage_types

class TestParseIgnClasses:
    def test_default_sol(self):
        """'sol' → code 2 seul."""
        from lidar_pipeline.dtm import parse_ign_classes
        assert parse_ign_classes("sol") == [2]

    def test_names_sorted_dedup(self):
        """Noms acceptés (EN/FR), triés et dédupliqués."""
        from lidar_pipeline.dtm import parse_ign_classes
        assert parse_ign_classes("sol,unclassified") == [1, 2]
        assert parse_ign_classes("unclassified,sol") == [1, 2]
        assert parse_ign_classes("non-classe") == [1]
        assert parse_ign_classes("sol,2") == [2]

    def test_numeric_codes(self):
        """Codes LAS directs, triés."""
        from lidar_pipeline.dtm import parse_ign_classes
        assert parse_ign_classes("2,1") == [1, 2]
        assert parse_ign_classes("66") == [66]

    def test_invalid_raises(self):
        """Nom inconnu, code hors bornes ou liste vide → ValueError."""
        from lidar_pipeline.dtm import parse_ign_classes
        with pytest.raises(ValueError):
            parse_ign_classes("foo")
        with pytest.raises(ValueError):
            parse_ign_classes("300")
        with pytest.raises(ValueError):
            parse_ign_classes("")

    def test_method_label(self):
        """'ign' seul pour le sol, combinaison encodée sinon (invalidation cache)."""
        from lidar_pipeline.dtm import ign_method_label
        assert ign_method_label([2]) == "ign"
        assert ign_method_label([1, 2]) == "ign_1_2"
        assert ign_method_label([2, 1]) == "ign_1_2"


class TestIGNPipelineMultiClasses:
    def test_multi_class_limits(self):
        """Plusieurs classes → plages OU logiques sur Classification."""
        from lidar_pipeline.dtm import _create_ground_pipeline
        result = _create_ground_pipeline("/input/a.laz", "/output/a_ground.las",
                                         'ign', ign_codes=[1, 2])
        pipeline = json.loads(result)
        range_stages = [s for s in pipeline["pipeline"]
                        if isinstance(s, dict) and s.get("type") == "filters.range"]
        limits = [str(s.get("limits", "")) for s in range_stages]
        assert any("Classification[1:1]" in l and "Classification[2:2]" in l
                   for l in limits)

    def test_default_sol_only(self):
        """Sans ign_codes, la voie IGN reste sol seul (2) — rétrocompatible."""
        from lidar_pipeline.dtm import _create_ground_pipeline
        result = _create_ground_pipeline("/input/a.laz", "/output/a_ground.las", 'ign')
        pipeline = json.loads(result)
        range_stages = [s for s in pipeline["pipeline"]
                        if isinstance(s, dict) and s.get("type") == "filters.range"]
        limits = [str(s.get("limits", "")) for s in range_stages]
        assert any("Classification[2:2]" in l and "Classification[1:1]" not in l
                   for l in limits)


class TestClassifyGroundIgnClasses:
    @patch('lidar_pipeline.dtm.subprocess')
    def test_ign_classes_encoded_in_filenames(self, mock_subprocess):
        """--ign-classes sol,unclassified → fichiers ign_1_2 + filtre multi-classes."""
        import tempfile
        from lidar_pipeline.dtm import classify_ground

        mock_subprocess.run.return_value = MagicMock(returncode=0)

        with tempfile.TemporaryDirectory() as tmpdir:
            tmpdir = Path(tmpdir)
            classify_ground(Path("/data/input/test.laz"), tmpdir,
                            method='ign', ign_classes="sol,unclassified")

            pipeline_file = tmpdir / "pipeline_ign_1_2.json"
            assert pipeline_file.exists()
            pipeline = json.loads(pipeline_file.read_text())
            limits = [str(s.get("limits", "")) for s in pipeline["pipeline"]
                      if isinstance(s, dict) and s.get("type") == "filters.range"]
            assert any("Classification[1:1]" in l and "Classification[2:2]" in l
                       for l in limits)

    @patch('lidar_pipeline.dtm.subprocess')
    def test_ign_default_label_unchanged(self, mock_subprocess):
        """--ign-classes sol (défaut) → noms 'ign' inchangés (cache préservé)."""
        import tempfile
        from lidar_pipeline.dtm import classify_ground

        mock_subprocess.run.return_value = MagicMock(returncode=0)

        with tempfile.TemporaryDirectory() as tmpdir:
            tmpdir = Path(tmpdir)
            classify_ground(Path("/data/input/test.laz"), tmpdir, method='ign')

            assert (tmpdir / "pipeline_ign.json").exists()
            assert not (tmpdir / "pipeline_ign_1_2.json").exists()


class TestPureDtm:
    """Mode pur (classification IGN) : aucune retouche, trous en nodata."""

    def _write_las(self, path, points):
        import laspy
        hdr = laspy.LasHeader(version='1.2', point_format=0)
        las = laspy.LasData(hdr)
        las.x = [p[0] for p in points]
        las.y = [p[1] for p in points]
        las.z = [p[2] for p in points]
        las.write(str(path))
        return path

    def _make_clouds(self, tmp_output_dir):
        """Grille 2x2 (res=1.0). Sol sur 3 cellules (z=10), trou en (1,1).
        Le nuage complet a un retour plus bas (z=7) dans le trou."""
        corners = [(0.05, 0.05, 10.0), (1.95, 0.05, 10.0), (0.05, 1.95, 10.0)]
        ground = [(0.5, 0.5, 10.0), (1.5, 0.5, 10.0), (0.5, 1.5, 10.0)] + corners
        source = list(ground) + [(1.5, 1.5, 7.0)]
        las_file = self._write_las(tmp_output_dir / "ground_pure.las", ground)
        source_laz = self._write_las(tmp_output_dir / "source_pure.las", source)
        return las_file, source_laz

    def _dtm_array(self, tmp_output_dir, pure):
        from lidar_pipeline.dtm import create_dtm_fast
        import rasterio
        las_file, source_laz = self._make_clouds(tmp_output_dir)
        out = create_dtm_fast(las_file, "tile_pure", tmp_output_dir, 1.0,
                              force=True, source_laz=source_laz, pure=pure)
        assert out is not None
        with rasterio.open(str(out)) as src:
            return src.read(1).astype("float64")

    def test_pure_fills_holes_without_floor(self, tmp_output_dir):
        """pur=True : trous comblés par interpolation, sans plancher à 7.

        Le comblement est actif dans tous les modes (comportement
        historique) ; « pur » ne désactive que l'abaissement au retour
        le plus bas.
        """
        arr = self._dtm_array(tmp_output_dir, pure=True)
        assert int(np.isnan(arr).sum()) == 0
        vals = sorted(float(v) for v in arr.flatten())
        assert vals == [10.0, 10.0, 10.0, 10.0]

    def test_not_pure_fills_holes(self, tmp_output_dir):
        """pur=False : le trou est comblé (comportement historique conservé)."""
        arr = self._dtm_array(tmp_output_dir, pure=False)
        assert int(np.isnan(arr).sum()) == 0
        vals = sorted(float(v) for v in arr.flatten())
        assert len(vals) == 4
        assert vals[-1] == 10.0


class TestStripLidarExt:
    def test_copc_laz(self):
        from lidar_pipeline.dtm import _strip_lidar_ext
        assert _strip_lidar_ext("LHD_FXX_1000_6881_PTS_LAMB93_IGN69.copc.laz") == "LHD_FXX_1000_6881_PTS_LAMB93_IGN69"

    def test_laz(self):
        from lidar_pipeline.dtm import _strip_lidar_ext
        assert _strip_lidar_ext("file.laz") == "file"

    def test_las(self):
        from lidar_pipeline.dtm import _strip_lidar_ext
        assert _strip_lidar_ext("file.las") == "file"

    def test_path_object(self):
        from lidar_pipeline.dtm import _strip_lidar_ext
        from pathlib import Path
        assert _strip_lidar_ext(Path("/data/input/file.copc.laz")) == "file"


class TestStripVerticalOffsets:
    """Vertical de-bias of flight-line point sources (strip alignment)."""

    def _synthetic(self, offsets, n=1_500_000, extent=300.0, seed=0):
        """Interleaved sources over one tile; each carries a known Z bias."""
        rng = np.random.default_rng(seed)
        x = rng.uniform(0, extent, n)
        y = rng.uniform(0, extent, n)
        z = 100 + 0.02 * x - 0.01 * y + rng.normal(0, 0.01, n)
        psid = rng.integers(0, len(offsets), n)
        return x, y, z + np.asarray(offsets)[psid], psid.astype(np.uint16)

    def test_recovers_known_offsets(self):
        """±5 cm biases are recovered; the aligned source stays untouched."""
        from lidar_pipeline.dtm import _strip_vertical_offsets
        x, y, z, psid = self._synthetic((0.0, 0.05, -0.05))
        offs = _strip_vertical_offsets(x, y, z, psid)
        assert abs(offs.get(1, 0.0) - 0.05) < 0.01
        assert abs(offs.get(2, 0.0) + 0.05) < 0.01
        assert 0 not in offs  # biais ~0 < seuil : pas de correction

    def test_small_bias_below_threshold_ignored(self):
        """Biases under the 0.5 cm threshold trigger no correction."""
        from lidar_pipeline.dtm import _strip_vertical_offsets
        x, y, z, psid = self._synthetic((0.0, 0.003, -0.003))
        assert _strip_vertical_offsets(x, y, z, psid) == {}

    def test_single_source_returns_empty(self):
        """A single point source cannot be compared: no offsets."""
        from lidar_pipeline.dtm import _strip_vertical_offsets
        x, y, z, psid = self._synthetic((0.05,))
        assert _strip_vertical_offsets(x, y, z, psid) == {}

    def test_no_shared_cells_returns_empty(self):
        """Sources covering disjoint areas (no overlap) are not corrected."""
        from lidar_pipeline.dtm import _strip_vertical_offsets
        rng = np.random.default_rng(1)
        n = 200_000
        x = np.concatenate([rng.uniform(0, 100, n), rng.uniform(200, 300, n)])
        y = rng.uniform(0, 300, 2 * n)
        z = 100 + rng.normal(0, 0.01, 2 * n)
        psid = np.concatenate([np.zeros(n, np.uint16), np.ones(n, np.uint16)])
        z[psid == 1] += 0.10
        assert _strip_vertical_offsets(x, y, z, psid) == {}


class TestStripJitterOffsets:
    """Gigue verticale intra-faisceau par fenêtres de temps GPS."""

    def _synthetic(self, bias_fn, n=800_000, extent=300.0, duration=20.0, seed=0):
        """Deux faisceaux entrelacés ; le n° 1 porte un biais dépendant du temps."""
        rng = np.random.default_rng(seed)
        x = rng.uniform(0, extent, n)
        y = rng.uniform(0, extent, n)
        clean = 100 + 0.02 * x - 0.01 * y + rng.normal(0, 0.01, n)
        t = rng.uniform(0, duration, n)
        psid = rng.integers(0, 2, n).astype(np.uint16)
        return x, y, clean + bias_fn(t, psid), psid, t, clean

    def test_recovers_time_varying_offset(self):
        """Une oscillation lente ±4 cm du faisceau 1 est retirée du terrain vrai.

        En recouvrement à deux, la référence hors-faisceau attribue une série
        à chaque faisceau (chacun absorbe sa part) : on vérifie le résidu
        contre le terrain synthétique propre, fenêtre par fenêtre.
        """
        from lidar_pipeline.dtm import _strip_jitter_offsets, _apply_strip_jitter
        w = 2 * np.pi / 8.0
        x, y, z, psid, t, clean = self._synthetic(
            lambda tt, p: 0.04 * np.sin(w * tt) * (p == 1))
        jitter = _strip_jitter_offsets(x, y, z, psid, t)
        assert set(jitter) == {0, 1}
        resid = z - _apply_strip_jitter(psid, t, jitter) - clean
        assert np.sqrt(np.mean(resid ** 2)) < 0.012
        for lo in np.arange(0, 20.0, 2.0):
            m = (t >= lo) & (t < lo + 2.0)
            assert abs(resid[m].mean()) < 0.012, f"fenêtre {lo:.0f}-{lo + 2:.0f} s"

    def test_tracks_step_offset(self):
        """Un échelon −3 cm sur la seconde moitié du vol est suivi."""
        from lidar_pipeline.dtm import _strip_jitter_offsets, _apply_strip_jitter
        x, y, z, psid, t, clean = self._synthetic(
            lambda tt, p: np.where(tt >= 10.0, -0.03, 0.0) * (p == 1))
        jitter = _strip_jitter_offsets(x, y, z, psid, t)
        resid = z - _apply_strip_jitter(psid, t, jitter) - clean
        assert np.sqrt(np.mean(resid ** 2)) < 0.012
        for lo in (3.0, 6.0, 13.0, 16.0):  # loin de la transition lissée
            m = (t >= lo) & (t < lo + 2.0)
            assert abs(resid[m].mean()) < 0.012, f"fenêtre {lo:.0f}-{lo + 2:.0f} s"

    def test_apply_interpolates_linearly(self):
        """Interpolation entre centres de fenêtres ; 0 hors faisceau connu."""
        from lidar_pipeline.dtm import _apply_strip_jitter
        jitter = {7: (np.array([10.0, 11.0]), np.array([0.0, 0.1]))}
        psid = np.array([7, 7, 7, 3], dtype=np.uint16)
        t = np.array([10.0, 10.5, 15.0, 10.5])
        np.testing.assert_allclose(
            _apply_strip_jitter(psid, t, jitter), [0.0, 0.05, 0.1, 0.0])

    def test_requires_two_sources_and_time(self):
        """Faisceau unique ou temps non fini : rien à corriger."""
        from lidar_pipeline.dtm import _strip_jitter_offsets
        x, y, z, psid, t, _clean = self._synthetic(lambda tt, p: 0.04 * np.sin(tt) * (p == 1))
        assert _strip_jitter_offsets(x, y, z, np.zeros_like(psid), t) == {}
        t_nan = t.copy()
        t_nan[0] = np.nan
        assert _strip_jitter_offsets(x, y, z, psid, t_nan) == {}


class TestStripAlignSidecar:
    def test_sidecar_roundtrip_and_threshold(self, tmp_path):
        """Sidecar records version/threshold/offsets and matches config."""
        from lidar_pipeline.dtm import (
            _write_strip_align_sidecar, STRIP_ALIGN_VERSION, STRIP_ALIGN_THRESHOLD)
        import json
        offsets = {1049: 0.026, 1147: -0.026}
        _write_strip_align_sidecar(tmp_path, "TILE", "_r0p2", offsets)
        data = json.loads((tmp_path / "TILE_dtm_r0p2_stripalign.json").read_text())
        assert data["version"] == STRIP_ALIGN_VERSION
        assert data["threshold"] == STRIP_ALIGN_THRESHOLD
        assert data["offsets"] == {"1049": 0.026, "1147": -0.026}  # clés JSON en chaînes
        assert data["jitter"] == {}

    def test_sidecar_records_jitter_series(self, tmp_path):
        """Le sidecar consigne les séries de gigue et leurs paramètres."""
        from lidar_pipeline.dtm import (
            _write_strip_align_sidecar, STRIP_JITTER_BIN, STRIP_JITTER_SMOOTH)
        import json
        jitter = {11: (np.array([0.05, 0.15]), np.array([0.012, -0.008]))}
        _write_strip_align_sidecar(tmp_path, "T", "", {}, jitter)
        data = json.loads((tmp_path / "T_dtm_stripalign.json").read_text())
        assert data["jitter_bin"] == STRIP_JITTER_BIN
        assert data["jitter_smooth"] == STRIP_JITTER_SMOOTH
        entry = data["jitter"]["11"]
        assert entry["bins"] == 2
        assert entry["series_m"] == [0.012, -0.008]
        assert entry["max_m"] == 0.012

    def test_pipeline_match_logic(self, tmp_path):
        """_strip_align_matches invalidates legacy DTMs and config changes."""
        from lidar_pipeline.pipeline import LidarArchaeoPipeline
        from lidar_pipeline.dtm import _write_strip_align_sidecar
        import json

        class P(LidarArchaeoPipeline):
            def __init__(self, out, strip_align):
                self.output_dir = out
                self.dtm_dir = out / "DTM"
                self.dtm_dir.mkdir(exist_ok=True)
                self.strip_align = strip_align

        p = P(tmp_path, strip_align=True)
        # DTM hérité sans sidecar : à régénérer
        assert not p._strip_align_matches("TILE", "_r0p2")
        # Sidecar conforme : valide
        _write_strip_align_sidecar(p.dtm_dir, "TILE", "_r0p2", {})
        assert p._strip_align_matches("TILE", "_r0p2")
        # Seuil différent : à régénérer
        bad = p.dtm_dir / "TILE_dtm_r0p2_stripalign.json"
        bad.write_text(json.dumps({"version": 1, "threshold": 0.02, "offsets": {}}))
        assert not p._strip_align_matches("TILE", "_r0p2")
        # Paramètres de gigue différents : à régénérer
        bad.write_text(json.dumps({"version": 2, "threshold": 0.005, "offsets": {},
                                   "jitter_bin": 0.5, "jitter_smooth": 5}))
        assert not p._strip_align_matches("TILE", "_r0p2")
        # Calage désactivé + DTM calé : à régénérer
        assert not P(tmp_path, strip_align=False)._strip_align_matches("TILE", "_r0p2")
        # Paramètres ligne à ligne différents : à régénérer
        _write_strip_align_sidecar(p.dtm_dir, "TILE", "_r0p2", {})
        data = json.loads(bad.read_text())
        data["line_window"] = 99
        bad.write_text(json.dumps(data))
        assert not p._strip_align_matches("TILE", "_r0p2")



class TestEdgeBuffer:
    """Raccord des bords : MNT étendu par les points sol des tuiles voisines."""

    BASENAME = "LHD_FXX_0638_6628_PTS_LAMB93_IGN69"
    # Grille LHD : (col, row) = coin nord-ouest → 0638_6628 couvre
    # X ∈ [638000, 639000], Y ∈ [6627000, 6628000] (bord nord = 6628 km).
    NOMINAL = (638000.0, 6627000.0, 639000.0, 6628000.0)  # dalle 1 km

    def _write_las(self, path, points, classification=None):
        import laspy
        hdr = laspy.LasHeader(version='1.2', point_format=0)
        las = laspy.LasData(hdr)
        las.x = [p[0] for p in points]
        las.y = [p[1] for p in points]
        las.z = [p[2] for p in points]
        if classification is not None:
            las.classification = classification
        las.write(str(path))
        return path

    def _write_clouds(self, root, res=50.0):
        """Tuile centrale à z=10 + voisine EST à z=20 (un point par maille res)."""
        import numpy as np
        input_dir = root / "input"
        input_dir.mkdir(exist_ok=True)
        min_x, min_y, max_x, max_y = self.NOMINAL
        xs = np.arange(min_x + res / 2, max_x, res)
        ys = np.arange(min_y + res / 2, max_y, res)
        gx, gy = np.meshgrid(xs, ys)
        main = list(zip(gx.ravel(), gy.ravel(), np.full(gx.size, 10.0)))
        nxs = xs + 1000.0
        ngx, ngy = np.meshgrid(nxs, ys)
        east = list(zip(ngx.ravel(), ngy.ravel(), np.full(ngx.size, 20.0)))
        ground = self._write_las(root / "ground.las", main)
        source = self._write_las(input_dir / f"{self.BASENAME}.copc.laz", main)
        self._write_las(input_dir / "LHD_FXX_0639_6628_PTS_LAMB93_IGN69.copc.laz",
                        east, classification=[2] * len(east))
        return ground, source, input_dir

    def test_neighbor_discovery(self, tmp_output_dir):
        from lidar_pipeline.dtm import _neighbor_laz_files
        (tmp_output_dir / f"{self.BASENAME}.copc.laz").touch()
        present = [(637, 6627), (639, 6629), (638, 6629)]
        for c, r in present:
            (tmp_output_dir / f"LHD_FXX_{c}_{r}_PTS_LAMB93_IGN69.copc.laz").touch()
        # Bruit non voisin : jamais retenu
        (tmp_output_dir / "LHD_FXX_0650_6700_PTS_LAMB93_IGN69.copc.laz").touch()
        found = _neighbor_laz_files(tmp_output_dir / f"{self.BASENAME}.copc.laz")
        assert {f.name for f in found} == {
            f"LHD_FXX_{c}_{r}_PTS_LAMB93_IGN69.copc.laz" for c, r in present}

    def test_neighbor_discovery_non_lhd(self, tmp_output_dir):
        from lidar_pipeline.dtm import _neighbor_laz_files
        src = tmp_output_dir / "nuage_arbitraire.laz"
        src.touch()
        assert _neighbor_laz_files(src) == []

    def test_neighbor_found_in_edge_subdir(self, tmp_output_dir):
        """Une voisine isolée dans edge_neighbors/ est trouvée ; la priorité
        reste à une dalle à plat dans input/."""
        from lidar_pipeline.dtm import _neighbor_laz_files, EDGE_NEIGHBORS_DIRNAME
        base = "LHD_FXX_0637_6627_PTS_LAMB93_IGN69"
        (tmp_output_dir / f"{base}.copc.laz").touch()
        edge = tmp_output_dir / EDGE_NEIGHBORS_DIRNAME
        edge.mkdir()
        # Voisine uniquement dans le sous-dossier de raccord
        (edge / "LHD_FXX_0638_6628_PTS_LAMB93_IGN69.copc.laz").touch()
        # Voisine présente aux deux endroits : la version input/ gagne
        (edge / "LHD_FXX_0636_6626_PTS_LAMB93_IGN69.copc.laz").touch()
        (tmp_output_dir / "LHD_FXX_0636_6626_PTS_LAMB93_IGN69.copc.laz").touch()
        found = _neighbor_laz_files(tmp_output_dir / f"{base}.copc.laz")
        by_name = {f.name: f for f in found}
        assert by_name["LHD_FXX_0638_6628_PTS_LAMB93_IGN69.copc.laz"].parent == edge
        assert by_name["LHD_FXX_0636_6626_PTS_LAMB93_IGN69.copc.laz"].parent == tmp_output_dir

    def test_buffered_dtm_extends_into_neighbor(self, tmp_output_dir):
        """MNT 24x24 (dalle 20x20 + bande 100 m), bande EST remplie à z=20 par la voisine."""
        from lidar_pipeline.dtm import create_dtm_fast, read_dtm_edge_buffer, EDGE_BUFFER_TAG
        import rasterio
        ground, source, _ = self._write_clouds(tmp_output_dir)
        dtm = create_dtm_fast(ground, self.BASENAME, tmp_output_dir, 50.0,
                              force=True, source_laz=source, strip_align=False,
                              edge_buffer=100.0, neighbor_classes=[2])
        assert dtm is not None
        with rasterio.open(str(dtm)) as src:
            assert (src.width, src.height) == (24, 24)
            assert abs(src.bounds.left - 637900.0) < 1e-6
            assert abs(src.bounds.top - 6628100.0) < 1e-6
            assert src.tags().get(EDGE_BUFFER_TAG) == "100"
            arr = src.read(1)
        assert arr[12, 12] == 10.0      # cœur : tuile centrale
        assert arr[12, 23] == 20.0      # bande EST : points de la voisine
        assert np.isnan(arr[0, 0])      # bande OUEST sans voisine : vide
        assert read_dtm_edge_buffer(dtm) == 100.0

    def test_unbuffered_dtm_has_no_tag(self, tmp_output_dir):
        from lidar_pipeline.dtm import create_dtm_fast, read_dtm_edge_buffer, EDGE_BUFFER_TAG
        import rasterio
        ground, source, _ = self._write_clouds(tmp_output_dir)
        dtm = create_dtm_fast(ground, self.BASENAME, tmp_output_dir, 50.0,
                              force=True, source_laz=source, strip_align=False)
        assert dtm is not None
        with rasterio.open(str(dtm)) as src:
            assert EDGE_BUFFER_TAG not in src.tags()
        assert read_dtm_edge_buffer(dtm) == 0.0

    def test_buffered_dtm_non_lhd_falls_back(self, tmp_output_dir):
        """Nom hors pattern LHD : pas de tuile nominale, bornes d'en-tête conservées."""
        from lidar_pipeline.dtm import create_dtm_fast, read_dtm_edge_buffer
        import rasterio
        ground = self._write_las(tmp_output_dir / "ground.las",
                                 [(0.5, 0.5, 10.0), (0.05, 0.05, 10.0), (1.95, 1.95, 10.0)])
        source = self._write_las(tmp_output_dir / "nuage.laz",
                                 [(0.5, 0.5, 10.0), (0.05, 0.05, 10.0), (1.95, 1.95, 10.0)])
        dtm = create_dtm_fast(ground, "nuage", tmp_output_dir, 1.0,
                              force=True, source_laz=source, strip_align=False,
                              edge_buffer=100.0)
        assert dtm is not None
        with rasterio.open(str(dtm)) as src:
            assert abs(src.bounds.left - 0.05) < 1e-6  # bornes de l'en-tête
        assert read_dtm_edge_buffer(dtm) == 0.0


def _synthetic_beam(n_lines=240, spacing=0.4, seed=0, offsets=None, tilts=None):
    """Faisceau synthétique : lignes de balayage (dents de scie de scan_angle)
    sur un terrain en pente traversé par un fossé ; offsets verticaux par ligne."""
    rng = np.random.default_rng(seed)
    xs, ys, zs, ts, angs, ids = [], [], [], [], [], []
    x_line = np.arange(0, 100, 0.12)
    for k in range(n_lines):
        y = k * spacing + rng.normal(0, 0.02, x_line.size)
        z = 100 + 0.03 * x_line + 0.05 * y
        z = z - 0.5 * (np.abs(x_line - 41) < 1.0)          # fossé perpendiculaire aux lignes
        u = np.linspace(-1, 1, x_line.size)
        z = z + offsets[k] + (0 if tilts is None else tilts[k]) * u + rng.normal(0, 0.01, x_line.size)
        xs.append(x_line); ys.append(y); zs.append(z); ids.append(np.full(x_line.size, k))
        ts.append(k * 0.0067 + np.linspace(0, 0.006, x_line.size))
        angs.append(np.linspace(-3300, 3300, x_line.size))
    return (np.concatenate(xs), np.concatenate(ys), np.concatenate(zs),
            np.concatenate(ts), np.concatenate(angs), np.concatenate(ids))


class TestScanLineAlignment:
    """3ᵉ passe du calage : décalage vertical entre lignes de balayage successives."""

    def test_line_ids_follow_sawtooth_not_vegetation_gaps(self):
        from lidar_pipeline.dtm import _scan_line_ids
        t = np.arange(50) * 0.0001
        ang = np.tile(np.linspace(-3000, 3000, 10), 5)
        keep = np.ones(50, bool); keep[13:17] = False        # trou de végétation dans la ligne 2
        ids = _scan_line_ids(t[keep], ang[keep])
        assert ids.max() == 4
        t2 = t.copy(); t2[30:] += 1.0                          # fin de passe : nouvelle ligne
        assert _scan_line_ids(t2, np.zeros(50)).max() == 1

    def test_group_median_matches_numpy(self):
        from lidar_pipeline.dtm import _group_median
        rng = np.random.default_rng(3)
        g = rng.integers(0, 20, 2000)
        v = rng.normal(size=2000); v[::17] = np.nan
        med = _group_median(v, g, 20, 1)
        for k in range(20):
            vals = v[(g == k) & np.isfinite(v)]
            assert med[k] == pytest.approx(np.median(vals))

    def test_removes_alternating_line_offsets_and_keeps_ditch(self):
        from lidar_pipeline.dtm import _scan_line_corrections_beam
        n = 240
        rng = np.random.default_rng(1)
        offsets = 0.015 * (-1.0) ** np.arange(n) + rng.normal(0, 0.006, n)
        x, y, z, t, ang, ids = _synthetic_beam(n, offsets=offsets)
        corr, per_line, _ = _scan_line_corrections_beam(x, y, z, t, ang)
        assert len(per_line) == n
        from scipy.ndimage import gaussian_filter1d
        hp = lambda v: v - gaussian_filter1d(v, 3, mode="nearest")
        before, after = hp(offsets)[10:-10].std(), hp(offsets - per_line)[10:-10].std()
        assert after < 0.2 * before, f"{before*1000:.1f} → {after*1000:.1f} mm"
        ditch = np.abs(x - 41) < 0.8
        depth = lambda zz: np.median(zz[~ditch & (np.abs(x - 41) < 4)]) - np.median(zz[ditch])
        assert depth(z - corr) == pytest.approx(depth(z), abs=0.01)

    def test_removes_alternating_line_tilts(self):
        """Roulis : lignes basculées alternativement (un bout haut, l'autre bas)."""
        from scipy.ndimage import gaussian_filter1d
        from lidar_pipeline.dtm import _scan_line_corrections_beam
        n = 240
        rng = np.random.default_rng(4)
        tilts = 0.02 * (-1.0) ** np.arange(n) + rng.normal(0, 0.008, n)
        x, y, z, t, ang, ids = _synthetic_beam(n, offsets=np.zeros(n), tilts=tilts)
        corr, _, per_tilt = _scan_line_corrections_beam(x, y, z, t, ang)
        hp = lambda v: v - gaussian_filter1d(v, 3, mode="nearest")
        before, after = hp(tilts)[10:-10].std(), hp(tilts - per_tilt)[10:-10].std()
        assert after < 0.2 * before, f"{before*1000:.1f} → {after*1000:.1f} mm"

    def test_joint_adjustment_fixes_all_scales_against_other_beam(self):
        """Deux faisceaux superposés : le faisceau fautif (dérive lente, roulis,
        ligne isolée à −8 cm) est recalé sur l'autre à toutes les échelles,
        sans dérive de l'altitude d'ensemble."""
        from lidar_pipeline.dtm import _joint_line_corrections
        n = 240
        k = np.arange(n)
        bad_off = 0.03 * np.sin(2 * np.pi * k / 120)          # dérive lente (non vue sur sa propre surface)
        bad_off[100] -= 0.08                                   # ligne isolée très décalée
        bad_tilt = 0.02 * np.cos(2 * np.pi * k / 60)           # roulis lent
        xa, ya, za, ta, aa, _ = _synthetic_beam(n, seed=1, offsets=bad_off, tilts=bad_tilt)
        xb, yb, zb, tb, ab, _ = _synthetic_beam(n, seed=2, offsets=np.zeros(n))
        tb = tb + 1000.0
        x = np.r_[xa, xb]; y = np.r_[ya, yb]; z = np.r_[za, zb]; t = np.r_[ta, tb]
        ang = np.r_[aa, ab]; psid = np.r_[np.full(len(za), 1), np.full(len(zb), 2)]
        corr, gl, touched, iters = _joint_line_corrections(x, y, z, psid, t, ang)
        truth = np.r_[bad_off[np.repeat(k, len(za) // n)] + bad_tilt[np.repeat(k, len(za) // n)]
                      * np.tile(np.linspace(-1, 1, len(za) // n), n), np.zeros(len(zb))]
        # Sans vérité terrain, l'écart est partagé entre les faisceaux : c'est
        # l'écart ENTRE faisceaux (points homologues, même géométrie) qui doit
        # disparaître.
        na = len(za)
        before = truth[:na] - truth[:na].mean()
        rel = (truth[:na] - corr[:na]) - (0.0 - corr[na:])
        after = rel - rel.mean()
        assert np.std(after) < 0.25 * np.std(before), \
            f"{np.std(before)*1000:.1f} → {np.std(after)*1000:.1f} mm"
        assert abs(np.mean(corr)) < 0.002                     # pas de dérive d'ensemble
        assert touched.mean() > 0.9 and iters <= 8

    def test_joint_adjustment_removes_static_angle_profile(self):
        """Étalonnage en arc selon l'angle (même pour toutes les lignes) :
        non linéaire, invisible pour décalage + inclinaison, retiré par le
        profil par faisceau et classe d'angle."""
        from lidar_pipeline.dtm import _joint_line_corrections
        n = 240
        xa, ya, za, ta, aa, _ = _synthetic_beam(n, seed=5, offsets=np.zeros(n))
        uu = aa / 3300.0
        prof = 0.02 * (uu ** 2 - np.mean(uu ** 2))
        za = za + prof
        xb, yb, zb, tb, ab, _ = _synthetic_beam(n, seed=6, offsets=np.zeros(n))
        x = np.r_[xa, xb]; y = np.r_[ya, yb]; z = np.r_[za, zb]; t = np.r_[ta, tb + 1000.0]
        ang = np.r_[aa, ab]; psid = np.r_[np.full(len(za), 1), np.full(len(zb), 2)]
        corr, _, _, _ = _joint_line_corrections(x, y, z, psid, t, ang)
        na = len(za)
        rel = (prof - corr[:na]) - (0.0 - corr[na:])
        assert np.std(rel - rel.mean()) < 0.3 * np.std(prof), \
            f"{np.std(prof)*1000:.1f} → {np.std(rel - rel.mean())*1000:.1f} mm"

    def test_clean_beam_left_untouched(self):
        from lidar_pipeline.dtm import _scan_line_corrections
        x, y, z, t, ang, ids = _synthetic_beam(240, offsets=np.zeros(240))
        corr, stats = _scan_line_corrections(x, y, z, np.full(len(z), 7), t, ang)
        assert stats == {} and not corr.any()


class TestIgnDirectExtraction:
    """Classification IGN : extraction directe par laspy (PDAL en secours)."""

    def test_keeps_requested_classes_and_valid_returns(self, tmp_path):
        import laspy
        from lidar_pipeline.dtm import _extract_ign_ground
        n = 1000
        rng = np.random.default_rng(0)
        hdr = laspy.LasHeader(point_format=6, version="1.4")
        hdr.scales = np.array([0.01, 0.01, 0.001]); hdr.offsets = np.array([1000.0, 6800000.0, 0.0])
        las = laspy.LasData(hdr)
        las.x = 1000 + rng.uniform(0, 100, n); las.y = 6800000 + rng.uniform(0, 100, n)
        las.z = rng.uniform(100, 110, n)
        cls = rng.choice([1, 2, 3, 6, 9], n); las.classification = cls
        rn = np.ones(n, dtype=np.uint8); rn[:10] = 0; las.return_number = rn
        las.number_of_returns = np.ones(n, dtype=np.uint8)
        las.point_source_id = np.full(n, 42, dtype=np.uint16)
        src = tmp_path / "t.las"; las.write(str(src))
        out = tmp_path / "g.las"
        assert _extract_ign_ground(src, out, [2])
        g = laspy.read(str(out))
        expected = (cls == 2) & (rn >= 1)
        assert len(g.points) == int(expected.sum())
        assert set(np.unique(np.asarray(g.classification))) == {2}
        assert g.header.point_format.id == 6 and np.all(np.asarray(g.point_source_id) == 42)
        assert _extract_ign_ground(src, tmp_path / "e.las", [1, 2])
        assert len(laspy.read(str(tmp_path / "e.las")).points) == int(((np.isin(cls, [1, 2])) & (rn >= 1)).sum())