# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
sourceval/geomatch.py — line-geometry conflation with a generalization dead-band.
================================================================================

Compares a CHARTED linestring against its SOURCE counterpart (quay walls,
shoreline construction, …) the way a cartographer judges it, not the way a
vertex diff would: chart geometry is *supposed* to differ from the source —
generalization is the cartographer's job — so agreement is measured on SHAPE,
never on vertex sets, and small differences live inside a deliberate DEAD-BAND
that stays silent.

Metrics (all densification-independent):
  * corridor coverage — the fraction of the charted line's length lying within
    ``CORRIDOR_X × tol`` of the source line;
  * out-of-corridor RUNS — contiguous charted stretches outside the corridor
    (a single displaced vertex is noise; a 40 m run is a finding);
  * Hausdorff distance — the worst deviation anywhere, symmetric, densified.

Verdict (the dead-band): a pair is flagged ONLY when the disagreement is
obvious — a qualifying run outside the corridor, or a worst deviation beyond
``HAUSDORFF_CAP_X × tol``. Within-tolerance and almost-similar shapes pass
silently. Thresholds are module constants, tunable like QC1's.

Pure Python (no GDAL/QGIS): geometry arrives as [(lon, lat), …] and is
projected to local metres with an equirectangular projection about the data's
own latitude — at quay scale its error against an ellipsoidal computation is
far below the tolerances used here (same argument as reconcile.distance_m).
"""

import math

# ── the dead-band (tunable, mirrors QC1's constant style) ────────────────────
PASS_FRACTION = 0.95      # ≥ this much charted length inside the corridor
CORRIDOR_X = 2.0          # corridor half-width = CORRIDOR_X × tol
HAUSDORFF_CAP_X = 3.0     # worst deviation beyond this × tol → obvious
MIN_RUN_M = 25.0          # a run outside the corridor at least this long …
MIN_RUN_FRAC = 0.10       # … or this share of the line, counts as obvious
MATCH_GATE_X = 10.0       # pairing gate: mean distance must beat this × tol

_EARTH_R = 6371008.8      # IUGG mean radius, metres (same as reconcile)


# ── local metric frame ───────────────────────────────────────────────────────
def project_local(coords, lon0, lat0):
    """[(lon, lat)] → [(x, y)] metres in an equirectangular frame at lat0."""
    k = math.cos(math.radians(lat0)) * _EARTH_R * math.pi / 180.0
    m = _EARTH_R * math.pi / 180.0
    return [((lon - lon0) * k, (lat - lat0) * m) for lon, lat in coords]


def unproject_local(x, y, lon0, lat0):
    """Inverse of project_local for reporting a worst-point back in lon/lat."""
    k = math.cos(math.radians(lat0)) * _EARTH_R * math.pi / 180.0
    m = _EARTH_R * math.pi / 180.0
    return (lon0 + (x / k if k else 0.0), lat0 + y / m)


# ── polyline primitives (metres space) ───────────────────────────────────────
def _seg_dist(px, py, ax, ay, bx, by):
    """Distance from point (px,py) to segment (a,b)."""
    dx, dy = bx - ax, by - ay
    seg2 = dx * dx + dy * dy
    if seg2 <= 0.0:
        ex, ey = px - ax, py - ay
    else:
        t = ((px - ax) * dx + (py - ay) * dy) / seg2
        if t < 0.0:
            t = 0.0
        elif t > 1.0:
            t = 1.0
        ex, ey = px - (ax + t * dx), py - (ay + t * dy)
    return math.hypot(ex, ey)


def dist_to_polyline(p, line):
    """Distance from point ``p`` to polyline ``line`` (≥2 vertices)."""
    px, py = p
    best = float("inf")
    for (ax, ay), (bx, by) in zip(line, line[1:]):
        d = _seg_dist(px, py, ax, ay, bx, by)
        if d < best:
            best = d
    return best


def polyline_length(line):
    return sum(math.hypot(bx - ax, by - ay)
               for (ax, ay), (bx, by) in zip(line, line[1:]))


def densify(line, step):
    """Points along ``line`` spaced ≤ ``step`` apart (vertices included)."""
    if len(line) < 2:
        return list(line)
    out = [line[0]]
    for (ax, ay), (bx, by) in zip(line, line[1:]):
        seg = math.hypot(bx - ax, by - ay)
        n = max(1, int(math.ceil(seg / step)))
        for i in range(1, n + 1):
            t = i / n
            out.append((ax + t * (bx - ax), ay + t * (by - ay)))
    return out


def mean_distance(a, b, step):
    """Mean distance from densified ``a`` to polyline ``b`` (pairing score)."""
    pts = densify(a, step)
    if not pts or len(b) < 2:
        return float("inf")
    return sum(dist_to_polyline(p, b) for p in pts) / len(pts)


def hausdorff(a, b, step):
    """Symmetric densified Hausdorff distance between two polylines."""
    if len(a) < 2 or len(b) < 2:
        return float("inf")
    d1 = max(dist_to_polyline(p, b) for p in densify(a, step))
    d2 = max(dist_to_polyline(p, a) for p in densify(b, step))
    return max(d1, d2)


def corridor_runs(chart, source, half_width, step):
    """Walk the densified charted line and measure agreement with the corridor
    of ``half_width`` metres about the source line.

    Returns ``(fraction_inside, runs)`` where each run is a contiguous
    out-of-corridor stretch: ``(length_m, worst_m, worst_xy)``."""
    pts = densify(chart, step)
    if len(pts) < 2 or len(source) < 2:
        return 0.0, []
    inside_len = total_len = 0.0
    runs = []
    cur_len, cur_worst, cur_worst_xy = 0.0, 0.0, None
    prev = pts[0]
    prev_out = dist_to_polyline(prev, source) > half_width
    for p in pts[1:]:
        seg = math.hypot(p[0] - prev[0], p[1] - prev[1])
        total_len += seg
        d = dist_to_polyline(p, source)
        out = d > half_width
        if out:
            cur_len += seg
            if d > cur_worst:
                cur_worst, cur_worst_xy = d, p
        else:
            inside_len += seg
            if prev_out and cur_len > 0.0:
                runs.append((cur_len, cur_worst, cur_worst_xy))
                cur_len, cur_worst, cur_worst_xy = 0.0, 0.0, None
        prev, prev_out = p, out
    if cur_len > 0.0:
        runs.append((cur_len, cur_worst, cur_worst_xy))
    frac = (inside_len / total_len) if total_len > 0 else 0.0
    return frac, runs


# ── the verdict ──────────────────────────────────────────────────────────────
def compare_lines(chart_ll, source_ll, tol_m):
    """Compare a charted line against its source counterpart (both
    [(lon, lat), …]) with tolerance ``tol_m`` metres.

    Returns a dict:
      flag        True only for an OBVIOUS difference (see module docstring)
      fraction    share of charted length within the corridor (2×tol)
      hausdorff_m worst deviation anywhere (metres)
      worst_run   (length_m, worst_m, (lon, lat)) of the worst run, or None
      length_m    charted line length
    """
    tol = max(0.01, float(tol_m))
    lon0 = chart_ll[0][0]
    lat0 = chart_ll[0][1]
    chart = project_local(chart_ll, lon0, lat0)
    source = project_local(source_ll, lon0, lat0)
    step = max(0.5, tol / 2.0)
    corridor = CORRIDOR_X * tol
    frac, runs = corridor_runs(chart, source, corridor, step)
    hd = hausdorff(chart, source, step)
    length = polyline_length(chart)

    worst_run = None
    if runs:
        r = max(runs, key=lambda t: t[0])
        wx, wy = r[2] if r[2] else (chart[0][0], chart[0][1])
        worst_run = (r[0], r[1], unproject_local(wx, wy, lon0, lat0))

    run_bar = max(MIN_RUN_M, MIN_RUN_FRAC * length)
    obvious_run = (frac < PASS_FRACTION and worst_run is not None
                   and worst_run[0] >= run_bar)
    obvious_worst = hd > HAUSDORFF_CAP_X * tol
    return {
        "flag": bool(obvious_run or obvious_worst),
        "fraction": frac,
        "hausdorff_m": hd,
        "worst_run": worst_run,
        "length_m": length,
    }


# ── pairing ──────────────────────────────────────────────────────────────────
def _bbox(coords, pad):
    xs = [c[0] for c in coords]
    ys = [c[1] for c in coords]
    return (min(xs) - pad, min(ys) - pad, max(xs) + pad, max(ys) + pad)


def _bbox_overlap(a, b):
    return not (a[2] < b[0] or b[2] < a[0] or a[3] < b[1] or b[3] < a[1])


def match_lines(chart_lines, source_lines, tol_m):
    """Pair charted lines with source lines: bbox prefilter, then best mean
    distance, gated at ``MATCH_GATE_X × tol`` so far-apart lines never pair.

    ``chart_lines`` / ``source_lines`` are lists of [(lon, lat), …].
    Returns ``(pairs, unmatched_chart, unmatched_source)`` — pairs as
    (chart_index, source_index), the rest as index lists."""
    tol = max(0.01, float(tol_m))
    gate = MATCH_GATE_X * tol
    if not chart_lines or not source_lines:
        return [], list(range(len(chart_lines))), list(range(len(source_lines)))
    lon0, lat0 = chart_lines[0][0]
    ch = [project_local(c, lon0, lat0) for c in chart_lines]
    so = [project_local(s, lon0, lat0) for s in source_lines]
    step = max(0.5, tol / 2.0)
    ch_bb = [_bbox(c, gate) for c in ch]
    so_bb = [_bbox(s, gate) for s in so]
    pairs = []
    used_source = set()
    for i, c in enumerate(ch):
        best, best_j = float("inf"), -1
        for j, s in enumerate(so):
            if j in used_source or not _bbox_overlap(ch_bb[i], so_bb[j]):
                continue
            d = mean_distance(c, s, step)
            if d < best:
                best, best_j = d, j
        if best_j >= 0 and best <= gate:
            pairs.append((i, best_j))
            used_source.add(best_j)
    matched_chart = {i for i, _ in pairs}
    unmatched_chart = [i for i in range(len(ch)) if i not in matched_chart]
    unmatched_source = [j for j in range(len(so)) if j not in used_source]
    return pairs, unmatched_chart, unmatched_source
