
from functools import singledispatch
import numbers
from typing import Optional, Union, List

import numpy as np

from geogst.core.geometries.shapes.points import Point
from geogst.core.geometries.shapes.segments import Segment, shortest_segment_or_point3d, point_or_segment3d
from geogst.core.geometries.shapes.lines import Ln
from geogst.core.geometries.shapes.multilines import MultiLine
from geogst.core.geometries.shapes.polygons import Polygon
from geogst.core.mathematics.defaults import PRACTICAL_MIN_DIST
from geogst.core.geometries.shapes.shape import Shape
from geogst.core.geometries.shapes.segments import UnionPtSegment2D


@singledispatch
def intersect(
    shape1: Shape,
    shape2: Shape
) -> Optional[Union[List[UnionPtSegment2D], NotImplemented]]:

    return NotImplemented


def _bboxes_overlap(
    bbox1: Optional[tuple],
    bbox2: Optional[tuple],
) -> bool:
    """
    Cheap necessary-condition check for two axis-aligned 2D bounding boxes
    (each as xmin, ymin, xmax, ymax) to possibly intersect. Used to skip
    the much costlier exact segment-segment intersection math for pairs
    that plainly cannot intersect.
    """

    if bbox1 is None or bbox2 is None:
        return True

    xmin1, ymin1, xmax1, ymax1 = bbox1
    xmin2, ymin2, xmax2, ymax2 = bbox2

    return xmin1 <= xmax2 and xmax1 >= xmin2 and ymin1 <= ymax2 and ymax1 >= ymin2


@intersect.register(Segment)
def _(
    shape1: Segment,
    shape2: Shape
) -> Union[List[UnionPtSegment2D], type(NotImplemented)]:
    """
    Calculates the possible intersections between the provided segment and another shape.

    :param shape1: the first shape, a 2D segment
    :param shape2: the second generic Shape instance
    :return: the optional intersections, as a PointsSegments2D instance, or NotImplemented
    """

    shape1_bbox = shape1.bounding_box()

    if isinstance(shape2, Ln):

        if not _bboxes_overlap(shape1_bbox, shape2.bounding_box()):
            return []

        intersections = [
            curr_segment.intersect_segments2d(shape1)
            for curr_segment in shape2.segments()
            if curr_segment is not None and _bboxes_overlap(shape1_bbox, curr_segment.bounding_box())
        ]

    elif isinstance(shape2, MultiLine):

        if not _bboxes_overlap(shape1_bbox, shape2.bounding_box()):
            return []

        intersections = [
            curr_segment.intersect_segments2d(shape1)
            for line in shape2.lines()
            for curr_segment in line.segments()
            if curr_segment is not None and _bboxes_overlap(shape1_bbox, curr_segment.bounding_box())
        ]

    else:
        print(f"Error: second shape is a {type(shape2)} instance")
        return NotImplemented

    return list(filter(lambda val: val is not None, intersections))


class PreparedLine:
    """
    A Ln with its bounding box and per-segment bounding boxes computed once,
    so that testing it against many segments (e.g. many profile segments)
    does not repeat the same bbox_/segments() computation over and over.
    """

    __slots__ = ("bbox", "segments_with_bbox")

    def __init__(self, bbox, segments_with_bbox):

        self.bbox = bbox
        self.segments_with_bbox = segments_with_bbox


def prepare_line_for_intersection(
    line: Ln,
) -> PreparedLine:
    """
    Precomputes the bounding box and per-segment bounding boxes of a line,
    for reuse across many calls to intersect_segment_with_prepared_line.
    """

    bbox = line.bounding_box()
    segments = line.segments() or []

    segments_with_bbox = [
        (curr_segment, curr_segment.bounding_box())
        for curr_segment in segments
        if curr_segment is not None
    ]

    return PreparedLine(bbox, segments_with_bbox)


def intersect_segment_with_prepared_line(
    segment: Segment,
    segment_bbox: Optional[tuple],
    prepared_line: PreparedLine,
) -> List[UnionPtSegment2D]:
    """
    Same result as intersect(segment, line), but using a PreparedLine
    (precomputed bounding boxes) instead of recomputing them, and taking
    the caller-provided segment bounding box instead of recomputing it
    on every call.

    :param segment: the segment to test (e.g. a profile segment).
    :param segment_bbox: the precomputed bounding box of segment.
    :param prepared_line: the line to test against, already prepared.
    :return: the possible intersections.
    """

    if not _bboxes_overlap(segment_bbox, prepared_line.bbox):
        return []

    intersections = [
        curr_segment.intersect_segments2d(segment)
        for curr_segment, curr_segment_bbox in prepared_line.segments_with_bbox
        if _bboxes_overlap(segment_bbox, curr_segment_bbox)
    ]

    return list(filter(lambda val: val is not None, intersections))


class PreparedPolygon:
    """
    A Polygon with everything that does not depend on the intersecting
    segment computed once: the polygon bounding box, the boundary edges
    (outer ring plus inner holes) as coordinate arrays, and the ring edge
    arrays used by the winding-number point-in-polygon test.

    Intersecting one polygon against many profile segments (and, with
    parallel profiles, the same segments over and over) otherwise rebuilds
    a Segment object per polygon edge, and re-walks all polygon vertices
    in pure Python, on every single (profile-segment, polygon) pair.

    The boundary edges are kept as arrays, and turned into Segment instances
    only for the few edges that a given profile segment can actually meet.
    """

    __slots__ = ("bbox", "boundary_starts", "boundary_ends", "boundary_edges", "outer_ring", "inner_rings")

    def __init__(self, bbox, boundary_starts, boundary_ends, boundary_edges, outer_ring, inner_rings):

        self.bbox = bbox
        self.boundary_starts = boundary_starts
        self.boundary_ends = boundary_ends
        self.boundary_edges = boundary_edges
        self.outer_ring = outer_ring
        self.inner_rings = inner_rings

    def boundary_segment(self,
                         ndx: numbers.Integral,
                         ) -> Segment:
        """
        Returns the boundary edge of given index as a Segment instance,
        identical to the one that Ln.segments() would build for the ring.
        """

        return Segment(
            Point(self.boundary_starts[ndx]),
            Point(self.boundary_ends[ndx]),
        )


def _ring_edge_arrays(
    ring: Optional[Ln],
):
    """
    Decomposes a (closed) ring into the per-edge coordinate arrays used by
    the vectorized winding-number test: edge start x/y, edge end y, and
    edge x/y deltas.

    :param ring: the polygon ring, already closed by the Polygon constructor.
    :return: the edge arrays, or None when the ring has no usable geometry.
    """

    if ring is None:
        return None

    coords = ring.coords

    if coords is None or coords.shape[0] < 2:
        return None

    xy = np.asarray(coords[:, :2], dtype=float)

    vx = xy[:-1, 0]
    vy = xy[:-1, 1]
    wy = xy[1:, 1]

    return vx, vy, wy, xy[1:, 0] - vx, wy - vy


def prepare_polygon_for_intersection(
    polygon: Polygon,
) -> PreparedPolygon:
    """
    Precomputes the bounding box, boundary edge arrays and ring edge arrays
    of a polygon, for reuse across many calls to
    intersect_segment_with_prepared_polygon / point_in_prepared_polygon.
    """

    outer_line = polygon.outer_line()

    bbox = outer_line.bounding_box() if outer_line is not None else None

    # the rings are closed by the Polygon constructor, so that consecutive
    # coordinate pairs are exactly the edges that Ln.segments() would build

    ring_coords = [
        ring.coords for ring in [outer_line] + list(polygon.inner_lines() or [])
        if ring is not None and ring.coords is not None and ring.coords.shape[0] >= 2
    ]

    if ring_coords:
        boundary_starts = np.vstack([coords[:-1] for coords in ring_coords])
        boundary_ends = np.vstack([coords[1:] for coords in ring_coords])
    else:
        boundary_starts = None
        boundary_ends = None

    inner_rings = [
        edge_arrays for edge_arrays in
        [_ring_edge_arrays(inner_line) for inner_line in polygon.inner_lines() or []]
        if edge_arrays is not None
    ]

    return PreparedPolygon(
        bbox,
        boundary_starts,
        boundary_ends,
        _boundary_edge_arrays(boundary_starts, boundary_ends),
        _ring_edge_arrays(outer_line),
        inner_rings,
    )


def _boundary_edge_arrays(
    boundary_starts: Optional[np.ndarray],
    boundary_ends: Optional[np.ndarray],
):
    """
    Packs the boundary edges into the arrays used by the vectorized
    segment-straddling pre-check: edge start x/y, edge x/y deltas and the
    reciprocal of the edge lengths (zero for degenerate, zero-length edges,
    which are then never pruned).
    """

    if boundary_starts is None or boundary_starts.shape[0] == 0:
        return None

    starts_x = np.asarray(boundary_starts[:, 0], dtype=float)
    starts_y = np.asarray(boundary_starts[:, 1], dtype=float)

    deltas_x = np.asarray(boundary_ends[:, 0], dtype=float) - starts_x
    deltas_y = np.asarray(boundary_ends[:, 1], dtype=float) - starts_y

    lengths = np.hypot(deltas_x, deltas_y)

    inverse_lengths = np.where(lengths > 0.0, 1.0 / np.where(lengths > 0.0, lengths, 1.0), 0.0)

    return starts_x, starts_y, deltas_x, deltas_y, inverse_lengths


def _winding_number(
    x: numbers.Real,
    y: numbers.Real,
    ring_edge_arrays,
) -> numbers.Integral:
    """
    Vectorized counterpart of 'calculate_wn': the winding number of a point
    with respect to a closed ring, computed with numpy over all ring edges
    at once instead of one Python loop iteration per edge.
    """

    vx, vy, wy, dx, dy = ring_edge_arrays

    is_left = dx * (y - vy) - (x - vx) * dy

    upward_crossings = (vy <= y) & (wy > y) & (is_left > 0.0)
    downward_crossings = (vy > y) & (wy <= y) & (is_left < 0.0)

    return int(np.count_nonzero(upward_crossings)) - int(np.count_nonzero(downward_crossings))


def point_in_prepared_polygon(
    x: numbers.Real,
    y: numbers.Real,
    prepared_polygon: PreparedPolygon,
) -> bool:
    """
    Same result as Polygon.point_in_polygon, but on a PreparedPolygon and
    with a bounding box pre-check (a point outside the polygon bounding box
    has, provably, a null winding number).
    """

    if prepared_polygon.outer_ring is None:
        return False

    if prepared_polygon.bbox is not None:

        xmin, ymin, xmax, ymax = prepared_polygon.bbox

        if x < xmin or x > xmax or y < ymin or y > ymax:
            return False

    for inner_ring in prepared_polygon.inner_rings:
        if _winding_number(x, y, inner_ring) != 0:
            return False

    return _winding_number(x, y, prepared_polygon.outer_ring) != 0


def select_segments_inside_prepared_polygon(
    segments: List[Segment],
    prepared_polygon: PreparedPolygon,
) -> List[Segment]:
    """
    Same result as Polygon.select_intersecting_segments, i.e. the segments
    whose midpoint lies inside the polygon, but on a PreparedPolygon and
    reading the segment coordinates directly, without building a Point
    instance per tested midpoint.
    """

    inner_segments = []

    for curr_segment in segments:

        coords = curr_segment.coords

        # midpoint coordinates, as Segment.point_at_factor(0.5) computes them

        if point_in_prepared_polygon(
            coords[0, 0] + (coords[1, 0] - coords[0, 0]) * 0.5,
            coords[0, 1] + (coords[1, 1] - coords[0, 1]) * 0.5,
            prepared_polygon,
        ):
            inner_segments.append(curr_segment)

    return inner_segments


def _straddling_boundary_segments(
    segment: Segment,
    prepared_polygon: PreparedPolygon,
    tol: numbers.Real = PRACTICAL_MIN_DIST,
) -> List[Segment]:
    """
    Returns the polygon boundary segments that can possibly intersect the
    provided segment, i.e. those for which each segment straddles (within
    tolerance) the infinite line of the other one - a necessary condition
    for a real intersection, tested with numpy over all boundary edges at
    once.

    A bounding box pre-check is useless here, because a profile segment is
    typically long and oblique, so its bounding box contains whole polygons
    and every one of their edges; the point-to-line distances below prune
    instead on the actual segment direction.

    :param segment: the segment to test (e.g. a profile segment).
    :param prepared_polygon: the polygon to test against, already prepared.
    :param tol: the distance below which a point counts as lying on a line.
        Kept at the tolerance used by the exact intersection routine, so no
        boundary edge that it would report as intersecting gets pruned.
    :return: the boundary segments to submit to the exact intersection math.
    """

    boundary_edges = prepared_polygon.boundary_edges

    if boundary_edges is None:
        return []

    edges_x, edges_y, edges_dx, edges_dy, edges_inverse_length = boundary_edges

    start_x, start_y = segment.start_pt.x, segment.start_pt.y
    end_x, end_y = segment.end_pt.x, segment.end_pt.y

    segment_dx, segment_dy = end_x - start_x, end_y - start_y

    segment_length = np.hypot(segment_dx, segment_dy)
    segment_inverse_length = 1.0 / segment_length if segment_length > 0.0 else 0.0

    # signed distances of the two boundary edge ends from the segment line

    edge_start_offsets = (segment_dx * (edges_y - start_y) - (edges_x - start_x) * segment_dy) * segment_inverse_length
    edge_end_offsets = (segment_dx * (edges_y + edges_dy - start_y) - (edges_x + edges_dx - start_x) * segment_dy) * segment_inverse_length

    apart = ((edge_start_offsets > tol) & (edge_end_offsets > tol)) | \
            ((edge_start_offsets < -tol) & (edge_end_offsets < -tol))

    # signed distances of the two segment ends from the boundary edge lines

    segment_start_offsets = (edges_dx * (start_y - edges_y) - (start_x - edges_x) * edges_dy) * edges_inverse_length
    segment_end_offsets = (edges_dx * (end_y - edges_y) - (end_x - edges_x) * edges_dy) * edges_inverse_length

    apart |= ((segment_start_offsets > tol) & (segment_end_offsets > tol)) | \
             ((segment_start_offsets < -tol) & (segment_end_offsets < -tol))

    return [prepared_polygon.boundary_segment(ndx) for ndx in np.flatnonzero(~apart)]


def intersect_segment_with_prepared_polygon(
    segment: Segment,
    segment_bbox: Optional[tuple],
    prepared_polygon: PreparedPolygon,
) -> Optional[List[UnionPtSegment2D]]:
    """
    Intersects a segment with the boundaries (outer ring plus inner holes)
    of a prepared polygon, skipping the exact segment-segment math for the
    boundary edges that plainly cannot intersect.

    :param segment: the segment to test (e.g. a profile segment).
    :param segment_bbox: the precomputed bounding box of segment.
    :param prepared_polygon: the polygon to test against, already prepared.
    :return: the boundary intersections, or None when the segment lies
        entirely outside the polygon bounding box (so that the caller can
        skip the whole polygon: neither boundary crossings nor containment
        are then possible).
    """

    if prepared_polygon.outer_ring is None:
        return None

    if not _bboxes_overlap(segment_bbox, prepared_polygon.bbox):
        return None

    intersections = [
        curr_segment.intersect_segments2d(segment)
        for curr_segment in _straddling_boundary_segments(segment, prepared_polygon)
    ]

    return list(filter(lambda val: val is not None, intersections))


class PreparedPolygons:
    """
    A category's polygons, prepared for intersection and indexed by their
    bounding boxes, so that a profile segment can be tested against only
    the few polygons whose bounding box it overlaps, instead of against
    all the (possibly many thousands of) polygons of the layer.
    """

    __slots__ = ("polygons", "num_source_polygons", "_bboxes")

    def __init__(self,
                 polygons: List[Polygon]):

        self.num_source_polygons = len(polygons)

        prepared = [prepare_polygon_for_intersection(polygon) for polygon in polygons]

        self.polygons = [
            prepared_polygon for prepared_polygon in prepared
            if prepared_polygon.outer_ring is not None and prepared_polygon.bbox is not None
        ]

        if self.polygons:
            self._bboxes = np.array([prepared_polygon.bbox for prepared_polygon in self.polygons], dtype=float)
        else:
            self._bboxes = None

    def candidates(self,
                   bbox: Optional[tuple],
                   ) -> List[PreparedPolygon]:
        """
        Returns the prepared polygons whose bounding box overlaps the
        provided one, i.e. the only ones that can possibly intersect or
        contain a shape with that bounding box.
        """

        if self._bboxes is None:
            return []

        if bbox is None:
            return self.polygons

        xmin, ymin, xmax, ymax = bbox

        bboxes = self._bboxes

        overlapping = (bboxes[:, 0] <= xmax) & (bboxes[:, 2] >= xmin) & \
                      (bboxes[:, 1] <= ymax) & (bboxes[:, 3] >= ymin)

        return [self.polygons[ndx] for ndx in np.flatnonzero(overlapping)]


def intersect_segments3d(
    segment1: Segment,
    segment2: Segment,
    tol: numbers.Real = PRACTICAL_MIN_DIST
) -> Optional[Union[Point, Segment]]:
    """
    Determines the optional point or segment intersection between the segment pair.

    :param segment1: the first segment
    :param segment2: the second segment
    :param tol: the distance tolerance for collapsing a intersection segment into a point
    :return: the optional point or segment intersection between the segment pair.

    Examples:
      >>> s2 = Segment(Point(0,0,0), Point(1,0,0))
      >>> s1 = Segment(Point(0,0,0), Point(1,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0. 0. 0.]), end_pt=Point([1. 0. 0.]))
      >>> s1 = Segment(Point(-2,0,0), Point(-1,0,0))
      >>> intersect_segments3d(s1, s2)
      None
      >>> s1 = Segment(Point(-2,0,0), Point(0,0,0))
      >>> intersect_segments3d(s1, s2)
      Point([0. 0. 0.])
      >>> s1 = Segment(Point(-2,0,0), Point(0.5,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0. 0. 0.]), end_pt=Point([0.5 0.  0. ]))
      >>> s1 = Segment(Point(-2,0,0), Point(1,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0. 0. 0.]), end_pt=Point([1. 0. 0.]))
      >>> s1 = Segment(Point(-2,0,0), Point(2,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0. 0. 0.]), end_pt=Point([1. 0. 0.]))
      >>> s1 = Segment(Point(0,0,0), Point(0.5,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0. 0. 0.]), end_pt=Point([0.5 0.  0. ]))
      >>> s1 = Segment(Point(0.25,0,0), Point(0.75,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0.25 0.   0.  ]), end_pt=Point([0.75 0.   0.  ]))
      >>> s1 = Segment(Point(0.25,0,0), Point(1,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0.25 0.   0.  ]), end_pt=Point([1. 0. 0.]))
      >>> s1 = Segment(Point(0.25,0,0), Point(1.25,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0.25 0.   0.  ]), end_pt=Point([1. 0. 0.]))
      >>> s1 = Segment(Point(0,0,0), Point(1.25,0,0))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0. 0. 0.]), end_pt=Point([1. 0. 0.]))
      >>> s1 = Segment(Point(1,0,0), Point(1.25,0,0))
      >>> intersect_segments3d(s1, s2)
      Point([1. 0. 0.])
      >>> s2 = Segment(Point(0,0,0), Point(1,1,1))
      >>> s1 = Segment(Point(0.25,0.25,0.25), Point(0.75,0.75,0.75))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0.25 0.25 0.25]), end_pt=Point([0.75 0.75 0.75]))
      >>> s1 = Segment(Point(0.25,0.25,0.25), Point(1.75,1.75,1.75))
      >>> intersect_segments3d(s1, s2)
      Segment(start_pt=Point([0.25 0.25 0.25]), end_pt=Point([1. 1. 1.]))
      >>> s1 = Segment(Point(0.25,0.25,0.25), Point(1.75,0,1.75))
      >>> intersect_segments3d(s1, s2)
      Point([0.25 0.25 0.25])
      >>> s1 = Segment(Point(0.25,1,0.25), Point(0.75,0.75,0.75))
      >>> intersect_segments3d(s1, s2)
      Point([0.75 0.75 0.75])
      >>> s2 = Segment(Point(-1,-1,-1), Point(1,1,1))
      >>> s1 = Segment(Point(-1,1,1), Point(1,-1,-1))
      >>> intersect_segments3d(s1, s2)
      Point([0. 0. 0.])
    """

    s1_startpt_inside = segment1.segment_start_in(segment2)
    s2_startpt_inside = segment2.segment_start_in(segment1)

    s1_endpt_inside = segment1.segment_end_in(segment2)
    s2_endpt_inside = segment2.segment_end_in(segment1)

    elements = [s1_startpt_inside, s2_startpt_inside, s1_endpt_inside, s2_endpt_inside]

    if all(elements):
        return segment1.clone()

    if s1_startpt_inside and s1_endpt_inside:
        return segment1.clone()

    if s2_startpt_inside and s2_endpt_inside:
        return segment2.clone()

    if s1_startpt_inside and s2_startpt_inside:
        return point_or_segment3d(
            segment1.start_pt,
            segment2.start_pt,
            tol=tol
        )

    if s1_startpt_inside and s2_endpt_inside:
        return point_or_segment3d(
            segment1.start_pt,
            segment2.end_pt,
            tol=tol
        )

    if s1_endpt_inside and s2_startpt_inside:
        return point_or_segment3d(
            segment2.start_pt,
            segment1.end_pt,
            tol=tol
        )

    if s1_endpt_inside and s2_endpt_inside:
        return point_or_segment3d(
            segment1.end_pt,
            segment2.end_pt,
            tol=tol
        )

    if s1_startpt_inside:
        return segment1.start_pt.clone()

    if s1_endpt_inside:
        return segment1.end_pt.clone()

    if s2_startpt_inside:
        return segment2.start_pt.clone()

    if s2_endpt_inside:
        return segment2.end_pt.clone()

    shortest_segm_or_pt = shortest_segment_or_point3d(
        segment1,
        segment2,
        tol=tol
    )

    if not shortest_segm_or_pt:
        return None

    if not isinstance(shortest_segm_or_pt, Point):
        return None

    inters_pt = shortest_segm_or_pt

    if not segment1.contains_pt(inters_pt):
        return None

    if not segment2.contains_pt(inters_pt):
        return None

    return inters_pt
