how to fit a plane to my point cloud data and rotate it so that the plane is parallel to the x-y plane.

Hi All,

I have a set of approximately 1 million (x,y,z) points collected through photogrammetric tecniques of seabeds. The seafloor is characterized by a slope. I would like to:

1. fit a plane to my data;
2. rotate my 3d point cloud so that the fitted plane is parallel to xy-plane;
3. export the new 3d point cloud into .txt or .xyz file

i fit the plane in this way:

%%Loading data
%data = load('Coralligeno_outcrop.txt');
data = load('./input.dat');
x = data(:,1);
y = data(:,2);
z = data(:,3);
clear data;

%%Finds the mean plan
% Combine X,Y,Z into one array
XYZ = [x,y,z];

% column means. This will allow us to do the fit properly with no
% constant term needed. It is necessary to do it this way, as
% otherwise the fit would not be properly scale independent
cm = mean(XYZ,1);

% subtract off the column means
XYZ0 = bsxfun(@minus,XYZ,cm);

% The "regression" as a planar fit is now accomplished by SVD.
% This presumes errors in all three variables. In fact, it makes
% presumptions that the noise variance is the same for all the
% variables. Be very careful, as this fact is built into the model.
% If your goal is merely to fit z(x,y), where x and y were known
% and only z had errors in it, then this is the wrong way
% to do the fit.
[U,S,V] = svd(XYZ0,0);

% The singular values are ordered in decreasing order for svd.
% The vector corresponding to the zero singular value tells us
% the direction of the normal vector to the plane. Note that if
% your data actually fell on a straight line, this will be a problem
% as then there are two vectors normal to your data, so no plane fit.
% LOOK at the values on the diagonal of S. If the last one is NOT
% essentially small compared to the others, then you have a problem
% here. (If it is numerically zero, then the points fell exactly in
% a plane, with no noise.)
diag(S)

% Assuming that S(3,3) was small compared to S(1,1), AND that S(2,2)
% is significantly larger than S(3,3), then we are ok to proceed.
% See that if the second and third singular values are roughly equal,
% this would indicate a points on a line, not a plane.
% You do need to check these numbers, as they will be indicative of a
% potential problem.
% Finally, the magnitude of S(3,3) would be a measure of the noise
% in your regression. It is a measure of the deviations from your
% fitted plane.

% The normal vector is given by the third singular vector, so the
% third (well, last in general) column of V. I'll call the normal
% vector P to be consistent with the question notation.
P = V(:,3);

% The equation of the plane for ANY point on the plane [x,y,z]
% is given by
%
%    dot(P,[x,y,z] - cm) == 0
%
% Essentially this means we subtract off the column mean from our
% point, and then take the dot product with the normal vector. That
% must yield zero for a point on the plane. We can also think of it
% in a different way, that if a point were to lie OFF the plane,
% then this dot product would see some projection along the normal
% vector.
%
% So if your goal is now to predict Z, as a function of X and Y,
% we simply expand that dot product.
% 
%   dot(P,[x,y,z]) - dot(P,cm) = 0
%
%   P(1)*X + P(2)*Y + P(3)*Z - dot(P,cm) = 0
%
% or simply (assuming P(3), the coefficient of Z) is not zero...
Z_mp = (dot(P,cm) - P(1)*x - P(2)*y)/P(3);

i don't know how to do the translation

Any help or direction to appropriate answers would be very welcome.

Thanks in advance!

Ubi

Answer 1

[Matt J]

i don't know how to do the translation

I assume you mean the rotation. You can use my AxelRot tool ( Download ).

 u=cross(P/norm(P),[0,0,1]);
 deg=asind(norm(u));

 [XYZnew, R, t] = AxelRot(XYZ', deg, u, [0;0;0]);

[Ubaldo Pantaleo]

Hi Matt, Thanks for your reply! I tried to apply your suggestion to my case but may be i did some mistakes.

Essentially i loaded the xyz data and combine them into one array: data = load('./input.dat');

x = data(:,1);
y = data(:,2);
z = data(:,3);
clear data;

XYZold = [x,y,z];

Then i found the best-fit plane:

cm_old = mean(XYZold,1);

XYZold0 = bsxfun(@minus,XYZold,cm_old);
[U,S,V] = svd(XYZold0,0);
% The singular values are ordered in decreasing order for svd.
% The vector corresponding to the zero singular value tells us
% the direction of the normal vector to the plane. Note that if
% your data actually fell on a straight line, this will be a problem
% as then there are two vectors normal to your data, so no plane fit.

diag(S)

P = V(:,3);

After i applied your rotation:

u=cross(P/norm(P),[0,0,1]);
 deg=asind(norm(u));

[XYZnew, R, t] = AxelRot(XYZold', deg, u, [0;0;0]);

and found the best-fit plane of the new dataset (XYZnew):

X=(XYZnew(:,1));
Y=(XYZnew(:,2));
Z=(XYZnew(:,3));

>> cm = mean(XYZnew,1);

% subtract off the column means
XYZnew0 = bsxfun(@minus,XYZnew,cm);

% The "regression" as a planar fit is now accomplished by SVD.
% This presumes errors in all three variables. In fact, it makes
% presumptions that the noise variance is the same for all the
% variables. Be very careful, as this fact is built into the model.
% If your goal is merely to fit z(x,y), where x and y were known
% and only z had errors in it, then this is the wrong way
% to do the fit.
[F,G,H] = svd(XYZnew0,0);

% The singular values are ordered in decreasing order for svd.
% The vector corresponding to the zero singular value tells us
% the direction of the normal vector to the plane. Note that if
% your data actually fell on a straight line, this will be a problem
% as then there are two vectors normal to your data, so no plane fit.
% LOOK at the values on the diagonal of G. If the last one is NOT
% essentially small compared to the others, then you have a problem
% here. (If it is numerically zero, then the points fell exactly in
% a plane, with no noise.)
diag(G)

% Assuming that G(3,3) was small compared to G(1,1), AND that G(2,2)
% is significantly larger than G(3,3), then we are ok to proceed.
% See that if the second and third singular values are roughly equal,
% this would indicate a points on a line, not a plane.
% You do need to check these numbers, as they will be indicative of a
% potential problem.
% Finally, the magnitude of G(3,3) would be a measure of the noise
% in your regression. It is a measure of the deviations from your
% fitted plane.

% The normal vector is given by the third singular vector, so the
% third (well, last in general) column of H. I'll call the normal
% vector J to be consistent with the question notation.
J = H(:,3);

Z_mp_new = (dot(J,cm) - J(1)*X - J(2)*y)/J(3);

If i plot the data of the new fitted-plane, i should obtain a plane parallel to xy-plane:

    figure(100)
    scatter(X,Y,2,Z_mp_new)

but it doesn't seem so.

Can you provide me any suggestions?!

Thank you again!

Ubaldo

[Matt J]

As a side note, please font your code separately from your text as I have done for you now, for the sake of readability.

As for the results you're getting, try calculating the rotation angle this way, instead.

deg=acosd(dot(P/norm(P),[0 0 1]));

Beyond that, I would need to see your point cloud data. You can attach it here in a .mat file.

[Ubaldo Pantaleo]

It sounds good!

i'm trying, in the meanwhile an email address where can i send you the .mat file (113 MB)?!

Let you know soon.

Thank you again Matt!

[Ubaldo Pantaleo]

Hi Matt, i did some test but even if numerically it seems that work, when i plot the data the plane is not parallel to xy-plane.

scatter(X,Y,2,Z_mp_new) or scatter3(X,Y,Z_mp_new)

Have you got idea about the reason?! Here i attacched a sample of my xyz data. I really appreciate if you could please help me with this problem. Thank you! Ubaldo

[Matt J]

Here's what I did, and the plot I see looks very reasonable. Also, when I perform a plane fit to the rotated data, I get a normal vector of nnew=[0 0 1].' to very high precision

    XYZ=[x,y,z];

    %%plane fit
    c=mean(XYZ,1);
    Pc=bsxfun(@minus,XYZ,c);    
    [~,~,V]=svd(Pc,0);
    n=V(:,end);  %normal
    n=n*sign(dot(n,[0,0,1])); %orientation convention

    %rotate/align with xy-plane
    u=cross(n,[0,0,1]);
    deg=acosd(dot(n,[0,0,1]));
    XYZnew=AxelRot(XYZ.',deg,u,[0,0,0]).';

    %redo plane fit
    cnew=mean(XYZnew,1);
    Pcnew=bsxfun(@minus,XYZnew,cnew);    
    [~,~,V]=svd(Pcnew,0);

    nnew=V(:,end),  %normal to rotated data

    xyz=XYZnew(1:100:end,:);
    scatter3(xyz(:,1),  xyz(:,2), xyz(:,3))

[Robert Stöger]

Hi all,

I was wanted to use thise code Matt provided to correct for the slope in my point clouds. I have added two pictures, one where you can see the point cloud befor I applied the code and one after.

In the first picture you can clearly see, that in this Point cloud there is a slope along the X-Axis. So fitting this point cloud the the XY-Plane seems to be a good solution to me.

Unfortunately I get a strange output using this code. It just flips the point cloud 90 degrees, so that the Z-Axis now is the former Y-Axis and the Y-Axis now is the former Z-Axis (maybe it's better understandable when you take a look at the pictures).

I do not really understand what each parameter from the code does and hence I do not know how to change it to get the output I want. I did some "try and error" but that did not result in something satisfying.

Any suggestion are warmly welcome :) And please ask If I did not point out well, what I want to achieve.

Best regards, Robert

ORIGINAL POINTCLOUD

POINTCLOUD AFTER APLICATION OF MATT'S CODE

[Matt J]

What does the code give for deg, u, n, and nnew?

Also, how do you extract the particular xyz data that in your diagram defines the plane to be rotated? In any case, if you attach the extracted data here in a .mat file, we can run the code on it ourselves.

[Robert Stöger]

Hi Matt,

sorry for the late reply, I thought I will get an email notification, thats why I did not check here.

Anyways, I found a solution for my problem. It turned out, that your code works perfectly fine if my pointcloud itself is more or less a plane, or flat. But since I have those errected parts it somehow gets confused about the angle and the direction it should turn when I use the raw data as it is. I did not figured out why that is, but what I did I extract the lowest plane from my point cloud with the function "pcfitplane". I have attached a picture of the extracted plane. To this plane I applied your code to get the angle and the other values (I have to admitt I still do not fully understand your code) and apply them to the raw data. This gives me the output I need.

Best regards, Robert

[Tasha Tardieu]

Hi Robert,

I have a very similar problem to yours where I have 3D data that I want to make horizontal. After you found the lowest plane of your point cloud, how did you use AxelRot and apply it to all of your points? It seems that I am still getting the same image before and after rotation. I have attached 2D images of the points for reference with respect to the z axis.

<<

>>

Answer 2

[Preetham Manjunatha]

Fit a oriented bounding box and obtain the orientation of that bounding box and invert the rotation matrix of the Euler angles. This will produce a re-orientated point cloud that is parallel to XY, YZ and XZ planes. If you need to rotate again, Eulaer angles [phi theta psi] can be used, for example [90 0 0] or [-90 0 0] and others. The program which does this can be downloaded at 3DspaceTile.

This program actually divides the 3D space, one can use the divisions as [1 1 1] in XYZ directions. Same program can reorient the point clouds to be parallel to XY, YZ and XZ planes. For the program requirements, please visit the GitHub page.

Hope this helps someone in future.