How to correctly calculate the angle between two planes?

7 visualizaciones (últimos 30 días)
Ali
Ali el 16 de Jul. de 2020
Comentada: Ali el 17 de Jul. de 2020
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);

Respuesta aceptada

Ali
Ali el 16 de Jul. de 2020
Editada: Ali el 16 de Jul. de 2020
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.

Más respuestas (1)

Matt J
Matt J el 16 de Jul. de 2020
Editada: Matt J el 16 de Jul. de 2020
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.
  12 comentarios
Matt J
Matt J el 16 de Jul. de 2020
Editada: Matt J el 16 de Jul. de 2020
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
Ali el 17 de Jul. de 2020
"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. :-)

Iniciar sesión para comentar.

Categorías

Más información sobre Point Cloud Processing en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by