Equation de plan à partir de points [Résolu]

[ocelote]

Bonjour,

j'ai suivi ce tutoriel pour déterminer un plan a partir de plusieurs points.

J'ai essayé de suivre au mieux les indications, mais mon resultat est peu probant :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    X = Pts(:,1)
    Y = Pts(:,2)
    Z = Pts(:,3)
    U = ones(length(X), 1);
    M = [U X Y X.^2 X.*Y Y.^2]
    K = M\Z;
    Z;
    Z_modele = M*K;
    residus = Z-Z_modele; 
    variance_d_origine = var(Z);
    variance_expliquee = var(Z_modele);
    variance_residuelle = var(residus);
    coefficient_correlation_r2 = variance_expliquee/variance_d_origine

Ou X contient tous les X de mes points (x1,x2,...xn), Y et Z egalement sur le meme principe.

Je souhaiterais dessiner mon plan calculé mais je ne sais pas comment m'y prendre.

Merci

[davcha]

Une autre façon de faire consiste à calculer les valeurs singulières de la matrice formée par tes points, centrés en zéro.

En gros :

Code :

1
2
3
4
5
6
7
M=Pts - mean(Pts)
[u s v]=svd(M)
n=v(:,3) % normale à l'hyperplan.
n/=n(3)
[X Y]=meshgrid(-10:10)
Z=-n(1)*X-n(2)*Y
mesh(X,Y,Z)

[ocelote]

Merci pour ton aide, mais il y a un probleme a la ligne : n/=n(3) et je n'ai pas suffisamment compris pour corriger moi meme

[davcha]

n/=n(3) c'est équivalent à n=n/n(3).

[ocelote]

je vois, mais j'ai une erreur de dimension en utilisant mean(). Je dois soustraire la valeur moyenne de X, Y et Z pour chaque point ?

Par exemple : X(1) - moy(X) ?

[davcha]

oui.

[ocelote]

Merci beaucoup, cela fonctionne "presque". J'obtiens en effet un plan avec la bonne inclinaison, mais j'aurais aimé qu'ils soit proche des points, il se positionne pas toujours correctement, j'ajoute + mean(Z)

Code :

1
2
 
    Z=-n(1)*X-n(2)*Y + mean(Z)

[ocelote]

Cette méthode donne la bonne inclinaison du plan, mais la position en Z est erronée.

Mes modifications precedentes n'y ont rien changé

[davcha]

ok, une autre façon de faire...

Code :

1
2
3
4
5
6
7
8
9
10
11
12
M = [Pts ones(size(Pts,1),1)];
[u s v] = svd(M);
v=v(:,4);
# v(1)*x + v(2)*y + v(3)*z + v(4) = 0
v/=v(3);
v=[v(1:2) v(4)];

# nouveaux points...
x = [randn(1000,2) ones(1000,1)];
x(:,3) = x*v;

# tous les x sont sur ton plan.

[ocelote]

Au risque de jouer les rabajoie, j'obtiens une erreur de dimension a

Code :

v=[v(1:2) v(4)];

et la

Code :

x(:,3) = x*v;

Merci quand meme pour ton temps

[davcha]

v=[v(1:2);v(4)];

pardon.

[ocelote]

J'ai essayé de modifier mon code de cette maniere, mais je crois faire fausse route. Le probleme est que je n'ai qu une vague idee de ce que je poste la...

Code :

1
2
3
4
5
6
7
8
9
10
11
12
M = [Pts ones(size(Pts,1),1)];
    [u s v] = svd(M);
    v=v(:,4);
    % # v(1)*x + v(2)*y + v(3)*z + v(4) = 0
    v=v/v(3);
    v=[v(1:2);v(4)];
    [X Y]=meshgrid(0:5:300);
    Z=-v(1)*X-v(2)*Y + mean(Z);
    hold on
    o = mesh(X,Y,Z,'FaceLighting','gouraud','LineWidth',0.3);
    colormap hot;

[davcha]

La ligne 8 est pas bonne.

Quand tu fais la SVD de la matrice M, la 4e colonne des vecteurs propres à droite correspond à l'équation du plan, comme indiqué dans la ligne 4.

Donc, pour obtenir Z, la 3e coordonnée, tu fais : -(v(1)*x+v(2)*y+v(4)) / v(3)
Donc, pour simplifier le calcul, tu pars de ton vecteur v original (celui de la ligne 3), tu divises chaque composante par v(3). Ensuite tu crées un nouveau vecteur V=-[v(1:2);v(4)]. A partir de là, tout point [x y 1] multiplié par le vecteur V te donnera la 3e coordonnée, Z.

Dans le cas des meshgrid, tu dois donc faire ceci :

Z=-(v(1)*x+v(2)*y+v(4))/v(3)

[ocelote]

Merci beaucoup, j'ai enfin trouvé !

j'avoue m'etre aidé de ce poste, mais c'est juste parce que j'ai du mal.

Je donne quand meme la solution car il y avait un pti souci avec le vecteur v(4):

Code :

1
2
3
4
5
6
7
8
9
10
11
M = [Pts ones(size(Pts,1),1)];
    [u s v] = svd(M);
    v=v(:,4);
    % # v(1)*x + v(2)*y + v(3)*z + v(4) = 0
    v=v/v(3);
    v=-[v(1:2);v(4)];
    [X Y]=meshgrid(0:5:300);
    Z=v(1)*X+v(2)*Y+v(3)
    hold on
    o = mesh(X,Y,Z);

Une autre méthode (celle du lien) plus simple pour moi sans svd() :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
    x = Pts(:,1)';
    y = Pts(:,2)';
    z = Pts(:,3)';
    M = [x; y; ones(1,size(Pts,1))];
    c = z/M;
    [X Y]=meshgrid(0:5:300);
    Z=c(1)*X+c(2)*Y+c(3);
    hold on
    plot3(x,y,z,'r*')
    hold on
    o = mesh(X,Y,Z)

Conclusion : je suis vraiment bidon en maths.

Merci beaucoup pour ton aide !