Discussion :

[KeV31]

Bonjour,

Je cherche l'équation d'un plan moyen défini par N points.

Pour ce faire, la méthode la plus appropriée semble être la minimisation de la somme des carrés des distances des points à ce plan.

Ce lien propose la solution suivante :

Soit z=a.x+b.y+c l'équation du plan .

la solution est donné par la résolution du système :
(Sxx).a+(Sxy).b+(Sx).c = (Sxz)
(Sxy).a+(Syy).b+(Sy).c = (Syz)
(Sx).a+(Sy).b+n.c = (Sz)

Les coefficients étant calculés préalablement par :
Sxx = Somme des (xk)² pour k=1 à k=n
Sxy = Somme des (xk)(yk)
Sx = Somme des (xk)
Sxz = Somme de (xk)(zk)
Syy = Somme de (yk)²
Syz = Somme de (yk)(zk)
Szz = Somme de (zk)²
Sz = Somme de (zk)

Malheureusement, cet algorithme ne semble pas donner des résultats justes (j'ai testé en le codant en java). Est-ce quelqu’un pourrait certifier qu'il est correct, y apporter des modifications ou me proposer une méthode différente.
Merci.

[souviron34]

En général il y a 2 types de solutions pour ce problème :

• tu peux calculer les coefficients pour chacun des points :
  moindre carrés généralisés à N (nombre de coefficients) paramètres
• tu ne connais pas la forme des coefficients pour chacun des points, mais la forme théorique de l'équation :
  algo style simplex (résolution d'un système d'équations linéaires avec contrainte)

[pseudocode]

1. Calculer

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
a00=Somme[(xi-xc)^2]
    i=1..n
 
a11=Somme[(yi-yc)^2]
    i=1..n
 
a22=Somme[(zi-zc)^2]
    i=1..n
 
a01=a10=Somme[(xi-xc)(yi-yc)]
        i=1..n
 
a02=a20=Somme[(xi-xc)(zi-zc)]
        i=1..n
 
a12=a21=Somme[(yi-yc)(zi-zc)]
        i=1..n

avec xc, yc et zc les isobarycentres des xi, yi et zi


2. Rechercher les valeurs et vecteurs propres de la matrice

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
| a00 a01 a02 |
| a10 a11 a12 |
| a20 a21 a22 |

le vecteur propre (A,B,C) associée a la plus petite valeur propre est le vecteur normal au plan.

Donc le plan est d'equation Ax+By+Cz+D = 0

3. Résoudre D avec les isobarycentres xc, yc et zc.

[Zavonen]

L'équation générale d'un plan dans l'espace est
ux+vy+wz+t=0
A priori 4 inconnues, mais on peut homogéneiser, supposer par exemple que u=1.
Ce faisant on élimine seulement les plans paralléles à un plan de coordonnées.
Donc déterminer un plan revient à déterminer 3 inconnues.
P(v,w,t) plan d'équation: x+vy+wz+t=0
Pour chaque point d'une famille indexée Ai 1<=i<=n, on peut déterminer l'intersection Bi de P avec la normale à P passant par Ai, les coordonnées sont des fonctions linéaires de v,w,t.
Puis on peut calculer la distance AiBi au carré (fonction polynomiale de degré 2 en les 3 variables v,w,t).
Le résultat est (de mémoire,je crois) (xi+vyi+wzi+t)^2
On peut enfin calculer la somme S2(v,w,t) de ces carrés qui reste du même type (polynomiale degré 2).
Le problème est donc de chercher le minimum d'une fonction de 3 variables, on annule les dérivées partielles.
On tombe donc sur un système de 3 équations à 3 inconnues, mais il est malheureusement de degré 2.
On peut utiliser une équation pour calculer une inconnue en fonction des deux autres, la dérivée partielle en v est, par exemple, du premier degré.
v=f(w,t)
On tombe alors sur un système de deux équations à deux inconnues de degré deux, on a alors à calculer l'intersection de deux coniques. On a (le plus souvent) deux points
(w1,t1) (w2,t2)
On a alors les possibilités:

w1 ,t1, f(w1,t1)
w2,t2, f(w2,t2)

On calcule S2 avec ces 2 possibilités et on retient ce qui donne le minimum effectif.

[KeV31]

L'algorithme que je présentais dans mon premier post était bon, les erreurs que j'obtenais sur le résultat provenaient d'une erreur dans mon code.

Je laisse ci-dessous mon code en java qui fonctionne correctement, au cas ou ça puisse aider qelqu'un.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
double sommex = 0;
        double sommey = 0;
        double sommez = 0;
        double sommexy = 0;
        double sommexz = 0;
        double sommeyz = 0;
        double sommexx = 0;
        double sommeyy = 0;
        double sommezz = 0;
       //tabPoints est un tableau qui contient les N points
        for (int k=0; k<nbPoints ; k++){
            // Sommexx = Somme des (xk)² pour k=1 à k=n 
            sommexx = sommexx + (tabPoints[k].x)*(tabPoints[k].x);
            //Sxy = Somme des (xk) (yk)
            sommexy = sommexy + ( tabPoints[k].x)*(tabPoints[k].y);
            //Sx = Somme des (xk)
            sommex = sommex + (tabPoints[k].x);
            //Sxz = Somme de (xk)(zk)
            sommexz = sommexz + (tabPoints[k].x)*(tabPoints[k].z);
            //Syy = Somme de (yk)²
            sommeyy = sommeyy + (tabPoints[k].y)*(tabPoints[k].y);
            //Syz = Somme de (yk)(zk)
            sommeyz = sommeyz + (tabPoints[k].y)*(tabPoints[k].z); 
            //Szz = Somme de (zk)²
            sommezz = sommezz + (tabPoints[k].z)*(tabPoints[k].z); 
            //Sy = Somme de (yk)
            sommey = sommey + (tabPoints[k].y); 
            //Sz = Somme de (zk) 
            sommez = sommez + (tabPoints[k].z);
        }
 
 
        //Matrice (M)= Sxx, Sxy, Sx
        //             Sxy, Syy, Sy
        //             Sx, Sy, n
        Matrix3d matrice = new Matrix3d();
        matrice.m00 = sommexx;
        matrice.m01 = sommexy;
        matrice.m02 = sommex;
        matrice.m10 = sommexy;
        matrice.m11 = sommeyy;
        matrice.m12 = sommey;
        matrice.m20 = sommex;
        matrice.m21 = sommey;
        matrice.m22 = nbPoints;
 
 
        //Vecteur (V)=Sxz,Syz,Sz
        Vector3d vecteurV =  new Vector3d();
        vecteurV.x = sommexz;
        vecteurV.y = sommeyz;
        vecteurV.z = sommez;
 
        //VecteurS est le vecteur solution du système
        //On a (Matrice)*(vecteurS)=(VecteurV)
        //Donc (vecteurS) = ((Matrice)^(-1))*(VecteurV)
        Vector3d vecteurS =  new Vector3d();
 
        matrice.invert();
        matrice.transform(vecteurV);
        vecteurS = vecteurV;

[bernard6]

Bonjour,

J'ai implémenté cela sous Matlab mais je ne suis pas sur que c'est correct..

Pourriez vous me donner par exemple les coordonnées de 4 points et le résultat (a, b, c) après ajustement d'un plan dessus histoire de voir si j'ai juste...

Voici mon script :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
  % Construction de la matrice M
sxx = sum (power(PointsSegment(:,1),2))
sxy = sum (PointsSegment(:,1).*PointsSegment(:,2))
sx = sum (PointsSegment(:,1))
syy = sum (power(PointsSegment(:,2),2))
sy = sum (PointsSegment(:,2))
n = PO
MatriceM = [sxx sxy sx; sxy syy sy; sx sy n]
 
% Construction de la matrice V
sxz = sum (PointsSegment(:,1).*PointsSegment(:,3))
syz = sum (PointsSegment(:,2).*PointsSegment(:,3))
sz = sum (PointsSegment(:,3))
MatriceV = [sxz; syz; sz]
 
% Recherche de la matrice X (a b c)
MatriceX = inv(MatriceM) * MatriceV

Merci d'avance !

[Krystof13]

Bonjour,
J'ai le même problème mais en plus je voudrait connaître la dispersion de mes points autour du plan calculé, je crois qu'on appelle ça le résidu mais je ne sais pas comment le calculer?
Est-ce que que quelqu'un pourrait m'indiquer comment faire?
Merci d'avance

[gaspard_bourgeois]

Comme évoqué plus la haut le calcul de la matrice de covariance permet de quantifier les écarts des points au plan. lire cet article https://fr.wikipedia.org/wiki/Covariance
La VSD permet d'obtenir la direction du vecteur qui possède la distance au plan la plus faible, le vecteur normal. Lire cet article https://fr.wikipedia.org/wiki/D%C3%A...on_statistique

En code python ça peut donner ça:

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
import numpy as np
def calculNormal(vectorList):
    """ Calcul le plan moyen passant par un nuage de point"""
    vectorList = np.array(vectorList)
    covariance = np.cov(vectorList.transpose()) # calcul de la matrice de covariance
    u = np.linalg.svd(covariance)[0] # calcul des vecteurs propres ordonnées par valeur propre décroissante
    vecNormal = u[:,-1] # vecteur propre associé à la valeur propre représentant la distance des points au plan la plus faible
    vecNormalNorme = vecNormal / np.linalg.norm(vecNormal) # calcul du vecteur normal normée
    barycentre = np.average(vectorList, axis=0) # calcul du barycentre dans le repère orthogonal d'origine
    d = - np.vdot(vecNormalNorme, barycentre) # calcul de la distance du barycentre au plan
    a, b, c = vecNormalNorme # extraction des coordonnées du vecteur directeur

Code à tester.