Source code for sect.core.delaunay.triangulation

from __future__ import annotations

from collections import deque
from functools import partial
from itertools import (accumulate,
                       chain,
                       repeat)
from typing import (Callable,
                    Dict,
                    Iterable,
                    List,
                    Optional,
                    Sequence,
                    Set,
                    Type)

from decision.partition import coin_change
from ground.base import (Context,
                         Location,
                         Orientation,
                         Relation)
from ground.hints import (Contour,
                          Point,
                          Polygon,
                          Segment)
from reprit.base import generate_repr

from sect.core.utils import (flatten,
                             pairwise)
from .events_queue import EventsQueue
from .hints import (PointInCircleLocator,
                    SegmentsRelater)
from .quad_edge import (QuadEdge,
                        edge_to_neighbours,
                        edges_with_opposites)
from .utils import (ceil_log2,
                    complete_vertices,
                    contour_to_oriented_edges_endpoints,
                    normalize_contour_vertices,
                    to_distinct,
                    to_endpoints)


[docs]class Triangulation: """Represents triangulation."""
[docs] @classmethod def constrained_delaunay(cls, polygon: Polygon, *, extra_constraints: Sequence[Segment] = (), extra_points: Sequence[Point] = (), context: Context) -> Triangulation: """ Constructs constrained Delaunay triangulation of given polygon (with potentially extra points and constraints). Based on * divide-and-conquer algorithm by L. Guibas & J. Stolfi for generating Delaunay triangulation, * algorithm by S. W. Sloan for adding constraints to Delaunay triangulation, * clipping algorithm by F. Martinez et al. for deleting in-hole triangles. Time complexity: ``O(vertices_count * log vertices_count)`` for convex polygons without extra constraints, ``O(vertices_count ** 2)`` otherwise Memory complexity: ``O(vertices_count)`` where .. code-block:: python vertices_count = (len(polygon.border.vertices) + sum(len(hole.vertices) for hole in polygon.holes) + len(extra_points) + len(extra_constraints)) Reference: http://www.sccg.sk/~samuelcik/dgs/quad_edge.pdf https://www.newcastle.edu.au/__data/assets/pdf_file/0019/22519/23_A-fast-algortithm-for-generating-constrained-Delaunay-triangulations.pdf https://doi.org/10.1016/j.advengsoft.2013.04.004 http://www4.ujaen.es/~fmartin/bool_op.html :param polygon: target polygon. :param extra_points: additional points to be presented in the triangulation. :param extra_constraints: additional constraints to be presented in the triangulation. :param context: geometric context. :returns: triangulation of the border, holes & extra points considering constraints. """ border, holes = polygon.border, polygon.holes if extra_points: border, holes, extra_points = complete_vertices( border, holes, extra_points, context ) result = cls.delaunay(list(chain(border.vertices, flatten(hole.vertices for hole in holes), extra_points)), context=context) border_edges = context.contour_segments(border) constrain(result, chain(border_edges, flatten(map(context.contour_segments, holes)), extra_constraints)) bound(result, border_edges) cut(result, holes) result._triangular_holes_vertices.update( frozenset(hole.vertices) for hole in holes if len(hole.vertices) == 3 ) return result
[docs] @classmethod def delaunay(cls, points: Sequence[Point], *, context: Context) -> Triangulation: """ Constructs Delaunay triangulation of given points. Based on divide-and-conquer algorithm by L. Guibas & J. Stolfi. Time complexity: ``O(len(points) * log len(points))`` Memory complexity: ``O(len(points))`` Reference: http://www.sccg.sk/~samuelcik/dgs/quad_edge.pdf :param points: 3 or more points to triangulate. :param context: geometric context. :returns: triangulation of the points. """ points = sorted(to_distinct(points)) lengths = coin_change(len(points), base_cases) result = [cls._initialize_triangulation(points[start:stop], context) for start, stop in pairwise(accumulate((0,) + lengths))] for _ in repeat(None, ceil_log2(len(result))): parts_to_merge_count = len(result) // 2 * 2 result = ([merge(result[offset], result[offset + 1]) for offset in range(0, parts_to_merge_count, 2)] + result[parts_to_merge_count:]) return result[0]
__slots__ = ('context', 'left_side', 'right_side', '_triangular_holes_vertices') def __init__(self, left_side: QuadEdge, right_side: QuadEdge, context: Context) -> None: self.context, self.left_side, self.right_side = (context, left_side, right_side) self._triangular_holes_vertices = set() __repr__ = generate_repr(__init__)
[docs] def delete(self, edge: QuadEdge) -> None: """Deletes given edge from the triangulation.""" if edge is self.right_side or edge.opposite is self.right_side: self.right_side = self.right_side.right_from_end.opposite if edge is self.left_side or edge.opposite is self.left_side: self.left_side = self.left_side.left_from_start edge.delete()
[docs] def triangles(self) -> List[Contour]: """Returns triangles of the triangulation.""" vertices_sets = to_distinct( frozenset((edge.start, edge.end, edge.left_from_start.end)) for edge in to_edges(self) if (edge.left_from_start.end == edge.opposite.right_from_start.end and (edge.orientation_of(edge.left_from_start.end) is Orientation.COUNTERCLOCKWISE)) ) contour_cls, orienteer = (self.context.contour_cls, self.context.angle_orientation) return [contour_cls(normalize_contour_vertices(list(vertices), orienteer)) for vertices in vertices_sets if vertices not in self._triangular_holes_vertices]
@classmethod def _initialize_triangulation(cls, points: Sequence[Point], context: Context) -> Triangulation: return base_cases[len(points)](cls, points, context)
def bound(triangulation: Triangulation, border_edges: Sequence[Segment]) -> None: border_endpoints = {to_endpoints(edge) for edge in border_edges} non_boundary = {edge for edge in to_unique_boundary_edges(triangulation) if to_endpoints(edge) not in border_endpoints} while non_boundary: edge = non_boundary.pop() candidates = edge_to_neighbours(edge) triangulation.delete(edge) non_boundary.update(candidate for candidate in candidates if to_endpoints(candidate) not in border_endpoints) def connect(base_edge: QuadEdge, point_in_circle_locator: PointInCircleLocator) -> None: while True: left_candidate, right_candidate = ( to_left_candidate(base_edge, point_in_circle_locator), to_right_candidate(base_edge, point_in_circle_locator) ) if left_candidate is right_candidate is None: break base_edge = ( right_candidate.connect(base_edge.opposite) if (left_candidate is None or right_candidate is not None and (point_in_circle_locator(right_candidate.end, left_candidate.end, base_edge.end, base_edge.start) is Location.INTERIOR)) else base_edge.opposite.connect(left_candidate.opposite) ) def constrain(triangulation: Triangulation, constraints: Iterable[Segment]) -> None: endpoints = {to_endpoints(edge) for edge in to_edges(triangulation)} inner_edges = to_unique_inner_edges(triangulation) point_in_circle_locator, segments_relater = ( triangulation.context.locate_point_in_point_point_point_circle, triangulation.context.segments_relation ) for constraint in constraints: constraint_endpoints = to_endpoints(constraint) if constraint_endpoints in endpoints: continue crossings = detect_crossings(inner_edges, constraint, segments_relater) inner_edges.difference_update(crossings) endpoints.difference_update(to_endpoints(edge) for edge in crossings) new_edges = resolve_crossings(crossings, constraint, segments_relater) set_criterion({edge for edge in new_edges if to_endpoints(edge) != constraint_endpoints}, point_in_circle_locator) endpoints.update(to_endpoints(edge) for edge in new_edges) inner_edges.update(new_edges) def cut(triangulation: Triangulation, holes: Sequence[Contour]) -> None: if not holes: return events_queue = EventsQueue(triangulation.context) for edge in to_unique_inner_edges(triangulation): events_queue.register_edge(edge, from_first=True, is_counterclockwise_contour=True) orienteer = triangulation.context.angle_orientation for hole in holes: for endpoints in contour_to_oriented_edges_endpoints( hole, clockwise=True, orienteer=orienteer ): events_queue.register_segment(endpoints, from_first=False, is_counterclockwise_contour=False) for event in events_queue.sweep(): if event.from_first and event.inside: triangulation.delete(event.edge) def detect_crossings(inner_edges: Iterable[QuadEdge], constraint: Segment, segments_relater: SegmentsRelater) -> List[QuadEdge]: return [edge for edge in inner_edges if segments_relater(edge, constraint) is Relation.CROSS] def edge_should_be_swapped(edge: QuadEdge, point_in_circle_locator: PointInCircleLocator ) -> bool: return (is_convex_quadrilateral_diagonal(edge) and (point_in_circle_locator(edge.right_from_start.end, edge.start, edge.end, edge.left_from_start.end) is Location.INTERIOR or (point_in_circle_locator(edge.left_from_start.end, edge.end, edge.start, edge.right_from_start.end) is Location.INTERIOR))) def find_base_edge(first: Triangulation, second: Triangulation) -> QuadEdge: while True: if (first.right_side.orientation_of(second.left_side.start) is Orientation.COUNTERCLOCKWISE): first.right_side = first.right_side.left_from_end elif (second.left_side.orientation_of(first.right_side.start) is Orientation.CLOCKWISE): second.left_side = second.left_side.right_from_end else: break base_edge = second.left_side.opposite.connect(first.right_side) if first.right_side.start == first.left_side.start: first.left_side = base_edge.opposite if second.left_side.start == second.right_side.start: second.right_side = base_edge return base_edge def is_convex_quadrilateral_diagonal(edge: QuadEdge) -> bool: return (edge.right_from_start.orientation_of(edge.end) is Orientation.COUNTERCLOCKWISE is edge.right_from_end.opposite.orientation_of( edge.left_from_start.end ) is edge.left_from_end.orientation_of(edge.start) is edge.left_from_start.opposite.orientation_of( edge.right_from_start.end )) def merge(first: Triangulation, second: Triangulation) -> Triangulation: connect(find_base_edge(first, second), first.context.locate_point_in_point_point_point_circle) return type(first)(first.left_side, second.right_side, first.context) def resolve_crossings(crossings: List[QuadEdge], constraint: Segment, segments_relater: SegmentsRelater) -> List[QuadEdge]: result = [] crossings = deque(crossings, maxlen=len(crossings)) while crossings: edge = crossings.popleft() if is_convex_quadrilateral_diagonal(edge): edge.swap() if segments_relater(edge, constraint) is Relation.CROSS: crossings.append(edge) else: result.append(edge) else: crossings.append(edge) return result def set_criterion(target_edges: Set[QuadEdge], point_in_circle_locator: PointInCircleLocator) -> None: while True: edges_to_swap = { edge for edge in target_edges if edge_should_be_swapped(edge, point_in_circle_locator) } if not edges_to_swap: break for edge in edges_to_swap: edge.swap() target_edges.difference_update(edges_to_swap) def to_boundary_edges(triangulation: Triangulation) -> Iterable[QuadEdge]: return edges_with_opposites(to_unique_boundary_edges(triangulation)) def to_edges(triangulation: Triangulation) -> Iterable[QuadEdge]: return edges_with_opposites(to_unique_edges(triangulation)) def to_left_candidate(base_edge: QuadEdge, point_in_circle_locator: PointInCircleLocator ) -> Optional[QuadEdge]: result = base_edge.opposite.left_from_start if base_edge.orientation_of(result.end) is not Orientation.CLOCKWISE: return None while (point_in_circle_locator(result.left_from_start.end, base_edge.end, base_edge.start, result.end) is Location.INTERIOR and (base_edge.orientation_of(result.left_from_start.end) is Orientation.CLOCKWISE)): next_candidate = result.left_from_start result.delete() result = next_candidate return result def to_right_candidate(base_edge: QuadEdge, point_in_circle_locator: PointInCircleLocator ) -> Optional[QuadEdge]: result = base_edge.right_from_start if (base_edge.orientation_of(result.end) is not Orientation.CLOCKWISE): return None while (point_in_circle_locator(result.right_from_start.end, base_edge.end, base_edge.start, result.end) is Location.INTERIOR and (base_edge.orientation_of(result.right_from_start.end) is Orientation.CLOCKWISE)): next_candidate = result.right_from_start result.delete() result = next_candidate return result def to_unique_boundary_edges(triangulation: Triangulation ) -> Iterable[QuadEdge]: start = triangulation.left_side edge = start while True: yield edge if edge.right_from_end is start: break edge = edge.right_from_end def to_unique_edges(triangulation: Triangulation) -> Iterable[QuadEdge]: visited_edges = set() is_visited, visit_multiple = (visited_edges.__contains__, visited_edges.update) queue = [triangulation.left_side, triangulation.right_side] while queue: edge = queue.pop() if is_visited(edge): continue yield edge visit_multiple((edge, edge.opposite)) queue.extend((edge.left_from_start, edge.left_from_end, edge.right_from_start, edge.right_from_end)) def to_unique_inner_edges(triangulation: Triangulation) -> Set[QuadEdge]: return (set(to_unique_edges(triangulation)) .difference(to_boundary_edges(triangulation))) BaseCase = Callable[ [Type[Triangulation], Sequence[Point], Context], Triangulation ] base_cases = {} # type: Dict[int, BaseCase] register_base_case = partial(partial, base_cases.setdefault) @register_base_case(2) def triangulate_two_points(cls: Type[Triangulation], sorted_points: Sequence[Point], context: Context) -> Triangulation: first_edge = QuadEdge.from_endpoints(*sorted_points, context=context) return cls(first_edge, first_edge.opposite, context) @register_base_case(3) def triangulate_three_points(cls: Type[Triangulation], sorted_points: Sequence[Point], context: Context) -> Triangulation: left_point, mid_point, right_point = sorted_points first_edge, second_edge = (QuadEdge.from_endpoints(left_point, mid_point, context=context), QuadEdge.from_endpoints(mid_point, right_point, context=context)) first_edge.opposite.splice(second_edge) orientation = first_edge.orientation_of(right_point) if orientation is Orientation.COUNTERCLOCKWISE: second_edge.connect(first_edge) return cls(first_edge, second_edge.opposite, context) elif orientation is Orientation.CLOCKWISE: third_edge = second_edge.connect(first_edge) return cls(third_edge.opposite, third_edge, context) else: return cls(first_edge, second_edge.opposite, context)