I am finding it very difficult to follow fortune's algorithm, I have gone through all the resources available on the net and fairly understood the theory behind it. But when I am implementing it on my own the little details which were left out in the sources have been a real pain for me. The intersection points, which I understood are the center of a circle with two given points on it and a horizontal tangent with equation y = c but working through the equations I am unable to get the co-ordinates of the center(the intersection point). Can someone please help me on how to find the co-ordinates of the intersection points.
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]))