Centre d’un cercle passant par 3 points

Dans le cadre d’un développement, j’ai eu besoin de calculer en 3D le centre d’un cercle passant par 3 points.

 

1- L’astuce

Nous remarquerons que s’ils ne sont pas alignés, les 3 points (M1,M2,M3) forment un triangle (voir figure) : le centre n’est autre que le point d’intersection des médiatrices (en 2D) ou des plans médiateurs (en 3D).

 

cercle_3pts 

(intersection des médiatrices – 2D)

 

2- Un peu de géométrie analytique

Nous allons directement travailler en 3D. Si O est le centre du cercle, alors (XO, YO, ZO) sont ses coordonnées. Ce système à 3 inconnues requière 3 équations qui ne sont ni plus ni moins que les équations des 3 plans qui se croisent en O :

  1. Plan n°1 : plan passant par (M1,M2,M3)
  2. Plan n°2 : plan médiateur au segment [M1M2]
  3. Plan n°3 : plan médiateur au segment [M2M3]

 

 

Soit I_1 et I_2 les points milieux tels que

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).

 

Règle : Le point P(x,y,z) appartient au plan médiateur ssi \left( \overrightarrow{I_i P_i} \right) \cdot \left( \overrightarrow{M_i M_{i+1}} \right) = 0

 

a) Premier plan médiateur au segment [M1M2]

\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

 

En développant le produit scalaire et en posant :
\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.,

 

nous obtenons :

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

 

 

b) Second plan médiateur au segment [M2M3]

Il est obtenu par permutation circulaire de l’équation précédente :

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

avec
\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) troisième plan passant par les points M_1, M_2 et M_3

En passant par le produit vectoriel de :
\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)

 

Nous aboutissons à

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

avec
\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.

 

Sachant que le plan passe par les 3 points (M_1 est choisi ici), alors D_3 = A_3 x_1 + B_3 y_1 + C_3 z_1.

 

Sous forme matricielle, le système devient :

(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*}

 

c’est-à-dire que les coordonnées de O sont le résultat du produit matriciel :

(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 : A noter que si les 3 points ne sont pas alignés, alors
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

 

Les programmes Scilab réalisant ce calcul (ici circle_3points_scilab) est visible ci-dessous (un simple produit scalaire vérifie initialement que les points ne sont pas alignés – une tolérance de 1% est ajoutée pour limiter les effets des erreurs d’arrondi) :

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

 

Les coordonnées ci-dessous proviennent d’un simple test où des sphères sont générées sous Openscad (ici sphere_openscad), exportées au format STL et lu avec Gmsh pour récupérer les coordonnées de 3 points quelconques en surface (le plan de coupe passe obligatoirement par le centre de la sphère) : 

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");
		}
	}
}

Les écarts sont dus aux tolérances