Incorrect angle detected between two planes

Question

I want to calculate the angle between 2 planes, Reference plane and Plane1. When I feed the X,Y,Z co-ordinates of pointCloud to the function plane_fit.m (by Kevin Mattheus Moerman), I get the output coefficients:

(MATLAB)

reference_plane_coeff: [-0.13766204 -0.070385590 130.69409]

Plane1_coeff: [0.0044337390 -0.0013548643 95.890228]

Next, I find the intersection of both planes, separately on the XZ plane and get a line equation; ref_line_XZ and plane1_line_XZ respectively. For this, I make the second coefficient 0. (Is this right?)

Aref = reference_plane_coeff(1);
Cref = reference_plane_coeff(3);
ref_line_XZ = [Aref Cref];
Arun = Plane1_coeff(1);
Crun = Plane1_coeff(3);
plane1_line_XZ = [Arun Crun];
angle_XZ = acosd(  dot(ref_line_XZ,plane1_line_XZ ) / (norm(ref_line_XZ) * norm(plane1_line_XZ ))        )

I get the angle_XZ value as 0.0685 degrees.

When I plot these planes on a graph and view it, the angle seems to be much more than 0.0685 degrees. I'm talking about the angle made by the two lines upon intersection of both planes with the XZ plane.

What am I doing wrong?

Also, when I tried to find angle between its normals, using:

angle_beta_deg = acosd(  dot(reference_plane_coeff,Plane1_coeff) / (norm(reference_plane_coeff) *     norm(Plane1_coeff))    )

I got the angle as 0.0713.

On visual inspection of both planes' plots and manually calculating from the plot, angle_XZ should be around 9 degrees.

I tried with another code, Plane Fitting and Normal Calculation, by Dan Couture.

With the same X, Y and Z co-ordinates as input, I got the coefficients (of normal) as

n_reference = [0.1512 -0.0327 0.9880]

n_Plane1 = [0.0005 0.0156 0.9999]

Using the same equation:

angle_XZ = acosd(  dot(ref_line_XZ,plane1_line_XZ ) / (norm(ref_line_XZ) * norm(plane1_line_XZ ))        )

I get angle as 8.6731 degrees.

But why is there difference in the Z-axis? The function plane_fit by Kevin Moerman appears to fit plane properly but does not give the correct angle. Maybe I am not using the output of this function in the right way.

I want to use the one giving highest accuracy.

Could you suggest any other more accurate function?

Or explain why there is a difference in my angles?

I can't visualize the svd method of plane fitting as I am not clear about it. I have an intuition of the least squares regression method.

I want a method which tells me the orientation of my checkerboard (reference) and a flat (almost flat.....similar to a wheel) circular object (Plane1) as accurately as possible, with the X, Y and Z axes.

Link for the NX3 matrix (X_Y_X) on which plane is to be fit:

Reference plane: www.dropbox.com/s/pxtdl1ywuje7og9/X_Y_Z%28ref%29.mat?dl=0

Plane1: www.dropbox.com/s/yfyout8ppz9ai4o/X_Y_Z%28Plane1%29.mat?dl=0

This is the code I'm using to plot:

x_temp = X_Y_Z(:,1);
y_temp = X_Y_Z(:,2);
z_temp = X_Y_Z(:,3);

%^^^^^^^^^^DanCouture's^^^^^^^^^^^%
n = fitNormal(X_Y_Z);       
A = n(1);
B = n(2);
C = n(3);
%^^^^^^^^^^^^^^^^^^^^^%
%*****************Moerman's************************%
% [A,B,C]=plane_fit(x_temp,y_temp,z_temp);
%*****************************************%
%&&&&&&&&&&&Adrien's&&&&&&&&&&&&%
% [n,V,p] = affine_fit(X_Y_Z);
% A = n(1);
% B = n(2);
% C = n(3);
%&&&&&&&&&&&&&&&&&&&&&&&%

[X,Y]=meshgrid(linspace(min(x_temp),max(x_temp),20),linspace(min(y_temp),max(y_temp),20));

Z=(A*X)+(B*Y)+C;
figure;
showPointCloud(pointCloudDisp, J1, 'VerticalAxis', 'Y',...
    'VerticalAxisDir', 'Down' );
hold on; grid on;
surf(X,Y,Z,'FaceColor','r'); alpha(0.5);
xlabel('X'); 
ylabel('Y');
zlabel('Z'); `

When I give C = -1, instead of n(3), it too doesnt work.

Answer 1

The question of adjusting a line or a plane to a set of points depends on how you are looking at it.

Basically, you can fit by minimizing the "vertical" distance between the line/plane and the cloud of points or by minimizing the normal distance (perpendicular to the plane).

In the first case you have a linear system to solve, in the second case it becomes an eigenvalue problem to solve.

the code of Moerman gives A,B,C such that:

Z=(A*X)+(B*Y)+C;

Therefore the normal vector is [A B -1] and is not normalized. I do not think you used the output of Moerman's code correctly.

An alternative to Moerman's code is Plane fit which I submitted to the Matlab file exchange and which has the advantage of giving directly the normalized vector normal to the plane.

The code by Dan Couture (as I understand it in 2 min) seems to try to fit the data by finding plane that minimizes the x distance then a plane that minimizes the y distance then a plane that minimizes the z distance and then selects the best of the 3 fits. This approach is more oriented towards finding the plane that provides the best linear explanation of a coordinate as a function of two others. This approach is quite popular because of its simplicity, but you should actually first decide what you want to do and then choose the appropriate code.

Hope this helps.

Adrien.

update: I just added an example in my Matlab FileExchange contribution. This should help you: (http://www.mathworks.com/matlabcentral/fileexchange/43305-plane-fit "Plane fit") .