# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
validation/checks_geometry.py
-----------------------------
Spatial logical-consistency checks. These use the GeomOps operators and need
real OGR geometries, so they run inside QGIS (or any GDAL-enabled process); the
headless dict-based tests skip them because the test doubles carry no geometry.

The pattern for a relational check is: build a spatial index of the "container"
class once, then test each "subject" feature against only the candidates its
bounding box overlaps. We keep it simple and readable here (nested scan with an
early bbox reject); swap in QgsSpatialIndex/STRtree in the loader if a cell is
large enough to need it.
"""

from .engine import check, Severity, Finding


def _geoms(ds, klass):
    return [f for f in ds.features(klass) if f.geom is not None]


def _bbox_overlap(a, b):
    """Cheap envelope reject before the expensive predicate. OGR envelope is
    (minX, maxX, minY, maxY)."""
    ax0, ax1, ay0, ay1 = a.GetEnvelope()
    bx0, bx1, by0, by1 = b.GetEnvelope()
    return not (ax1 < bx0 or bx1 < ax0 or ay1 < by0 or by1 < ay0)


def _envelope_overlap_pairs(envs):
    """Yield index pairs (i, j) whose envelopes overlap, via an x-axis sweep so
    spatially separated geometries are never compared and each envelope is read
    once. ``envs`` is a list of (minX, maxX, minY, maxY). Replaces the O(n^2)
    all-pairs bbox scan (which also recomputed every envelope on each compare)
    with O(n log n + P), P = number of overlapping pairs. The set of pairs is
    identical to the brute-force bbox scan, so callers find exactly the same
    crossings."""
    import bisect
    n = len(envs)
    order = sorted(range(n), key=lambda k: envs[k][0])   # by minX
    active = []          # active indices, kept sorted by their maxX
    active_max = []      # parallel maxX list for bisect eviction
    for idx in order:
        x0, x1, y0, y1 = envs[idx]
        cut = bisect.bisect_left(active_max, x0)          # evict maxX < x0
        if cut:
            del active[:cut]
            del active_max[:cut]
        for a in active:                                  # x already overlaps
            _ax0, _ax1, ay0, ay1 = envs[a]
            if not (y1 < ay0 or ay1 < y0):                # y overlaps too
                yield a, idx
        pos = bisect.bisect_left(active_max, x1)
        active.insert(pos, idx)
        active_max.insert(pos, x1)


def _segment_bboxes(geom):
    """Per-segment bounding boxes (x0, y0, x1, y1) of an OGR (multi)linestring.

    The DEPCNT-cross broad phase works at *segment* granularity because a depth
    contour's whole envelope is a poor discriminator on dense bathymetry —
    nested/parallel contours' envelopes all overlap — whereas their individual
    segments do not. OGR-touching, so it is exercised by the in-QGIS path, not
    the headless pair-generator tests.
    """
    out = []
    try:
        gc = geom.GetGeometryCount()
    except Exception:
        gc = 0
    parts = [geom.GetGeometryRef(k) for k in range(gc)] if gc else [geom]
    for part in parts:
        try:
            pts = part.GetPoints()
        except Exception:
            pts = None
        if not pts:
            continue
        for a, b in zip(pts, pts[1:]):
            x0, x1 = (a[0], b[0]) if a[0] <= b[0] else (b[0], a[0])
            y0, y1 = (a[1], b[1]) if a[1] <= b[1] else (b[1], a[1])
            out.append((x0, y0, x1, y1))
    return out


def _segment_overlap_pairs(seg_bboxes_per_item, cell=None):
    """Yield index pairs (i, j), i < j, of items that have at least one pair of
    segments with overlapping bounding boxes — the cheap broad phase before the
    exact, expensive GEOS predicate.

    ``seg_bboxes_per_item[i]`` is a list of ``(x0, y0, x1, y1)`` segment boxes
    for item i. Every segment box is bucketed into the uniform-grid cells it
    touches (cell auto-sized to the median segment extent); two segments are
    compared only when they share a cell and reported when their boxes overlap.
    Pure Python and grid-accelerated.

    Correctness: this is a SUPERSET of every truly intersecting pair — two
    polylines that cross share the crossing point on a segment each, and both
    those boxes cover that point's cell, so the pair is always emitted. It is in
    fact exactly the set of "items with >=1 overlapping segment-box pair", which
    is why it can replace the whole-envelope sweep without changing which
    crossings the exact predicate then confirms. Pairs are yielded sorted, for a
    deterministic report order.
    """
    segs = []
    extents = []
    for owner, boxes in enumerate(seg_bboxes_per_item):
        for (x0, y0, x1, y1) in boxes:
            segs.append((owner, x0, y0, x1, y1))
            extents.append((x1 - x0) + (y1 - y0))
    if not segs:
        return
    if cell is None:
        extents.sort()
        cell = extents[len(extents) // 2] or extents[-1] or 1.0
    if cell <= 0:
        cell = 1.0
    inv = 1.0 / cell
    grid = {}
    pairs = set()
    for owner, x0, y0, x1, y1 in segs:
        cx0, cy0 = int(x0 * inv), int(y0 * inv)
        cx1, cy1 = int(x1 * inv), int(y1 * inv)
        for cx in range(cx0, cx1 + 1):
            for cy in range(cy0, cy1 + 1):
                bucket = grid.setdefault((cx, cy), [])
                for (o2, ox0, oy0, ox1, oy1) in bucket:
                    if o2 == owner:
                        continue
                    if x1 < ox0 or ox1 < x0 or y1 < oy0 or oy1 < y0:
                        continue                       # segment boxes disjoint
                    pairs.add((owner, o2) if owner < o2 else (o2, owner))
                bucket.append((owner, x0, y0, x1, y1))
    for pair in sorted(pairs):
        yield pair


@check(56, Severity.ERROR, "Logical consistency",
       "BUAARE is not located on LNDARE.",
       "Amend BUAARE so that it is covered by a LNDARE.", "Logical consistency")
def buaare_on_land(ds, cat, geom):
    land = _geoms(ds, "LNDARE")
    for f in ds.features("BUAARE"):
        if f.geom is None:
            continue
        if not any(_bbox_overlap(f.geom, la.geom) and geom.covered_by(f.geom, la.geom)
                   for la in land):
            yield Finding(_meta(56), f.cell, f.klass, f.fid)


# S-58 #54a — the classes a FORSTC/LNDMRK/SILTNK must sit on (S-58 Ed: the
# "suitable supporting object" set).
_54_COVER = ("BRIDGE", "COALNE", "DAMCON", "FLODOC", "LNDARE", "OFSPLF",
             "PILPNT", "PONTON", "PYLONS", "SLCONS", "UWTROC")


@check("54a", Severity.CRITICAL, "Logical consistency",
       "FORSTC, LNDMRK or SILTNK not covered by a suitable supporting object.",
       "Amend object to ensure it is situated on a suitable object.",
       "Logical consistency")
def forstc_lndmrk_siltnk_covered(ds, cat, geom):
    """S-58 #54a: each FORSTC, LNDMRK or SILTNK feature must be COVERED_BY one of
    the supporting classes in _54_COVER (bridge, coastline, land area, pontoon,
    pylon, …). Mirrors #56's COVERED_BY-some-single-feature semantics."""
    covers = [c for cls in _54_COVER for c in _geoms(ds, cls)]
    for kls in ("FORSTC", "LNDMRK", "SILTNK"):
        for f in ds.features(kls):
            if f.geom is None:
                continue
            if not any(_bbox_overlap(f.geom, c.geom)
                       and geom.covered_by(f.geom, c.geom) for c in covers):
                yield Finding(_meta("54a"), f.cell, f.klass, f.fid)


@check(60, Severity.WARNING, "Logical consistency",
       "CBLSUB is within or crosses a LNDARE.",
       "Ensure CBLSUB is not covered by a LNDARE.", "Logical consistency")
def cblsub_on_land(ds, cat, geom):
    land = [la for la in _geoms(ds, "LNDARE") if la.primitive == "A"]
    for f in ds.features("CBLSUB"):
        if f.geom is None:
            continue
        for la in land:
            if _bbox_overlap(f.geom, la.geom) and (
                    geom.within(f.geom, la.geom) or geom.crosses(f.geom, la.geom)):
                yield Finding(_meta(60), f.cell, f.klass, f.fid)
                break


@check("1500a", Severity.WARNING, "Object catalogue",
       "CBLARE object overlaps a LNDARE object.",
       "Amend objects to remove overlap.", "Logical consistency")
def cblare_overlaps_land(ds, cat, geom):
    land = [la for la in _geoms(ds, "LNDARE") if la.primitive == "A"]
    for f in ds.features("CBLARE"):
        if f.geom is None:
            continue
        for la in land:
            if _bbox_overlap(f.geom, la.geom) and (
                    geom.within(f.geom, la.geom) or geom.overlaps(f.geom, la.geom)):
                yield Finding(_meta("1500a"), f.cell, f.klass, f.fid)
                break


@check(77, Severity.CRITICAL, "Logical consistency",
       "DEPCNT objects cross.",
       "Amend DEPCNT objects so they do not cross.", "Logical consistency")
def depcnt_cross(ds, cat, geom):
    contours = [c for c in _geoms(ds, "DEPCNT") if c.geom is not None]
    # Broad phase at SEGMENT granularity. A depth contour's whole envelope is a
    # poor filter on dense bathymetry — contours every 0.1 m are nested, so their
    # envelopes all overlap and the old whole-envelope sweep let ~every pair
    # through to the costly GEOS crosses(). Segment-box co-location instead
    # rejects nested/offset-parallel contours cheaply (two parallel segments
    # offset in y have disjoint boxes) while never dropping a real crossing — the
    # crossing point lies on a segment of each, whose boxes share its cell — so
    # GEOS still confirms each crossing and the findings are identical.
    seg_boxes = [_segment_bboxes(c.geom) for c in contours]
    for i, j in _segment_overlap_pairs(seg_boxes):
        if geom.crosses(contours[i].geom, contours[j].geom):
            lo, hi = (i, j) if i < j else (j, i)          # stable report order
            yield Finding(_meta(77), contours[lo].cell, "DEPCNT",
                          contours[lo].fid,
                          detail=f"crosses fid {contours[hi].fid}")


@check(1786, Severity.WARNING, "Logical consistency",
       "Area object with WATLEV = 2 (always dry) is not covered by a LNDARE.",
       "Amend WATLEV value or ensure object is on land.", "Logical consistency")
def watlev2_on_land(ds, cat, geom):
    land = [la for la in _geoms(ds, "LNDARE") if la.primitive == "A"]
    for f in ds.features():
        if f.primitive != "A" or f.geom is None:
            continue
        if str(f.attrs.get("WATLEV")) != "2":
            continue
        if not any(_bbox_overlap(f.geom, la.geom) and geom.covered_by(f.geom, la.geom)
                   for la in land):
            yield Finding(_meta(1786), f.cell, f.klass, f.fid, detail="WATLEV=2")


@check(1573, Severity.ERROR, "Object catalogue",
       "DRYDOC is not covered by a LNDARE.",
       "Amend LNDARE or DRYDOC as required.", "4.6.6.1")
def drydoc_on_land(ds, cat, geom):
    land = [la for la in _geoms(ds, "LNDARE") if la.primitive == "A"]
    for f in ds.features("DRYDOC"):
        if f.geom is None:
            continue
        if not any(_bbox_overlap(f.geom, la.geom) and geom.covered_by(f.geom, la.geom)
                   for la in land):
            yield Finding(_meta(1573), f.cell, f.klass, f.fid)


# ── ENCoreGIS quality check: redundant (collinear) vertex bloat ──────────────
# ENC production increasingly starts from GIS sources whose geometries are
# densified: a straight quay that needs 2 shape points arrives as 100 collinear
# vertices, a rectangle needs 4 but carries hundreds. The product still passes
# conformance, but it bloats cells, slows every consumer, and buries the real
# shape points. QC1 measures, per line/area feature, how many vertices a
# Douglas-Peucker pass at a tolerance far below chart accuracy would discard —
# vertices that add NO shape — and reports the feature when the bloat is
# substantial. This is an ENCoreGIS QUALITY check, not an S-58 conformance
# requirement: findings are advice to simplify at source, never pass/fail.
_QC1_EPS_DEG = 1e-6          # ≈ 0.1 m N-S — far below chart accuracy anywhere
_QC1_MIN_REMOVABLE = 20      # ignore trifles: at least this many dead vertices
_QC1_MIN_RATIO = 0.5         # and more than half the feature adds no shape


def _dp_keep_count(pts, eps):
    """Vertices Douglas-Peucker KEEPS at tolerance ``eps`` (same units as the
    coordinates). Iterative (no recursion limit); endpoints always kept.
    Duplicate consecutive vertices count as removable automatically."""
    n = len(pts)
    if n <= 2:
        return n
    eps2 = eps * eps
    keep = 2                              # the two endpoints
    stack = [(0, n - 1)]
    while stack:
        i0, i1 = stack.pop()
        if i1 - i0 < 2:
            continue
        ax, ay = pts[i0]
        bx, by = pts[i1]
        dx, dy = bx - ax, by - ay
        seg2 = dx * dx + dy * dy
        dmax = -1.0
        imax = -1
        for i in range(i0 + 1, i1):
            px, py = pts[i]
            if seg2 <= 0.0:               # zero-length anchor: plain distance
                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)
            d2 = ex * ex + ey * ey
            if d2 > dmax:
                dmax = d2
                imax = i
        if dmax > eps2:
            keep += 1
            stack.append((i0, imax))
            stack.append((imax, i1))
    return keep


def _coord_runs(g):
    """Yield the vertex runs ([(x, y), …] per linestring/ring) of an OGR
    geometry, recursing through polygons and multi-geometries."""
    name = g.GetGeometryName()
    if name in ("LINESTRING", "LINEARRING"):
        pts = g.GetPoints()
        if pts:
            yield [(p[0], p[1]) for p in pts]
    else:
        for i in range(g.GetGeometryCount()):
            yield from _coord_runs(g.GetGeometryRef(i))


@check("QC1", Severity.WARNING, "Geometry quality",
       "Geometry is bloated with redundant (collinear) vertices.",
       "Simplify the geometry at source — an ENC only needs the "
       "shape-defining vertices (a straight line needs 2, a rectangle 4).",
       "ENCoreGIS quality check — not an S-58 conformance requirement",
       name="geometry-quality/redundant-collinear-vertices",
       validated=False)
def redundant_vertices(ds, cat, geom):
    for f in ds.features():
        if f.geom is None or f.primitive not in ("L", "A"):
            continue
        total = kept = 0
        try:
            for run in _coord_runs(f.geom):
                total += len(run)
                kept += (len(run) if len(run) < 3
                         else _dp_keep_count(run, _QC1_EPS_DEG))
        except Exception:                               # noqa: BLE001
            continue
        removable = total - kept
        if (removable >= _QC1_MIN_REMOVABLE and total
                and removable / total >= _QC1_MIN_RATIO):
            pct = int(round(100.0 * removable / total))
            yield Finding(
                _meta("QC1"), f.cell, f.klass, f.fid,
                detail=(f"{removable} of {total} vertices ({pct}%) add no "
                        f"shape within {_QC1_EPS_DEG}° (≈ 0.1 m) — "
                        "GIS-source densification"))


def _meta(number):
    from .engine import REGISTRY
    for m in REGISTRY:
        if m.number == number:
            return m
    from .engine import CheckMeta
    return CheckMeta(number, Severity.ERROR, "", "", "", "", None)
