Find Arc/Circle equation given three points in space (3D)

Question

Given 3 points in space (3D): A = (x1, y1, z1), B = (x2, y2, z2) C = (x3, y3, z3); then how to find the center and radius of the circle (arc) that passes through these three points, i.e. find circle equation? Using Python and Numpy here is my initial code

import numpy as np
A = np.array([x1, y1, z1])
B = np.array([x2, y2, z2])
C = np.array([x3, y3, z3])

#Find vectors connecting the three points and the length of each vector
AB = B - A 
BC = C - B 
AC = C - A 

# Triangle Lengths
a = np.linalg.norm(AB)
b = np.linalg.norm(BC)
c = np.linalg.norm(AC)

From the Circumradius definition, the radius can be found using:

R = (a * b * c) / np.sqrt(2.0 * a**2 * b**2 +
                          2.0 * b**2 * c**2 +
                          2.0 * c**2 * a**2 -
                          a**4 - b**4 - c**4)

However, I am having problems finding the Cartesian coordinates of the center. One possible solution is to use the "Barycentric Coordinates" of the triangle points to find the trilinear coordinates of the circumcenter (Circumcenter).

First (using this source) we find the circumcenter barcyntric coordinates:

#barcyntric coordinates of center
b1 = a**2 * (b**2 + c**2 - a**2)
b2 = b**2 * (c**2 + a**2 - b**2)
b3 = c**2 * (a**2 + b**2 - c**2)

Then the Cartesian coordinates of the center (P) would be:

Px = (b1 * A[0]) + (b2 * B[0]) + (b3 * C[0])
Py = (b1 * A[1]) + (b2 * B[1]) + (b3 * C[1])
Pz = (b1 * A[2]) + (b2 * B[2]) + (b3 * C[2])   

However, the barcyntric coordinates values above do not seem to be correct. When solved with an example of known values, the radius is correct, but the coordinates of the center are not.

Example: For these three points:

A = np.array([2.0, 1.5, 0.0])
B = np.array([6.0, 4.5, 0.0])
C = np.array([11.75, 6.25, 0.0]) 

The radius and center coordinates are:

R = 15.899002930062595
P = [13.4207317073, -9.56097560967, 0]

Any ideas on how to find the center coordinates?

Answer 1

There are two issues with your code.

The first is in the naming convention. For all the formulas you are using to hold, the side of length a has to be the one opposite the point A, and similarly for b and B and c and C. You can solve that by computing them as:

a = np.linalg.norm(C - B)
b = np.linalg.norm(C - A)
c = np.linalg.norm(B - A)

The second has to do with the note in your source for the barycentric coordinates of the circumcenter: not necessarily homogeneous. That is, they need not be normalized in any way, and the formula you are using to compute the Cartesian coordinates from the barycentric ones is only valid when they add up to one.

Fortunately, you only need to divide the resulting Cartesian coordinates by b1 + b2 + b3 to get the result you are after. Streamlining a little bit your code for efficiency, I get the results you expect:

>>> A = np.array([2.0, 1.5, 0.0])
>>> B = np.array([6.0, 4.5, 0.0])
>>> C = np.array([11.75, 6.25, 0.0])
>>> a = np.linalg.norm(C - B)
>>> b = np.linalg.norm(C - A)
>>> c = np.linalg.norm(B - A)
>>> s = (a + b + c) / 2
>>> R = a*b*c / 4 / np.sqrt(s * (s - a) * (s - b) * (s - c))
>>> b1 = a*a * (b*b + c*c - a*a)
>>> b2 = b*b * (a*a + c*c - b*b)
>>> b3 = c*c * (a*a + b*b - c*c)
>>> P = np.column_stack((A, B, C)).dot(np.hstack((b1, b2, b3)))
>>> P /= b1 + b2 + b3
>>> R
15.899002930062531
>>> P
array([ 13.42073171,  -9.56097561,   0.        ])

Answer 2

As an extension to the original problem: Now suppose we have an arc of known length (e.g. 1 unit) extending from point A as defined above (away from point B). How to find the 3D coordinates of the end point (N)? The new point N lies on the same circle passing through points A, B & C.

This can be solved by first finding the angle between the two vectors PA & PN:

# L = Arc Length 
theta = L / R

Now all we need to do is rotate the vector PA (Radius) by this angle theta in the correct direction. In order to do that, we need the 3D rotation matrix. For that we use the Euler–Rodrigues formula:

def rotation_matrix_3d(axis, theta):
    axis = axis / np.linalg.norm(axis)
    a = np.cos(theta / 2.0)
    b, c, d = axis * np.sin(theta / 2.0)   
    rot = np.array([[a*a+b*b-c*c-d*d,     2*(b*c-a*d),     2*(b*d+a*c)],
                    [    2*(b*c+a*d), a*a+c*c-b*b-d*d,     2*(c*d-a*b)],
                    [    2*(b*d-a*c),     2*(c*d+a*b), a*a+d*d-b*b-c*c]])
    return rot

Next, we need to find the rotation axis to rotate PA around. This can be achieved by finding the axis normal to the plane of the circle passing through A, B & C. The coordinates of point N can then be found from the center of the circle.

PA = P - A
PB = P - B
xx = np.cross(PB, PA)
r3d = rotation_matrix_3d(xx, theta)
PN = np.dot(r3d, PA)
N = P - PN

as an example, we want to find coordinates of point N that are 3 degrees away from point A, the coordinates would be

N = (1.43676498,  0.8871264,   0.) 

Which is correct! (manually verified using CAD program)