"""DTM generation from classified LiDAR point clouds.

Handles ground classification (IGN supplier pre-classification extracted
directly with laspy, PDAL as a fallback; or PDAL SMRF / CSF), vertical
alignment of the flight strips, and DTM rasterisation (GPU-first via
gpu.bin_mean_2d, scipy binned_statistic_2d as fallback). Gaps between points
are filled only within the envelope of the measured points; larger holes
without LiDAR data stay nodata.
"""

import json
import logging
import subprocess
import time
from pathlib import Path

import numpy as np
import rasterio
from rasterio.transform import from_bounds
from scipy.stats import binned_statistic_2d

from .gpu import bin_mean_2d

logger = logging.getLogger("lidar")

# Usable LAS classes of the LiDAR HD pre-classification (names → codes)
IGN_CLASS_NAMES = {
    "sol": 2,
    "unclassified": 1,
    "non-classe": 1,
    "eau": 9,
    "virtuel": 66,
    "pont": 17,
    "sursol": 64,
}


# Vertical alignment of flight strips (strip alignment). A LiDAR HD tile is
# covered by several acquisition passes, each carried by one or two
# PointSourceIds (strips). Passes are well aligned horizontally, but some
# strips carry a vertical bias of a few cm (measured up to ~5 cm on
# LHD_FXX_1000_6882: the two strips of one pass were ±2.5 cm apart). At
# 0.2 m/px these offsets create steps and noise along the seams of overlap
# areas. The robust offset of each strip is measured on the ground points of
# the tile itself (PointSourceIds change between campaigns, nothing is
# hard-coded) and subtracted before rasterisation. Offsets drift along a
# flight line (sign measured to flip between neighbouring tiles), so the
# computation is per tile, never global.
#
# 2nd pass: intra-strip jitter. The vertical offset can also vary along ONE
# pass (sensor vibration / high-frequency trajectory noise): successive scan
# lines of a strip then appear randomly shifted vertically. Each strip is cut
# into GPS time windows, the robust offset of each window is measured against
# the median surface of the other strips (same 1 m cell), the series is
# smoothed (rolling median) and then interpolated at each point's GPS time.
# Requires the gps_time dimension (silently skipped otherwise).
#
# 3rd pass: line-by-line offset. The 0.1 s windows (smoothed over 0.5 s)
# group ~15 scan lines (~150 lines/s) and only work in overlaps. Yet two
# SUCCESSIVE lines of one pass can differ by 1 to 2 cm (alternating pattern,
# measured on LHD_FXX_0999_6882): fine stripes perpendicular to the flight
# direction over the whole DTM. Each strip is cut into lines (sawtooth jumps
# of scan_angle), the robust offset of each line (offset AND tilt along the
# line: roll tilts the lines) is measured against the surface of ITS OWN
# strip (0.5 m cell, 1.5 m box, local plane fitted at the point's actual
# position, 3 iterations so the line does not bias its own reference), and
# only the line-to-line component is removed (series − Gaussian smoothing
# σ 3 lines; a rolling median would follow an alternating pattern instead of
# erasing it): slow variations are left to the previous passes. Works without
# overlap; requires gps_time and scan_angle.
#
# Joint adjustment (takes precedence): a strip's own surface absorbs any error
# wider than its reference box; where several strips overlap, each line
# (offset + tilt) is therefore re-aligned against the consensus of the OTHER
# strips, at every scale (slow roll of a whole pass, isolated strongly
# shifted lines), with damped, iterated steps. The correction against the
# strip's own surface is now only used for lines without overlap.
STRIP_ALIGN_VERSION = 3
STRIP_ALIGN_THRESHOLD = 0.005  # m: minimum offset to correct a strip (0.5 cm)
STRIP_ALIGN_CELL = 1.0         # m: strip comparison cell
STRIP_ALIGN_MIN_SHARED = 500   # minimum shared ground cells to validate an offset
STRIP_JITTER_BIN = 0.1         # s: length of a GPS time window (jitter)
STRIP_JITTER_SMOOTH = 5        # windows: width of the rolling median
STRIP_JITTER_MIN_CELLS = 40    # minimum shared ground cells to validate a window
STRIP_JITTER_MAX = 0.10        # m: maximum amplitude of a jitter correction (safeguard)
STRIP_LINE_CELL = 0.5          # m: cell of a strip's reference surface
STRIP_LINE_BOX = 3             # cells: reference smoothing (1.5 m box)
STRIP_LINE_WINDOW = 3          # lines: σ of the removed Gaussian smoothing (keeps line-to-line)
STRIP_LINE_ITERS = 3           # iterations (attenuation by the line itself < 1 %)
STRIP_LINE_MAX = 0.05          # m: maximum correction of a line (safeguard)
STRIP_LINE_MIN_POINTS = 100    # minimum ground points to measure a line
STRIP_LINE_GAP = 0.05          # s: time gap that breaks a line (end of pass)
STRIP_LINE_MIN_RMS = 0.002     # m: strip left untouched below this level
STRIP_LINE_MODEL = "conjoint+decalage+inclinaison+profil-angle"  # model (recorded, invalidates the cache)
_SCAN_ANGLE_UNIT = 0.006        # ° per scan_angle unit (LAS 1.4)
STRIP_ANGLE_BIN = 1.0          # °: step of the per-strip correction profile by angle
STRIP_ANGLE_MIN_POINTS = 500   # minimum overlap points per angle bin (else neighbour value)
STRIP_LINE_SUBSAMPLE = 3       # estimate on 1 point in 3 (correction applied to all)
STRIP_JOINT_CELL = 1.0         # m: cell of the other strips' consensus
STRIP_JOINT_ITERS = 8          # maximum iterations of the joint adjustment
STRIP_JOINT_DAMPING = 0.5      # damped step: two strips converge without crossing
STRIP_JOINT_TOL = 0.001        # m: stop when the mean step falls below 1 mm
STRIP_JOINT_MIN_POINTS = 30    # minimum overlap points to re-align a line
STRIP_JOINT_MAX = 0.15         # m: maximum correction of a point (safeguard)

# Per-file memo of offsets: classification is shared between resolutions,
# the same ground LAS is rasterised at 0.5 m and then 0.2 m.
_STRIP_OFFSETS_CACHE = {}
_STRIP_JITTER_CACHE = {}
_STRIP_LINES_CACHE = {}


def _strip_surface_grid(x, y, z, inv, n_sources, cell=STRIP_ALIGN_CELL):
    """Per-strip ground surfaces on a 1 m grid (mean of the low points).

    For each strip and each cell: mean of the points lying within 0.5 m of
    the strip's minimum in that cell (robust to residual vegetation), cells
    with ≥ 2 points only.

    Returns:
        (allcells, grid, key): sorted cells, grid[strip, cell] = mean Z (NaN
        if absent) and the cell key of EVERY point; None if no surface is
        populated.
    """
    x0 = np.floor(np.min(x) / cell) * cell
    y0 = np.floor(np.min(y) / cell) * cell
    xi = ((x - x0) / cell).astype(np.int64)
    yi = ((y - y0) / cell).astype(np.int64)
    ny = int(yi.max()) + 1
    key = xi * ny + yi

    def _surface(k):
        m = inv == k
        kk, zz = key[m], z[m]
        order = np.argsort(kk, kind='stable')
        k_s, z_s = kk[order], zz[order]
        starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]])
        mins = np.minimum.reduceat(z_s, starts)
        ukey = k_s[starts]
        thr = mins[np.searchsorted(ukey, kk)]
        sel = np.flatnonzero(zz <= thr + 0.5)
        k2, z2 = kk[sel], zz[sel]
        cnt = np.bincount(k2)
        sums = np.bincount(k2, weights=z2)
        v = np.flatnonzero(cnt >= 2)
        return v, sums[v] / cnt[v]

    surfaces = [_surface(k) for k in range(n_sources)]
    populated = [c for c, _ in surfaces if len(c)]
    if not populated:
        return None
    allcells = np.unique(np.concatenate(populated))
    grid = np.full((n_sources, len(allcells)), np.nan)
    for k, (c, zs) in enumerate(surfaces):
        grid[k, np.searchsorted(allcells, c)] = zs
    return allcells, grid, key


def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL,
                            threshold=STRIP_ALIGN_THRESHOLD,
                            min_shared=STRIP_ALIGN_MIN_SHARED):
    """Measure the relative vertical offsets between the strips of a tile.

    Method: ground surface per strip (mean of the points within 0.5 m of the
    minimum of each 1 m cell, robust to residual vegetation), then the offset
    of each strip = median of its deviation from the median reference
    surface (iterated 3 times), on cells covered by at least two strips only.
    Corrects nothing but the vertical: the horizontal alignment measured on
    the data is excellent (≤ 1 cm).

    Args:
        x, y, z, psid: coordinates and PointSourceId of the ground points.
        cell: comparison cell size (m).
        threshold: threshold (m) below which an offset is ignored.
        min_shared: minimum shared cells to validate a strip.

    Returns:
        dict {psid: offset} of the offsets to SUBTRACT (z - offset), holding
        only |offset| >= threshold; empty if there is nothing to correct
        (single strip, no overlap, tile already aligned).
    """
    us, inv = np.unique(psid, return_inverse=True)
    if len(us) < 2:
        return {}
    built = _strip_surface_grid(x, y, z, inv, len(us), cell)
    if built is None:
        return {}
    allcells, grid, _key = built
    comparable = np.sum(~np.isnan(grid), axis=0) >= 2

    offsets = np.zeros(len(us))
    for _ in range(3):
        ref = np.nanmedian(grid + offsets[:, None], axis=0)
        for k in range(len(us)):
            m = ~np.isnan(grid[k]) & comparable
            if int(m.sum()) >= min_shared:
                offsets[k] = np.median(grid[k][m] - ref[m])

    return {int(p): round(float(offsets[k]), 3)
            for k, p in enumerate(us) if abs(offsets[k]) >= threshold}


def _rolling_median(values, window):
    """Rolling median (window truncated at the edges) over a small vector."""
    n = len(values)
    if window <= 1 or n == 0:
        return np.array(values, dtype=np.float64)
    half = window // 2
    return np.array([np.median(values[max(0, i - half):i + half + 1])
                     for i in range(n)])


def _strip_jitter_offsets(x, y, z, psid, t, cell=STRIP_ALIGN_CELL,
                          bin_seconds=STRIP_JITTER_BIN,
                          smooth=STRIP_JITTER_SMOOTH,
                          min_cells=STRIP_JITTER_MIN_CELLS,
                          max_corr=STRIP_JITTER_MAX,
                          threshold=STRIP_ALIGN_THRESHOLD):
    """Measure intra-strip vertical jitter over GPS time windows.

    The constant alignment removes one offset per strip, but the vertical
    offset can also vary along a single pass (sensor vibration,
    high-frequency trajectory noise). Each strip is cut into GPS time
    windows; the robust offset of each window is measured against the median
    surface of the OTHER strips (1 m cell, same low-point selection as the
    constant alignment; with a two-strip overlap, a reference including the
    tested strip would reveal only half of the offset), then the series is
    smoothed (rolling median) so as not to follow measurement noise. z must
    already be corrected for the constant offsets.

    Args:
        x, y, z, psid: coordinates, Z (constant-aligned) and PointSourceId
            of the ground points.
        t: GPS time (s) of each point.
        cell: comparison cell size (m).
        bin_seconds: length of a time window (s).
        smooth: width of the rolling median (windows).
        min_cells: minimum shared cells to validate a window.
        max_corr: maximum amplitude of a correction (safeguard, m).
        threshold: series amplitude below which a strip is considered
            stable (no correction).

    Returns:
        dict {psid: (times, corrections)} of the corrections to SUBTRACT,
        linearly interpolable at each point's GPS time; empty if nothing is
        measurable (single strip, no time, no overlap).
    """
    us, inv = np.unique(psid, return_inverse=True)
    if len(us) < 2:
        return {}
    t = np.asarray(t, dtype=np.float64)
    if t.size != z.size or not np.isfinite(t).all():
        return {}
    built = _strip_surface_grid(x, y, z, inv, len(us), cell)
    if built is None:
        return {}
    allcells, grid, key = built
    covered = np.sum(~np.isnan(grid), axis=0) >= 2

    # Reference of a strip = median of the other strips.
    ref = np.full(grid.shape, np.nan)
    if len(us) == 2:
        ref[0], ref[1] = grid[1], grid[0]
    else:
        import warnings
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", RuntimeWarning)  # all-NaN slices
            for k in range(len(us)):
                ref[k] = np.nanmedian(np.delete(grid, k, axis=0), axis=0)

    # Vertical residual of each point against its strip's reference.
    idx = np.minimum(np.searchsorted(allcells, key), len(allcells) - 1)
    ref_point = ref[inv, idx]
    hit = (allcells[idx] == key) & covered[idx] & np.isfinite(ref_point)
    residual = np.full(np.asarray(z).shape, np.nan)
    residual[hit] = np.asarray(z)[hit] - ref_point[hit]

    # Per-strip time origin: passes over one tile can be hours apart, so the
    # windows stay dense around the actual flight (and a GPS week rollover
    # does not create huge empty ranges).
    origin = np.full(len(us), np.inf)
    np.minimum.at(origin, inv, t)
    tb = np.floor((t - origin[inv]) / bin_seconds).astype(np.int64)
    n_bins = int(tb.max()) + 1
    group = inv * n_bins + tb

    # Robust median per (strip, window) group: finite residuals are sorted to
    # the head of each segment, NaNs (no reference) are ignored.
    order = np.lexsort((np.where(np.isfinite(residual), residual, np.inf), group))
    g_s, r_s = group[order], residual[order]
    starts = np.flatnonzero(np.r_[True, g_s[1:] != g_s[:-1]])
    ends = np.r_[starts[1:], len(g_s)]

    gids, times_c, med = [], [], []
    for s, e in zip(starts, ends):
        finite = np.isfinite(r_s[s:e])
        if int(finite.sum()) < min_cells:
            continue
        gids.append(int(g_s[s]))
        med.append(float(np.median(r_s[s:e][finite])))
    if not gids:
        return {}

    gids = np.asarray(gids, dtype=np.int64)
    med = np.asarray(med, dtype=np.float64)
    result = {}
    for k in range(len(us)):
        selk = gids // n_bins == k
        if int(selk.sum()) < 3:  # series too short: unreliable correction
            continue
        b = gids[selk] % n_bins
        cs = np.clip(_rolling_median(med[selk], smooth), -max_corr, max_corr)
        if float(np.max(np.abs(cs))) < threshold:
            continue
        result[int(us[k])] = (origin[k] + (b + 0.5) * bin_seconds, cs)
    return result


def _module_of(arr):
    """Array module (numpy or cupy) of an array."""
    mod = type(arr).__module__
    if mod.startswith("cupy"):
        import cupy
        return cupy
    return np


def _array_module(use_gpu):
    """cupy if the GPU is active and requested, numpy otherwise."""
    if use_gpu:
        from . import gpu as _gpu
        if _gpu.is_gpu_active() and _gpu._cp is not None:
            return _gpu._cp
    return np


def _scan_line_ids(t, angle, gap=STRIP_LINE_GAP):
    """Scan-line id of the points of ONE strip, sorted by time.

    A new line starts at every scan_angle jump larger than half the
    amplitude in the direction opposite to the sweep (sawtooth flyback) or
    at every time gap > gap (end of pass). Holes left by removed vegetation
    inside a line do not break it.
    """
    m = _module_of(angle)
    angle = m.asarray(angle, dtype=m.float64)
    if len(angle) == 0:
        return m.zeros(0, dtype=m.int64)
    amp = float(m.percentile(angle, 99) - m.percentile(angle, 1))
    d = m.diff(angle)
    moving = d[d != 0]
    sweep = float(m.sign(m.median(moving))) if len(moving) else 1.0
    # Sawtooth flyback: a large jump OPPOSITE to the sweep direction. A
    # vegetation hole also makes the angle jump, but in the sweep direction:
    # it does not break the line.
    new = m.concatenate([m.ones(1, dtype=bool), (d * sweep < -0.5 * amp) | (m.diff(t) > gap)])
    return m.cumsum(new) - 1


def _group_median(values, groups, n_groups, min_count):
    """Median of the finite values per group (vectorised), NaN below min_count."""
    finite = np.isfinite(values)
    order = np.lexsort((np.where(finite, values, np.inf), groups))
    g = groups[order]
    v = values[order]
    nf = np.bincount(groups[finite], minlength=n_groups)
    start = np.searchsorted(g, np.arange(n_groups))
    med = np.full(n_groups, np.nan)
    ok = nf >= max(1, min_count)
    lo = start[ok] + (nf[ok] - 1) // 2
    hi = start[ok] + nf[ok] // 2
    med[ok] = 0.5 * (v[lo] + v[hi])
    return med


def _line_fit(r, w, line, u, n_lines, min_points):
    """Weighted least squares per line: r ≈ a + b·u (sums via bincount).

    Returns:
        (a, b, ok); b is 0 where the spread of u does not allow a tilt to be
        estimated (a alone, weighted mean).
    """
    m = _module_of(r)
    rw = m.where(w, r, 0.0)
    wf = w.astype(m.float64)
    s0 = m.bincount(line, weights=wf, minlength=n_lines)
    s1 = m.bincount(line, weights=wf * u, minlength=n_lines)
    s2 = m.bincount(line, weights=wf * u * u, minlength=n_lines)
    sr = m.bincount(line, weights=rw, minlength=n_lines)
    sur = m.bincount(line, weights=rw * u, minlength=n_lines)
    det = s0 * s2 - s1 * s1
    ok = s0 >= min_points
    tilt = ok & (det > 1e-2 * m.maximum(s0, 1) ** 2)
    a = m.where(ok, sr / m.maximum(s0, 1), 0.0)
    b = m.zeros(n_lines)
    safe = m.where(tilt, det, 1.0)
    a = m.where(tilt, (s2 * sr - s1 * sur) / safe, a)
    b = m.where(tilt, (s0 * sur - s1 * sr) / safe, b)
    return a, b, ok


def _robust_mask(r, floor=0.03, k=5.0):
    """Finite residuals below k MAD (with a floor): rejects low vegetation,
    misclassified points and hole edges without sorting each line."""
    m = _module_of(r)
    f = m.isfinite(r)
    if not bool(f.any()):
        return f
    mad = 1.4826 * float(m.median(m.abs(r[f])))
    return f & (m.abs(r) < max(floor, k * mad))


def _scan_line_corrections_beam(x, y, z, t, angle, cell=STRIP_LINE_CELL,
                                box=STRIP_LINE_BOX, window=STRIP_LINE_WINDOW,
                                iters=STRIP_LINE_ITERS, max_corr=STRIP_LINE_MAX,
                                min_points=STRIP_LINE_MIN_POINTS,
                                subsample=STRIP_LINE_SUBSAMPLE):
    """Line-by-line corrections of a strip against ITS OWN surface (to SUBTRACT).

    Each line is modelled by an offset AND a tilt along the line (a + b·u,
    u = normalised scan_angle): a roll error tilts the line, one end higher
    than the other. Only the line-to-line component is removed (series −
    Gaussian of σ window lines). Estimated on 1 point in subsample, with
    vectorised trimmed least squares (5 MAD).

    Returns:
        (corr per point in input order, offsets per line, tilts per line at
        the swath edge).
    """
    from scipy.ndimage import gaussian_filter1d, uniform_filter
    n = len(z)
    o = np.argsort(t, kind="stable")
    line_all = np.empty(n, dtype=np.int64)
    line_all[o] = _scan_line_ids(t[o], angle[o])
    nl = int(line_all.max()) + 1 if n else 0
    tot_a, tot_b = np.zeros(nl), np.zeros(nl)
    if nl < 10 * window:
        return np.zeros(n), tot_a, tot_b
    u_all = np.asarray(angle, dtype=np.float64) / max(np.percentile(np.abs(angle), 99), 1e-9)
    sel = np.arange(n) % max(1, int(subsample)) == 0
    xs, ys, zs = x[sel], y[sel], np.asarray(z, dtype=np.float64)[sel]
    line, u = line_all[sel], u_all[sel]
    # Reference = local plane: (x, y, z) means per box of cells, slope of the
    # smoothed surface; evaluated at the point's actual position (a cell's
    # mean is not at its centre: on a slope, interpolating at the centre
    # creates a bias that depends on where the line lies).
    x0, y0 = xs.min(), ys.min()
    ix = np.floor((xs - x0) / cell).astype(np.int64)
    iy = np.floor((ys - y0) / cell).astype(np.int64)
    W, H = int(ix.max()) + 1, int(iy.max()) + 1
    flat = iy * W + ix
    count_b = uniform_filter(np.bincount(flat, minlength=W * H).reshape(H, W).astype(np.float64), box)
    valid = count_b > 0
    inv_c = np.where(valid, 1.0 / np.maximum(count_b, 1e-12), np.nan)
    mean_x = uniform_filter(np.bincount(flat, weights=xs - x0, minlength=W * H).reshape(H, W), box) * inv_c
    mean_y = uniform_filter(np.bincount(flat, weights=ys - y0, minlength=W * H).reshape(H, W), box) * inv_c
    dxp = (xs - x0) - mean_x.ravel()[flat]
    dyp = (ys - y0) - mean_y.ravel()[flat]
    min_pts = max(10, min_points // max(1, int(subsample)))
    z_work = zs.copy()
    for _ in range(iters):
        mean_z = uniform_filter(np.bincount(flat, weights=z_work, minlength=W * H).reshape(H, W), box) * inv_c
        gz_y, gz_x = np.gradient(np.where(valid, mean_z, np.nanmean(mean_z)), cell)
        r = z_work - (mean_z.ravel()[flat] + gz_x.ravel()[flat] * dxp + gz_y.ravel()[flat] * dyp)
        a, b, ok = _line_fit(r, _robust_mask(r), line, u, nl, min_pts)
        if ok.sum() < 10 * window:
            break
        idx = np.flatnonzero(ok)
        a = np.interp(np.arange(nl), idx, a[ok])
        b = np.interp(np.arange(nl), idx, b[ok])
        # Only the line-to-line component is removed (a rolling median would
        # follow an alternating pattern instead of erasing it)
        da = a - gaussian_filter1d(a, window, mode="nearest")
        db = b - gaussian_filter1d(b, window, mode="nearest")
        da[~ok] = 0.0
        db[~ok] = 0.0
        tot_a += da
        tot_b += db
        z_work = zs - np.clip(tot_a[line] + tot_b[line] * u, -max_corr, max_corr)
    corr = np.clip(tot_a[line_all] + tot_b[line_all] * u_all, -max_corr, max_corr)
    return corr, tot_a, tot_b


def _joint_line_corrections(x, y, z, psid, t, angle, cell=STRIP_JOINT_CELL,
                            iters=STRIP_JOINT_ITERS, damping=STRIP_JOINT_DAMPING,
                            tol=STRIP_JOINT_TOL, min_points=STRIP_JOINT_MIN_POINTS,
                            subsample=STRIP_LINE_SUBSAMPLE, max_corr=STRIP_JOINT_MAX,
                            use_gpu=True):
    """Joint adjustment of the lines of all strips against the consensus of
    the OTHER strips (offset + tilt per line), plus a per-strip profile by
    scan-angle bin.

    At each iteration, a point's residual is measured against the mean of
    the other strips in its cell (brought to its position by the surface
    slope), each line is fitted by trimmed least squares, and a damping
    fraction of the step is applied to all lines at once; the mean elevation
    is re-centred (no overall drift). Runs on GPU (CuPy) when available,
    falling back to numpy on any error.

    Returns:
        (corr to SUBTRACT per point, global line id per point, re-aligned
        lines (bool per line), number of iterations), as numpy arrays.
    """
    m = _array_module(use_gpu)
    if m is not np:
        try:
            return _joint_line_corrections_impl(m, x, y, z, psid, t, angle, cell, iters,
                                                damping, tol, min_points, subsample, max_corr)
        except Exception as e:
            logger.warning(f"    GPU joint adjustment failed ({e}), falling back to CPU")
    return _joint_line_corrections_impl(np, x, y, z, psid, t, angle, cell, iters,
                                        damping, tol, min_points, subsample, max_corr)


def _joint_line_corrections_impl(m, x, y, z, psid, t, angle, cell, iters, damping,
                                 tol, min_points, subsample, max_corr):
    def host(a):
        return a.get() if m is not np else a
    n = len(z)
    psid_h = np.asarray(psid)
    beams, bidx_h = np.unique(psid_h, return_inverse=True)
    # Angle bin per strip: calibration profile shared by all lines of a strip
    # (non-linear deviation across the swath).
    abin_h = np.round(np.asarray(angle, dtype=np.float64) * _SCAN_ANGLE_UNIT / STRIP_ANGLE_BIN).astype(np.int64)
    amin = int(abin_h.min()) if n else 0
    n_ab = int(abin_h.max()) - amin + 1 if n else 1
    abin_h = bidx_h * n_ab + (abin_h - amin)
    x, y = m.asarray(x, dtype=m.float64), m.asarray(y, dtype=m.float64)
    z, t = m.asarray(z, dtype=m.float64), m.asarray(t, dtype=m.float64)
    angle = m.asarray(angle, dtype=m.float64)
    bidx = m.asarray(bidx_h)
    gl = m.zeros(n, dtype=m.int64)
    u = m.zeros(n)
    base = 0
    for i in range(len(beams)):
        idx = m.flatnonzero(bidx == i)
        o = m.argsort(t[idx])
        gl[idx[o]] = _scan_line_ids(t[idx][o], angle[idx][o]) + base
        base = int(gl[idx].max()) + 1
        u[idx] = angle[idx] / max(float(m.percentile(m.abs(angle[idx]), 99)), 1e-9)
    n_lines = base
    touched = m.zeros(n_lines, dtype=bool)
    if len(beams) < 2 or n == 0:
        return np.zeros(n), host(gl), host(touched), 0
    abin = m.asarray(abin_h)
    n_bins = len(beams) * n_ab
    sel = m.arange(0, n, max(1, int(subsample)))
    xs, ys, zs = x[sel], y[sel], z[sel]
    gs, us, bs, cs = gl[sel], u[sel], bidx[sel], abin[sel]
    x0, y0 = float(xs.min()), float(ys.min())
    ix = m.floor((xs - x0) / cell).astype(m.int64)
    iy = m.floor((ys - y0) / cell).astype(m.int64)
    W, H = int(ix.max()) + 1, int(iy.max()) + 1
    nk = W * H
    key = iy * W + ix
    bkey = bs * nk + key
    dxc = (xs - x0) - (ix + 0.5) * cell
    dyc = (ys - y0) - (iy + 0.5) * cell
    n_all = m.bincount(key, minlength=nk)
    n_own = m.bincount(bkey, minlength=len(beams) * nk)
    n_other = n_all[key] - n_own[bkey]
    overlap = n_other >= 2
    min_pts = max(10, min_points // max(1, int(subsample)))
    A = m.zeros(n_lines)
    B = m.zeros(n_lines)
    C = m.zeros(n_bins)
    min_bin = max(30, STRIP_ANGLE_MIN_POINTS // max(1, int(subsample)))
    it = 0
    for it in range(1, iters + 1):
        zc = zs - m.clip(A[gs] + B[gs] * us + C[cs], -max_corr, max_corr)
        s_all = m.bincount(key, weights=zc, minlength=nk)
        s_own = m.bincount(bkey, weights=zc, minlength=len(beams) * nk)
        mean = m.where(n_all > 0, s_all / m.maximum(n_all, 1), m.nan).reshape(H, W)
        gy, gx = m.gradient(m.where(m.isfinite(mean), mean, m.nanmean(mean)), cell)
        ref = ((s_all[key] - s_own[bkey]) / m.maximum(n_other, 1)
               + gx.ravel()[key] * dxc + gy.ravel()[key] * dyc)
        r = m.where(overlap, zc - ref, m.nan)
        w = m.zeros(len(r), dtype=bool)
        for i in range(len(beams)):
            mb = bs == i
            w[mb] = _robust_mask(r[mb])
        a, b, ok = _line_fit(r, w, gs, us, n_lines, min_pts)
        touched |= ok
        # Profile per strip and angle bin: trimmed mean of the residual left
        # after the line step; its mean (carried by the lines) is removed
        # strip by strip. Its slope is kept: when the swath only partly
        # crosses the tile, the line tilt is undetermined and only the
        # profile can carry it.
        r2 = m.where(w, r - (a[gs] + b[gs] * us), 0.0)
        wf = w.astype(m.float64)
        cnt = m.bincount(cs, weights=wf, minlength=n_bins)
        prof = m.where(cnt >= min_bin, m.bincount(cs, weights=r2, minlength=n_bins) / m.maximum(cnt, 1), 0.0)
        su = m.bincount(cs, weights=wf * us, minlength=n_bins)
        ub = m.where(cnt > 0, su / m.maximum(cnt, 1), 0.0)
        for i in range(len(beams)):
            sl = slice(i * n_ab, (i + 1) * n_ab)
            k = cnt[sl] >= min_bin
            if int(k.sum()) >= 3:
                wk = cnt[sl][k]
                uk, pk = ub[sl][k], prof[sl][k]
                um, pm = float((wk * uk).sum() / wk.sum()), float((wk * pk).sum() / wk.sum())
                fitted = prof[sl] - pm
                # Bin too sparse (swath edge, incomplete bin): value of the
                # nearest valid bin, not zero; the deviation is precisely
                # strongest at the edge.
                idx = m.arange(n_ab, dtype=m.float64)
                prof[sl] = m.interp(idx, idx[k], fitted[k])
            else:
                prof[sl] = 0.0
        A += damping * a
        B += damping * b
        C += damping * prof
        A -= float(m.mean(A[gs] + B[gs] * us + C[cs]))      # no overall drift
        step_lines = float(m.mean(m.abs(a[ok]))) * damping if bool(ok.any()) else 0.0
        step_prof = float(m.max(m.abs(prof))) * damping
        if max(step_lines, step_prof) < tol:
            break
    corr = m.clip(A[gl] + B[gl] * u + C[abin], -max_corr, max_corr)
    return host(corr), host(gl), host(touched), it


def _rolling_median_fast(values, window):
    """Centred rolling median (edges repeated), vectorised."""
    from numpy.lib.stride_tricks import sliding_window_view
    half = window // 2
    return np.median(sliding_window_view(np.pad(values, half, mode="edge"), window), axis=1)


def _scan_line_corrections(x, y, z, psid, t, angle, min_rms=STRIP_LINE_MIN_RMS):
    """Line-by-line corrections of all strips of a tile.

    1. Joint adjustment against the other strips (lines in overlaps).
    2. Lines without overlap: line-to-line component against the surface of
       their own strip (only when more than 5 % of the strip's points are in
       such lines).
    z must already be aligned (constant offsets and, if computed, time
    jitter).

    Returns:
        (corrections to SUBTRACT per point, {psid: (lines, rms, max)} of the
        corrected strips).
    """
    psid = np.asarray(psid)
    z = np.asarray(z, dtype=np.float64)
    corr, gl, touched, _ = _joint_line_corrections(x, y, z, psid, t, angle)
    stats = {}
    for p in np.unique(psid):
        m = psid == p
        if int(m.sum()) < 50 * STRIP_LINE_MIN_POINTS:
            continue
        alone = ~touched[gl[m]]
        if alone.mean() > 0.05:
            c_self, per_line, _ = _scan_line_corrections_beam(
                x[m], y[m], z[m] - corr[m], t[m], angle[m])
            c_beam = corr[m] + np.where(alone, c_self, 0.0)
        else:
            c_beam = corr[m]
        rms = float(np.sqrt(np.mean(c_beam ** 2))) if m.any() else 0.0
        if rms < min_rms:
            corr[m] = 0.0
            continue
        corr[m] = c_beam
        n_lines = int(len(np.unique(gl[m])))
        stats[int(p)] = (n_lines, rms, float(np.max(np.abs(c_beam))))
    return corr, stats


def _scan_angle(las):
    """scan_angle (LAS 1.4) or scan_angle_rank (LAS ≤ 1.3), None if absent."""
    for name in ("scan_angle", "scan_angle_rank"):
        try:
            return np.asarray(getattr(las, name), dtype=np.float64)
        except AttributeError:
            continue
    return None


def _scan_lines_for_file(las_file, las, z_aligned, t):
    """Line-by-line corrections of a ground LAS, memoised by (path, mtime)."""
    try:
        p = Path(las_file)
        cache_key = (str(p), p.stat().st_mtime_ns)
    except OSError:
        cache_key = (str(las_file), 0)
    if cache_key in _STRIP_LINES_CACHE:
        return _STRIP_LINES_CACHE[cache_key]
    angle = _scan_angle(las)
    result = (np.zeros(len(z_aligned)), {})
    if angle is not None and len(angle) == len(z_aligned):
        result = _scan_line_corrections(
            np.asarray(las.x, dtype=np.float64), np.asarray(las.y, dtype=np.float64),
            z_aligned, np.asarray(las.point_source_id), t, angle)
    _STRIP_LINES_CACHE[cache_key] = result
    return result


def _apply_strip_jitter(psid, t, jitter):
    """Jitter corrections interpolated at each point's GPS time.

    Args:
        psid, t: PointSourceId and GPS time (s) of each point.
        jitter: dict {psid: (times, corrections)} from _strip_jitter_offsets.

    Returns:
        array of the corrections to SUBTRACT (0 for points without a series).
    """
    corr = np.zeros(len(t))
    if not jitter:
        return corr
    psid = np.asarray(psid)
    t = np.asarray(t, dtype=np.float64)
    for p, (times, cs) in jitter.items():
        m = psid == p
        if m.any():
            corr[m] = np.interp(t[m], times, cs)
    return corr


def _strip_jitter_for_file(las_file, las, offsets):
    """Time jitter of a ground LAS (after the measured constant offsets), memoised."""
    try:
        p = Path(las_file)
        cache_key = (str(p), p.stat().st_mtime_ns)
    except OSError:
        cache_key = (str(las_file), 0)
    if cache_key in _STRIP_JITTER_CACHE:
        return _STRIP_JITTER_CACHE[cache_key]
    jitter = {}
    try:
        t = np.asarray(las.gps_time, dtype=np.float64)
        psid = np.asarray(las.point_source_id)
    except AttributeError:
        t, psid = None, None  # missing dimensions: no measurable jitter
    if (t is not None and psid is not None
            and len(t) == len(las.points) == len(psid)):
        z = np.asarray(las.z, dtype=np.float64)
        if offsets:
            lut = np.zeros(65536)
            for p_, off in offsets.items():
                lut[int(p_) & 0xFFFF] = off
            z = z - lut[np.asarray(las.point_source_id, dtype=np.int64)]
        jitter = _strip_jitter_offsets(
            np.asarray(las.x, dtype=np.float64),
            np.asarray(las.y, dtype=np.float64),
            z, psid, t)
    _STRIP_JITTER_CACHE[cache_key] = jitter
    return jitter


def _strip_offsets_for_file(las_file, las):
    """Alignment offsets of a ground LAS, memoised by (path, mtime)."""
    try:
        p = Path(las_file)
        cache_key = (str(p), p.stat().st_mtime_ns)
    except OSError:
        cache_key = (str(las_file), 0)
    if cache_key in _STRIP_OFFSETS_CACHE:
        return _STRIP_OFFSETS_CACHE[cache_key]
    try:
        psid = np.asarray(las.point_source_id)
    except AttributeError:
        psid = None  # missing dimension (third-party producer): no alignment
    if psid is not None and len(psid) == len(las.points):
        offsets = _strip_vertical_offsets(
            np.asarray(las.x, dtype=np.float64),
            np.asarray(las.y, dtype=np.float64),
            np.asarray(las.z, dtype=np.float64),
            psid)
    else:
        offsets = {}
    _STRIP_OFFSETS_CACHE[cache_key] = offsets
    return offsets


def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets,
                               jitter=None, lines=None):
    """Record the applied alignments (version, threshold, offsets, jitter, lines).

    The sidecar doubles as a cache tracker: a DTM without a sidecar, or one
    produced with a different version/threshold/jitter or line parameters,
    is regenerated. It is written even when no correction was applied, so a
    tile already known to be well aligned is not measured again.
    """
    jitter_payload = {}
    for p, (times, cs) in (jitter or {}).items():
        cs = np.asarray(cs, dtype=np.float64)
        jitter_payload[str(p)] = {
            "bins": int(len(times)),
            "rms_m": round(float(np.sqrt(np.mean(cs ** 2))), 4),
            "max_m": round(float(np.max(np.abs(cs))), 4),
            "series_m": [round(float(c), 4) for c in cs],
        }
    payload = {
        "version": STRIP_ALIGN_VERSION,
        "threshold": STRIP_ALIGN_THRESHOLD,
        "offsets": offsets,
        "jitter_bin": STRIP_JITTER_BIN,
        "jitter_smooth": STRIP_JITTER_SMOOTH,
        "jitter": jitter_payload,
        "line_window": STRIP_LINE_WINDOW,
        "line_cell": STRIP_LINE_CELL,
        "line_model": STRIP_LINE_MODEL,
        "lines": {str(p): {"lines": n, "rms_m": round(r, 4), "max_m": round(m, 4)}
                  for p, (n, r, m) in (lines or {}).items()},
    }
    try:
        sidecar = Path(dtm_dir) / f"{basename}_dtm{output_suffix}_stripalign.json"
        sidecar.write_text(json.dumps(payload, ensure_ascii=False),
                           encoding="utf-8")
    except Exception as e:
        logger.warning(f"  Could not write the strip alignment sidecar: {e}")


def _strip_lidar_ext(path):
    """Extract base name from a LAZ/LAS file (mirrors pipeline._file_basename)."""
    name = Path(path).name
    for ext in ('.copc.laz', '.copc.las', '.laz', '.las'):
        if name.lower().endswith(ext):
            return name[:-len(ext)]
    return Path(path).stem


def parse_ign_classes(spec):
    """Convert a list of IGN classes (names or codes) into sorted LAS codes.

    Args:
        spec: Comma-separated string, e.g. "sol,unclassified" or "2,1".

    Returns:
        Sorted list of unique LAS codes.

    Raises:
        ValueError: If an item is neither a known name nor a LAS code 0-255,
            or if the list is empty.
    """
    codes = set()
    for token in str(spec).split(","):
        token = token.strip().lower()
        if not token:
            continue
        if token in IGN_CLASS_NAMES:
            codes.add(IGN_CLASS_NAMES[token])
        else:
            try:
                code = int(token)
            except ValueError:
                raise ValueError(
                    f"Unknown IGN class: '{token}' "
                    f"(names: {', '.join(sorted(IGN_CLASS_NAMES))} or LAS code 0-255)")
            if not 0 <= code <= 255:
                raise ValueError(f"LAS code out of range (0-255): {code}")
            codes.add(code)
    if not codes:
        raise ValueError("No IGN class given")
    return sorted(codes)


def ign_method_label(codes):
    """Method label encoding the IGN classes (e.g. 'ign_1_2').

    'ign' alone = ground only (code 2), backward compatible with existing
    classification files. Any other combination is encoded in the name to
    invalidate the cache and trigger reclassification.
    """
    codes = sorted(codes)
    if codes == [2]:
        return "ign"
    return "ign_" + "_".join(str(c) for c in codes)


def _create_ground_pipeline(input_laz, output_las, method, ign_codes=None):
    """Create a PDAL pipeline JSON for ground classification.

    All methods include a ReturnNumber/NumberOfReturns >= 1 filter to handle
    LiDAR HD files that may contain points with invalid return numbers.

    Pre-processing steps (PDAL recommended workflow):
    1. Reset Classification to 0
    2. ELM (Extended Local Minimum) — mark low outliers as noise (Classification=7)
    3. Statistical outlier removal
    4. Ground classification (SMRF or CSF)
    5. Extract ground points (Classification=2)

    Args:
        input_laz: Path to input LAZ/LAS file.
        output_las: Path to output classified LAS file.
        method: Ground classification method ('ign', 'smrf' or 'csf').
        ign_codes: LAS class codes to extract with the 'ign' method
            (default: [2] = sol). Multiple ranges on Classification are
            logically ORed by filters.range (documented PDAL semantics).

    Returns:
        JSON string of the PDAL pipeline.
    """
    # Common ReturnNumber filter for LiDAR HD compatibility
    return_filter = {
        "type": "filters.range",
        "limits": "ReturnNumber[1:],NumberOfReturns[1:]"
    }

    # Classification filter (ground points only)
    ground_filter = {
        "type": "filters.range",
        "limits": "Classification[2:2]"
    }

    # IGN LiDAR HD: the file is pre-classified by the supplier. The
    # classification is reused as is (pure mode): the extracted classes are
    # configurable, ground only (2) by default, but e.g. unclassified (1) can
    # be added to fill holes without reprocessing. Multiple ranges on
    # Classification are ORed by filters.range (documented PDAL semantics).
    if method == 'ign':
        codes = sorted(ign_codes) if ign_codes else [2]
        ign_filter = {
            "type": "filters.range",
            "limits": ",".join(f"Classification[{c}:{c}]" for c in codes)
        }
        pipeline = {
            "pipeline": [
                str(input_laz),
                return_filter,
                ign_filter,
                {
                    "type": "writers.las",
                    "filename": str(output_las),
                    "extra_dims": "all"
                }
            ]
        }
        return json.dumps(pipeline)

    # Reset Classification to 0 before preprocessing
    reset_classification = {
        "type": "filters.assign",
        "assignment": "Classification[:]=0"
    }

    # ELM (Extended Local Minimum) — mark low outliers as noise
    # Parameters tuned for rocky limestone terrain with low vegetation:
    # - cell=5.0m: fine resolution to capture rocky relief
    # - threshold=2.0m: high threshold to avoid marking rock outcrops as noise
    elm_filter = {
        "type": "filters.elm",
        "cell": 5.0,
        "threshold": 2.0
    }

    # Statistical outlier removal
    outlier_filter = {
        "type": "filters.outlier",
        "method": "statistical",
        "mean_k": 8,
        "multiplier": 3.0
    }

    # Method-specific ground classification filter
    if method == 'smrf':
        ground_step = {
            "type": "filters.smrf",
            "ignore": "Classification[7:7]",
            "slope": 1.0,
            "window": 16.0,
            "threshold": 0.5,
            "scalar": 1.25
        }
    elif method == 'csf':
        # resolution 1.0 m: a 0.5 m cloth (= 4 M particles per km²) makes
        # classification ~4× slower with no visible gain on the DTM (the
        # final DTM resolution comes from rasterisation, not from the cloth).
        ground_step = {
            "type": "filters.csf",
            "resolution": 1.0,
            "rigidness": 3,
            "smooth": True,
            "threshold": 0.5
        }
    else:
        raise ValueError(f"Unknown classification method: {method}")

    pipeline = {
        "pipeline": [
            str(input_laz),
            return_filter,
            reset_classification,
            elm_filter,
            outlier_filter,
            ground_step,
            ground_filter,
            {
                "type": "writers.las",
                "filename": str(output_las),
                "extra_dims": "all"
            }
        ]
    }
    return json.dumps(pipeline)


def create_smrf_pipeline(input_laz, output_las):
    """Create a PDAL pipeline JSON for SMRF ground classification."""
    return _create_ground_pipeline(input_laz, output_las, 'smrf')


def create_ign_pipeline(input_laz, output_las):
    """Create a PDAL pipeline JSON using the IGN supplier pre-classification."""
    return _create_ground_pipeline(input_laz, output_las, 'ign')


def create_csf_pipeline(input_laz, output_las):
    """Create a PDAL pipeline JSON for CSF ground classification."""
    return _create_ground_pipeline(input_laz, output_las, 'csf')


def validate_laz(laz_file):
    """Integrity check for a LAZ/LAS file.

    Verifies that both the header AND point data are readable. Some corrupted
    COPC files have valid headers but inaccessible point data (LazrsError:
    failed to fill whole buffer). Such files must be re-downloaded.

    Returns:
        True if file is readable and contains accessible points, False otherwise.
    """
    import laspy
    try:
        with laspy.open(str(laz_file)) as f:
            header = f.header
            point_count = header.point_count
            if point_count == 0:
                logger.error(f"  ✗ Empty file (0 points): {laz_file.name}")
                logger.error(f"    → Download it again from https://ign.fr/lidar-hd")
                return False
            # Verify point data is actually accessible (not just header metadata)
            try:
                for _ in f.chunk_iterator(1000):
                    break  # Read just one chunk to confirm data integrity
            except Exception as chunk_err:
                logger.error(f"  ✗ Point data unreadable (corrupted file?): {laz_file.name}")
                logger.error(f"    Error: {chunk_err}")
                logger.error(f"    → Download it again from https://ign.fr/lidar-hd")
                return False
        return True
    except Exception:
        pass

    # Fallback: try PDAL (handles COPC v1.1 that laspy can't read)
    try:
        result = subprocess.run(
            ["pdal", "info", str(laz_file), "--summary"],
            capture_output=True, text=True, timeout=30
        )
        if result.returncode == 0:
            # Check point count from PDAL info output
            import json as _json
            try:
                info = _json.loads(result.stdout)
                count = info.get('summary', {}).get('num_points', 0)
                if count == 0:
                    logger.error(f"  ✗ Empty file (0 points per PDAL): {laz_file.name}")
                    logger.error(f"    → Download it again from https://ign.fr/lidar-hd")
                    return False
            except Exception:
                pass  # Can't parse — assume valid
            # Verify PDAL can actually read point data (not just header)
            try:
                test_result = subprocess.run(
                    ["pdal", "info", str(laz_file), "--point", "1"],
                    capture_output=True, text=True, timeout=60
                )
                if test_result.returncode != 0:
                    logger.error(f"  ✗ Point data unreadable (PDAL): {laz_file.name}")
                    logger.error(f"    → Download it again from https://ign.fr/lidar-hd")
                    return False
            except (subprocess.TimeoutExpired, FileNotFoundError):
                pass  # Timeout — assume valid, will fail later if corrupted
            return True
        logger.error(f"  ✗ Unreadable file: {laz_file.name}")
        logger.error(f"    PDAL: {result.stderr.strip()[:200]}")
    except (subprocess.TimeoutExpired, FileNotFoundError):
        logger.error(f"  ✗ Could not check the file: {laz_file.name}")
    logger.error(f"    → Download it again from https://ign.fr/lidar-hd")
    return False


def _read_with_pdal(laz_file):
    """Read a LAZ/LAS file via PDAL when laspy fails (e.g. COPC v1.1).

    Returns a laspy.LasData object, or None on failure.
    """
    import subprocess
    import tempfile
    import os

    tmp_path = None
    try:
        # Convert COPC to LAS via PDAL, then read with laspy
        with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp:
            tmp_path = tmp.name

        pipeline = json.dumps({
            "pipeline": [
                str(laz_file),
                {"type": "writers.las", "filename": tmp_path}
            ]
        })

        result = subprocess.run(
            ["pdal", "pipeline", "--stdin"],
            input=pipeline, capture_output=True, text=True, timeout=300
        )

        if result.returncode != 0:
            logger.warning(f"  PDAL conversion failed: {result.stderr[:200]}")
            return None

        import laspy
        las = laspy.read(tmp_path)
        if len(las.points) == 0:
            logger.warning(f"  PDAL: conversion succeeded but produced 0 points")
            return None
        return las

    except Exception as e:
        logger.warning(f"  PDAL fallback failed: {e}")
        return None
    finally:
        # Always delete the temporary LAS, even on an exception (subprocess
        # timeout, failed laspy.read...): otherwise ~GB are leaked.
        if tmp_path:
            try:
                os.unlink(tmp_path)
            except Exception:
                pass


def detect_ground_method(laz_file):
    """Detect the best ground classification method based on point cloud statistics.

    Auto-selects between IGN, SMRF and CSF:
    - IGN: supplier pre-classification, chosen when >= 20 % of the points
      are class 2 (ground)
    - SMRF: fast, robust for most natural terrain (PDAL recommended default)
    - CSF: cloth simulation, better for complex/urban terrain

    Falls back to SMRF if the file cannot be read or attributes are missing.

    Args:
        laz_file: Path to input LAZ/LAS file.

    Returns:
        String: 'ign', 'smrf' or 'csf'
    """
    import laspy

    # Try laspy first, then PDAL for COPC files
    las = None
    try:
        las = laspy.read(str(laz_file))
    except Exception as e:
        logger.warning(f"  laspy: {e}")
        logger.info(f"  → Reading via PDAL for auto-detection...")
        las = _read_with_pdal(laz_file)

    if las is None:
        logger.info(f"  → Method: SMRF (default, file unreadable)")
        return 'smrf'
    _LAST_READ.clear()
    _LAST_READ[str(laz_file)] = las       # reused by the IGN extraction

    total_points = len(las.points)
    if total_points == 0:
        logger.warning(f"  Empty point cloud (0 points), default method: SMRF")
        return 'smrf'

    # IGN LiDAR HD: the delivered data is pre-classified by the supplier
    # (class 2 = ground). It is the fastest base (~10 s) and the reference.
    # Holes left by the pre-classification (dense forest / relief) are then
    # filled by create_dtm_fast (bounded gap filling + interpolation), so it
    # is preferred as soon as a reasonable share of the points is classified
    # as ground, rather than filtering again.
    try:
        cls = np.asarray(las.classification, dtype=np.int32)
        ground_ratio = float(np.mean(cls == 2))
    except Exception:
        ground_ratio = 0.0
    if ground_ratio >= 0.2:
        logger.info(f"  → Method: IGN (supplier pre-classification, "
                    f"{ground_ratio * 100:.1f}% class 2 points)")
        return 'ign'

    z = np.array(las.z)

    # Height variance (always available)
    z_std = float(np.std(z))
    z_range = float(np.max(z) - np.min(z))

    # Try to get NumberOfReturns (may not exist in all point formats)
    single_return_ratio = 0.0
    try:
        num_returns = np.array(las.NumberOfReturns)
        single_return_count = int(np.sum(num_returns == 1))
        single_return_ratio = single_return_count / total_points if total_points > 0 else 0
    except AttributeError:
        logger.debug("  NumberOfReturns unavailable, using Z variance only")

    logger.info(f"  Point cloud analysis: {total_points:,} points, "
                f"single_return_ratio={single_return_ratio:.2f}, "
                f"z_std={z_std:.1f}m, z_range={z_range:.1f}m")

    # Decision logic:
    # - High single-return ratio (>0.6) → urban (buildings, roads) → CSF (cloth simulation)
    # - High elevation variance (>30m) → complex/mountainous terrain → CSF
    # - Default → SMRF (fast, robust for most natural terrain)
    if single_return_ratio > 0.6:
        method = 'csf'
        reason = f"single-return ratio={single_return_ratio:.2f} > 0.6 → urban area"
    elif z_std > 30:
        method = 'csf'
        reason = f"z_std={z_std:.1f}m > 30m → complex terrain"
    else:
        method = 'smrf'
        reason = f"standard natural terrain"

    logger.info(f"  → Method: {method.upper()} ({reason})")
    return method


_LAST_READ = {}


def _extract_ign_ground(laz_file, output_las, codes):
    """Direct extraction (laspy) of the selected IGN classes into a LAS.

    Same result as the PDAL pipeline (ReturnNumber ≥ 1, NumberOfReturns ≥ 1
    and class filters) but with no re-read or conversion: ~4.9 s instead of
    ~13.5 s per tile, and the read done by auto-detection is reused. Point
    format, scales and offsets of the source file are preserved.

    Returns:
        True if the file was written with at least one point.
    """
    import laspy
    las = _LAST_READ.pop(str(laz_file), None)
    if las is None:
        las = laspy.read(str(laz_file))
    keep = ((np.asarray(las.return_number) >= 1)
            & (np.asarray(las.number_of_returns) >= 1)
            & np.isin(np.asarray(las.classification), np.asarray(sorted(codes))))
    if not keep.any():
        return False
    header = laspy.LasHeader(point_format=las.header.point_format,
                             version=las.header.version)
    header.scales = las.header.scales
    header.offsets = las.header.offsets
    out = laspy.LasData(header)
    out.points = las.points[keep]
    out.write(str(output_las))
    return True


def classify_ground(laz_file, temp_dir, method='auto', force=False, ign_classes="sol"):
    """Classify ground points using PDAL ground classification filter.

    Args:
        laz_file: Path to input LAZ/LAS file.
        temp_dir: Directory for temporary files (pipeline.json, ground.las).
        method: Ground classification method ('auto', 'ign', 'smrf' or 'csf').
        force: If True, reclassify even if output file already exists.
        ign_classes: LAS classes extracted by the IGN method (comma-separated
            names or codes, e.g. "sol,unclassified"). Ignored by the other
            methods.

    Returns:
        Path to classified ground LAS file, or None on failure.
    """
    import laspy

    # Auto-detect method if requested
    if method == 'auto':
        method = detect_ground_method(laz_file)
        logger.info(f"  Ground classification: {method.upper()} (auto)")
    else:
        logger.info(f"  Ground classification: {method.upper()} (forced)")

    # The IGN classes are encoded in the file name (e.g. ign_1_2) so that a
    # change of classes invalidates the cache and triggers reclassification.
    ign_codes = parse_ign_classes(ign_classes) if method == 'ign' else None
    method_label = ign_method_label(ign_codes) if ign_codes else method

    laz_base = _strip_lidar_ext(laz_file)

    output_las = temp_dir / f"{laz_base}_ground_{method_label}.las"

    if output_las.exists() and not force:
        logger.info(f"  {method.upper()} classification already done, reusing the existing file")
        return output_las

    if force and output_las.exists():
        logger.info(f"  Forced reclassification, deleting {output_las.name}")
        output_las.unlink()

    # IGN pre-classification: direct extraction, PDAL as a fallback.
    if method == 'ign':
        try:
            if _extract_ign_ground(laz_file, output_las, ign_codes or [2]):
                logger.info(f"  ✓ IGN ground classification done (direct extraction)")
                return output_las
            logger.warning("  No point in the requested IGN classes, falling back to SMRF")
            output_las.unlink(missing_ok=True)
            return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)
        except Exception as e:
            logger.warning(f"  Direct IGN extraction failed ({e}), using the PDAL pipeline")
            output_las.unlink(missing_ok=True)
    _LAST_READ.clear()

    pipeline_json = _create_ground_pipeline(laz_file, output_las, method, ign_codes=ign_codes)
    pipeline_file = temp_dir / f"pipeline_{method_label}.json"

    with open(pipeline_file, 'w') as f:
        f.write(pipeline_json)

    try:
        subprocess.run(
            ["pdal", "pipeline", str(pipeline_file)],
            capture_output=True, check=True
        )
        # Verify that ground file has points
        if output_las.exists() and output_las.stat().st_size < 100:
            logger.error(f"  ✗ Empty ground file (size < 100 bytes)")
            output_las.unlink(missing_ok=True)
            # Fallback: if the method yields no ground point, retry with SMRF
            if method in ('csf', 'ign'):
                return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)
            return None
        logger.info(f"  ✓ {method.upper()} ground classification done")
        return output_las
    except subprocess.CalledProcessError as e:
        error_msg = e.stderr.decode() if e.stderr else str(e)
        logger.warning(f"  ✗ PDAL classification error ({method.upper()}): {error_msg}")

        # Fallback: if CSF or the pre-classification fails, retry with SMRF
        if method in ('csf', 'ign'):
            return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)

        # Try repairing file with laspy if PDAL fails on EVLR/VLR
        if 'VLR' in error_msg or 'Invalid' in error_msg:
            logger.info(f"  → Trying to repair the file with laspy...")
            repaired_las = temp_dir / f"{laz_base}_repaired.las"
            if _repair_laz_with_laspy(laz_file, repaired_las):
                # Retry PDAL pipeline with repaired file
                pipeline_json = _create_ground_pipeline(repaired_las, output_las, method)
                with open(pipeline_file, 'w') as f:
                    f.write(pipeline_json)
                try:
                    subprocess.run(
                        ["pdal", "pipeline", str(pipeline_file)],
                        capture_output=True, check=True
                    )
                    logger.info(f"  ✓ {method.upper()} ground classification done (repaired file)")
                    return output_las
                except subprocess.CalledProcessError as e2:
                    error_msg2 = e2.stderr.decode() if e2.stderr else str(e2)
                    logger.error(f"  ✗ Classification failed even after repair: {error_msg2}")
            else:
                logger.error(f"  ✗ Could not repair the file")
        return None


def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False, source='csf'):
    """Retry ground classification with SMRF when CSF/IGN fails.

    CSF (Cloth Simulation Filter) can fail on certain terrain types where
    SMRF (Simple Morphological Filter) succeeds, and a file without usable
    pre-classification produces an empty ground extract. This fallback ensures
    processing continues even when the selected method fails.

    Args:
        laz_file: Path to input LAZ/LAS file.
        temp_dir: Directory for temporary files.
        laz_base: Base name for the file.
        force: If True, reclassify even if output exists.
        source: Method that failed ('csf' or 'ign').

    Returns:
        Path to classified ground LAS file, or None on failure.
    """
    logger.info(f"  → Switching {source.upper()} → SMRF (fallback)")

    # Clean up failed output if it exists
    failed_output = temp_dir / f"{laz_base}_ground_{source}.las"
    if failed_output.exists():
        failed_output.unlink(missing_ok=True)

    output_las = temp_dir / f"{laz_base}_ground_smrf.las"

    if output_las.exists() and not force:
        logger.info(f"  SMRF classification already exists, reusing the file")
        return output_las

    pipeline_json = _create_ground_pipeline(laz_file, output_las, 'smrf')
    pipeline_file = temp_dir / "pipeline_smrf.json"

    with open(pipeline_file, 'w') as f:
        f.write(pipeline_json)

    try:
        subprocess.run(
            ["pdal", "pipeline", str(pipeline_file)],
            capture_output=True, check=True
        )
        if output_las.exists() and output_las.stat().st_size < 100:
            logger.error(f"  ✗ Empty SMRF ground file (size < 100 bytes)")
            output_las.unlink(missing_ok=True)
            return None
        logger.info(f"  ✓ SMRF ground classification done (fallback from {source.upper()})")
        return output_las
    except subprocess.CalledProcessError as e:
        error_msg = e.stderr.decode() if e.stderr else str(e)
        logger.error(f"  ✗ SMRF classification failed (fallback): {error_msg}")
        return None


def _repair_laz_with_laspy(input_laz, output_las):
    """Try to repair a corrupt LAZ file by re-reading with laspy and saving as LAS.

    Works around PDAL errors like 'Invalid Extended VLR size' by stripping
    problematic VLR/EVLR metadata during re-save.

    Args:
        input_laz: Path to corrupt LAZ/LAS file.
        output_las: Path for repaired LAS output.

    Returns:
        True if repair succeeded, False otherwise.
    """
    import laspy
    try:
        las = laspy.read(str(input_laz))
        las.write(str(output_las))
        logger.info(f"  ✓ File repaired via laspy ({len(las.points):,} points)")
        return True
    except Exception as e:
        logger.warning(f"  ✗ laspy repair failed: {e}")
        return False


# ------------------------------------------------------------
# Filling the gaps between points (small holes)
# ------------------------------------------------------------
# At 0.2 m, ~80 % of the pixels receive no point: the DTM is filled between
# the measured points. A "fixed-distance" fill (every pixel within 1 m of a
# point) grew each isolated point into a flat 2 m patch and extrapolated a
# 1 m band into every hole (copied plateaus, false slopes: pink patches and
# coloured fringes in the oriented relief). The fill is now bounded to the
# ENVELOPE of the points: a morphological closing whose radius follows the
# local point spacing (short in dense areas, long under sparse cover),
# followed by removal of isolated islands.
GAP_FILL_TAG = "LIDAR_GAP_FILL"   # GeoTIFF tag: gap-fill version
GAP_FILL_VERSION = 3               # 1 (absent) = fixed 1 m distance; 2 = envelope; 3 = + density file
GAP_RADIUS_K = 1.5                 # radius = K × local point spacing
GAP_RADII_M = (1.0, 1.5, 2.0, 3.0)  # radius steps (m); 1 m minimum: holes < 2 m filled
GAP_DENSITY_WINDOW_M = 5.0         # window for measuring the local density
GAP_MIN_ISLAND_M2 = 1.0            # smaller data islands are removed


def _morph_step(mask, square, erode):
    """One 3×3 step (square or cross) of binary dilation/erosion using numpy
    slices (~10× faster than scipy.ndimage). At the border, the missing
    neighbour counts as the pixel itself."""
    op = np.logical_and if erode else np.logical_or
    out = mask.copy()
    op(out[1:], mask[:-1], out=out[1:])
    op(out[:-1], mask[1:], out=out[:-1])
    src = out.copy() if square else mask  # square: separable (columns then rows)
    op(out[:, 1:], src[:, :-1], out=out[:, 1:])
    op(out[:, :-1], src[:, 1:], out=out[:, :-1])
    return out


def _morph_disk(mask, radius_px, erode=False, first_step=1):
    """Dilation/erosion by an approximate disc: an octagon, alternating square
    (odd) and cross (even) steps. `first_step` continues a chain."""
    for step in range(first_step, first_step + radius_px):
        mask = _morph_step(mask, step % 2 == 1, erode)
    return mask


def _gap_fill_envelope(valid, resolution):
    """Mask of the pixels to fill or keep: envelope of the measured points.

    Morphological closing (dilation then erosion) with a variable radius:
    gaps BETWEEN neighbouring points are filled, nothing is extended outwards
    (an isolated point stays one pixel). The radius is GAP_RADIUS_K × the
    local point spacing, measured over GAP_DENSITY_WINDOW_M within the
    covered area (the edge of a hole does not drag the density down), rounded
    to the nearest GAP_RADII_M step. Islands smaller than GAP_MIN_ISLAND_M2
    are removed (sparse returns on water, in a courtyard...).
    """
    from scipy import ndimage as nd
    radii_px = sorted({max(1, int(round(r / resolution))) for r in GAP_RADII_M})
    r_max = radii_px[-1]
    # Mask extended by mirroring: the grid border neither erodes the data nor
    # fills empty bands. A single chain of dilations serves every step (disc
    # of radius r = the first r steps).
    pad = 2 * r_max
    dilated = {}
    grown = np.pad(valid, pad, mode="symmetric")
    for step in range(1, r_max + 1):
        grown = _morph_disk(grown, 1, first_step=step)
        if step in radii_px:
            dilated[step] = grown
    del grown

    # Local density: share of measured pixels within the covered area (within
    # the largest radius of a point), converted into a mean point spacing.
    # Computed on a ~1 m grid (5 m window: more than enough).
    covered = dilated[r_max][pad:-pad, pad:-pad]
    h, w = valid.shape
    f = max(1, int(round(1.0 / resolution)))
    hc, wc = -(-h // f), -(-w // f)

    def _block_mean(a):
        full = np.zeros((hc * f, wc * f), dtype=np.float32)
        full[:h, :w] = a
        return full.reshape(hc, f, wc, f).mean(axis=(1, 3))

    win = max(3, int(round(GAP_DENSITY_WINDOW_M / (resolution * f))) | 1)
    frac = nd.uniform_filter(_block_mean(valid), win, mode="nearest")
    frac_cov = nd.uniform_filter(_block_mean(covered), win, mode="nearest")
    with np.errstate(divide="ignore", invalid="ignore"):
        # wanted radius in pixels = K × spacing / resolution
        wanted = GAP_RADIUS_K / np.sqrt(frac / np.maximum(frac_cov, 1e-6))
    # Wanted radius upsampled bilinearly (otherwise step changes form 1 m
    # staircases along hole edges), then rounded to the nearest step.
    wanted = np.nan_to_num(np.clip(wanted, 0, 2 * r_max), nan=2 * r_max)
    if f > 1:
        wanted = nd.zoom(wanted.astype(np.float32), f, order=1, grid_mode=True,
                         mode="nearest")[:h, :w]
    mids = [(a + b) / 2 for a, b in zip(radii_px, radii_px[1:])]
    level = np.digitize(wanted, mids).astype(np.uint8)
    del covered, frac, frac_cov, wanted

    # Erosion per step, restricted to the bounding box of its pixels (large
    # radii only concern sparse areas)
    envelope = valid.copy()
    for i, r in enumerate(radii_px):
        sel = (level == i) & ~valid
        rows = np.flatnonzero(sel.any(axis=1))
        if rows.size == 0:
            continue
        cols = np.flatnonzero(sel.any(axis=0))
        y0, y1 = rows[0], rows[-1] + 1
        x0, x1 = cols[0], cols[-1] + 1
        m = 2 * r  # margin: erosion at the crop edge stays outside the box
        crop = dilated[r][y0 + pad - m:y1 + pad + m, x0 + pad - m:x1 + pad + m]
        closed = _morph_disk(crop, r, erode=True)[m:-m, m:-m]
        envelope[y0:y1, x0:x1] |= sel[y0:y1, x0:x1] & closed
    del dilated, level

    # Islands too small: removed (8-connectivity, area in m²)
    labels, n = nd.label(envelope, structure=np.ones((3, 3), dtype=bool))
    if n:
        area = np.bincount(labels.ravel()) * resolution * resolution
        small = area < GAP_MIN_ISLAND_M2
        small[0] = False
        if small.any():
            envelope &= ~small[labels]
    return envelope


def _fill_small_gaps(dtm, resolution):
    """Fill the gaps between points within the data envelope.

    Returns:
        Tuple (dtm, filled_count, removed_count): filled pixels, and measured
        pixels removed along with isolated islands.
    """
    valid = ~np.isnan(dtm)
    if valid.all() or not valid.any():
        return dtm, 0, 0
    envelope = _gap_fill_envelope(valid, resolution)
    from rasterio.fill import fillnodata
    # Every envelope pixel lies within the largest radius of a point
    max_px = max(1, int(round(max(GAP_RADII_M) / resolution)))
    # fillnodata writes into the array it is given: pass a copy
    filled = fillnodata(dtm.copy(), mask=valid, max_search_distance=max_px)
    to_fill = envelope & ~valid & ~np.isnan(filled)
    removed = valid & ~envelope
    out = dtm.copy()
    out[to_fill] = filled[to_fill]
    out[removed] = np.nan
    return out, int(to_fill.sum()), int(removed.sum())


# Ground point density ("densite_sol" layer): count of the points kept for
# the DTM in DENSITY_CELL_M cells, averaged over DENSITY_SMOOTH ×
# DENSITY_SMOOTH cells (pts/m² over 9 m²: less small-integer noise). Written
# next to the DTM (*_dtm*_density.tif): the visualisations only receive the
# DTM, no longer the points.
DENSITY_CELL_M = 1.0
DENSITY_SMOOTH = 3


def density_path(dtm_path):
    """Density companion file of a DTM."""
    dtm_path = Path(dtm_path)
    return dtm_path.with_name(f"{dtm_path.stem}_density.tif")


def _write_density(xs, ys, bounds, dtm_path):
    """Write the ground point density (pts/m², float32 GeoTIFF at 1 m)."""
    from scipy import ndimage as nd
    min_x, min_y, max_x, max_y = bounds
    cell = DENSITY_CELL_M
    nx = max(1, int(np.ceil((max_x - min_x) / cell - 1e-6)))
    ny = max(1, int(np.ceil((max_y - min_y) / cell - 1e-6)))
    top = max_y
    counts, _, _ = np.histogram2d(top - ys, xs - min_x, bins=(ny, nx),
                                  range=[[0, ny * cell], [0, nx * cell]])
    density = nd.uniform_filter(counts, DENSITY_SMOOTH, mode="nearest") / (cell * cell)
    out = density_path(dtm_path)
    with rasterio.open(
        out, 'w', driver='GTiff', height=ny, width=nx, count=1,
        dtype='float32', crs='EPSG:2154',
        transform=from_bounds(min_x, top - ny * cell, min_x + nx * cell, top, nx, ny),
        compress='deflate',
    ) as dst:
        dst.write(density.astype('float32'), 1)
    return out


def read_dtm_gap_fill(dtm_path):
    """Gap-fill version recorded in a DTM (1 if absent/unreadable)."""
    try:
        with rasterio.open(dtm_path) as src:
            return int(src.tags().get(GAP_FILL_TAG, 1) or 1)
    except Exception:
        return 1


def _interpolate_holes(dtm, downsample=8):
    """Fill remaining NaN holes with a terrain-aware surface interpolation.

    In complex / rocky terrain the ground under-classification leaves interior
    holes far too large for the bounded gap fill, which otherwise become flat
    nearest-neighbor patches in the downstream layers. This helper triangulates
    the valid cells on a downsampled grid (linear, nearest as a fallback for
    cells outside the data hull) and bilinearly upsamples the result, keeping
    the operation fast even for large rasters.

    Args:
        dtm: 2-D float array (may contain NaN holes).
        downsample: Coarsening factor for the interpolation grid.

    Not called by create_dtm_fast (large holes are deliberately left as
    nodata there); kept as a standalone helper.

    Returns:
        Tuple (filled_array, filled_count).
    """
    holes = np.isnan(dtm)
    if not holes.any():
        return dtm, 0
    valid = ~holes
    if not valid.any():
        return dtm, 0

    height, width = dtm.shape
    step = max(1, downsample)
    coarse = dtm[::step, ::step].astype(np.float64)
    c_valid = ~np.isnan(coarse)
    c_holes = np.isnan(coarse)
    if not c_holes.any() or not c_valid.any():
        return dtm, 0

    from scipy.interpolate import griddata
    from scipy.ndimage import map_coordinates

    cy, cx = np.where(c_valid)
    c_coords = np.column_stack([cx, cy]).astype(np.float64)
    c_vals = coarse[c_valid]
    hy, hx = np.where(c_holes)
    h_coords = np.column_stack([hx, hy]).astype(np.float64)

    interp = griddata(c_coords, c_vals, h_coords, method='linear')
    bad = np.isnan(interp)
    if bad.any():
        interp[bad] = griddata(c_coords, c_vals, h_coords[bad], method='nearest')
    coarse_filled = coarse.copy()
    coarse_filled[c_holes] = interp

    # Coarse cell i represents fine column/row i*step, so fine index c maps to
    # coarse coordinate c/step (no half-cell offset).
    rows = np.arange(height) / step
    cols = np.arange(width) / step
    grid_y, grid_x = np.meshgrid(rows, cols, indexing='ij')
    upsampled = map_coordinates(coarse_filled, [grid_y, grid_x], order=1)

    filled = dtm.copy()
    filled[holes] = upsampled[holes]
    return filled, int(holes.sum())


# ------------------------------------------------------------
# Edge matching (edge buffer with the adjacent tiles)
# ------------------------------------------------------------
# Large-kernel visualisations (openness/SVF: radii up to 100 m, LRM: 15 m)
# truncate their window at the tile edge: renders then show a band of
# artefacts at every tile boundary. With edge_buffer > 0, the DTM is
# rasterised on the nominal 1 km tile EXTENDED by a band of `edge_buffer`
# metres filled with the ground points of the 8 neighbouring LAZ tiles
# (streamed PDAL read: crop + class filter). The visualisations thus see the
# real terrain beyond the edge, and the final images are cropped back to the
# exact tile (see rendering.py).
EDGE_BUFFER_TAG = "LIDAR_EDGE_BUFFER"  # GeoTIFF tag: buffer used (m)

# Subfolder of input/ where the neighbouring tiles downloaded only for edge
# matching are kept apart: they are not part of the corpus of tiles to render
# (global scans only list input/ itself, non-recursively).
EDGE_NEIGHBORS_DIRNAME = "edge_neighbors"

# Offsets of a tile's 8 neighbours (col, row) in km
_NEIGHBOR_OFFSETS = [(-1, -1), (0, -1), (1, -1), (-1, 0),

                     (1, 0), (-1, 1), (0, 1), (1, 1)]


def _tile_coords(name):
    """(col, row) coordinates in km of an LHD file name, or None."""
    from .index import parse_basename_coords
    return parse_basename_coords(Path(name).name)


def _neighbor_laz_files(source_laz):
    """List the (up to) 8 LAZ/LAS files adjacent to `source_laz`.

    Looks in the source's folder, then in its edge_neighbors/ subfolder: the
    neighbouring tiles downloaded only for edge matching are kept there;
    otherwise they would pile up directly in input/ and every global pass
    would add an extra ring of tiles to render.
    """
    coords = _tile_coords(source_laz)
    if coords is None:
        return []
    col, row = coords
    directory = Path(source_laz).parent
    search_dirs = [directory, directory / EDGE_NEIGHBORS_DIRNAME]
    neighbors = []
    for dcol, drow in _NEIGHBOR_OFFSETS:
        nc, nr = col + dcol, row + drow
        # LHD names: col/row in km zero-padded to 4 digits (e.g. 0638_6628);
        # the unpadded variant is also accepted for unusual tiles.
        names = {f"LHD_FXX_{nc:04d}_{nr:04d}", f"LHD_FXX_{nc}_{nr}"}
        matches = []
        for search_dir in search_dirs:
            for name in names:
                matches += list(search_dir.glob(f"{name}_*.las"))
                matches += list(search_dir.glob(f"{name}_*.laz"))
        if matches:
            neighbors.append(sorted(matches)[0])
        else:
            logger.debug(f"    Missing neighbour: LHD_FXX_{nc:04d}_{nr:04d} (empty edge band)")
    return neighbors


def _neighbor_ground_points(source_laz, bounds, classes):
    """Ground points of the neighbouring tiles within `bounds` (edge matching).

    Streamed PDAL read per neighbour: crop to the extended footprint, then
    class filter (same codes as the DTM). Best-effort: an unreadable or
    missing neighbour is skipped and the corresponding band stays empty.

    Args:
        source_laz: LAZ of the processed tile (used to locate the neighbours).
        bounds: (min_x, min_y, max_x, max_y) of the extended footprint.
        classes: LAS codes to extract (e.g. [2] = ground).

    Returns:
        Concatenated (xs, ys, zs) (empty arrays if there is no neighbour).
    """
    import tempfile

    neighbors = _neighbor_laz_files(source_laz)
    if not neighbors:
        return (np.empty(0),) * 3

    min_x, min_y, max_x, max_y = bounds
    codes = sorted(set(int(c) for c in classes)) or [2]
    limits = ",".join(f"Classification[{c}:{c}]" for c in codes)

    xs, ys, zs = [], [], []
    found = 0
    for neighbor in neighbors:
        tmp_path = None
        try:
            with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp:
                tmp_path = tmp.name
            pipeline = json.dumps({
                "pipeline": [
                    {"type": "readers.las", "filename": str(neighbor)},
                    {"type": "filters.crop",
                     "bounds": f"([{min_x},{max_x}],[{min_y},{max_y}])"},
                    {"type": "filters.range", "limits": limits},
                    {"type": "writers.las", "filename": tmp_path},
                ]
            })
            result = subprocess.run(
                ["pdal", "pipeline", "--stdin"],
                input=pipeline, capture_output=True, text=True, timeout=300
            )
            if result.returncode != 0:
                raise RuntimeError(result.stderr[:200])
            import laspy
            las = laspy.read(tmp_path)
            if len(las.points) > 0:
                xs.append(np.asarray(las.x, dtype=np.float64))
                ys.append(np.asarray(las.y, dtype=np.float64))
                zs.append(np.asarray(las.z, dtype=np.float64))
                found += 1
        except Exception as e:
            logger.debug(f"    Neighbour {Path(neighbor).name} skipped: {e}")
        finally:
            if tmp_path:
                try:
                    Path(tmp_path).unlink(missing_ok=True)
                except Exception:
                    pass

    if not found:
        logger.warning("    No readable neighbour, edge band left empty")
        return (np.empty(0),) * 3
    logger.info(f"    Edge matching: {found} neighbour(s), "
                f"{sum(len(a) for a in xs):,} ground pts")
    return np.concatenate(xs), np.concatenate(ys), np.concatenate(zs)


def read_dtm_edge_buffer(dtm_path):
    """Edge buffer recorded in a DTM (m; 0 if absent/unreadable)."""
    try:
        with rasterio.open(dtm_path) as src:
            tags = src.tags()
        return float(tags.get(EDGE_BUFFER_TAG, 0.0) or 0.0)
    except Exception:
        return 0.0


def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
                    output_suffix="", source_laz=None,
                    pure=False, strip_align=True, edge_buffer=0.0,
                    neighbor_classes=None):
    """Create DTM using fast binning method with gap filling.

    Args:
        las_file: Path to classified ground LAS file.
        basename: Base name for output file.
        dtm_dir: Directory for output DTM GeoTIFF.
        resolution: Grid resolution in meters per pixel.
        force: If True, regenerate even if DTM already exists.
        output_suffix: Suffix for output filename (e.g. '_r0p2' for additional resolutions).
        source_laz: Optional: full LAZ of the processed tile (used to read
            the 8 neighbouring tiles for edge matching, see edge_buffer).
        pure: No effect (kept for compatibility). Gaps between points are
            filled within the data envelope (_fill_small_gaps); large holes
            stay nodata (rendered as no-data in the outputs).
        strip_align: If True (default), measure and correct the vertical
            offsets between flight strips (PointSourceId) before
            rasterisation: constant offsets ≥ STRIP_ALIGN_THRESHOLD (0.5 cm),
            then GPS time-window jitter (only when scan_angle is missing),
            then the line-by-line joint adjustment (gps_time and scan_angle
            required); everything is recorded in a *_dtm*_stripalign.json
            sidecar.
        edge_buffer: Edge buffer in metres (0 = disabled). The DTM then covers
            the nominal 1 km tile extended by this band, filled with the
            ground points of the 8 neighbouring LAZ tiles (source_laz
            required); the visualisations compute on the extended footprint,
            then the images are cropped back to the exact tile
            (rendering.py). The buffer is written to the LIDAR_EDGE_BUFFER
            GeoTIFF tag for cache invalidation.
        neighbor_classes: LAS codes extracted from the neighbours (default
            [2]; the pipeline passes the DTM's IGN classes). Neighbours are
            read with their supplier pre-classification, even when the
            central tile is classified with SMRF/CSF (context band, a few cm
            of deviation at worst).

    Returns:
        Path to output DTM GeoTIFF, or None on failure.
    """
    output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"

    if output_tif.exists() and not force:
        logger.info(f"  DTM already exists, reusing the file: {output_tif.name}")
        return output_tif

    import laspy

    logger.info("  → Generating DTM...")

    try:
        t_read = time.perf_counter()
        las = laspy.read(str(las_file))
        logger.info(f"    Read {len(las.points):,} points "
                    f"({time.perf_counter() - t_read:.1f}s)")
    except Exception as e:
        # laspy can't read COPC v1.1 — try PDAL conversion
        logger.warning(f"  laspy: {e}")
        logger.info(f"  → Converting via PDAL to read COPC...")
        las = _read_with_pdal(las_file)
        if las is None:
            logger.error(f"  ✗ Could not read {las_file.name}")
            return None

    if len(las.points) == 0:
        logger.error(f"  ✗ Empty file (0 points): {las_file.name}")
        return None

    # Vertical alignment of the flight strips before rasterisation:
    # best-effort, if the measurement fails we carry on unaligned (never abort).
    strip_offsets = {}
    strip_jitter = {}
    gps_time = None
    t_align = time.perf_counter()
    if strip_align:
        try:
            strip_offsets = _strip_offsets_for_file(las_file, las)
        except Exception as e:
            logger.warning(f"  Strip alignment measurement failed ({e}), DTM left unaligned")
            strip_offsets = {}
        if strip_offsets:
            logger.info("    Strip alignment: " + ", ".join(
                f"PSID {p} {off:+.3f} m" for p, off in sorted(strip_offsets.items())))
        else:
            logger.debug("    Strip alignment: no offset >= "
                         f"{STRIP_ALIGN_THRESHOLD * 100:.1f} cm, nothing to correct")
        # 2nd pass: intra-strip jitter (GPS time required, skipped otherwise).
        try:
            gps_time = np.asarray(las.gps_time, dtype=np.float64)
        except AttributeError:
            gps_time = None
        # The line-by-line joint adjustment (3rd pass) already re-aligns every
        # line against the other strips, at every scale: time-window jitter
        # is only computed when scan_angle is missing.
        lines_possible = _scan_angle(las) is not None
        if gps_time is not None and len(gps_time) == len(las.points) and not lines_possible:
            try:
                strip_jitter = _strip_jitter_for_file(las_file, las, strip_offsets)
            except Exception as e:
                logger.warning(f"  Intra-strip jitter measurement failed ({e}), jitter left uncorrected")
                strip_jitter = {}
        if strip_jitter:
            for p_, (jt, jc) in sorted(strip_jitter.items()):
                logger.info(f"    Jitter PSID {p_}: ±{np.max(np.abs(jc)) * 100:.1f} cm "
                            f"(rms {np.sqrt(np.mean(np.asarray(jc) ** 2)) * 100:.1f} cm, "
                            f"{len(jt)} windows of {STRIP_JITTER_BIN:g} s)")
        else:
            logger.debug("    Intra-strip jitter: nothing to correct")
    logger.info(f"    Strip alignment: {time.perf_counter() - t_align:.1f}s")

    try:

        min_x, max_x = float(las.header.min[0]), float(las.header.max[0])
        min_y, max_y = float(las.header.min[1]), float(las.header.max[1])

        width = int(np.ceil((max_x - min_x) / resolution))
        height = int(np.ceil((max_y - min_y) / resolution))

        # Edge matching: footprint = nominal 1 km tile (aligned on the
        # multi-tile grid) + a band of edge_buffer metres filled with the
        # neighbours' ground points. Otherwise: header bounds (historical).
        used_edge_buffer = 0.0
        if edge_buffer > 0:
            coords = _tile_coords(source_laz or las_file)
            if coords is not None:
                buffer_px = max(1, int(round(edge_buffer / resolution)))
                buffer_m = buffer_px * resolution
                col_km, row_km = coords
                # LHD grid: (col, row) = north-west corner in km →
                # X ∈ [col, col+1] km, Y ∈ [row-1, row] km (north edge = row).
                min_x = float(col_km) * 1000.0
                max_x = min_x + 1000.0
                max_y = float(row_km) * 1000.0
                min_y = max_y - 1000.0
                ext_bounds = (min_x - buffer_m, min_y - buffer_m,
                              max_x + buffer_m, max_y + buffer_m)
                width = int(round(1000.0 / resolution)) + 2 * buffer_px
                height = width
                min_x, min_y, max_x, max_y = ext_bounds
                used_edge_buffer = float(edge_buffer)
            else:
                logger.warning("    Edge matching impossible: not an LHD "
                               f"file name ({basename}), tile alone")

        logger.debug(f"    Bounds: X[{min_x:.1f}, {max_x:.1f}] Y[{min_y:.1f}, {max_y:.1f}]")
        logger.debug(f"    Grid: {width}x{height} pixels ({len(las.points):,} points)")
        logger.info(f"    Rasterising {width}x{height} ({len(las.points):,} points)...")

        xs = np.asarray(las.x, dtype=np.float64)
        ys = np.asarray(las.y, dtype=np.float64)
        zs = np.asarray(las.z, dtype=np.float64)
        if strip_offsets:
            lut = np.zeros(65536)
            for p, off in strip_offsets.items():
                lut[int(p) & 0xFFFF] = off
            zs = zs - lut[np.asarray(las.point_source_id, dtype=np.int64)]
        if strip_jitter and gps_time is not None:
            zs = zs - _apply_strip_jitter(las.point_source_id, gps_time, strip_jitter)
        # 3rd pass: line-by-line offset (after the first two).
        strip_lines = {}
        if strip_align and gps_time is not None and len(gps_time) == len(zs):
            t_lines = time.perf_counter()
            try:
                line_corr, strip_lines = _scan_lines_for_file(las_file, las, zs, gps_time)
                zs = zs - line_corr
            except Exception as e:
                logger.warning(f"  Line-by-line offset measurement failed ({e}), left uncorrected")
                strip_lines = {}
            for p_, (nl_, rms_, max_) in sorted(strip_lines.items()):
                logger.info(f"    Lines PSID {p_}: rms {rms_ * 100:.1f} cm, max {max_ * 100:.1f} cm "
                            f"({nl_} lines)")
            logger.info(f"    Line-by-line offset: {time.perf_counter() - t_lines:.1f}s")

        # Ground points of the neighbouring tiles in the edge band
        # (best-effort, not strip-aligned: context band, the final image is
        # cropped to the tile before delivery).
        if used_edge_buffer > 0 and source_laz is not None:
            t_neigh = time.perf_counter()
            nx, ny, nz = _neighbor_ground_points(
                source_laz, (min_x, min_y, max_x, max_y),
                neighbor_classes if neighbor_classes is not None else [2])
            logger.info(f"    Neighbours (edge matching): {len(nx):,} points "
                        f"({time.perf_counter() - t_neigh:.1f}s)")
            if len(nx):
                xs = np.concatenate([xs, nx])
                ys = np.concatenate([ys, ny])
                zs = np.concatenate([zs, nz])

        # Density of the kept ground points (densite_sol layer), best-effort
        try:
            _write_density(xs, ys, (min_x, min_y, max_x, max_y),
                           dtm_dir / f"{basename}_dtm{output_suffix}.tif")
        except Exception as e:
            logger.warning(f"    Point density not written ({e})")

        t_raster = time.perf_counter()
        dtm = bin_mean_2d(xs, ys, zs, width, height,
                          (min_x, max_x), (min_y, max_y))
        if dtm is not None:
            # Output already (height, width): only the Y flip is left
            dtm = dtm[::-1, :]  # north at the top
            logger.info(f"    ✓ Rasterised {width}x{height} GPU "
                        f"({time.perf_counter() - t_raster:.1f}s)")
        else:
            stat = binned_statistic_2d(
                xs, ys, zs,
                statistic='mean',
                bins=[width, height],
                range=[[min_x, max_x], [min_y, max_y]]
            )
            dtm = stat.statistic.T
            dtm = dtm[::-1, :]  # Flip Y so north is at top
            logger.info(f"    ✓ Rasterised {width}x{height} CPU "
                        f"({time.perf_counter() - t_raster:.1f}s)")

        # Only small holes close to the data are filled; large holes stay
        # nodata in the renders. Deliberately NO lowest-return floor: under
        # dense canopy that return is vegetation, which would imprint the
        # trees into the DTM.

        # Gaps between points filled within the data envelope (radius adapted
        # to the local density), isolated islands removed
        nan_count = np.count_nonzero(np.isnan(dtm))
        if nan_count > 0:
            total = dtm.size
            nan_pct = 100.0 * nan_count / total
            logger.info(f"    {nan_count:,} pixels without data ({nan_pct:.1f}%)")
            t_fill = time.perf_counter()
            dtm, filled_count, removed_count = _fill_small_gaps(dtm, resolution)
            logger.info(f"    {filled_count:,} pixels filled between points, "
                        f"{removed_count:,} isolated-island pixels removed "
                        f"({time.perf_counter() - t_fill:.1f}s)")

        # Save as GeoTIFF
        output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"
        transform = from_bounds(min_x, min_y, max_x, max_y, width, height)

        t_write = time.perf_counter()
        with rasterio.open(
            output_tif, 'w',
            driver='GTiff', height=height, width=width,
            count=1, dtype='float32',
            crs='EPSG:2154', transform=transform,
            nodata=float('nan'),
            compress='lzw'
        ) as dst:
            dst.write(dtm.astype('float32'), 1)
            # Gap-fill version: a DTM from an earlier version is regenerated
            # (patches and fringes of the old gap fill).
            dst.update_tags(**{GAP_FILL_TAG: str(GAP_FILL_VERSION)})
            if used_edge_buffer > 0:
                # Edge buffer written into the file: changing --edge-buffer
                # ⇒ automatic invalidation of the DTM cache.
                dst.update_tags(**{EDGE_BUFFER_TAG: f"{used_edge_buffer:g}"})
        logger.info(f"    GeoTIFF write: {time.perf_counter() - t_write:.1f}s")

        if strip_align:
            _write_strip_align_sidecar(dtm_dir, basename, output_suffix,
                                       strip_offsets, strip_jitter, strip_lines)
        logger.info(f"    ✓ DTM created: {output_tif.name}")
        return output_tif

    except Exception as e:
        logger.error(f"  ✗ DTM error: {e}", exc_info=True)
        return None