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

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>>> R15.899002930062531>>> Parray([ 13.42073171,  -9.56097561,   0.        ])

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 - APB = P - Bxx = 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)