trouver la matrice de rotation à appliquer à un plan

**oieretxe** · 26/05/2014, 18h26

Bonjour,

je suis en train d'essayer de retrouver une matrice de rotation, je m'explique :

J'ai comme données un ensemble de points "expérimentaux" et je réussi à approximé un plan passant par ces points. J'aimerais maintenant pouvoir trouver la matrice de rotation à appliquer à ce plan afin qu'il soit parallèle au plan XZ, et appliquer cette matrice de rotation à l'ensemble de mes points afin de simplifier la suite de mon travail. Je vous mets le début de mon code pour être plus clair :

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
% On vide la mémoire
close all
clear all
clc
 
%% Plan_rotation.m
% Version 1 : trouver la matrice de rotation pour aligner le plan approximé
 
%% DONNEES
 
    % Nombre de points nominaux des anneaux
N = 360;
 
    % Centre de l'anneau 1
x_1 = 0;
y_1 = 0;
z_1 = 0;
 
    % Rayon max de l'anneau 1
R_max_1 = 15;
 
    % Rayon min de l'anneau 1
R_min_1 = 12;
 
    % Définition de l'anneau nominal 1 et 2
XAnneau_nom_1 = [];
YAnneau_nom_1 = [];
ZAnneau_nom_1 = [];
 
for VTheta = 0:(2*pi)/N:2*pi
    p_1 = rand(1);
    X_1 = x_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * cos(VTheta);
    Y_1 = y_1 * ones(size(VTheta));
    Z_1 = z_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * sin(VTheta);
    XAnneau_nom_1 = [XAnneau_nom_1 X_1];
    YAnneau_nom_1 = [YAnneau_nom_1 Y_1];
    ZAnneau_nom_1 = [ZAnneau_nom_1 Z_1];
 
end
 
    % Définition des anneaux expérimentaux 1 et 2
Ordre = 1 / 10000;
Anneau_exp_1 = [XAnneau_nom_1; YAnneau_nom_1; ZAnneau_nom_1];
 
for i_1 = 1 : N + 1
    p_1 = rand(1);
    Anneau_exp_1(2, i_1) = Anneau_exp_1(2, i_1) + 2 * Ordre * p_1 - Ordre;
end
 
X_1 = (Anneau_exp_1(1, :).');
Y_1 = (Anneau_exp_1(2, :).');
Z_1 = (Anneau_exp_1(3, :).');
b_1 = 10;
i_1 = size(X_1);
 
%% RESOLUTION
 
    % Plan approximé
M_1 = [X_1/b_1 Z_1/b_1 ones(i_1)/b_1];
Res_1 = M_1 \ (-Y_1);
a_1 = Res_1(1); c_1 = Res_1(2) ;d_1 = Res_1(3);
 
%% AFFICHAGE
 
    % Affichage du plan
figure('units','normalized','outerposition',[0 0 1 1])
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '*b'); hold on
 
[Xp_1, Zp_1] = meshgrid(-R_max_1 + x_1 : 0.5 : R_max_1 + x_1 , -R_max_1 + z_1 : 0.5 : R_max_1 + z_1); 
 
Yp_1 = @(Xp_1,Zp_1) (-d_1 - a_1 * Xp_1 - c_1 * Zp_1) / b_1;
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); hold off
title('Plan approximé de la tolérance A')
xlabel ('X')
ylabel('Y')
zlabel ('Z')
grid on

L'équation d'un plan étant : a*X + b*Y + c*Z + d = 0, un plan parallèle au plan XZ serait de la forme b*Y + d = 0, mais je ne sais pas comment faire le lien entre les 2. Je pensais utiliser une équation du genre : plan_quelconque * mat_rotation = plan_parallèle, et le résoudre en utilisant la division par la gauche \.
Vais-je dans la bonne direction ?

Merci d'avance !

**ERCO503** · 26/05/2014, 23h10

Bonjour,

Personnellement, j'aime bien passer par la définition d'un référentiel, comme dans ma réponse pour cette discussion : http://www.developpez.net/forums/d14...e/#post7826897

Par exemple :

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
xx = Xp_1;
yy =  Yp_1(Xp_1, Zp_1);
zz =  Zp_1;
oo = [xx(1);yy(1);zz(1)];
 
axex = ([xx(1,end);yy(1,end);zz(1,end)]-oo)/norm([xx(1,end);yy(1,end);zz(1,end)]-oo);
tmpz = ([xx(end,1);yy(end,1);zz(end,1)]-oo)/norm([xx(end,1);yy(end,1);zz(end,1)]-oo);
axey = cross(tmpz,axex);
axez = cross(axex,axey);
 
rotx = @(degrees) [1,0,0;0,cosd(degrees),-sind(degrees);0,sind(degrees),cosd(degrees)];
pts = rotx(-90)*[axex';axey';axez']*bsxfun(@minus,[Anneau_exp_1(1,:);Anneau_exp_1(2,:);Anneau_exp_1(3,:)],oo);
 
figure;
plot3(pts(1,:),pts(2,:),pts(3,:),'.');axis equal;hold on;
quiver3(0,0,0,axex(1),axex(2),axex(3),0,'r');
quiver3(0,0,0,axey(1),axey(2),axey(3),0,'g');
quiver3(0,0,0,axez(1),axez(2),axez(3),0,'b');
hold off;
xlabel ('X')
ylabel('Y')
zlabel ('Z')
axis equal;
grid on

**Paul-Ceria** · 27/05/2014, 00h36

Je t'ai mis ma ptite boucle à la fin de ton code, il te suffit juste de transformer la boucle en indexation vectorielle si tu veux que ça aille plus vite le jour où tu utilises beaucoup plus de points, trop la flemme de réfléchir là

J'ai utilisé le produit vectoriel afin de redresser tous les points de ton plan sur le vecteur [x 0 0]

Code m :

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
% On vide la mémoire
close all
clear all
clc
 
%% Plan_rotation.m
% Version 1 : trouver la matrice de rotation pour aligner le plan approximé
 
%% DONNEES
 
    % Nombre de points nominaux des anneaux
N = 360;
 
    % Centre de l'anneau 1
x_1 = 0;
y_1 = 0;
z_1 = 0;
 
    % Rayon max de l'anneau 1
R_max_1 = 15;
 
    % Rayon min de l'anneau 1
R_min_1 = 12;
 
    % Définition de l'anneau nominal 1 et 2
XAnneau_nom_1 = [];
YAnneau_nom_1 = [];
ZAnneau_nom_1 = [];
 
for VTheta = 0:(2*pi)/N:2*pi
    p_1 = rand(1);
    X_1 = x_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * cos(VTheta);
    Y_1 = y_1 * ones(size(VTheta));
    Z_1 = z_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * sin(VTheta);
    XAnneau_nom_1 = [XAnneau_nom_1 X_1];
    YAnneau_nom_1 = [YAnneau_nom_1 Y_1];
    ZAnneau_nom_1 = [ZAnneau_nom_1 Z_1];
 
end
 
    % Définition des anneaux expérimentaux 1 et 2
Ordre = 1 / 10000;
Anneau_exp_1 = [XAnneau_nom_1; YAnneau_nom_1; ZAnneau_nom_1];
 
for i_1 = 1 : N + 1
    p_1 = rand(1);
    Anneau_exp_1(2, i_1) = Anneau_exp_1(2, i_1) + 2 * Ordre * p_1 - Ordre;
end
 
X_1 = (Anneau_exp_1(1, :).');
Y_1 = (Anneau_exp_1(2, :).');
Z_1 = (Anneau_exp_1(3, :).');
b_1 = 10;
i_1 = size(X_1);
 
%% RESOLUTION
 
    % Plan approximé
M_1 = [X_1/b_1 Z_1/b_1 ones(i_1)/b_1];
Res_1 = M_1 \ (-Y_1);
a_1 = Res_1(1); c_1 = Res_1(2) ;d_1 = Res_1(3);
 
%% AFFICHAGE
% Affichage du plan
figure('units','normalized','outerposition',[0 0 1 1])
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '*b'); hold on
[Xp_1, Zp_1] = meshgrid(-R_max_1 + x_1 : 0.5 : R_max_1 + x_1 , -R_max_1 + z_1 : 0.5 : R_max_1 + z_1); 
hold on
Yp_1 = @(Xp_1,Zp_1) (-d_1 - a_1 * Xp_1 - c_1 * Zp_1) / b_1;
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); 
title('Plan approximé de la tolérance A')
xlabel ('X')
ylabel('Y')
zlabel ('Z')
grid on
 
%% A toi de jouer ;)
Xf=Xp_1;
Yf=Yp_1(Xp_1, Zp_1);
Zf=Zp_1;
for i=1:size(Xf,1)
D = [Xf(i)-Xf(i,:) ;Yf(i)-Yf(i,:); Zf(i)-Zf(i,:)]';
Num = cross(repmat([Xf(1),0,0],size(Xf,2),1),D);
Xf1(i,:)=Xf(i,:)-Num(:,3)';
end
plot3(Xf1,zeros(size(Xf1)),Zf,'.r')
%surf(Xf1,zeros(size(Xf1)),Zf)
hold off

**oieretxe** · 28/05/2014, 04h30

Merci beaucoup pour vos réponses !

Malheureusement je ne suis pas super doué en mathématiques, du coup pour le passer en Matlab je galère encore plus ^^'

J'ai testé vos codes, ils marches à merveilles, mais ce que j'aurais aimé obtenir c'est une matrice rotation(3, 3), et s'il faut une translation, qui permet au plan obtenu de devenir parallèle à XZ (avec Y=0 par exemple), car ensuite je voudrais appliquer cette rotation à tous mes points expérimentaux. Je sais pas si j'arrive à vous faire comprendre ce que j'aimerais obtenir...

Merci encore !

**oieretxe** · 28/05/2014, 13h31

Ou sinon m'expliquer un peu votre raisonnement afin que je puisse chercher de mon coté et l'appliquer afin d'obtenir ce que je recherche.

Merci en tout cas !

**ERCO503** · 28/05/2014, 15h15

Voilà mon code, en mieux, et avec des explications en commentaires.

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
99
100
101
102
103
104
105
106
107
108
109
110
% On vide la mémoire
close all
clear all
clc
 
%% Plan_rotation.m
% Version 1 : trouver la matrice de rotation pour aligner le plan approximé
 
%% DONNEES
 
    % Nombre de points nominaux des anneaux
N = 360;
 
    % Centre de l'anneau 1
x_1 = 0;
y_1 = 0;
z_1 = 0;
 
    % Rayon max de l'anneau 1
R_max_1 = 15;
 
    % Rayon min de l'anneau 1
R_min_1 = 12;
 
    % Définition de l'anneau nominal 1 et 2
XAnneau_nom_1 = [];
YAnneau_nom_1 = [];
ZAnneau_nom_1 = [];
 
for VTheta = 0:(2*pi)/N:2*pi
    p_1 = rand(1);
    X_1 = x_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * cos(VTheta);
    Y_1 = y_1 * ones(size(VTheta));
    Z_1 = z_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * sin(VTheta);
    XAnneau_nom_1 = [XAnneau_nom_1 X_1];
    YAnneau_nom_1 = [YAnneau_nom_1 Y_1];
    ZAnneau_nom_1 = [ZAnneau_nom_1 Z_1];
 
end
 
    % Définition des anneaux expérimentaux 1 et 2
Ordre = 1 / 10000;
Anneau_exp_1 = [XAnneau_nom_1; YAnneau_nom_1; ZAnneau_nom_1];
 
for i_1 = 1 : N + 1
    p_1 = rand(1);
    Anneau_exp_1(2, i_1) = Anneau_exp_1(2, i_1) + 2 * Ordre * p_1 - Ordre;
end
 
X_1 = (Anneau_exp_1(1, :).');
Y_1 = (Anneau_exp_1(2, :).');
Z_1 = (Anneau_exp_1(3, :).');
b_1 = 10;
i_1 = size(X_1);
 
%% RESOLUTION
 
    % Plan approximé
M_1 = [X_1/b_1 Z_1/b_1 ones(i_1)/b_1];
Res_1 = M_1 \ (-Y_1);
a_1 = Res_1(1); c_1 = Res_1(2) ;d_1 = Res_1(3);
 
%% AFFICHAGE
 
    % Affichage du plan
figure('units','normalized','outerposition',[0 0 1 1])
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '*b'); hold on
 
[Xp_1, Zp_1] = meshgrid(-R_max_1 + x_1 : 0.5 : R_max_1 + x_1 , -R_max_1 + z_1 : 0.5 : R_max_1 + z_1); 
 
Yp_1 = @(Xp_1,Zp_1) (-d_1 - a_1 * Xp_1 - c_1 * Zp_1) / b_1;
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); alpha(0.5);
title('Plan approximé de la tolérance A')
xlabel ('X')
ylabel('Y')
zlabel ('Z')
grid on
 
xx = Xp_1; % Données du plan en X
yy =  Yp_1(Xp_1, Zp_1); % Données du plan en Y
zz =  Zp_1; % Données du plan en Z
oo = [xx(1);yy(1);zz(1)]; % Origine du plan
 
v1 = [xx(1,end);yy(1,end);zz(1,end)]-oo; % Vecteur X du plan
v2 = [xx(end,1);yy(end,1);zz(end,1)]-oo; % Vecteur Y du plan
v3 = cross(v1,v2); % Vecteur Z du plan
axex = v1/norm(v1); % Référentiel axe X du plan
tmpy = v2/norm(v2); % Référentiel temporaire axe Y du plan (au où le plan n'aurait pas des coins orthogonaux)
axez = cross(axex,tmpy); % Référentiel axe Z du plan : Produit vectoriel des vecteurs X et Y
axey = cross(axez,axex); % Référentiel axe Y du plan : Produit vectoriel des vecteurs Z et X
 
quiver3(oo(1),oo(2),oo(3),v1(1),v1(2),v1(3),0,'c','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v2(1),v2(2),v2(3),0,'m','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v3(1)/100,v3(2)/100,v3(3)/100,0,'y','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),axex(1),axex(2),axex(3),0,'r','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axey(1),axey(2),axey(3),0,'g','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axez(1),axez(2),axez(3),0,'b','LineWidth',3);hold off;axis equal;
 
% [axex';axey';axez'] est la matrice de rotation 3x3 à appliquer au plan
% Il faut ramener les données à l'origine avant la rotation autour de l'origine, donc on soustrait l'origine "oo" du plan
pts = [axex';axey';axez']*bsxfun(@minus,[Anneau_exp_1(1,:);Anneau_exp_1(2,:);Anneau_exp_1(3,:)],oo);
 
figure;
plot3(pts(1,:),pts(2,:),pts(3,:),'.');axis equal;hold on;
hold off;
xlabel ('X')
ylabel('Y')
zlabel ('Z')
axis equal;
grid on

**oieretxe** · 28/05/2014, 18h11

Merci beaucoup pour avoir répondu aussi rapidement !

Par contre je n'ai toujours pas l'impression d'obtenir ce que je recherche. En fait j'ai l'impression que grâce à ton code j'obtiens mes points expérimentaux qui ont un z quasi constant, alors qu'avant c'était leur y qui l'était. Je vais essayer d'être plus clair dans ce que j'aimerais obtenir :
- à l'aide des points expérimentaux j'obtiens les paramètres d'un plan approximé.
- ce plan est quasiment parallèle au plan XZ comme on peut le voir dans la figure 1 du programme ci-après.
- je voudrais trouver la légère rotation à appliquer à ce plan pour qu'il soit parallèle au plan XZ, donc avec une valeur y constante
- enfin j'applique rotation à tous mes points

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
% On vide la mémoire
close all
clear all
clc
 
%% Plan_rotation.m
% Version 1 : trouver la matrice de rotation pour aligner le plan approximé
 
%% DONNEES
 
    % Nombre de points nominaux des anneaux
N = 360;
 
    % Centre de l'anneau 1
x_1 = 0;
y_1 = 0;
z_1 = 0;
 
    % Rayon max de l'anneau 1
R_max_1 = 15;
 
    % Rayon min de l'anneau 1
R_min_1 = 12;
 
    % Définition de l'anneau nominal 1 et 2
XAnneau_nom_1 = [];
YAnneau_nom_1 = [];
ZAnneau_nom_1 = [];
 
for VTheta = 0:(2*pi)/N:2*pi
    p_1 = rand(1);
    X_1 = x_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * cos(VTheta);
    Y_1 = y_1 * ones(size(VTheta));
    Z_1 = z_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * sin(VTheta);
    XAnneau_nom_1 = [XAnneau_nom_1 X_1];
    YAnneau_nom_1 = [YAnneau_nom_1 Y_1];
    ZAnneau_nom_1 = [ZAnneau_nom_1 Z_1];
 
end
 
    % Définition des anneaux expérimentaux 1 et 2
Ordre = 1 / 10000;
Anneau_exp_1 = [XAnneau_nom_1; YAnneau_nom_1; ZAnneau_nom_1];
 
for i_1 = 1 : N + 1
    p_1 = rand(1);
    Anneau_exp_1(2, i_1) = Anneau_exp_1(2, i_1) + 2 * Ordre * p_1 - Ordre;
end
 
X_1 = (Anneau_exp_1(1, :).');
Y_1 = (Anneau_exp_1(2, :).');
Z_1 = (Anneau_exp_1(3, :).');
b_1 = 10;
i_1 = size(X_1);
 
%% RESOLUTION
 
    % Plan approximé
M_1 = [X_1/b_1 Z_1/b_1 ones(i_1)/b_1];
Res_1 = M_1 \ (-Y_1);
a_1 = Res_1(1); c_1 = Res_1(2) ;d_1 = Res_1(3);
 
%% AFFICHAGE
 
    % Affichage du plan
figure('units','normalized','outerposition',[0 0 1 1])
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '*b'); hold on
 
[Xp_1, Zp_1] = meshgrid(-R_max_1 + x_1 : 0.5 : R_max_1 + x_1 , -R_max_1 + z_1 : 0.5 : R_max_1 + z_1); 
 
Yp_1 = @(Xp_1,Zp_1) (-d_1 - a_1 * Xp_1 - c_1 * Zp_1) / b_1;
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); alpha(0.5);
title('Plan approximé de la tolérance A')
xlabel ('X')
ylabel('Y')
zlabel ('Z')
grid on

Je vais essayer de poursuivre de mon coté avec ton code, j'espère réussir à trouver.

**ERCO503** · 28/05/2014, 18h39

Définis le référentiel autrement :

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
99
100
101
102
103
104
105
106
107
108
109
110
111
% On vide la mémoire
close all
clear all
clc
 
%% Plan_rotation.m
% Version 1 : trouver la matrice de rotation pour aligner le plan approximé
 
%% DONNEES
 
    % Nombre de points nominaux des anneaux
N = 360;
 
    % Centre de l'anneau 1
x_1 = 0;
y_1 = 0;
z_1 = 0;
 
    % Rayon max de l'anneau 1
R_max_1 = 15;
 
    % Rayon min de l'anneau 1
R_min_1 = 12;
 
    % Définition de l'anneau nominal 1 et 2
XAnneau_nom_1 = [];
YAnneau_nom_1 = [];
ZAnneau_nom_1 = [];
 
for VTheta = 0:(2*pi)/N:2*pi
    p_1 = rand(1);
    X_1 = x_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * cos(VTheta);
    Y_1 = y_1 * ones(size(VTheta));
    Z_1 = z_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * sin(VTheta);
    XAnneau_nom_1 = [XAnneau_nom_1 X_1];
    YAnneau_nom_1 = [YAnneau_nom_1 Y_1];
    ZAnneau_nom_1 = [ZAnneau_nom_1 Z_1];
 
end
 
    % Définition des anneaux expérimentaux 1 et 2
Ordre = 1 / 10000;
Anneau_exp_1 = [XAnneau_nom_1; YAnneau_nom_1; ZAnneau_nom_1];
 
for i_1 = 1 : N + 1
    p_1 = rand(1);
    Anneau_exp_1(2, i_1) = Anneau_exp_1(2, i_1) + 2 * Ordre * p_1 - Ordre;
end
 
X_1 = (Anneau_exp_1(1, :).');
Y_1 = (Anneau_exp_1(2, :).');
Z_1 = (Anneau_exp_1(3, :).');
b_1 = 10;
i_1 = size(X_1);
 
%% RESOLUTION
 
    % Plan approximé
M_1 = [X_1/b_1 Z_1/b_1 ones(i_1)/b_1];
Res_1 = M_1 \ (-Y_1);
a_1 = Res_1(1); c_1 = Res_1(2) ;d_1 = Res_1(3);
 
%% AFFICHAGE
 
    % Affichage du plan
figure('units','normalized','outerposition',[0 0 1 1])
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '*b'); hold on
 
[Xp_1, Zp_1] = meshgrid(-R_max_1 + x_1 : 0.5 : R_max_1 + x_1 , -R_max_1 + z_1 : 0.5 : R_max_1 + z_1); 
 
Yp_1 = @(Xp_1,Zp_1) (-d_1 - a_1 * Xp_1 - c_1 * Zp_1) / b_1;
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); alpha(0.5);
title('Plan approximé de la tolérance A')
xlabel ('X')
ylabel('Y')
zlabel ('Z')
grid on
 
xx = Xp_1; % Données du plan en X
yy =  Yp_1(Xp_1, Zp_1); % Données du plan en Y
zz =  Zp_1; % Données du plan en Z
oo = [xx(1);yy(1);zz(1)]; % Origine du plan
 
v1 = [xx(1,end);yy(1,end);zz(1,end)]-oo; % Vecteur X du plan
v2 = [xx(end,1);yy(end,1);zz(end,1)]-oo; % Vecteur Z du plan
v3 = cross(v2,v1); % Vecteur Y du plan
axex = v1/norm(v1); % Référentiel axe X du plan
tmpz = v2/norm(v2); % Référentiel temporaire axe Z du plan (au où le plan n'aurait pas des coins orthogonaux)
axey = cross(tmpz,axex); % Référentiel axe Y du plan : Produit vectoriel des vecteurs Z et X
axez = cross(axex,axey); % Référentiel axe Z du plan : Produit vectoriel des vecteurs X et Y
 
quiver3(oo(1),oo(2),oo(3),v1(1),v1(2),v1(3),0,'c','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v3(1)/100,v3(2)/100,v3(3)/100,0,'m','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v2(1),v2(2),v2(3),0,'y','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),axex(1),axex(2),axex(3),0,'r','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axey(1),axey(2),axey(3),0,'g','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axez(1),axez(2),axez(3),0,'b','LineWidth',3);hold off;axis equal;
 
% [axex';axey';axez'] est la matrice de rotation 3x3 à appliquer au plan
% Il faut ramener les données à l'origine avant la rotation autour de l'origine, donc on soustrait l'origine "oo" du plan
pts = [axex';axey';axez']*bsxfun(@minus,[Anneau_exp_1(1,:);Anneau_exp_1(2,:);Anneau_exp_1(3,:)],oo);
 
figure;
plot3(pts(1,:),pts(2,:),pts(3,:),'.');
xlabel ('X');
ylabel('Y');
zlabel ('Z');
grid on;
 
format long;
eye(3)-[axex';axey';axez']

PS : Ton plan est vraiment très près d'être parallèle, la différence entre la matrice identité et la matrice de rotation est de l'ordre de 10^-6.

**oieretxe** · 28/05/2014, 19h14

Merci beaucoup, c'est exactement ce que je cherchais !

Je sais qu'il est extrêmement proche, mais pour le moment je travail avec des points "expérimentaux" générés aléatoirement, et mon projet porte sur des pièces de très grande précision, c'est pour ça que je ne pouvais pas me permettre d'approximé le plan en disant qu'il était parallèle.

Vraiment merci à tout le monde!

Ce forum est vraiment super, on reçoit de l'aide rapidement, et les personnes prennent le temps d'expliquer la solution, c'est vraiment plaisant ! Continuez comme ça ;-)

**oieretxe** · 28/05/2014, 20h28

J'ai réussi à tracer ce que je voulais, mais pourquoi les nouveaux points sont-ils autant éloignés des anciens ? je suppose que c'est le fait de soustraire l'origine "oo" du plan, mais pourquoi doit-on passer par là ? N'est-il pas possible d'appliquer seulement une rotation ?

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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
% On vide la mémoire
close all
clear all
clc
 
%% Plan_rotation.m
% Version 1 : trouver la matrice de rotation pour aligner le plan approximé
 
%% DONNEES
 
    % Nombre de points nominaux des anneaux
N = 360;
 
    % Centre de l'anneau 1
x_1 = 0;
y_1 = 0;
z_1 = 0;
 
    % Rayon max de l'anneau 1
R_max_1 = 15;
 
    % Rayon min de l'anneau 1
R_min_1 = 12;
 
    % Définition de l'anneau nominal 1 et 2
XAnneau_nom_1 = [];
YAnneau_nom_1 = [];
ZAnneau_nom_1 = [];
 
for VTheta = 0:(2*pi)/N:2*pi
    p_1 = rand(1);
    X_1 = x_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * cos(VTheta);
    Y_1 = y_1 * ones(size(VTheta));
    Z_1 = z_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * sin(VTheta);
    XAnneau_nom_1 = [XAnneau_nom_1 X_1];
    YAnneau_nom_1 = [YAnneau_nom_1 Y_1];
    ZAnneau_nom_1 = [ZAnneau_nom_1 Z_1];
 
end
 
    % Définition des anneaux expérimentaux 1 et 2
Ordre = 1 / 10000;
Anneau_exp_1 = [XAnneau_nom_1; YAnneau_nom_1; ZAnneau_nom_1];
 
for i_1 = 1 : N + 1
    p_1 = rand(1);
    Anneau_exp_1(2, i_1) = Anneau_exp_1(2, i_1) + 2 * Ordre * p_1 - Ordre;
end
 
X_1 = (Anneau_exp_1(1, :).');
Y_1 = (Anneau_exp_1(2, :).');
Z_1 = (Anneau_exp_1(3, :).');
b_1 = 10;
i_1 = size(X_1);
 
%% RESOLUTION
 
    % Plan approximé
M_1 = [X_1/b_1 Z_1/b_1 ones(i_1)/b_1];
Res_1 = M_1 \ (-Y_1);
a_1 = Res_1(1); c_1 = Res_1(2) ;d_1 = Res_1(3);
 
%% AFFICHAGE
 
    % Affichage du plan
figure('units','normalized','outerposition',[0 0 1 1])
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '.b'); hold on
 
[Xp_1, Zp_1] = meshgrid(-R_max_1 + x_1 : 0.5 : R_max_1 + x_1 , -R_max_1 + z_1 : 0.5 : R_max_1 + z_1); 
 
Yp_1 = @(Xp_1,Zp_1) (-d_1 - a_1 * Xp_1 - c_1 * Zp_1) / b_1;
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); alpha(0.5);
title('Plan approximé de la tolérance A')
xlabel ('X')
ylabel('Y')
zlabel ('Z')
grid on
 
xx = Xp_1; % Données du plan en X
yy =  Yp_1(Xp_1, Zp_1); % Données du plan en Y
zz =  Zp_1; % Données du plan en Z
oo = [xx(1);yy(1);zz(1)]; % Origine du plan
 
v1 = [xx(1,end);yy(1,end);zz(1,end)]-oo; % Vecteur X du plan
v2 = [xx(end,1);yy(end,1);zz(end,1)]-oo; % Vecteur Z du plan
v3 = cross(v2,v1); % Vecteur Y du plan
axex = v1/norm(v1); % Référentiel axe X du plan
tmpz = v2/norm(v2); % Référentiel temporaire axe Z du plan (au où le plan n'aurait pas des coins orthogonaux)
axey = cross(tmpz,axex); % Référentiel axe Y du plan : Produit vectoriel des vecteurs Z et X
axez = cross(axex,axey); % Référentiel axe Z du plan : Produit vectoriel des vecteurs X et Y
 
quiver3(oo(1),oo(2),oo(3),v1(1),v1(2),v1(3),0,'c','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v3(1)/100,v3(2)/100,v3(3)/100,0,'m','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v2(1),v2(2),v2(3),0,'y','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),axex(1),axex(2),axex(3),0,'r','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axey(1),axey(2),axey(3),0,'g','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axez(1),axez(2),axez(3),0,'b','LineWidth',3);hold off
 
% [axex';axey';axez'] est la matrice de rotation 3x3 à appliquer au plan
% Il faut ramener les données à l'origine avant la rotation autour de l'origine, donc on soustrait l'origine "oo" du plan
rotation = [axex';axey';axez'];
Anneau_exp_2 = rotation*bsxfun(@minus,[Anneau_exp_1(1,:);Anneau_exp_1(2,:);Anneau_exp_1(3,:)],oo);
 
Yp_1_1 = Yp_1(Xp_1, Zp_1);
Plan_2 = rotation*bsxfun(@minus,[Xp_1(:).'; Yp_1_1(:).'; Zp_1(:).'],oo);
 
figure;
plot3(Anneau_exp_2(1,:),Anneau_exp_2(2,:),Anneau_exp_2(3,:),'.m');hold on;
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '.b');
plot3(Plan_2(1,:), Plan_2(2,:), Plan_2(3,:), '.m');
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); alpha(0.5);
xlabel ('X');
ylabel('Y');
zlabel ('Z');
grid on;

**ERCO503** · 28/05/2014, 21h18

Après avoir soustrait par l'origine et effectué la rotation, tu peux additionner l'origine pour revenir à l'endroit où tu étais.

J'ai aussi modifié l'origine, probablement que le centre moyen de tes points t'intéresse plus (deux possibilités).

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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
% On vide la mémoire
close all
clear all
clc
 
%% Plan_rotation.m
% Version 1 : trouver la matrice de rotation pour aligner le plan approximé
 
%% DONNEES
 
    % Nombre de points nominaux des anneaux
N = 360;
 
    % Centre de l'anneau 1
x_1 = 0;
y_1 = 0;
z_1 = 0;
 
    % Rayon max de l'anneau 1
R_max_1 = 15;
 
    % Rayon min de l'anneau 1
R_min_1 = 12;
 
    % Définition de l'anneau nominal 1 et 2
XAnneau_nom_1 = [];
YAnneau_nom_1 = [];
ZAnneau_nom_1 = [];
 
for VTheta = 0:(2*pi)/N:2*pi
    p_1 = rand(1);
    X_1 = x_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * cos(VTheta);
    Y_1 = y_1 * ones(size(VTheta));
    Z_1 = z_1 + (R_min_1 + (R_max_1-R_min_1) * p_1) * sin(VTheta);
    XAnneau_nom_1 = [XAnneau_nom_1 X_1];
    YAnneau_nom_1 = [YAnneau_nom_1 Y_1];
    ZAnneau_nom_1 = [ZAnneau_nom_1 Z_1];
 
end
 
    % Définition des anneaux expérimentaux 1 et 2
Ordre = 1 / 10000;
Anneau_exp_1 = [XAnneau_nom_1; YAnneau_nom_1; ZAnneau_nom_1];
 
for i_1 = 1 : N + 1
    p_1 = rand(1);
    Anneau_exp_1(2, i_1) = Anneau_exp_1(2, i_1) + 2 * Ordre * p_1 - Ordre;
end
 
X_1 = (Anneau_exp_1(1, :).');
Y_1 = (Anneau_exp_1(2, :).');
Z_1 = (Anneau_exp_1(3, :).');
b_1 = 10;
i_1 = size(X_1);
 
%% RESOLUTION
 
    % Plan approximé
M_1 = [X_1/b_1 Z_1/b_1 ones(i_1)/b_1];
Res_1 = M_1 \ (-Y_1);
a_1 = Res_1(1); c_1 = Res_1(2) ;d_1 = Res_1(3);
 
%% AFFICHAGE
 
    % Affichage du plan
figure('units','normalized','outerposition',[0 0 1 1])
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '.b'); hold on
 
[Xp_1, Zp_1] = meshgrid(-R_max_1 + x_1 : 0.5 : R_max_1 + x_1 , -R_max_1 + z_1 : 0.5 : R_max_1 + z_1); 
 
Yp_1 = @(Xp_1,Zp_1) (-d_1 - a_1 * Xp_1 - c_1 * Zp_1) / b_1;
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); alpha(0.5);
title('Plan approximé de la tolérance A')
xlabel ('X')
ylabel('Y')
zlabel ('Z')
grid on
 
xx = Xp_1; % Données du plan en X
yy =  Yp_1(Xp_1, Zp_1); % Données du plan en Y
zz =  Zp_1; % Données du plan en Z
oo = [xx(1);yy(1);zz(1)]; % Origine du plan
 
v1 = [xx(1,end);yy(1,end);zz(1,end)]-oo; % Vecteur X du plan
v2 = [xx(end,1);yy(end,1);zz(end,1)]-oo; % Vecteur Z du plan
v3 = cross(v2,v1); % Vecteur Y du plan
axex = v1/norm(v1); % Référentiel axe X du plan
tmpz = v2/norm(v2); % Référentiel temporaire axe Z du plan (au où le plan n'aurait pas des coins orthogonaux)
axey = cross(tmpz,axex); % Référentiel axe Y du plan : Produit vectoriel des vecteurs Z et X
axez = cross(axex,axey); % Référentiel axe Z du plan : Produit vectoriel des vecteurs X et Y
 
quiver3(oo(1),oo(2),oo(3),v1(1),v1(2),v1(3),0,'c','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v3(1)/100,v3(2)/100,v3(3)/100,0,'m','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),v2(1),v2(2),v2(3),0,'y','LineWidth',2);
quiver3(oo(1),oo(2),oo(3),axex(1),axex(2),axex(3),0,'r','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axey(1),axey(2),axey(3),0,'g','LineWidth',3);
quiver3(oo(1),oo(2),oo(3),axez(1),axez(2),axez(3),0,'b','LineWidth',3);hold off
 
% Autres origines possibles
oo = mean(Anneau_exp_1,2); % Au lieu d'utiliser le coin, utiliser la valeur moyenne des points pour une rotation autour du centre
% oo = mean([xx(:)';yy(:)';zz(:)'],2);  % Au lieu d'utiliser le coin, utiliser la valeur moyenne des données du plan pour une rotation autour du centre
 
% [axex';axey';axez'] est la matrice de rotation 3x3 à appliquer au plan
% Il faut ramener les données à l'origine avant la rotation autour de l'origine, donc on soustrait l'origine "oo" du plan
rotation = [axex';axey';axez'];
Anneau_exp_2 = bsxfun(@plus,rotation*bsxfun(@minus,[Anneau_exp_1(1,:);Anneau_exp_1(2,:);Anneau_exp_1(3,:)],oo),oo);
 
Yp_1_1 = Yp_1(Xp_1, Zp_1);
Plan_2 = bsxfun(@plus,rotation*bsxfun(@minus,[Xp_1(:).'; Yp_1_1(:).'; Zp_1(:).'],oo),oo);
 
figure;
plot3(Anneau_exp_2(1,:),Anneau_exp_2(2,:),Anneau_exp_2(3,:),'.m');hold on;
plot3(Anneau_exp_1(1,:), Anneau_exp_1(2,:), Anneau_exp_1(3,:), '.b');
plot3(Plan_2(1,:), Plan_2(2,:), Plan_2(3,:), '.m');
surf(Xp_1 ,Yp_1(Xp_1, Zp_1), Zp_1); alpha(0.5);
xlabel ('X');
ylabel('Y');
zlabel ('Z');
grid on;

**oieretxe** · 28/05/2014, 23h31

Parfait ça marche à merveille !

Un très grand merci

, je vais pouvoir avancer maintenant !

**oieretxe** · 24/07/2014, 17h11

Re-bonjour,

je reviens sur ce sujet après tout ce temps car j'ai de nouveau un blocage. Grâce à ton code ERCO503, je réussi parfaitement à appliquer la matrice de rotation à mon plan afin que ces valeurs en "y" soit constante. J'ai légèrement modifié mon programme et maintenant le plan que je trouve doit être repositionné afin que ces valeurs en "x" soit constantes, et là je n'y arrive pas... J'ai essayé de définir le référentiel autrement mais je n'y arrive pas, je fais ça un peu au "pif" ^^ si tu peux m'expliquer le raisonnement à suivre, c'est avec grand plaisir que je t'écouterais attentivement !

Merci d'avance

**ERCO503** · 24/07/2014, 18h49

Il faudrait voir tes modifications.

trouver la matrice de rotation à appliquer à un plan

MATLAB

Discussions similaires

Partager

Partager