How to calculate the angle between two sets of scatter plots?

6 次查看（过去 30 天）
Joe Pashley
Joe Pashley 2020-9-10
编辑： Bruno Luong 2020-9-23
Hello,
I have two sets of measurement data in the format of 5 XYZ data points. These points are taken by a coordiante measurment machine (CMM).
The CMM is measuring a square face which is being manouvered by a rotary stage through azimuth and elevation angles.
I would like to calculate the plane of best fit of both of these datasets, and then calculate the angles between these two planes. Ideally plotting both planes and the scatter data so the angles can be visualised.
Im finding this complicated as the plane is not only adjsuted in one angle but both an azimuth and elevation angle.
I have tried the planefit function by Adrien Leygue and using the cftool, but struggled to find an end solution.
Example data below:
o_data = 3.11009098091043 0.0161338808382554 -3.50871929189684
34.1536562864604 0.446152781388255 -3.50922770378684
34.1557207645004 0.220196060488255 -30.0608994319768
2.77355697927043 -0.211958502351745 -29.5764317994968
6.0245043967304 0.0441464660382553 -20.2055638418268
m_data = 5.8062542992593 14.1534300854169 -20.0235623701306
30.4060090931593 9.11456283141687 -2.59783606773058
30.4016325470993 -18.8937153876031 -10.9017247126206
7.4045188446593 -14.0475206206231 -27.1514875576406
17.5591295448693 -4.02446512820313 -16.3715286064806
Any advice on how to achieve this would be greatly appreciated.
Thanks

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2020-9-11
编辑：Bruno Luong 2020-9-11
clear
o_data = [3.11009098091043 0.0161338808382554 -3.50871929189684;
34.1536562864604 0.446152781388255 -3.50922770378684;
34.1557207645004 0.220196060488255 -30.0608994319768;
2.77355697927043 -0.211958502351745 -29.5764317994968;
6.0245043967304 0.0441464660382553 -20.2055638418268 ];
m_data = [5.8062542992593 14.1534300854169 -20.0235623701306;
30.4060090931593 9.11456283141687 -2.59783606773058;
30.4016325470993 -18.8937153876031 -10.9017247126206;
7.4045188446593 -14.0475206206231 -27.1514875576406;
17.5591295448693 -4.02446512820313 -16.3715286064806 ];
data = {o_data; m_data};
n = length(data);
% Best plane fit to data
xyzc = zeros(3,1,n);
N = zeros(3,1,n);
X = zeros(3,2,n);
for k=1:n
xyz = data{k};
meanxyc = mean(xyz,1);
xyzc(:,:,k) = meanxyc';
[~,~,V] = svd(xyz-meanxyc);
N(:,:,k) = V(:,3);
X(:,:,k) = V(:,1:2);
end
% Compute intersection point c and intersection direction C
C = cross(N(:,:,1),N(:,:,2));
C = C/norm(C);
T = zeros(3,1,n);
for k=1:n
Yk = C'*X(:,:,k);
Tk = X(:,:,k)*[Yk(2);-Yk(1)];
T(:,:,k) = Tk/norm(Tk);
end
M = [T(:,:,1),-T(:,:,2)];
b = xyzc(:,:,2)-xyzc(:,:,1);
x = M \ b;
c = mean(xyzc + T.*reshape(x,1,1,2),3);
% Compute angle
% NOTE: Angle between two N vectors gives the same numerical result, meaning
% abs(theta) = asin(norm(cross(N(:,:,1),N(:,:,2)))
% but they are not directly linked to the plane itself, so difficult to illustrate on graph
W1 = cross(C,T(:,:,1));
theta = asin(dot(T(:,:,2),W1));
% Graphic
AllP = cell2mat(data);
close all
hold on
Cs = c + 20*C.*[-1,1];
plot3(Cs(1,:), Cs(2,:), Cs(3,:), '-.b')
for k=1:n
xyz = data{k};
h = plot3(xyz(:,1), xyz(:,2), xyz(:,3), 'o', 'MarkerSize', 5, 'LineWidth',2);
color = get(h, 'Color');
P = (AllP-xyzc(:,:,k)')*X(:,:,k);
[Pmin,Pmax] = bounds(P,1);
Q = [Pmin;
Pmax];
R = [Q([1 2 2 1 1],1), Q([1 1 2 2 1],2)];
R = xyzc(:,:,k)' + R*X(:,:,k)';
fill3(R(:,1), R(:,2), R(:,3), color, 'FaceAlpha', 0.3);
Tk = 20*T(:,:,k);
l = [c, c+Tk];
plot3(l(1,:), l(2,:), l(3,:), 'Color', color,'Linewidth',2);
end
% Plot arc
tt = linspace(0,theta,33);
arc = c + 10*[T(:,:,1),W1]*[cos(tt); sin(tt)];
plot3(arc(1,:),arc(2,:),arc(3,:),'k','Linewidth',1);
a = arc(:,ceil(end/2));
text(a(1),a(2),a(3),sprintf('%2.1f deg',rad2deg(abs(theta))))
view(3)
axis equal
  2 个评论
Joe Pashley
Joe Pashley 2020-9-23
编辑：Joe Pashley 2020-9-23
Thank you so much for your work on this.
Im going to take my time to look through your code and try and understand it step by step.
The only other thing I was trying to achieve was to split the compound angle into its seperate azimuth and elevation angles. For example: The green plane was originally in the same position as the red plane, but has been rotated clockwise (azimuth angle) and then tilted upwards (elevation angle). This results in the 76.9 deg value between the two planes. However how much has it rotated in azimuth and elevation. The values should be very close too: 57 deg azimuth, 73 deg elevation.
Thanks
Joe
Bruno Luong
Bruno Luong 2020-9-23
编辑：Bruno Luong 2020-9-23
A generic rotation has 3 degrees of freedom not 2. And furthermore we don't know if one/both of the planes rotate about their normal. If they rotate the information is not retrieve after fitting the plane thought points.
To convert to whatever convention you use follow the conversion of Wikipedia page
If the points can be matched by the rotation, then this is completely different problem, and you really fail to tell us this crucial information.

请先登录，再进行评论。

更多回答（1 个）

KSSV
KSSV 2020-9-10
o = [ 3.11009098091043 0.0161338808382554 -3.50871929189684
34.1536562864604 0.446152781388255 -3.50922770378684
34.1557207645004 0.220196060488255 -30.0608994319768
2.77355697927043 -0.211958502351745 -29.5764317994968
6.0245043967304 0.0441464660382553 -20.2055638418268 ] ;
m = [5.8062542992593 14.1534300854169 -20.0235623701306
30.4060090931593 9.11456283141687 -2.59783606773058
30.4016325470993 -18.8937153876031 -10.9017247126206
7.4045188446593 -14.0475206206231 -27.1514875576406
17.5591295448693 -4.02446512820313 -16.3715286064806] ;
% Equation of first plane A1*x+B1*y+C1*z = 0
P1 = o(1,:) ;
P2 = o(2,:) ;
P3 = o(3,:) ;
v1=cross(P1-P2,P1-P3) ;
% Equation of second plane A2*x+B2*y+C2*z = 0
Q1 = m(1,:) ;
Q2 = m(2,:) ;
Q3 = m(3,:) ;
v2=cross(Q1-Q2,Q1-Q3) ;
% find angle
ang = atan2(norm(cross(v1,v2)),dot(v1,v2));

类别

在 Help Center 和 File Exchange 中查找有关 Delaunay Triangulation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by