Circle centre calculation passing through 3 points

For a case study, I had to calculate in 3D the centre of a circle passing through 3 points.

 

1- The trick

We can note that if they aren’t aligned, the 3 points (M1,M2,M3) form a triangle (see figure): then the centre is the intersecting of the perpendicular bisectors (2D) or perpendicular bisector planes (3D).

 

cercle_3pts 

(Perpendicular bisector intersection – 2D)

 

2- A bit of analytic geometry

Let us working directly in 3D. O is the centre of the circle so that its coordinates are (XO, YO, ZO). The mathematical system has 3 unknowns and it needs 3 equations that are nothing more than the equations of the 3 planes intersecting in the point O:

  1. Plane number 1 : plane passing by (M1,M2,M3)
  2. Plane number 2 : perpendicular bisector plane to [M1M2] segment
  3. Plane number 3 : perpendicular bisector plane to [M2M3] segment

 

 

Let I_1 and I_2 the middle points so that

I_1 = \frac{1}{2} \left( \begin{array}{c} x_1 + x_2 \\ y_1 + y_2 \\ z_1 + z_2 \end{array} \right) ;  I_2 = \frac{1}{2} \left( \begin{array}{c} x_3 + x_2 \\ y_3 + y_2 \\ z_3 + z_2\end{array} \right).

 

Rule : Point P(x,y,z) belongs to a perpendicular  bisector plane if \left( \overrightarrow{I_i P_i} \right) \cdot \left( \overrightarrow{M_i M_{i+1}} \right) = 0

 

a) First perpendicular bisector plane to [M1M2] segment

\frac{1}{2} \left[ \begin{array}{c}2x - \left( x_1 + x_2 \right) \\2y - \left( y_1 + y_2 \right) \\2z - \left( z_1 + z_2 \right)\end{array} \right] \cdot\left[\begin{array}{c}\left( x_2 - x_1 \right) \\\left( y_2 - y_1 \right) \\\left( z_2 - z_1 \right)\end{array} \right] = 0

 

By developing the scalar product and by equating:
\left\lbrace\begin{array}{l}A_1 = 2 \left( x_2 - x_1 \right) \\B_1 = 2 \left( y_2 - y_1 \right) \\C_1 = 2 \left( z_2 - z_1 \right) \\D_1 = \left[ x_2^2 + y_2^2 + z_2^2 - x_1^2 - y_1^2 - z_1^2 \right]\end{array} \right.,

 

We obtain:

(1)   \begin{equation*}A_1 x + B_1 y + C_1 z = D_1\end{equation*}

 

 

b) Second perpendicular bisector plane to [M2M3] segment

It is obtained by circular permutation of the previous equation:

(2)   \begin{equation*}A_2 x + B_2 y + C_2 z = D_2\end{equation*}

with
\left\lbrace\begin{array}{l}A_2 = 2 \left( x_3 - x_2 \right) \\B_2 = 2 \left( y_3 - y_2 \right) \\C_2 = 2 \left( z_3 - z_2 \right) \\D_2 = \left[ x_3^2 + y_3^2 + z_3^2 - x_2^2 - y_2^2 - z_2^2 \right]\end{array} \right.

 

 

c) Third plane passing through the points M_1, M_2 and M_3

Using the vector product of:
\overrightarrow{M_1 M_2} = \left( \begin{array}{c} x_2 - x_1 \\ y_2 - y_1 \\ z_2 - z_1 \end{array} \right)  et \overrightarrow{M_1 M_3} = \left( \begin{array}{c} x_3 - x_1 \\ y_3 - y_1 \\ z_3 - z_1 \end{array} \right)

 

We finally have :

(3)   \begin{equation*}A_3 x + B_3 y + C_3 z = D_3\end{equation*}

with
\left\lbrace \begin{array}{l}A_3 = \left( y_2 - y_1 \right) \left( z_3 - z_1 \right) - \left( z_2 - z_1 \right) \left( y_3 - y_1 \right) \\B_3 = \left( z_2 - z_1 \right) \left( x_3 - x_1 \right) - \left( x_2 - x_1 \right) \left( z_3 - z_1 \right) \\C_3 = \left( x_2 - x_1 \right) \left( y_3 - y_1 \right) - \left( y_2 - y_1 \right) \left( x_3 - x_1 \right) \\\end{array} \right.

 

We know that this plane passes through the 3 points (M_1 has been chosen here), then D_3 = A_3 x_1 + B_3 y_1 + C_3 z_1.

 

Under the matrix form, the system becomes:

(4)   \begin{equation*}\left[\begin{array}{ccc}A_1 & B_1 & C_1 \\A_2 & B_2 & C_2 \\A_3 & B_3 & C_3\end{array} \right]\cdot\left[\begin{array}{c}x \\y \\z\end{array} \right]=\left[\begin{array}{c}D_1 \\D_2 \\D_3\end{array} \right]\end{equation*}

 

i.e. the O point coordinates are the results of the matrix product:

(5)   \begin{equation*}\left[\begin{array}{c}x \\y \\z\end{array} \right]=\left[\begin{array}{ccc}A_1 & B_1 & C_1 \\A_2 & B_2 & C_2 \\A_3 & B_3 & C_3\end{array} \right]^{-1}\cdot\left[\begin{array}{c}D_1 \\D_2 \\D_3\end{array} \right]\end{equation*}

 

NB : note that if the 3 points are aligned, then the determinant is equal to zero ; here: 
det \left[\begin{array}{ccc}A_1 & B_1 & C_1 \\A_2 & B_2 & C_2 \\A_3 & B_3 & C_3\end{array} \right]\ne 0

 

A Scilab code was developed to perform this calculation (here circle_3points_scilab); it can be seen hereafter (a basic scalar product verifies first that the points are not aligned – a tolerance of 1% has been introduced in order to limit errors of rounding effects):

mode(0)

function G = centre_calculus(M1,M2,M3)
    A1 = 2*(M2(1) - M1(1));
    A2 = 2*(M3(1) - M2(1));
    A3 = (M2(2)-M1(2))*(M3(3)-M1(3))-(M2(3)-M1(3))*(M3(2)-M1(2));

    B1 = 2*(M2(2) - M1(2));
    B2 = 2*(M3(2) - M2(2));
    B3 = (M2(3)-M1(3))*(M3(1)-M1(1))-(M2(1)-M1(1))*(M3(3)-M1(3));

    C1 = 2*(M2(3) - M1(3));
    C2 = 2*(M3(3) - M2(3));
    C3 = (M2(1)-M1(1))*(M3(2)-M1(2))-(M2(2)-M1(2))*(M3(1)-M1(1));

    D1 = M2(1)^2 + M2(2)^2 + M2(3)^2 - M1(1)^2 - M1(2)^2 - M1(3)^2;
    D2 = M3(1)^2 + M3(2)^2 + M3(3)^2 - M2(1)^2 - M2(2)^2 - M2(3)^2;
    D3 = A3*M1(1)+B3*M1(2)+C3*M1(3);

    M = [A1 B1 C1 ; A2 B2 C2 ; A3 B3 C3]
    det_ = det(M)
    G = inv(M) * [D1 ; D2 ; D3]
endfunction

tolerance = 0.01  // 1% 

function straight = scalar_product(A,B,C,tol_)
    scal_status = %f;
    AB = [B(1) - A(1) ; B(2) - A(2) ; B(3) - A(3)];
    BC = [C(1) - A(1) ; C(2) - A(2) ; C(3) - A(3)];
    norm_AB = sqrt(AB(1)^2 + AB(2)^2 + AB(3)^2);
    norm_BC = sqrt(BC(1)^2 + BC(2)^2 + BC(3)^2);
    
    scal = (AB(1)*BC(1) + AB(2)*BC(2) + AB(3)*BC(3))/(norm_AB * norm_BC);
    
    if ( ((scal >= (-1 - tolerance)) & (scal <= (-1 + tolerance))) | ((scal >= (1 - tolerance)) & (scal <= (1 + tolerance)))) then
        straight = %t;
    else
        straight = %f;
    end
endfunction

// Sphere #1
Case1_reference = [10;25;-5]

M1 = [11.0453 ; 15.0548 ; -5]
M2 = [6.90983 ; 34.5106 ; -5]
M3 = [19.1355 ; 20.9326 ; -5]

if (scalar_product(M1,M2,M3,tolerance) == %f) then
    Centre_Case1 = centre_calculus(M1,M2,M3)
else
    printf("The 3 points are on a straight line !!!\n") ; 
end

diff1_ = Centre_Case1 - Case1_reference

// Sphere #2
Case2_reference = [-10;-25;5]
N1 = [-18.5845 ; -20.9326 ; 8.124510000000001]
N2 = [-16.2878 ; -32.4314 ; 7.28856]
N3 = [-5.30154 ; -33.6603 ; 3.2899]
if (scalar_product(M1,M2,M3,tolerance) == %f) then
    Centre_Case2 = centre_calculus(N1,N2,N3)
else
    printf("The 3 points are on a straight line !!!\n")
end

diff2_ = Centre_Case2 - Case2_reference


// Sphere #3
Case3_reference = [50;-100;-50]
P1 = [69.563 ; -103.364 ; -52.4441]
P2 = [40 ; -85.98739999999999 ; -39.8193]
P3 = [66.1803 ; -90.4894 ; -43.0902]
if (scalar_product(M1,M2,M3,tolerance) == %f) then
    Centre_Case3 = centre_calculus(P1,P2,P3)
else
    printf("The 3 points are on a straight line !!!\n")
end

diff3_ = Centre_Case3 - Case3_reference

 

The previous coordinates come from a basic test where the spheres are generated using Openscad (here sphere_openscad), exported in STL format and read with Gmsh in order to get the coordinates of 3 random nodes in the surface (the cut plane necessarely passes by the centre of the sphere):

sphere

// case #1
radius_=10;
diameter_ = 2*radius_;
translate([10,25,-5])
{
difference()
	{
	sphere(r=radius_, centre="true");
	translate([-radius_,-radius_,0]) cube([diameter_,diameter_,diameter_], center="true");
	}
}


// case #2
translate([-10,-25,5])
{
	rotate([0,20,0])
	{
	difference()
		{
		sphere(r=radius_, centre="true");
		translate([-radius_,-radius_,0]) cube([diameter_,diameter_,diameter_], center="true");
		}
	}
}


// case #3
translate([50,-100,-50])
{
	rotate([36,0,0])
	{
	difference()
		{
		sphere(r=2*radius_, centre="true");
		translate([-2*radius_,-2*radius_,0]) cube([2*diameter_,2*diameter_,2*diameter_], center="true");
		}
	}
}

The gaps come from the mesher tolerances and to the rounding errors as well.