用于计算 Voronoi 图的 Fortune 算法中的交点

Question

我发现很难遵循 fortune 的算法，我浏览了网上所有可用的资源并且相当了解其背后的理论。但是当我自己实现它时，源代码中遗漏的小细节对我来说真的很痛苦。据我所知，交点是一个圆心，上面有两个给定点和一个水平切线，方程 y = c 但是通过方程我无法得到中心的坐标（交点） .有人可以帮助我如何找到交点的坐标。

Answer 1

Is this what you are asking about?

INPUT:
two given points:
[x1, y1]
[x2, y2]
and a real number: 
c
where y = c is a horizontal line 

CALCULATION of the coordinates of the centre of the circle
that passes through the two given points [x1, y1] and [x2, y2] 
and is tangent to the horizontal line y = c:

    k = (x2 - x1) / (y2 - y1)
    x_P = x1 + k * (c - y1)
    t = sqrt((x_P - x1)^2 + (c - y1)^2) * math.sqrt((x_P - x2)^2 + (c - y2)^2)
    t = sqrt(t)
    x_center = x_P - t * k / abs(k)
    y_center = (y1 + y2) / 2  - k * (x_center - (x1 + x2)/2)

Geometric Construction. Let us have two points A = [x1, y1] and B = [x2, y2] and let y = c bet he horizontal line, where y1 < c and y2 < c. We are looking for the center O = [x_center, y_center] of the circle C that passes through the points A and B and is tangent to the line y = c.

Write the equation of the line AB that passes through the points A and B. This allows us to find the coordinates of the intersection point P = [x_P, y_P] between the lines AB and y = c. Without loss of generality, assume that point B is between points A and P.
Our next goal is to find the coordinates of the point T at which the circle C is tangent to the line y = c. To do that, we will first find be the distance t = PT between the points P and T. After that, we simply have to move distance t from the point P along the line y = c, which means we find the x coordinate x_T of T as either x_T = x_P - t or x_T = x_P + t (there are two solutions). The y coordinate of T is y_T = c.
Consider the two triangles PAT and PTB. They share a common angle angle TPA = angle BPT. Furthermore, because y = c is tangent to the circle C circumscribed around the triangle ABT, we can deduce that angle PTA = angle PBT. Therefore triangles PAT and PTB are similar.
The similarity between PAT and PTB yields the similarity identity PT / PA = PB / PT which can be rewritten as t^2 = PT^2 = PA * PB. Thus, we can find that t = sqrt(PA * PB) (this is what t is in the code).
Finally, we find that x_T = x_P - t or x_T = x_P + t, whichever is closer to the segment AB. The center O of the circle C is the intersection of the vertical line x = x_T and the orthogonal bisector of segment AB (the line through the midpoint of AB and perpendicular to AB)

Remark: There is a degenerate case when the line AB is parallel to the line y = c. Then the picture is more symmetric and one can use a slightly different method, which I have included in the code below

Test code in python:

import numpy as np
import math


def find_center(A, B, c):
    x1 = A[0]
    y1 = A[1]
    x2 = B[0]
    y2 = B[1]
    if abs(y1 - y2) < 1e-7:
        radius = ((2*c - y1 - y2)**2 + (x2 - x1)**2) / (8*(c - y1))
        x_center = (x1 + x2) / 2
        y_center = c - radius 
    else:
        k = (x2 - x1) / (y2 - y1)
        x_P = x1 + k * (c - y1)
        t = math.sqrt((x_P - x1)**2 + (c - y1)**2) * math.sqrt((x_P - x2)**2 + (c - y2)**2)
        t = math.sqrt(t)
        x_center = x_P - t * k / abs(k)
        y_center = (y1 + y2) / 2  - k * (x_center - (x1 + x2)/2)     
    return np.array([x_center, y_center])

A = np.array([0, 0])
B = np.array([3, 1.5])
c = 5

O = find_center(A, B, c)

print('print the distances between the calculated center and the two points and the horizontal line: ')
print((np.linalg.norm(A - O), np.linalg.norm(B - O),  c - O[1]))

Edit. Here is another method, which also takes care of the degenerate case when y1 = y2, i.e. the points A and B form a line parallel to y = c.

Construction. All points from the plane that are at an equal distance from the line y = c and the point A = [x1, y1] form a parabola with directrix y = c and focus A. So simply write the equation of this parabola:

distance(X, A)^2 = distance(X, {y=c})^2

or as a quadratic equation:

(x - x1)^2  +  (y - y1)^2 = (c - y)^2

which simplifies to 

y = ( x^2 - 2*x1* x + x1^2 + y1^2 - c^2 ) / (y1 - c)

The same applies to point B = [x2, y2] and line y = c, where the equation of the parabola with focus B and directrix a = c is

y = ( x^2 - 2*x2* x + x2^2 + y2^2 - c^2 ) / (y2 - c)

Therefore, the points that are at equal distance from A, B and y = c should be the two intersection points of the two parabolas. In terms of equations, you simply have to solve the equation

 ( x^2 - 2*x1* x + x1^2 + y1^2 - c^2 ) / (y1 - c) =  ( x^2 - 2*x2* x + x2^2 + y2^2 - c^2 ) / (y2 - c)

for the unknown x coordinate of O. Then the y coordinate is found by plugging the solution in one of the parabola equations above.

Here is some python code:

import numpy as np
import math

def sqe(a):
    def sqe(a):
    if abs(a[0]) < 1e-7:
        return - a[2] / a[1], - a[2] / a[1]
    D = math.sqrt(a[1]**2 - 4*a[0]*a[2])
    return ( - a[1] - D ) /(2*a[0]), ( - a[1] + D ) / (2*a[0])

def find_center(A, B, c):
    x1 = A[0]
    y1 = A[1]
    x2 = B[0]
    y2 = B[1]     
    a1 = np.array([1, -2*x1, x1**2 + y1**2 - c**2])
    a2 = np.array([1, -2*x2, x2**2 + y2**2 - c**2])
    x_center = sqe((y1-c)*a2 - (y2-c)*a1)[1]
    y_center = (a1[0]*x_center**2 + a1[1]*x_center + a1[2])  / (2*y1-2*c)      
    return np.array([x_center, y_center])

A = np.array([0, 0])
B = np.array([3, 0.5])
c = 5

O = find_center(A, B, c)
print('print the distances between the calculated center and the two points and the horizontal line: ')
print((np.linalg.norm(A - O), np.linalg.norm(B - O),  c - O[1]))

用于计算 Voronoi 图的 Fortune 算法中的交点

Intersection points in Fortune's algoritm for computing Voronoi diagram

algorithm

math

geometry

voronoi

computational-geometry