Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Mathématiques Discussion :

Équation d'un plan défini par N points


Sujet :

Mathématiques

  1. #1
    Membre à l'essai
    Équation d'un plan défini par N points
    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.

  2. #2
    Expert éminent sénior
    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)
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  3. #3
    Rédacteur

    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.
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  4. #4
    Rédacteur

    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.
    Ce qu'on trouve est plus important que ce qu'on cherche.
    Maths de base pour les nuls (et les autres...)

  5. #5
    Membre à l'essai
    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;

  6. #6
    Membre du Club
    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 !

  7. #7
    Nouveau Candidat au Club
    Calcul du résidu
    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

  8. #8
    Nouveau Candidat au Club
    Solution Python simple à l'aide de la bibliothèque Numpy
    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.