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

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

• Are all of the 2000 data points supposed to belong (with some error) to the same circle? If so, using just 3 of those points will probably not do a robust job of estimating the true circle. (What you may want instead is a method which uses all 2000 points to get an approximation.) – Semiclassical Aug 29 '16 at 20:21
• The points that form a helical spring are on the surface of a cylinder. You are interested in the axis (and possibly radius) of that cylinder. So, use cylinder fitting to fit your $N \approx 2000$ samples to a cylinder; the diameter of the cylinder and the axis of the cylinder will tell you the diameter and axis of your helical spring, too. (For conical springs, fit to a cone, and so on.) – Nominal Animal Aug 29 '16 at 22:24
• These 2000 points will basically have an ever-changing radius. All of the center points of these radii will make up the center line of the overall spring. – Help Me Please Aug 30 '16 at 11:40
• @HelpMePlease: for any three points on the helix, the center of the circle isn't necessarily on the axis, your approach is flawed. In particular, the center of curvature (circle by three very close points) is on another helix. – Yves Daoust Aug 31 '16 at 13:19
• You can make a circle out of three points in a plane going through those points. But for a helical spring no two points of the spring lie on the plane normal to the spring's axis, so (as @YvesDaoust said) no three points determine a circle, whose center coincides with the spring's axis. Also, probably no such circle has a radius equal to the radius of the spring. You need to find a cylinder of your helix first, then retrieve its axis and radius. – CiaPan Aug 31 '16 at 13:41

## 4 Answers

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.

• This looks promising, but I'm still not sure how to solve for x, y, and z given what you have (I'm only a college sophomore and haven't learned any of this stuff yet.) Would you be able to break it down even further to help give me a better understanding. Thanks!!! – Help Me Please Aug 31 '16 at 12:23
• @HelpMePlease Ready – Yuri Negometyanov Aug 31 '16 at 13:15

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}$$

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.

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$

• it's not enough, this gives a line solution. The system needs to add the intersection of this line with the plane of the 3 points – user354674 Aug 29 '16 at 19:57
• The three equations you find can be interpreted as the equations of the three median planes (for example at equal distance from $P_1$ and $P_2$). The third equation is thus redundant, and should be replaced, as @igael says, by the equation of plane $P_1P_2P_3$. – Jean Marie Aug 29 '16 at 22:15
• That system of equation determine a single point, that if you traveled nomal to the plane you would intersect at the circumcenter. – Doug M Aug 30 '16 at 18:46
• I'm trying to workout your updated equations and noticed a couple of what I believe to be typos. In N = ... you have X3-Xz. I'm figuring you mean X3-X1. However, I can't figure out what you wanted to say in Nx. The first two numbers are both Y1*Z2. These will cancel each other out which is why i think it's a typo. Please fix or tell me I'm wrong – Help Me Please Aug 31 '16 at 12:04
• @HelpMePlease I have updated $N$ and $N_x$ – Doug M Aug 31 '16 at 16:11