scipy.spatial.Voronoi: How to know where a ray crosses a given line?

Question

Good day everyone,

I have the following code segment:

import numpy as np
from random import randint
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d

NUM_OF_POINTS = 20

points = []
for i in range (0, NUM_OF_POINTS):
    points.append([randint(0, 500), randint(0, 500)]) 

points = np.array(points)
vor = Voronoi(points)
voronoi_plot_2d(vor)
plt.show()

That produces Voronoi plots such as this:

My goal is to find where the 'rays' (the lines that go out of the plot, dashed or solid) intersect with a given line (for example x = 500). How may I go about doing this?

I have already tried using ridge_vertices list in the Voronoi object, however, these 'rays' are associated with only a single vertex in the list, so I cannot figure out the line equation.

Edit:

My ultimate goal with this is, given the borders of the plane, to find the points that intersect with these borders for a given edge cell. For example, given the edge cell in the upper left, and the borders y = -50 and x = 525, I would find the points marked with the red X's.

So if you have any insights to this, they would be appreciated.

Thank you.

Answer 1

1. Corners are trivial because you know the x and the y coordinates.

2. The edges in a Voronoi diagram are equidistant to the centres of the two cells that are separated by that edge, which naturally includes the endpoint of a "ray" (in your terminology). Let the two centres be the points (x1, y_1) and (x_2, y_2), and the intersection of the ray with the border be (x*, y*), then the following holds:

(1) (x*-x_1)^2 + (y*-y_1)^2 = d^2

(2) (x*-x_2)^2 + (y*-y_2)^2 = d^2

You know either x* or y* because they are defined by the border. Then you have two equations, and two unknowns (x* or y* and d). Assuming you know y*, then you arrive at the following solution for x*:

x* = ((y*-y_1)^2 - (y*-y_2)^2 + x_1^2 - x_2^2) / (2 * (x_1 - x_2))

---

Now how do you figure out which pairs of points (x_1, y_1) and (x_2, y_2) to pick?

As a first pass, I would brute force it:

(1) Iterate over all combinations of points (n * (n-1) / 2 per border, so not that many), find x* or y*, respectively. That gives you a list (x_1, y_1), (x_2, y_2), (x*, y*) of potential solutions.

(2) For each candidate (x*, y*) pair, I would find the 2 nearest neighbours in your original set of data points (efficiently via scipy.spatial.KDTree). If these points are not (x_1, y_1) and (x_2, y_2), discard the candidate solution (x*, y*).

Finding nearest neighbours in a KD-tree is O(n log n) (IIRC), so you are still O(n^2) for the whole procedure.