How to find the center of a circle and its radius in a 3D space given 3 points

Question

I'm working on writing a code in Octave (C++) for a helical spring. I need to figure out the center line of this spring in an effort to find any trends between different spring platforms as the geometry changes throughout the manufacturing processes.

I have a file with the center line of the wire, which is made up of 2000 or so points. Since I learned that you need 3 points to make a circle, I'm figuring i will be able to approach it that way (but haven't figured out how) so....

Given 3 points in 3 Dimensional space, how could I find the center and the radius of a circle in a code friendly manor.

Thanks

---

Answer 1

Find (x,y,z) such that the distance to your three points are equal.

$(x-x_1)^2+ (y-y_1)^2 + (z - z_1)^2 = \\(x-x_2)^2+ (y-y_2)^2 + (z - z_2)^2 = \\(x-x_3)^2+ (y-y_3)^2 + (z - z_3)^2$

Multiply that out and we can turn that into a system of linear equations that should be easy enough to solve.

$2x_1 x + 2 y_1 y + 2 z_1 z = x_1^2 + y_1^2 + z_1^2\\ 2x_2 x + 2 y_2 y + 2 z_2 z = x_2^2 + y_2^2 + z_2^2\\ 2x_3 x + 2 y_3 y + 2 z_3 z = x_3^2 + y_3^2 + z_3^2$

Update

That system of equation determine a single point, that if you traveled nomal to the plane you would intersect at the circumcenter.

Lets call this point, $(x,y,z)$

I was trying to avoid finding the normal line, but here is the normal vector.

$N=(x_2-x_1,y_2-y_1, z_2-z_1)\times(x_3-x_1,y_3-y_1, z_3-z_1)\\ N_x = y_1z_2 - y_1z_3 + y_2z_3 - y_2 z_1 + y_3z_1-y_3z_2\\ N_y = x_1z_3 - x_1z_2 + x_2z_1 - x_2 z_3 + x_3z_2-x_3z_1\\ N_z = x_1y_2 - x_1y_3 + x_2y_3 - x_2 y_1 + x_3y_1-x_3y_2\\ $

Circumcenter = $(x,y,z) - \dfrac {N \cdot (x,y,z) - N\cdot (x_1, y_1,z_1)}{N_x^2+N_y^2+N_z^2} N$

Answer 2

Let $$M=(M_x, M_y, M_z),\quad N=(N_x, N_y, N_z),\quad P=(P_x, P_y, P_z)$$ are the given points of the circle,
$$C=(x,y,z)$$ is the center, then $$(x-M_x)^2+(y-M_y)^2+(z-M_z)^2 = (x-N_x)^2+(y-N_y)^2+(z-N_z)^2,$$ $$(x-M_x)^2+(y-M_y)^2+(z-M_z)^2 = (x-P_x)^2+(y-P_y)^2+(z-P_z)^2,$$ $$\begin{vmatrix} x-M_x &y-M_y&z-M_z\\ N_x-M_x &N_y-M_y&N_z-M_z\\ P_x-M_x &P_y-M_y&P_z-M_z\\ \end{vmatrix}=0$$ (coplanarity condition) gives the linear system in $x,y,z$, and $$R^2=LHS(1,2)=RHS(1,2).$$

Details

Let us simplify the system. $$\begin{cases} 2(N_x-M_x)x+2(N_y-M_y)y+2(N_z-M_z)z = N_x^2+N_y^2+N_z^2-M_x^2-M_y^2-M_z^2\\ 2(P_x-M_x)x+2(P_y-M_y)y+2(P_z-M_z)z = P_x^2+P_y^2+P_z^2-M_x^2-M_y^2-M_z^2\\ A_xx+A_yy+A_zz=A_xM_x+A_yM_y+A_zM_z, \end{cases}$$ where $$A_x=\begin{vmatrix}N_y-M_y&N_z-M_z\\P_y-M_y&P_z-M_z\end{vmatrix}=(N_y-M_y)(P_z-M_z)-(N_z-M_z)(P_y-M_y),$$ $$A_y=\begin{vmatrix}N_z-M_z&N_x-M_x\\P_z-M_z&P_x-M_x\end{vmatrix}=(N_z-M_z)(P_x-M_x)-(N_x-M_x)(P_z-M_z),$$ $$A_z=\begin{vmatrix}N_x-M_x&N_y-M_y\\P_x-M_x&P_y-M_y\end{vmatrix}=(N_x-M_x)(P_y-M_y)-(N_y-M_y)(P_x-M_x).$$

Solve the resulting system can be by Cramer's rule, Gaussian elimination or using Wolfram Alpha.

Answer 3

Let $A = (A_1, A_2, A_3)^T, \,\, B = (B_1, B_2, B_3)^T, \,\, C = (C_1, C_2, C_3)^T$ be the three points and $X = (x,y,z)^T$ be the center we are looking for. The point $X$ lies on the two planes passing through the midpoints of the segments $AB$ and $AC$. In fact $X$ also lies on the plane passing through the midpoint of segment $BC$, the three planes intersecting at the common line passing trough the center of the circle through $A, B, C$ and orthogonal to the plane defined by the three points. Also, $X$ lies on the plane through $A, B, C$. Thus, $X$is the unique intersection point of three planes, which constitute a system of three linear equations with three unknown variables $X = (x, y, z)^T$. Let us write the equations of the three planes: \begin{align} \overrightarrow{AB} \cdot \left(X - \frac{1}{2}(B+A)\right) &= 0\\ \overrightarrow{AC} \cdot \left(X - \frac{1}{2}(C+A)\right) &= 0\\ \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big) \cdot \left(X - A\right) &= 0 \end{align}
which turn into
\begin{align} (B - A)^T \left(X - \frac{1}{2}(B+A)\right) &= 0\\ (C - A)^T \left(X - \frac{1}{2}(C+A)\right) &= 0\\ \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big) \cdot \left(X - A\right) &= 0 \end{align}

and as a linear system it looks like this: \begin{align} 2(B - A)^T X &= (B - A)^T (B+A) \\ 2 (C - A)^T X &= (C - A)^T (C+A) \\ \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)^T X &= \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)^T A \end{align}
So if we put the vectors multiplying $X$ in as rows in the $3 \times 3$ matrix \begin{align} M= \left[ {\begin{array}{c} 2 (B - A)^T \\ 2 (C - A)^T \\ \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)^T \end{array} } \right] = \left[ {\begin{array}{c c c} 2 B_1 - 2 A_1 & 2 B_2 - 2 A_2 & 2 B_3 - 2 A_3\\ 2 C_1 - 2 A_1 & 2 C_2 - 2 A_2 & 2 C_3 - 2 A_3 \\ \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)_1 & \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)_2 & \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)_3 \end{array} } \right] \end{align} and the numbers in the right hand side of the equations into a $3 \times 1$ row vector \begin{align} K&= \left[ {\begin{array}{c} (B - A)^T (B+A) \\ (C - A)^T (C+A) \\ \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)^T A \end{array} } \right] = \left[ {\begin{array}{c} B^TB - A^TA\\ C^TC - A^TA \\ \det \Big|\overrightarrow{AB} \,\, \overrightarrow{AC} \,\, A \Big| \end{array} } \right] \\ &= \left[ {\begin{array}{c} B_1^2 + B_2^2 + B_3^2 - A_1^2 - A_2^2 - A_3^2\\ C_1^2 + C_2^2 + C_3^2 - A_1^2 - A_2^2 - A_3^2\\ \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)_1 A_1 + \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)_2 A_2 + \Big(\overrightarrow{AB} \times \overrightarrow{AC}\Big)_3 A_3 \end{array} } \right] \end{align} Then the center of the circle is $$X = M^{-1}K$$ and the radius is $$R = \sqrt{(X - A)^T(X-A)} = \sqrt{(x-A_1)^2 + (y-A_2)^2 + (z-A_3)^2}$$

Answer 4

Let the points be $a,b,c$. Form the unit vector $\vec u=\vec{ab}/\|\vec{ab}\|$. Then $\vec w=\vec u\times\vec {ac}/\|\vec u\times\vec {ac}\|$ and $v=\vec w\times\vec u$. Now $\vec u,\vec v,\vec w$ form a triorthogonal frame. We can take $a$ for its origin.

The coordinates of the three points in this frame are

$$a:(0,0),\\b:(\vec{ab}\cdot\vec u,0)=(u_b,0),\\c:(\vec{ac}\cdot\vec u,\vec{ac}\cdot\vec v)=(u_c,v_c).$$

The center is on the the bissectrix of $ab$, the vertical at $u=u_b/2$, and on the bissectrix of $ac$, of equation $u_cu+v_cv=(u_c^2+v_c^2)/2$. Compute $v$.

Now you have the coordinates of the center in the reduced frame and the radius follows. Converting to the original coordinates is immediate.