Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

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?

like image 882
Nader Avatar asked Dec 01 '13 16:12

Nader


People also ask

How do you find the center of a circle with 3 points in 3D?

1) Find a plane from the 3 points and create a 2D coordinate system on that plane. 2) Convert all 3 points to that 2D coordinate system. 3) Find the circle center using the link above. 4) Convert the circle center (in 2D) back to 3D.

Can you find the 3D circle centre with 3 points?

IDS's earlier post says "will find the 3D circle centre given either any 3 points on the circle, or two tangent points and the intersection of the tangents." rb1957 - no, 3 points define a triangle on a plane (unless they are co-linear), and you can only draw one circle through the vertices of a triangle.

What is the equation for a circle in R3?

A circle in R 3 doesn't exactly have "an equation" in the sense that a circle in the plane or a sphere in space has an equation. Intuitively, "an equation" cuts down the dimension by one, but to get a circle in space you have to lower the dimension by two.

How do you find the equation of a circle?

The task is to find the equation of the circle and then print the center and the radius of the circle. Equation of circle in general form is x² + y² + 2gx + 2fy + c = 0 and in radius form is (x – h)² + (y -k)² = r², where (h, k) is the center of the circle and r is the radius. Examples: Output: The equation of the circle is x 2 + y 2 = 1.

How to find the centre and radius of a circle?

The task is to find the equation of the circle and then print the centre and the radius of the circle. Equation of circle in general form is x² + y² + 2gx + 2fy + c = 0 and in radius form is (x – h)² + (y -k)² = r², where (h, k) is the centre of the circle and r is the radius. The equation of the circle is x 2 + y 2 = 1.


2 Answers

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.        ])
like image 130
Jaime Avatar answered Sep 25 '22 12:09

Jaime


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)

like image 36
Nader Avatar answered Sep 26 '22 12:09

Nader