Z-angle of a fit-plane having 4 points

13 views (last 30 days)
Hi,
how can I calculate the z-angle of the plane that best fits 4 points having their coordinates?
Then, I would like to plot it.
Moreover, how can I calculate the z-angle of the plane passing through 3 points having their coordinates?
Then, I would like to plot it.

Accepted Answer

William Rose
William Rose on 9 Sep 2021
Here is code to find the unit normal to a plane defined by three points, and then the angle with the z-axis.. In the example below, the points are chosen at random. The gamma=min(..) statement assures that the answer will be in the range 0 <= gamma <= pi/2.
p0=rand(3,1); p1=rand(3,1); p2=rand(3,1); %three random points in the unit cube
n=cross(p1-p0,p2-p0); %cross product of 2 vectors that lie in the plane will be normal to the plane
n=n/norm(n); %make the normal vector have length=1
gamma=min(acos(n(3)),pi-acos(n(3))); %angle with z-axis
fprintf('Angle with Z-axis=%.3f radians\n',gamma)
Angle with Z-axis=0.890 radians
I asked whether you know the x,y,z coordinates exactly, because the answer affects you choice of what to minimize to get the "best fit". A. If you know the x and y values exactly, andt there is measurement error or of other uncertainty along z, then the best-fit is the fit that minimizes the sum squared error in the z-direction. This is ordinary least squares, or standard multiple linear regression. B. If you have measurement errors or other uncertinty on all three axes, then the best fit is the fit that minimizes the orthogonal sum squared error from the plane, also known as orthogonal or total least squares fitting.
A. Ordinary least squares=standard multiple regresion: If you have access to the Statistics toolbox, use fitlm(). If you don;t, you can do it the old fashioned way. I will put code for both in a separate post.
B. Orthogonal (total) least squares: use pca() to find the principal components (PCs). There will be thee PCs. The lthird (last) PC, which accounts for the smallest portion of the variance, is the normal to the best fit plane. I will put code in a separate post.
Once you have a unit normal vector n, from method A or B, then the angle with the Z axis is given as before:
gamma=min(acos(n(3)),pi-acos(n(3))); %angle with z-axis
Code for A and B in separate post.

More Answers (5)

William Rose
William Rose on 9 Sep 2021
@Gaetano Pavone, YOu also requested help for how to plot the plane, the points, the normal vector. See attached. This is for the case of a plane through three points. If you rest your pointr over the live plot, icons appear at the top right of the plot. Click the icon that looks like a cube with an arrow around it. Then you may click and drag to rotate the plot in three dimensions.
%fitPlane0.m WCRose 20210908
%Fit a plane through 3 points.
%For Gaetano Pavone.
clear;
p0=rand(3,1); p1=rand(3,1); p2=rand(3,1); %three random points in unit cube
n=cross(p1-p0,p2-p0); %normal to the plane
n=n/norm(n); %unit normal
gamma=min(acos(n(3)),pi-acos(n(3))); %angle from z-axis (radians)
%Plot best fit plane
P=[p0';p1';p2']; %matrix of points; each point is a row, columns 1,2,3=x,y,z
xv=linspace(min(P(:,1)),max(P(:,1)),11); %vector of x values for mesh
yv=linspace(min(P(:,2)),max(P(:,2)),11); %vector of y values for mesh
[X,Y]=meshgrid(xv,yv);
Z=-(n(1)*(X-p0(1))+n(2)*(Y-p0(2)))/n(3)+p0(3); %z-values for mesh
figure;
mesh(X,Y,Z,'FaceAlpha',0.4); %draw mesh with semi-transparent faces
xlabel('X'); ylabel('Y'); zlabel('Z');
%Plot the points as black circles
hold on;
scatter3(P(:,1),P(:,2),P(:,3),'MarkerFaceColor','k');
axis equal;
%Plot the normal vector (red) and z-axis vector (black)
arrowbase=(min(P)+max(P))/2; %base in center of the mesh
nup=n*sign(n(3)); %normal vector that alway points up
quiver3(arrowbase(1),arrowbase(2),arrowbase(3),nup(1),nup(2),nup(3),0,'r');
quiver3(arrowbase(1),arrowbase(2),arrowbase(3),0,0,1,0,'k');

William Rose
William Rose on 8 Sep 2021
Edited: William Rose on 8 Sep 2021
Is each data point a set of 3 coordinates? If the answer is yes, then do you know the x and y values exactly? In that case, you would minimize the sum squared errors between the measured z values and the corresponding z-values on the plane. Or is there uncertainty or measurement error in all 3 coordinates? In that case, you should minimize the sum of the normal distances between each point and the plane. The solution for the plane is different for these cases.
You asked for "the z-angle of the plane". Do you mean the angle between the normal to the plane and the z-axis?
  1 Comment
Gaetano Pavone
Gaetano Pavone on 8 Sep 2021
I have for each of the 4 points a set of exact 3 coordinates, and I need the angle between the normal to the plane and the z-axis

Sign in to comment.


Matt J
Matt J on 8 Sep 2021
Edited: Matt J on 8 Sep 2021
You can fit a plane to the points using planarFit() in this downloadable package:
fitObj=planarFit(xyz);
Calculating the angle is then simply,
zAngle=acos(fitObj.normal(3));
  1 Comment
Matt J
Matt J on 9 Sep 2021
Then, I would like to plot it.
To plot, use the plot() method of the fit object.
plot(fitObj)
The procedure remains the same whether for 4 or 3 points.

Sign in to comment.


William Rose
William Rose on 10 Sep 2021
Here's a script (attached) that fits a plane to 4 or more points. It fits by ordinary least squares, which means it minimizes the sum squared error in the z-direction. (A total least squares fit minimizes the sum squared perpendicular-to-the-plane errors.) This script does it two ways, and they give the same exact results. Method 1 uses fitlm(), which is in the Statistics & Machine Learning Toolbox. Method 2 does not require the Statistics & M.L. toolbox. Method 2 just uses matrix math. The script displays the unit normal vector, the z-axis intercept, the root mean square error in the z-direction, and the root mean square error in the perpendicular direction, for ech of the two methods.
You can adapt the plotting code I posted earlier to display these results.

William Rose
William Rose on 10 Sep 2021
In my previous post, I supplied a script (fitPlane3a.m) that finds the plane that is the best fit to 4 or more points. The plane minimizes the vertical (z-drection) squared distances of the points from the plane, i.e. ordinary least squares (OLS) minimization. The script did it two ways: with and without using the Statistics and Machine Learning toolbox.
Now I am posting an updated version: fitPlane3.m. This finds an additional, different plane: the plane that minimizes the orthogonal squared distances of the points from the plane, i.e. total least squares (TLS) minimization. The script does it two ways: with and without using the Statistics and Machine Learning toolbox. The first TLS method uses pca(), principal component analysis, which is in the Stats toolbox. The second TLS method uses eigenvector analysis of the covariance matrix, which does not require any toolbox.
The first part of the script creates m data points for fitting: matrix X and vector z. X is m-by-2, z is m-by-1. Replace these with the data you wish to fit.
An example of console output is below.
>> fitPlane3
Points to fit, each point is a column.
11.0000 11.0000 10.0000 9.0000 9.0000 10.0000
9.0000 10.0000 11.0000 11.0000 10.0000 9.0000
-14.9112 -16.1103 -15.9203 -15.0552 -13.8149 -14.1331
Method nHat zIntercept RMSEz RMSEtot
OLS 1: [0.5916, 0.5707, 0.5694] 5.42219 0.11193 0.06373
OLS 2: [0.5916, 0.5707, 0.5694] 5.42219 0.11193 0.06373
TLS 1: [0.5939, 0.5731, 0.5647] 5.67192 0.11239 0.06347
TLS 2: [0.5939, 0.5731, 0.5647] 5.67192 0.11239 0.06347
The console output shows:
A. The two OLS methods give identical results.
B. The two TLS methods give identical results.
C. The plane found by OLS is slightly different than the plane found by TLS.
D. The OLS plane has smaller RMS error in the z-direction than the TLS plane.
E. The TLS plane has smaller total (orthogonal) RMS error than the OLS plane.
  1 Comment
William Rose
William Rose on 10 Sep 2021
In the console output above, nHat is the unit normal. nHat and zIntercept completely define a plane. RMSEz and RMSEtot are measures of goodness of fit, where smaller root mean square error is better. RMSEz is the average error in the vertical direction. RMSEtot is the average perpendicular distance of each point from the plane.

Sign in to comment.

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!