财富算法中用于计算沃罗诺伊图的交点

Question

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

Answer 1

这是你要问的吗？

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)     

几何结构。让我们有两点A = [x1, y1]和B = [x2, y2]，y = c打赌他的水平线，其中y1 < c和y2 < c。我们正在寻找圆C的中心O = [x_center, y_center]，该圆穿过A点和B点，并且与线y = c相切。

1. 写下通过点的线AB的方程A和B。这使我们能够找到AB和y = c线之间的交点P = [x_P, y_P]坐标。在不失去一般性的情况下，假设点B位于点A和P点之间。

2. 我们的下一个目标是找到圆C与直线y = c相切的点T的坐标。为此，我们将首先找到点之间的距离t = PTP和T.之后，我们只需要将距离从点tP沿线y = c移动，这意味着我们发现T的x坐标x_T为x_T = x_P - t或x_T = x_P + t(有两种解决方案)。T的y坐标是y_T = c。

3. 考虑两个三角形PAT和PTB.他们有一个共同的角度angle TPA = angle BPT。此外，因为y = c与围绕三角形ABTC外接的圆相切，我们可以推断出angle PTA = angle PBT。因此，三角形PAT和PTB是相似的。

4. PAT和PTB之间的相似性产生了相似性恒等式PT / PA = PB / PT可以重写为t^2 = PT^2 = PA * PB。因此，我们可以找到t = sqrt(PA * PB)(这就是代码中t的内容)。

5. 最后，我们发现x_T = x_P - t或x_T = x_P + t，以更接近段AB者为准。圆C的中心O是垂直线x = x_T与段AB的正交平分线(穿过AB中点并垂直于AB的线)的交

   点

备注：当线AB平行于线y = c时，存在退化情况。然后图片更加对称，可以使用稍微不同的方法，我已经包含在下面的代码中

在 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]))

这是另一种方法，它也处理了y1 = y2时的退化情况，即点A和B形成一条平行于y = c的直线。

建设。平面上与直线y = c距离相等的所有点和该点A = [x1, y1]形成具有直线y = c和焦点A的抛物线。所以简单地写出这个抛物线的方程：

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)

这同样适用于点B = [x2, y2]和直线y = c，其中抛物线与焦点B和直线a = c的方程为

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

因此，与A、B和y = c距离相等的点应该是两条抛物线的两个交点。就方程而言，您只需求解方程

( 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)

对于未知的x坐标O.然后通过将解代入上面的抛物线方程之一来找到y坐标。

以下是一些python代码：

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]))