Issue with concave bounding box

[casparvitch]

Hiya, lovely package.

Everything works nicely for my use-case, except where I try to define a 'heavily' concave polygon bounding box. The issue is that the input Polygon appears to be clockwise ordered after input: '_order_points' and a convex polygon bounding box seems to be assumed elsewhere too.

This is quite easy to see with a slight adjustment to the Quick Start example:

import foronoi

points = [(2.5, 2.5), (4, 7.5), (7.5, 2.5), (6, 7.5), (4, 4), (3, 3), (6, 3)]

poly_nodes = [(2.5, 10), (5, 10), (10, 5), (10, 2.5), (5, 0), (2.5, 0), (0, 2.5), (0, 5)]
poly_nodes.insert(3, (4, 3))

v = foronoi.Voronoi(foronoi.Polygon(poly_nodes))

v.create_diagram(points=points)
foronoi.Visualizer(v).plot_sites().plot_edges(show_labels=False).plot_vertices().show()

Is it possible to generalize this to allow for a Polygon with pre-ordered vertices (in any manner)? Additionally, not all of the interior regions appear to be clipped to the boundary box.

Cheers!

[casparvitch]

I've had a bit of an explore, the issue is a little more complex than I expected.

Issue overcome:

• keep polygon vertices/points ordered, inserted new vertices to correct position in that ordering.

Issue created:

• Each cell/'source' point does not keep track of which vertices are on its boundary (e.g. [vertex.xy for vertex in site.vertices()] gives incorrect answer), even for concave bounding polygons.

Issues that remain:

• When clipping the diagram with a convex polygon, I have not worked out how to detect that an interior edge (i.e. not an edge in bounding polygon) is not completely inside the polygon. For example in the example given above the internal edge for P5 leaves the polygon and then re-enters.
• If I could detect this occurrence, I would need to decide what happens to the disconnected sections. For P5 above this looks simple, but is more complex for the example below. Other internal edges that would now be 'broken'/'hanging' would need to be continued to the boundary in some intelligent/self-consistent way. Of course an alternative is to build this restriction into the algorithm, but that sounds messy.

Example

import foronoi
from foronoi.contrib import ConcavePolygon

points = [(2.5, 2.5), (4, 7.5), (7.5, 2.5), (6, 7.5), (4, 4), (3, 3), (6, 2)]

poly_nodes = [(2.5, 10), (5, 10), (10, 5), (10, 2.5), (5, 0), (2.5, 0), (0, 2.5), (0, 5)]
poly_nodes.insert(3, (4, 3))

# v = foronoi.Voronoi(foronoi.Polygon(poly_nodes))
v = foronoi.Voronoi(ConcavePolygon(poly_nodes))

v.create_diagram(points=points)
foronoi.Visualizer(v).plot_sites().plot_edges(show_labels=False).plot_vertices().show()

Example output

New contrib

fornoi/contrib/concave_polygon.py:

from foronoi import Polygon

from foronoi.graph import Coordinate, Vertex, HalfEdge
from foronoi.graph.algebra import Algebra

import numpy as np

class ConcavePolygon(Polygon):
    def __init__(self, tuples):
        """
        A bounding polygon that will clip the Voronoi diagram.

        Parameters
        ----------
        tuples: list[(float, float)]
            x,y-coordinates of the polygon's vertices -> must be pre-ordered in chosen manner.
        """

        self._observers = []
        self._children = []
        self._root_sender = self
        self._child_sender = self
        points = [Coordinate(x, y) for x, y in tuples]
        self.points = points
        min_y = min([p.yd for p in self.points])
        min_x = min([p.xd for p in self.points])
        max_y = max([p.yd for p in self.points])
        max_x = max([p.xd for p in self.points])
        center = Coordinate((max_x + min_x) / 2, (max_y + min_y) / 2)
        self.min_y, self.min_x, self.max_y, self.max_x, self.center = (
            min_y,
            min_x,
            max_y,
            max_x,
            center,
        )

        self.polygon_vertices = [Vertex(point.xd, point.yd) for point in self.points]

    def _get_ordered_vertices(self, vertices):
        return [vertex for vertex in vertices if vertex.xd is not None]

    def _finish_edge(self, edge):
        # Sweep line position
        sweep_line = self.min_y - abs(self.max_y)

        # Start should be a breakpoint
        start = edge.get_origin(y=sweep_line, max_y=self.max_y)

        # End should be a vertex
        end = edge.twin.get_origin(y=sweep_line, max_y=self.max_y)

        # Get point of intersection
        point, prev_idx = self._get_intersection_point(end, start)

        # Create vertex
        v = Vertex(point.x, point.y) if point is not None else Vertex(None, None)
        v.connected_edges.append(edge)
        edge.origin = v
        if point is not None:
            self.polygon_vertices.insert(prev_idx + 1, v)
            self.points.insert(prev_idx + 1, Coordinate(*v.xy))

        return edge

    def _get_intersection_point(self, orig, end):
        p = self.points + [self.points[0]]
        intersection_points = []
        prev_idx_lst = []

        intersection = None  # point of intersection
        # index in self.points (i.e. polygon_vertices) this intersection is immediately above
        # i.e.: interesection occurs between self.points[prev_idx] and self.points[prev_idx + 1]
        # so new vertex should be inserted as self.polygon_vertices.insert(prev_idx + 1, NEW_VERTEX)
        prev_idx = None

        for i in range(0, len(p) - 1):
            intersection_point = Algebra.get_intersection(orig, end, p[i], p[i + 1])
            if intersection_point:
                intersection_points.append(intersection_point)
                prev_idx_lst.append(i)

        if not intersection_points:
            return None, None

        max_distance = Algebra.distance(orig, end)

        # Find the intersection point that is furthest away from the start
        if intersection_points:
            distances = [Algebra.distance(orig, p) for p in intersection_points]
            distances = [i for i in distances if i <= max_distance]
            if distances:
                idx = np.argmax(distances)
                intersection = intersection_points[idx]
                prev_idx = prev_idx_lst[idx]

        return intersection, prev_idx

foronoi/contrib/__init__.py:

from foronoi.contrib.bounding_circle import BoundingCircle
from foronoi.contrib.concave_polygon import ConcavePolygon