Approximation d'un plan

**oieretxe** · 15/05/2014, 23h14

Bonjour,

le but de mon programme sera au final de réussir à approximer un plan à partir de points expérimentaux formant un anneau dans l'espace. Pour le moment, je suis en train de réaliser la théorie, je travaille donc sur des points tirés d'un anneau "parfait" auxquels j'ai ajouté du décalage aléatoire. Pour résoudre le système d'équation, je suis parti sur l'équation d'un plan : aX + bY + cZ + d = 0, je suis donc parti de M = [-1 -X -Y] et Res = [d a b] en forçant c = 1. X, y et Z étant les coordonnées de mes points expérimentaux. J'obtiens des résultats avec Res = M \ Z, mais lorsque je trace le plan, il me parait totalement faux... J'en reviens donc à vous car je n'arrive pas à trouver mon erreur...

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
% On vide la mémoire
close all
clear all
clc
 
%% Tolerance_A.m
 
% Nombre de points nominaux de l'anneau
N = 360;
 
% Centre du cylindre
x = 0;
y = 0;
z = 0;
 
% Rayon max de l'anneau
R_max = 15;
 
% Rayon min de l'anneau
R_min = 12;
 
% Définition de l'anneau nominal
XAnneau_nom = [];
YAnneau_nom = [];
ZAnneau_nom = [];
for VTheta = 0:(2*pi)/N:2*pi
    p = rand(1);
    X = x  + (R_min + (R_max-R_min) * p) * cos(VTheta);
    Y = y * ones(size(VTheta));
    Z = z  + (R_min + (R_max-R_min) * p) * sin(VTheta);
    XAnneau_nom = [XAnneau_nom X];
    YAnneau_nom = [YAnneau_nom Y];
    ZAnneau_nom = [ZAnneau_nom Z];
end
 
Anneau_nom = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
 
% Affichage de l'anneau nominale
figure
plot3(Anneau_nom(1,:), Anneau_nom(2,:), Anneau_nom(3,:), '*r'); hold on
xlabel ('X')
ylabel('Y')
zlabel ('Z')
% Définition de l'anneau expérimental
Ordre = 1 / 10000;
Anneau_exp = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
for i = 1 : N + 1
    p = rand(1);
    Anneau_exp(2, i) = Anneau_exp(2, i) + 2 * Ordre * p - Ordre;
end
% Affichage de l'anneau expérimental
plot3(Anneau_exp(1,:), Anneau_exp(2,:), Anneau_exp(3,:), '*b');
grid on
 
% Approximation du plan
% Nous savons que l'équation d'un plan est de la forme : aX + bY + cZ + d = 0
% On pose c = 1, donc z = - (d + aX + bY)
% Nous allors résoudre cette équation
i = size(Anneau_exp, 2);
M = [-ones(i, 1) -(Anneau_exp(1, :).') -(Anneau_exp(2, :).')];
Res = M \ ((Anneau_exp(3, :).'));
 
a = Res(2); b = Res(3); c = 1 ; d = Res(1);
 
% Affichage du plan
min_x = min(Anneau_exp(1,:));
max_x = max(Anneau_exp(1,:));
min_z = min(Anneau_exp(3,:));
max_z = max(Anneau_exp(3,:));
 
xLim = (min_x : 0.1 : max_x);
zLim = (min_z : 0.1 : max_z);
[X,Z] = meshgrid(xLim,zLim);
Y = (a * X + c * Z + d)/ (-b);
surf(X, Y, Z)

Merci d'avance !

**le fab** · 16/05/2014, 08h53

salut

ton équation de plan sous forme matricielle est :

Code :

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

1
2
3
[X Y Z] * [a    [-d
           b     -d
           c] =  -d] = -d*ones(size(X))

ou X, Y et Z sont des vecteurs colonnes de taille identique

donc [a;b;c] = [X Y Z]\(-d*ones(size(X)));
plus de lecture par ici

Fabien

**oieretxe** · 16/05/2014, 16h17

Merci pour ta réponse.

Par contre je ne comprends pas trop ton équation :

Code :

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

1
2
3
[X Y Z] * [a    [-d
           b     -d
           c] =  -d] = -d*ones(size(X))

J'ai été sur le lien que tu m'as donné, et en suivant les étapes je trouve bien :
aX + bY + cZ + d = 0 (équation du plan)
Z = (-d - aX - bY)/c
Après on pause c = 1, car il y a une infinité de solution, il faut donc forcer un paramètre.
Donc Z = -d - aX - bY

On pose Res = [d a b].' et M = [-ones(i) -X -Y] avec i : le nombre de points.

Ce qui donne finalement : M * Res = Z

D'où Res = M\Z

Il y a un soucis ?

Tu m'as écris :

Code :

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

donc [a;b;c] = [X Y Z]\(-d*ones(size(X)));

Tu penses qu'il vaut mieux forcer la valeur de d plutôt que celle de c ? Cela change-t-il quelque chose ?

En tout cas merci de ta réponse !

**oieretxe** · 16/05/2014, 17h36

J'ai du nouveau ! Je trouve un plan cohérent avec ce code :

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
% On vide la mémoire
close all
clear all
clc
 
%% Tolerance_A.m
 
% Nombre de points nominaux de l'anneau
N = 360;
 
% Centre de l'anneau
x = 0;
y = 1;
z = 0;
 
% Rayon max de l'anneau
R_max = 15;
 
% Rayon min de l'anneau
R_min = 12;
 
% Définition de l'anneau nominal
XAnneau_nom = [];
YAnneau_nom = [];
ZAnneau_nom = [];
for VTheta = 0:(2*pi)/N:2*pi
    p = rand(1);
    X = x  + (R_min + (R_max-R_min) * p) * cos(VTheta);
    Y = y * ones(size(VTheta));
    Z = z  + (R_min + (R_max-R_min) * p) * sin(VTheta);
    XAnneau_nom = [XAnneau_nom X];
    YAnneau_nom = [YAnneau_nom Y];
    ZAnneau_nom = [ZAnneau_nom Z];
end
 
Anneau_nom = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
 
% Affichage de l'anneau nominale
figure
plot3(Anneau_nom(1,:), Anneau_nom(2,:), Anneau_nom(3,:), '*r'); hold on
xlabel ('X')
ylabel('Y')
zlabel ('Z')
% Définition de l'anneau expérimental
Ordre = 1 / 10000;
Anneau_exp = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
for i = 1 : N + 1
    p = rand(1);
    Anneau_exp(2, i) = Anneau_exp(2, i) + 2 * Ordre * p - Ordre;
end
% Affichage de l'anneau expérimental
plot3(Anneau_exp(1,:), Anneau_exp(2,:), Anneau_exp(3,:), '*b');
grid on
 
% Résolution du plan
% Nous savons que l'équation d'un plan est de la forme : aX + bY + cZ + d = 0
% On pose aX + bY + cZ = -d
% On fixe la valeur de d = 1
% Nous allors résoudre cette équation
d = 4;
X = (Anneau_exp(1, :).');
Y = (Anneau_exp(2, :).');
Z = (Anneau_exp(3, :).');
 
i = size(X);
M = [X Y Z];
 
Res = M \ (-d * ones(i));
 
a = Res(1); b = Res(2); c = Res(3) ;
 
% Affichage du plan
[Xp, Zp] = meshgrid(-R_max:0.5:R_max); 
Yp = @(Xp,Zp) (-d - a*Xp - c*Zp)/b;
mesh(Xp,Yp(Xp,Zp),Zp);

J'ai donc suivi ton conseil en résolvant aX + bY + cZ = -d.

Le problème est que lorsque le centre de mon anneau à y = 0, alors là le plan se décale complètement, je ne sais pas pourquoi...

Code :

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

1
2
3
4
% Centre de l'anneau
x = 0;
y = 0;
z = 0;

Encore merci !

**oieretxe** · 18/05/2014, 22h58

Personne ne saurait comment palier au problème qu'engendre le fait de centrer les points sur y = 0 svp ? ça fait 2 jours que je suis dessus et je n'arrive pas à trouver d'où vient le problème ...

**le fab** · 19/05/2014, 08h53

difficile de te répondre comme ca, donne des jeux de données (qui marche te qui marche "pas")

**oieretxe** · 19/05/2014, 15h11

Un code où le plan ne se positionne pas correctement :

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
% On vide la mémoire
close all
clear all
clc
 
%% Tolerance_A.m
% Version 2 : ajout d'une rotation pour tester l'approximation du plan
 
% Nombre de points nominaux de l'anneau
N = 360;
 
% Centre de l'anneau
x = 0;
y = 0;
z = 0;
 
% Rayon max de l'anneau
R_max = 15;
 
% Rayon min de l'anneau
R_min = 12;
 
% Définition de l'anneau nominal
XAnneau_nom = [];
YAnneau_nom = [];
ZAnneau_nom = [];
for VTheta = 0:(2*pi)/N:2*pi
    p = rand(1);
    X = x  + (R_min + (R_max-R_min) * p) * cos(VTheta);
    Y = y * ones(size(VTheta));
    Z = z  + (R_min + (R_max-R_min) * p) * sin(VTheta);
    XAnneau_nom = [XAnneau_nom X];
    YAnneau_nom = [YAnneau_nom Y];
    ZAnneau_nom = [ZAnneau_nom Z];
end
 
Anneau_nom = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
 
% Affichage de l'anneau nominale
figure
plot3(Anneau_nom(1,:), Anneau_nom(2,:), Anneau_nom(3,:), '*r'); hold on
xlabel ('X')
ylabel('Y')
zlabel ('Z')
% Définition de l'anneau expérimental
Ordre = 1 / 10000;
Anneau_exp = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
for i = 1 : N + 1
    p = rand(1);
    Anneau_exp(2, i) = Anneau_exp(2, i) + 2 * Ordre * p - Ordre;
end
 
% Valeurs de rotation suivant les axes (radian)
Alpha = 0; % rotation autour de l'axe x
Betha = 0; % rotation autour de l'axe y
Gamma = 0; % rotation autour de l'axe z
 
Rotation_x = [1, 0, 0;
            0, cos(Alpha), -sin(Alpha);
            0, sin(Alpha), cos(Alpha)];
Rotation_y = [cos(Betha), 0, sin(Betha);
               0, 1, 0;
            -sin(Betha), 0, cos(Betha)];
Rotation_z = [cos(Gamma), -sin(Gamma), 0;
            sin(Gamma), cos(Gamma), 0;
            0, 0, 1];
 
% Matrice de rotation
Rotation = Rotation_x*Rotation_y*Rotation_z;
 
% Ajout de la rotation et de la translation
Anneau_exp = Rotation * Anneau_exp;
 
% Affichage de l'anneau expérimental
plot3(Anneau_exp(1,:), Anneau_exp(2,:), Anneau_exp(3,:), '*b');
grid on
 
% Résolution du plan
% Nous savons que l'équation d'un plan est de la forme : aX + bY + cZ + d = 0
% On pose aX + bY + cZ = -d
% On fixe la valeur de d = 1
% Nous allors résoudre cette équation
d = 4;
X = (Anneau_exp(1, :).');
Y = (Anneau_exp(2, :).');
Z = (Anneau_exp(3, :).');
 
i = size(X);
M = [X Y Z];
 
Res = M \ (-d * ones(i));
 
a = Res(1); b = Res(2); c = Res(3) ;
 
% Affichage du plan
[Xp, Zp] = meshgrid(-R_max:0.5:R_max); 
Yp = @(Xp,Zp) (-d - a*Xp - c*Zp)/b;
mesh(Xp,Yp(Xp,Zp),Zp);

Et un code où tout semble être bon :

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
% On vide la mémoire
close all
clear all
clc

%% Tolerance_A.m
% Version 2 : ajout d'une rotation pour tester l'approximation du plan

% Nombre de points nominaux de l'anneau
N = 360;

% Centre de l'anneau
x = 0;
y = 1;
z = 0;

% Rayon max de l'anneau
R_max = 15;

% Rayon min de l'anneau
R_min = 12;

% Définition de l'anneau nominal
XAnneau_nom = [];
YAnneau_nom = [];
ZAnneau_nom = [];
for VTheta = 0:(2*pi)/N:2*pi
    p = rand(1);
    X = x  + (R_min + (R_max-R_min) * p) * cos(VTheta);
    Y = y * ones(size(VTheta));
    Z = z  + (R_min + (R_max-R_min) * p) * sin(VTheta);
    XAnneau_nom = [XAnneau_nom X];
    YAnneau_nom = [YAnneau_nom Y];
    ZAnneau_nom = [ZAnneau_nom Z];
end

Anneau_nom = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];

% Affichage de l'anneau nominale
figure
plot3(Anneau_nom(1,:), Anneau_nom(2,:), Anneau_nom(3,:), '*r'); hold on
xlabel ('X')
ylabel('Y')
zlabel ('Z')
% Définition de l'anneau expérimental
Ordre = 1 / 10000;
Anneau_exp = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
for i = 1 : N + 1
    p = rand(1);
    Anneau_exp(2, i) = Anneau_exp(2, i) + 2 * Ordre * p - Ordre;
end

% Valeurs de rotation suivant les axes (radian)
Alpha = 0; % rotation autour de l'axe x
Betha = 0; % rotation autour de l'axe y
Gamma = 0; % rotation autour de l'axe z

Rotation_x = [1, 0, 0;
            0, cos(Alpha), -sin(Alpha);
            0, sin(Alpha), cos(Alpha)];
Rotation_y = [cos(Betha), 0, sin(Betha);
               0, 1, 0;
            -sin(Betha), 0, cos(Betha)];
Rotation_z = [cos(Gamma), -sin(Gamma), 0;
            sin(Gamma), cos(Gamma), 0;
            0, 0, 1];
        
% Matrice de rotation
Rotation = Rotation_x*Rotation_y*Rotation_z;

% Ajout de la rotation et de la translation
Anneau_exp = Rotation * Anneau_exp;

% Affichage de l'anneau expérimental
plot3(Anneau_exp(1,:), Anneau_exp(2,:), Anneau_exp(3,:), '*b');
grid on

% Résolution du plan
% Nous savons que l'équation d'un plan est de la forme : aX + bY + cZ + d = 0
% On pose aX + bY + cZ = -d
% On fixe la valeur de d = 1
% Nous allors résoudre cette équation
d = 4;
X = (Anneau_exp(1, :).');
Y = (Anneau_exp(2, :).');
Z = (Anneau_exp(3, :).');

i = size(X);
M = [X Y Z];

Res = M \ (-d * ones(i));

a = Res(1); b = Res(2); c = Res(3) ;

% Affichage du plan
[Xp, Zp] = meshgrid(-R_max:0.5:R_max); 
Yp = @(Xp,Zp) (-d - a*Xp - c*Zp)/b;
mesh(Xp,Yp(Xp,Zp),Zp);

La seule chose qui change c'est la position du centre de l'anneau, nominal, avant l'ajout de décalage. Dans le cas où cet anneau est centré sur y = 0, le plan est complètement décalé et n'approxime pas les points "expérimentaux" correctement. Je ne vois pas pourquoi j'ai cette erreur.

Merci d'avance !

**le fab** · 19/05/2014, 17h17

excuses moi je t'ai complètement enduit d'erreur

je te laisse regarder ca :

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
% On vide la mémoire
close all
clear all
clc
 
%% Tolerance_A.m
% Version 2 : ajout d'une rotation pour tester l'approximation du plan
 
% Nombre de points nominaux de l'anneau
N = 360;
 
% Centre de l'anneau
x = 0;
y = 1;
z = 0;
 
% Rayon max de l'anneau
R_max = 15;
 
% Rayon min de l'anneau
R_min = 12;
 
% Définition de l'anneau nominal
XAnneau_nom = [];
YAnneau_nom = [];
ZAnneau_nom = [];
for VTheta = 0:(2*pi)/N:2*pi
    p = rand(1);
    X = x  + (R_min + (R_max-R_min) * p) * cos(VTheta);
    Y = y * ones(size(VTheta));
    Z = z  + (R_min + (R_max-R_min) * p) * sin(VTheta);
    XAnneau_nom = [XAnneau_nom X];
    YAnneau_nom = [YAnneau_nom Y];
    ZAnneau_nom = [ZAnneau_nom Z];
end
 
Anneau_nom = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
 
% Affichage de l'anneau nominale
figure
plot3(Anneau_nom(1,:), Anneau_nom(2,:), Anneau_nom(3,:), '*r'); hold on
xlabel ('X')
ylabel('Y')
zlabel ('Z')
% Définition de l'anneau expérimental
Ordre = 1 / 10000;
Anneau_exp = [XAnneau_nom; YAnneau_nom; ZAnneau_nom];
for i = 1 : N + 1
    p = rand(1);
    Anneau_exp(2, i) = Anneau_exp(2, i) + 2 * Ordre * p - Ordre;
end
 
% Valeurs de rotation suivant les axes (radian)
Alpha = 0; % rotation autour de l'axe x
Betha = 0; % rotation autour de l'axe y
Gamma = 0; % rotation autour de l'axe z
 
Rotation_x = [1, 0, 0;
            0, cos(Alpha), -sin(Alpha);
            0, sin(Alpha), cos(Alpha)];
Rotation_y = [cos(Betha), 0, sin(Betha);
               0, 1, 0;
            -sin(Betha), 0, cos(Betha)];
Rotation_z = [cos(Gamma), -sin(Gamma), 0;
            sin(Gamma), cos(Gamma), 0;
            0, 0, 1];
 
% Matrice de rotation
Rotation = Rotation_x*Rotation_y*Rotation_z;
 
% Ajout de la rotation et de la translation
Anneau_exp = Rotation * Anneau_exp;
 
% Affichage de l'anneau expérimental
plot3(Anneau_exp(1,:), Anneau_exp(2,:), Anneau_exp(3,:), '*b');
grid on
 
% Résolution du plan
% Nous savons que l'équation d'un plan est de la forme : aX + bY + cZ + d = 0
% On pose aX + bY + cZ = -d
% On fixe la valeur de d = 1
% Nous allors résoudre cette équation
X = (Anneau_exp(1, :).');
Y = (Anneau_exp(2, :).');
Z = (Anneau_exp(3, :).');
 
i = size(X);
M = [X Z ones(i)];
 
Res = M \ (-Y);
 
a = Res(1); b = 1; c = Res(2) ;d = Res(3);
 
% Affichage du plan
[Xp, Zp] = meshgrid(-R_max:0.5:R_max); 
Yp = @(Xp,Zp) (-d - a*Xp - c*Zp)/b;
mesh(Xp,Yp(Xp,Zp),Zp);

en gros :
il ne faut pas fixer d=1 !! (cela t'impose un plan ne passant pas par l'origine)
ensuite il faut écrire ton équation de plan sous la forme par exemple
- Z = A*X + B*Y + D ( = a/c*X + b/c *Y + d/c) ce qui suppose c~=0 donc ca ne couvre pas tous les plans possible !!
puis res = [X Y ones(size(X)]\(-Z);

ceci est un exemple car tu ne peux pas couvrir tous les cas comme ca (c=0, ca marche pas, faut changer de variable)

prend le temps de bien lire l'excellente contribution que j'ai mis en lien plus haut

il va falloir que tu détermine si ton plan est parallèle à un des plan XY, YZ, ZX puis choisir par quelle inconue (a,b,c,d) tu divise ton équation de plan

Fabien

**oieretxe** · 21/05/2014, 17h49

Merci beaucoup !

tu m'as complètement débloqué, ça marche à la perfection !

Je marque le sujet résolu !

Encore merci !

Approximation d'un plan

MATLAB

Discussions similaires

Partager

Partager