Find the arc / circle equation given three points in space (3D)

Question

Given 3 points in (3D) space: 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 the equation of a circle? Using Python and Numpy is my source 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 definition of Circumradius, 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 triangular points to search for the trilinear coordinates of a circle ( Circumcenter ).

First ( using this source ) we find the dashed coordinates of the circle:

 #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) will 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 above barcode values ​​do not seem correct. When resolving with an example of known values, the radius is correct, but there is no center coordinate.

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

Radius and coordinates of the center:

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

Any ideas on how to find the center coordinates?

Answer 1

There are two problems in the code.

The first is a naming convention. For all formulas that you use to hold, the side of length a must be the one opposite the point a , and similarly for b and b and c and c . You can solve this by calculating them as:

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

The second is related to a note in your source for the barycentric coordinates of the circle: not necessarily uniform . That is, they should not be normalized in any way, and the formula that you use to calculate the Cartesian coordinates from the barycentric ones is valid only 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. By streamlining a bit of your code to increase efficiency, I get the expected results:

 >>> 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 of the original problem: now suppose that we have an arc of known length (for example, 1 unit) extending from point A as defined above (away from point B). How to find the 3D coordinates of the endpoint (N)? The new point N lies on the same circle passing through points A, B and C.

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

 # L = Arc Length theta = L / R 

Now we only need to rotate the vector PA (radius) by this angle theta in the right direction. For this we need a 3D rotation matrix. To do this, 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*bc*cd*d, 2*(b*ca*d), 2*(b*d+a*c)], [ 2*(b*c+a*d), a*a+c*cb*bd*d, 2*(c*da*b)], [ 2*(b*da*c), 2*(c*d+a*b), a*a+d*db*bc*c]]) return rot 

Then 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 and 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 the coordinates of point N, which are 3 degrees from point A, the coordinates will be

 N = (1.43676498, 0.8871264, 0.) 

It is right! (checked manually using CAD software)