How to correctly calculate the angle between two planes?

I have 2 point cloud files from a buildings gable rooftop (.ply attached) and after fitting a plane to each and getting the normal vectors in plane model object I'd like to calculate the angle between these two, but the result is different depending on which point cloud is input first!

One time I get 2.68555450304481 and if I change the order of input I get 0.456276904837795 !!

What am I missing? Also if run the code a couple of more times I get some changes in fractions like 0.45xxx turns to 0.46xxx !!!

% input point cloud file gui
[FileName,PathName] = uigetfile({'*.pcd;*.ply;',...
    'Point Cloud Files (*.pcd,*.ply)';
    '*.pcd','Point Cloud library files (*.pcd)'; ...
    '*.ply','Polygon Mesh Point Cloud files (*.ply)';
    '*.*',  'All Files (*.*)'}, ...
    'Select a Point Cloud File');
ptCloud=pcread([PathName,FileName]);

[FileName,PathName] = uigetfile({'*.pcd;*.ply;',...
    'Point Cloud Files (*.pcd,*.ply)';
    '*.pcd','Point Cloud library files (*.pcd)'; ...
    '*.ply','Polygon Mesh Point Cloud files (*.ply)';
    '*.*',  'All Files (*.*)'}, ...
    'Select a Point Cloud File');
ptCloud2=pcread([PathName,FileName]);

maxDistance = 0.02;
[model1,inlierIndices1,outlierIndices1] = pcfitplane(ptCloud,maxDistance);
[model2,inlierIndices2,outlierIndices2] = pcfitplane(ptCloud2,maxDistance);

% calculate the angle between two planes
angle = atan2(norm(cross(model1.Normal,model2.Normal)),dot(model1.Normal,model2.Normal));

% another method but the same results!!!
A1 = model1.Parameters(1);B1 = model1.Parameters(2);C1 = model1.Parameters(3);
A2 = model2.Parameters(1);B2 = model2.Parameters(2);C2 = model2.Parameters(3);
dotproduct = (A1*A2) + (B1*B2) + (C1*C2);
angle2 = acos(dotproduct);

Answer 1

[Ali]

I found out the reason, I have two planes in diagonal shape / and \ so the difference is because the angle is calculated clockwise therefore the angle for the order ( \./ ) is 0.45 radians and 2.68 radians for the order of ( /.\) input.

I misunderstood the calculation and assumed no matter what the order, they somehow merge in XYZ space and always turn out like (/.\)!

Thanks everyone for their answers.

Answer 2

[Matt J]

From the documentation for pcfitplane:

This function uses the M-estimator SAmple Consensus (MSAC) algorithm to find the plane. The MSAC algorithm is a variant of the RANdom SAmple Consensus (RANSAC) algorithm.

The fact that the algorithm uses random sampling would explain why the result changes each time. The RANSAC algorithm is indeed a good choice if your point cloud data has significant outliers. If not, you can use my attached plane fitting routine instead, which will give non-stochastic results.

[Matt J]

Your code won't work on my data because it's a 3x1 matrix,

If your point cloud is just a 3x1 matrix, then it is not a "cloud". It's just a single point.

[Ali]

Nevertheless I found the solution, thank you for your time.

[Ali]

@Matt_J would you happen to know why am I getting negative values for the normal vector of some plane object? And what does it mean mathematically? For some it's [-a -b c], for some other [ a -b c], for some other [ -a b c] and finally there are some with all positive values [ a b c]. (I don't care about value c here, just a and b).

[Matt J]

If you attach a .mat file with the data as a 3xN matrix, I could examine it.

[Ali]

I don't know how to save my data as .mat file, perhaps if you give some directions on that, but I have attached the data in ply format along my code, after running the code you can find model1 and model2 in workspace and examine them.

[Matt J]

@Ali, to save one or several matlab variables to a .mat file, use the save command

https://www.mathworks.com/help/matlab/ref/save.html#btppzj8-1

Providing the data as .ply is not helpful for those who don't have the Computer Vision Toolbox.

[Ali]

@Bruno thank you for your explanation, I guess you haven't noticed the accepted answer yet. Right now I have another question that I asked in the comments and totally unrelated.

[Ali]

@Matt_J my bad! Didn't realize that, here you go, points in 3xN .mat.

[Matt J]

Ali,

The data looks puzzling, firstly because it looks like the units of distance you are using are micrometers. That seems like a very strange choice of units for an architectural building design.

>> min(pointTrasposed,[],2)

ans =

   1.0e+06 *

    1.5694
    5.1862
    0.0000

>> max(pointTrasposed,[],2)

ans =

   1.0e+06 *

    1.5694
    5.1863
    0.0000

Secondly, since the z-component is essentially zero throughout the data set, the plane normal should almost trivially be [0,0,1]. Since this is not realistic for any gabled roof, I am assuming there is a mistake in the units for the z-coordinate data. Thirdly, your .mat file only contained one 3xN data set, whereas I think you should have a separate data set for each plane.

[Ali]

@Matt_J I'm not sure how does a 3xN matrix works for XYZ points or if something went wrong during the process of transposing! But in the original point cloud object there is a property called 'Locations' that is an Nx3 matrix and each column corresponds to X Y Z coordinates, so puzzled with your comment I did a min on Z column and it returned the following:

>> min(points(:,3))
ans =
  13.861000061035156
  

As for the unit of distance in 'micrometers' I literally have no clue about it, if it's wrongly set I'd appreciate any extra tips on how to fix it.

Also the points in Nx3 matrix is attached, one for the merged segments that results in a complete roof point cloud, and two others as the two sides of gable.

Thank you for your time and effort. :)

[Matt J]

Never mind, I figured it out. Anyway, here is the code I used and, as the plots below show, planefit gives a very good fit to the roof planes. Keep in mind, of course, that the normals [a,b,c] are invariant to a sign change [-a,-b,-c] .

load points

P=pointTrasposed-mean(pointTrasposed,2);
Pleft=P(:,P(2,:)<0,:);
Pright=P(:,P(2,:)>0);
[xl,yl,zl]=deal(Pleft(1,:),Pleft(2,:),Pleft(3,:));
[xr,yr,zr]=deal(Pright(1,:),Pright(2,:),Pright(3,:));


[al,bl,cl,dl]=planefit([xl;yl;zl]); %left plane fit
[ar,br,cr,dr]=planefit([xr;yr;zr]); %right plane fit



hold on
sl=scatter3(xl,yl,zl,'MarkerFaceColor','r','MarkerEdgeColor','none');
sr=scatter3(xr,yr,zr,'MarkerFaceColor','b','MarkerEdgeColor','none');
fl=fimplicit3(@(x,y,z) al*x+bl*y+cl*z-dl,'FaceColor','r','FaceAlpha',0.2,'EdgeColor','none');
fr=fimplicit3(@(x,y,z) ar*x+br*y+cr*z-dr,'FaceColor','b','FaceAlpha',0.2,'EdgeColor','none');
xlabel x, ylabel y
hold off

[Ali]

"Keep in mind, of course, that the normals [a,b,c] are invariant to a sign change [-a,-b,-c]" Thank you @Matt_J, great work and I really appreciate your time and help, respect. :-)