# SPDX-FileCopyrightText: 2025 XLeitstelle Planen und Bauen <xleitstelle@gv.hamburg.de>
# SPDX-FileContributor: Tobias Kraft <tobias.kraft@gv.hamburg.de>
# SPDX-FileContributor: Michael Holzapfel <michael.holzapfel@geocledian.com>
# SPDX-FileContributor: Anton Jacobsson <anton.jacobsson@init.com>
#
# SPDX-License-Identifier: EUPL-1.2
import logging
from collections import Counter
from dataclasses import dataclass
from math import atan2, pi
from typing import Any, List

from qgis.core import (
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsGeometry,
    QgsLayerTreeGroup,
    QgsLayerTreeLayer,
    QgsLineString,
    QgsPointXY,
    QgsProject,
    QgsVectorLayer,
    QgsWkbTypes,
)

from xmas_plugin.util.metadata import PLUGIN_DIR_NAME

logger = logging.getLogger(PLUGIN_DIR_NAME)


def transform_to_project_crs(
    geom: QgsGeometry, crs: QgsCoordinateReferenceSystem
) -> QgsGeometry:
    """
    Clones and transforms a QgsGeometry from given CRS to projects current CRS.
    """
    _geom = QgsGeometry()
    _geom.fromWkb(geom.asWkb())
    if (project_crs := QgsProject.instance().crs()) != crs:
        transformer = QgsCoordinateTransform(crs, project_crs, QgsProject.instance())
        _geom.transform(transformer)
    return _geom


def transform_to_crs(
    geom: QgsGeometry,
    source_crs: QgsCoordinateReferenceSystem,
    target_crs: QgsCoordinateReferenceSystem,
) -> QgsGeometry:
    """
    Clones and transforms a QgsGeometry from a source to a target CRS.
    """
    _geom = QgsGeometry()
    _geom.fromWkb(geom.asWkb())
    if source_crs != target_crs:
        transformer = QgsCoordinateTransform(
            source_crs, target_crs, QgsProject.instance()
        )
        _geom.transform(transformer)
    return _geom


def ensure_ccw_geom(g: QgsGeometry) -> QgsGeometry:
    """Return a copy with OGC-style polygon ring winding (shell CCW, holes CW)."""
    if not g or g.isEmpty() or g.isNull():
        return g
    # Orientation normalization is a GEOS-boundary invariant, not a semantic change
    if QgsWkbTypes.geometryType(g.wkbType()) != QgsWkbTypes.PolygonGeometry:
        return g
    gg = QgsGeometry(g).forcePolygonCounterClockwise()
    # Return copy if exists otherwise original
    return gg if gg and not gg.isEmpty() else g


def segmentize_if_no_curves(geometry: QgsGeometry) -> QgsGeometry:
    """
    Normalize curve-capable geometries to a linear representation.

    Even if no curved segments are present, this ensures a consistent
    non-curve encoding for persistence.
    """
    ref = geometry.constGet()
    if ref and not ref.hasCurvedSegments():
        abstract_geom = ref.segmentize()
        if abstract_geom:
            geometry = QgsGeometry(abstract_geom)
    return geometry


def collapse_singlepart_multi(geometry: QgsGeometry) -> QgsGeometry:
    """Convert single-part multi geometries to their single-type equivalent."""
    if not geometry or geometry.isEmpty() or geometry.isNull():
        return geometry

    part_count = sum(1 for _ in geometry.constParts())
    if part_count != 1:
        return geometry

    g = QgsGeometry(geometry)
    g.convertToSingleType()
    return g


def prepare_geometry_for_persistence(geometry: QgsGeometry) -> QgsGeometry:
    """Prepares geometry for persistence."""
    geometry = segmentize_if_no_curves(geometry)
    geometry = ensure_ccw_geom(geometry)
    geometry = collapse_singlepart_multi(geometry)
    return geometry


# ----------------------------
# makeValid evaluation helpers
# ----------------------------


def _signed_area_xy(ring: list[QgsPointXY]) -> float:
    """Signed area (shoelace). Positive => CCW, negative => CW (in standard XY)."""
    if not ring or len(ring) < 3:
        return 0.0
    s = 0.0
    # ring in QGIS often closed (first==last). Works either way.
    n = len(ring)
    for i in range(n - 1):
        x1, y1 = ring[i].x(), ring[i].y()
        x2, y2 = ring[i + 1].x(), ring[i + 1].y()
        s += x1 * y2 - x2 * y1
    return 0.5 * s


def _ring_orientation(ring: list[QgsPointXY]) -> str:
    a = _signed_area_xy(ring)
    if a > 0:
        return "CCW"
    if a < 0:
        return "CW"
    return "DEGENERATE"


def _extract_polygon_rings_xy(geom: QgsGeometry) -> list[list[list[QgsPointXY]]]:
    """
    Returns parts -> rings -> points (XY).
    Works for linear Polygon/MultiPolygon only. For curved polygon types, this may be empty.
    """
    # asPolygon/asMultiPolygon return lists of QgsPointXY
    if geom.isMultipart():
        mp = (
            geom.asMultiPolygon()
        )  # list[polygon]; polygon=list[rings]; ring=list[points]
        return mp or []
    else:
        p = geom.asPolygon()
        return [p] if p else []


def _fmt_geom_errors(errs: list[Any], limit: int = 3) -> list[str]:
    out = []
    for e in errs[:limit]:
        try:
            out.append(str(e.what()))
        except Exception:
            out.append(str(e))
    return out


def _safe_hausdorff(a: QgsGeometry, b: QgsGeometry) -> float | None:
    """
    Hausdorff distance is great for 'shape change' detection.
    Some QGIS builds expose hausdorffDistance; if not, return None.
    """
    try:
        return float(a.hausdorffDistance(b))
    except Exception:
        return None


@dataclass
class MakeValidStats:
    wkb: int
    wkb_name: str
    geom_type: int
    is_multipart: bool
    is_empty: bool
    is_null: bool
    area: float
    length: float
    bbox_w: float
    bbox_h: float
    n_parts: int
    rings_per_part: list[int]  # for polygons
    ring_orientations: list[list[str]]  # parts -> rings -> "CW/CCW/DEGENERATE"
    valid_errs: list[Any]  # raw QGIS errors


def _stats(geom: QgsGeometry) -> MakeValidStats:
    wkb = geom.wkbType() if geom else QgsWkbTypes.Unknown
    wkb_name = QgsWkbTypes.displayString(wkb)
    gtype = QgsWkbTypes.geometryType(wkb)

    bb = geom.boundingBox()
    bbox_w = float(bb.width())
    bbox_h = float(bb.height())

    area = float(geom.area()) if geom else 0.0
    length = float(geom.length()) if geom else 0.0

    valid_errs = geom.validateGeometry() if geom else []

    rings_per_part: list[int] = []
    ring_orientations: list[list[str]] = []
    n_parts = 1

    if gtype == QgsWkbTypes.PolygonGeometry:
        parts = _extract_polygon_rings_xy(geom)
        n_parts = len(parts)
        for poly in parts:
            rings_per_part.append(len(poly))
            ring_orientations.append([_ring_orientation(r) for r in poly])
    else:
        n_parts = (
            1 if not geom.isMultipart() else 2
        )  # not super meaningful for non-polygons

    return MakeValidStats(
        wkb=wkb,
        wkb_name=wkb_name,
        geom_type=gtype,
        is_multipart=geom.isMultipart(),
        is_empty=geom.isEmpty(),
        is_null=geom.isNull(),
        area=area,
        length=length,
        bbox_w=bbox_w,
        bbox_h=bbox_h,
        n_parts=n_parts,
        rings_per_part=rings_per_part,
        ring_orientations=ring_orientations,
        valid_errs=valid_errs,
    )


def evaluate_make_valid_on_layers(
    layers: list, tol: float = 0.002, log_limit: int = 10
) -> dict:
    """
    Evaluates makeValid() effects without modifying layers.
    tol is in layer units (meters in EPSG:25832).
    Returns a dict summary you can inspect in debugger/log.
    """
    tol2 = tol * tol

    changed_topo = 0
    changed_type = 0
    changed_orientation = 0
    became_valid = 0
    became_invalid = 0

    worst_by_symdiff = []  # list of (symdiff_area, layer_name, fid)
    worst_by_hausdorff = []  # list of (hausdorff, layer_name, fid)

    per_layer = Counter()
    examples = []

    for layer in layers:
        for f in layer.getFeatures():
            g0 = f.geometry()
            if not g0 or g0.isEmpty():
                continue

            s0 = _stats(g0)

            # Make a safe copy and run makeValid
            try:
                g1 = QgsGeometry(g0)  # copy
                g1 = g1.makeValid()
            except Exception as e:
                # makeValid itself failed
                examples.append((layer.name(), f.id(), "makeValid raised", repr(e)))
                continue

            if not g1 or g1.isEmpty():
                examples.append((layer.name(), f.id(), "makeValid produced empty", ""))
                continue

            s1 = _stats(g1)

            # type change?
            if s0.wkb != s1.wkb:
                changed_type += 1
                per_layer[(layer.name(), "type_changed")] += 1

            # topology change
            topo_equal = True
            try:
                topo_equal = bool(g0.equals(g1))
            except Exception:
                # if equals not available or errors, fall back to WKT string compare (weaker)
                topo_equal = g0.asWkt() == g1.asWkt()

            if not topo_equal:
                changed_topo += 1
                per_layer[(layer.name(), "topo_changed")] += 1

            # validity delta
            v0 = len(s0.valid_errs) == 0
            v1 = len(s1.valid_errs) == 0

            if (not v0) and v1:
                became_valid += 1
                per_layer[(layer.name(), "became_valid")] += 1
            if v0 and (not v1):
                became_invalid += 1
                per_layer[(layer.name(), "became_invalid")] += 1

            # ring orientation change
            if s0.ring_orientations and s1.ring_orientations:
                if s0.ring_orientations != s1.ring_orientations:
                    changed_orientation += 1
                    per_layer[(layer.name(), "orientation_changed")] += 1

            # symmetric difference area
            symdiff_area = None
            try:
                sym = g0.symDifference(g1)
                symdiff_area = float(sym.area()) if sym else 0.0
            except Exception:
                symdiff_area = None

            if symdiff_area is not None:
                worst_by_symdiff.append((symdiff_area, layer.name(), f.id()))

            # Hausdorff: "shape drift"
            hd = _safe_hausdorff(g0, g1)
            if hd is not None:
                worst_by_hausdorff.append((hd, layer.name(), f.id()))

            # Heuristic: anything with non-trivial change or validity flip
            if (not topo_equal) or (v0 != v1) or (s0.wkb != s1.wkb):
                examples.append(
                    (
                        layer.name(),
                        f.id(),
                        f"topo_equal={topo_equal} type={s0.wkb_name}->{s1.wkb_name} valid={v0}->{v1}",
                        f"before_errs={_fmt_geom_errors(s0.valid_errs)} after_errs={_fmt_geom_errors(s1.valid_errs)} "
                        f"rings0={s0.rings_per_part} rings1={s1.rings_per_part} "
                        f"ori0={s0.ring_orientations} ori1={s1.ring_orientations}",
                    )
                )
                per_layer[(layer.name(), "flagged")] += 1

    worst_by_symdiff.sort(reverse=True, key=lambda x: x[0])
    worst_by_hausdorff.sort(reverse=True, key=lambda x: x[0])

    # Log summary
    summary = {
        "changed_topo": changed_topo,
        "changed_type": changed_type,
        "changed_orientation": changed_orientation,
        "became_valid": became_valid,
        "became_invalid": became_invalid,
        "worst_symdiff": worst_by_symdiff[:log_limit],
        "worst_hausdorff": worst_by_hausdorff[:log_limit],
        "per_layer": dict(per_layer),
        "examples": examples[: max(50, log_limit * 5)],
        "tol": tol,
        "tol2": tol2,
    }
    return summary


def get_all_plan_layers(
    plan_qgis_layer_id: str, only_xplan_layers: bool = True
) -> list[QgsVectorLayer]:
    # search all layers that belong to plan_layers group and make sure they all have the same CRS

    plan_layer = QgsProject.instance().mapLayer(plan_qgis_layer_id)
    try:
        plan_layergroup: QgsLayerTreeGroup = (
            QgsProject.instance().layerTreeRoot().findLayer(plan_layer).parent()
        )
    except AttributeError:
        return []
    if not isinstance(plan_layergroup, QgsLayerTreeGroup):
        raise RuntimeError(f"plan-layer {plan_layer} is not part of a group")
    xplan_id = plan_layer.customProperties().value(f"{PLUGIN_DIR_NAME}/plan_id", None)
    if not xplan_id:
        raise RuntimeError(f"xplan-id not found for layer {plan_layer}")

    plan_layers: List[QgsVectorLayer] = []
    crs_ids: list[str] = []
    seen: set[str] = set()

    for item in plan_layergroup.findLayers():
        if not isinstance(item, QgsLayerTreeLayer):
            continue
        layer = item.layer()
        if not isinstance(layer, QgsVectorLayer):
            continue

        if only_xplan_layers:
            layer_plan_id = layer.customProperty("xmas_plugin/plan_id", None)
            if not layer_plan_id or layer_plan_id != xplan_id:
                continue

        if not layer.crs().isValid():
            continue

        if layer.id() in seen:
            continue
        seen.add(layer.id())

        plan_layers.append(layer)
        crs_ids.append(layer.crs().authid())

    if not plan_layers:
        raise RuntimeError("no layers found for plan_id")

    if len(set(crs_ids)) != 1:
        raise RuntimeError("more than 1 CRS for plan layers")

    return plan_layers


def snap_to_boundary(
    geom: QgsGeometry, boundary: QgsGeometry, tol: float
) -> QgsGeometry:
    if not geom or geom.isEmpty():
        return geom
    if not boundary or boundary.isEmpty():
        return QgsGeometry(geom)

    # Dynamic import (as requested)
    try:
        from qgis.analysis import QgsGeometrySnapper
    except Exception as e:
        logger.debug("qgis.analysis import failed: %r", e)
        return QgsGeometry(geom)

    # IMPORTANT: snap lines to boundary linework, not polygon interior
    ref = boundary
    try:
        if QgsWkbTypes.geometryType(geom.wkbType()) == QgsWkbTypes.LineGeometry:
            ref = boundary.boundary()  # polygon -> multiline boundary
    except Exception:
        pass

    if not ref or ref.isEmpty():
        return QgsGeometry(geom)

    try:
        mode = QgsGeometrySnapper.PreferNodes
        if hasattr(QgsGeometrySnapper, "EndPointToVertex"):
            mode = QgsGeometrySnapper.EndPointToVertex  # helps “last vertex drift”
        return QgsGeometrySnapper.snapGeometry(geom, tol, [ref], mode)
    except Exception as e:
        logger.debug("Snapping didn't work: %r", e)
        return QgsGeometry(geom)


def get_layer_segments(layer):
    """
    collect all segments from a Line or Polygon layer.

    Returns a list of dicts:
      {
        "layer_id": <string>,
        "fid": <int>,
        "start_vertex_id": <int>,
        "end_vertex_id": <int>,
        "start_point": QgsPoinY,
        "end_point": QgsPoint,
        "segment_geometry": QgsLineString,
        "direction_angle": float
      }
    """

    if not isinstance(layer, QgsVectorLayer):
        print(TypeError("layer must be a QgsVectorLayer"))
        return {}

    geom_type = layer.geometryType()
    if geom_type not in (QgsWkbTypes.LineGeometry, QgsWkbTypes.PolygonGeometry):
        logger.debug(
            f"Layer must be a Line or Polygon layer {layer.id()}/{layer.name()}"
        )
        return {}

    result = []
    layer_id = layer.id()
    layer_name = layer.name()
    layer_type = layer.customProperty(f"{PLUGIN_DIR_NAME}/layer_type", None)

    for feat in layer.getFeatures():
        fid = feat.id()
        geom = feat.geometry()
        xplan_fid = feat.attributeMap().get("id", None)
        if geom is None or geom.isEmpty():
            continue

        vertices = list(geom.vertices())
        for v in range(len(vertices) - 1):
            p1 = vertices[v]
            p2 = vertices[v + 1]

            p1_vtxid = geom.vertexIdFromVertexNr(v)[1]
            p2_vtxid = geom.vertexIdFromVertexNr(v + 1)[1]

            if p1_vtxid.part != p2_vtxid.part or p1_vtxid.ring != p2_vtxid.ring:
                continue

            seg_geom = QgsLineString([p1, p2])
            # segment direction (0 to 2π clockwise from north)
            dx = p2.x() - p1.x()
            dy = p2.y() - p1.y()
            angle = (2 * pi - (atan2(dx, dy) + pi)) % (2 * pi)

            result.append(
                {
                    "layer_id": layer_id,
                    "layer_name": layer_name,
                    "layer_type": layer_type,
                    "fid": fid,
                    "xplan_id": xplan_fid,
                    "start_vertex_id": v,
                    "end_vertex_id": v + 1,
                    "start_point": p1,
                    "end_point": p2,
                    "segment_geometry": seg_geom,
                    "direction_angle": angle,
                }
            )
    return result
