9 views (last 30 days)

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
on 11 Sep 2020

Edited: Bruno Luong
on 11 Sep 2020

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

Bruno Luong
on 23 Sep 2020

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.

KSSV
on 10 Sep 2020

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));

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.