ajustement plan moindres carrés

Posté par bernard6 (invité)


Bonjour,

J'aimerais retrouver les paramètres a b c d'un plan qui s'ajuste au mieux (par les moindres carrés) à un ensemble de points de coordonnées x1 y1 z1... xN yN zN

J'ai trouvé comme solution sur le net le résultat suivant (voir pièce jointe)

Sauriez-vous me dire si cette solution est correcte parce que lorsque je l'implémente, j'obtiens une valeur de c beaucoup plus grande que celle attendue...

Merci d'avance !

Posté par JJa JJa

C'est bon, tu peux y aller sans hésitation !
Un petit coup de pub à l'occasion (au sujet des méthodes d'ajustement) :
Magazine QUADRATURE n°63, pp.33-40
http:///www.edpsciences.org/quadrature/

Posté par bernard6 (invité)

Ok Merci !

PS : votre lien internet ne fonctionne pas  

Posté par JJa JJa

Normalement, ça marche.
Si non, il y a probablement un bug de programmation. Dans ce cas, montre ton listing et on verra...

Posté par J-P J-P

JJ,

C'est le lien internet que tu as donné qui ne fonctionne pas.

Peut-être était ce celui-ci:

Posté par JJa JJa

OK, c'est bien cela. J'avais recopié une adresse trop ancienne. Merci pour la correction.

Posté par bernard6 (invité)

En fait, il s'avère que j'ai des résultats qui ne me satisfont pas vraiment...

par exemple, j'ai repris trois points et déterminé l'équation du plan passant par ces trois points (a b c) et les ai comparé avec les a b c déterminés par cette technique et il s'avère que ce ne sont pas les mêmes...

A tout hasard voici mon code matlab mais je l'ai vérifié en effectuant les opérations sous Excel et il semble correct ...

% 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

Posté par JJa JJa

En ne donnant que trois points, tu dois trouver exactement les paramètres du plan qui passe par les trois points. C'est garanti (sauf si les points sont alignés).
Il est clair que si tu ne trouves pas le résultat exact, il y a un bug de programmation ( dans l'emplacement des données ? ou peut-être dans la vérification elle-même ? )

Posté par JJa JJa

Si tu n'y arrives pas, montre ton exemple numérique (coordonnées des 3 points, matrice(M) et vecteur(V) que tu calcules). Je comparerai avec ce que donne mon programme (j'en ai un qui fait cela)

Posté par JJa JJa

Au fait, le lien que j'avais donné aurait du marcher. Je viens de vérifier.
il y avait une faute de frappe (un / en trop).
Mais, puisque le lien donné par J-P marche, n'en demandons pas plus !

Posté par bernard6 (invité)

Bonjour,

Prenons un exemple avec les 3 points :
x y z
1 2 3
6 2 2
3 3 3

Mon code pour retrouver les paramètres a b c du plan passant par ces trois points est :

a = (Y2 - Y1)*(Z3 - Z1) - (Z2 - Z1)*(Y3 - Y1)
b = (-1)*((X2 - X1)*(Z3 - Z1) - (Z2 -Z1)*(X3 - X1))
c = (X2 - X1)*(Y3 - Y1) - (Y2 - Y1)*(X3 - X1)

J'obtiens comme valeurs :
a 1
b -2
c 5


Mon code pour ajuster un plan à ces points par moindres carrés est :
% 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
a = MatriceX(1,1)
b = MatriceX(2,1)
c = MatriceX(3,1)

NB : j'ai vérifié ce code en faisant les calculs en parallèle sur Excel

J'obtiens comme valeurs :

a -0.2
b 0.4
c 2.4


Comme vous l'avez dit, ce résultat n'est pas normal car je devrais retrouver le plan passant par les trois points...

Je vous remercie d'avance pour l'aide que vous pourrez m'apporter !

Posté par JJa JJa

C'est le premier calcul de a, b, c qui est faux.
En effet, avec a=1, b=-2, c=5, pour le premier point (x=1, y=2, z=3) on trouve :
z = a*x+b*y+c = 1*1 -2*2 +5 = 2 au lieu de 3.
Je n'ai pas vérifié votre code correspondant. Vous devriez trouver aisément l'erreur.
Au contraire, le calcul avec la méthode des moindres carrés donne le bon résultat :
a=-0.2 , b=0.4 , c=2.4, et par exemple pour le premier point :
z = -0.2*1 +0.4*2 +2.4 = 3 ce qui est correct.
(j'avais vérifié votre code correspondant. Je n'avais pas trouvé d'erreur)

Posté par bernard6 (invité)

je ne comprends pas très bien la manière dont vous considérer le c...

Si on s'en tient à l'équation d'un plan :

ax + by +cz + d = 0

Ce serait plutôt le premier calcul qui serait correct et pas le deuxième...

Je dois sans doute avoir considéré a b c de façon différente pour les deux calculs, non?

A mon avis, l'erreur doit venir de là...

Maintentant le problème reste à exprimer a1 b1 c1 en fonctier de a2 b2 c2...

Enfin je pense...

Posté par bernard6 (invité)

voilà pour la conversion des paramètres :

a2 = -a1/c1
b2 = -b1/c1
c2 = -d1/c1

Merci pour ton intervention JJa !

Sans toi, j'aurais pas la distinction !

Bon we...

Posté par JJa JJa

Si vous changez de notations en cours de route, vous n'y arriverez jamais !!!
Vous appliquez des formules de calcul figurant dans votre message  posté le 05/04/2007 à 14:29. Il y est très clairement écrit :
z = a*x+b*y+c   (formule 1)
Maintenant vous utilisez une formule différente :
ax + by +cz + d = 0   (formule 2)
Ce ne sont plus les mêmes paramètres : par exemple, le b de la première n'est pas pareil que le b de la seconde.
La méthode des moindres carrés ayant été établie pour la formule 1, c'est dans cette formule que l'on doit mettre les (a, b, c) calculés et non pas dans une formule différente, qui donnera évidemment un résultat faux.
Bien entendu, il serait possible de partir de la formule 2 pour la méthode des moindres carré, mais alors, les équations seraient différentes de celles qui apparaissent dans le document que vous utilisez.
Alors, bon courage pour la suite, votre problème est quasiment résolu ...

Posté par JJa JJa

Lorsque j'ai posté mon message précédent, je n'avais pas encore le tiens de 13:52.
Donc, c'est clair maintenant.
Tout est bien qui fini bien.

Posté par mathjava mathjava

Bonsoir,

enfaite j'ai un problème pareil : je veux trouver un plan qui passe approximativement par un ensemble de points dans l'espace et je trouve la méthode des moindres carrées est la meilleure.
je dois implémenter mon code en java ce que je demande comment résoudre le produit des matrices pour trouver la matrice (a,b,c)

(matrice des sommes1 ) * (matrice a,b,c) = (matrice des sommes2)

est ce que
(matrice a,b,c)   =   inverse(matrice des sommes1 )*(matrice des sommes2) ??

merci