How to rotate a pointcloud data and align parallel with y-z axis?

I have a point cloud data with x,y and z values. It is little inclined in the y-axis little bit.I want to fit and well align with the axis. I have tried PCA based rotation but it causes over rotation. I have attached the image shows inclined data and excel file of xyz values. Please help me to rotate the data and align perfect with the axis. Thanks for your support.

Answer 1

[Matt J]

yz=readmatrix('CS_DATA.xls','Range',[1,2]);

[~,v1]=min(yz*[+1;+1]);
[~,v2]=min(yz*[-1;+1]);

[x1,y1,x2,y2]=deal(yz(v1,1), yz(v1,2),   yz(v2,1), yz(v2,2));

t=atan2d(y2-y1,x2-x1); %rotation angle

YZ=(yz-mean(yz))*[cosd(t), -sind(t);sind(t), cosd(t)] + mean(yz);


scatter(yz(:,1), yz(:,2)); hold on; scatter(YZ(:,1), YZ(:,2)); hold off
legend('Initial','Corrected', Location='east')
ylabel Z; xlabel Y

[William Rose]

@Mathanraj, The solution of @Matt J is great!

He finds the (x1,y1) and (x2,y2) points, which are the points at the bottom left and bottom right corners. He finds them as by finding the indices of the points whose projections onto the +45 and +135 degree lines are minimal.*

[~,v1]=min(yz*[+1;+1]);
[~,v2]=min(yz*[-1;+1]);
[x1,y1,x2,y2]=deal(yz(v1,1),yz(v1,2),yz(v2,1),yz(v2,2));

t=the angle of the line that connects the bottom corner points, relative to the +Y axis. Then he rotates all the points by angle (-t). Genius!

*If he had used the -45 degree line ([+1,-1]) instead of the +135 degree line ([-1,+1]), or if the beam had been centered below the +45 d line instead of above the +45 d line, then (x2,y2), found as above, would not have been the bottom right corner. And so he would have used the -45 d line for reference, if the beam had been below the +45 d line.

[Matt J]

@William Rose Very kind of you to say. I should mention, though, that the solution is not very noise-robust. If you have a noisy outlier down near one of the bottom corners, the method will fail.

[William Rose]

@Matt J, I admire solution that are simple and elegant and insightful, and yours is.

@Mathanraj, The figure below shows orthogonal views of the raw point cloud. The beam looks sheared in the X-Z plane (right panel). The real beam might not be sheared, but the sampled points are. Were these points taken with a 3D scanner? A partial view blockage can produce scan samples like this. The principal component approach will give inaccurate results in this case. Even if the samples do not show shearing, the principal component approach will give misleading results if the "random" samples are not random, for example if there are significantly more samples on the bottom right than the bottom left, etc.

[William Rose]

@Mathanraj, the shearing in the X-Z plane looks worse in the figure above than it really is, due to axis auto-scaling. The R and B axes don't look orthognal in the right panel, due to autoscaling, but they are orthogonal - they have to be.

[Mathanraj]

@Matt J Thanks for your support. The code you provided is simple to understand and works fine. I am aware that the outliers at the bottom corners affect the output. And, Sorry for the late reply.

[Mathanraj]

@William Rose Thanks for your detailed explanation and it really helps me to understand the concept of the Matt J's code. Yes, The data was collected by the scanner.

About the sheared look of the data, I used "axis equal" to rectify this issue. when I am using this command I can get the actual fitted look of the data without any shear. Is it okay to use this command?

And, sorry for the late reply.

[William Rose]

Yes it is fine and good to use axis equal.

Answer 2

[Constantino Carlos Reyes-Aldasoro]

This answer may solve your problem

https://uk.mathworks.com/matlabcentral/answers/123763-how-to-rotate-entire-3d-data-with-x-y-z-values-along-a-particular-axis-say-x-axis

Answer 3

[William Rose]

@Mathanraj,

xyz=xlsread('CS_DATA.xls');
C=cov(xyz);
[V,D]=eig(C);  % eigenvectors and eigenvalues of the covariance matrix
diag(C);

The eigenvectors are the directions of the axes of the best fit ellipse. The square roots of the eigenvalues are the lengths of the axes. Multiply by 2 so they'll be easier to see on the plot.

ellAxLen=2*diag(D).^0.5;

Display the axis lengths and the standard deviations along x,y,z. The similarity suggests that we are on the right track.

disp([std(xyz);ellAxLen'])

    0.0005    0.0499    0.0865
    0.0007    0.0997    0.1731

Plot the data

figure(1)
plot3(xyz(:,1),xyz(:,2),xyz(:,3),'k.'); 
xlabel('X'); ylabel('Y'); zlabel('Z'); grid on; hold on

Make matrices to represent the endpoints of the best fit axes:

ax1=[mean(xyz);mean(xyz)+ellAxLen(1)*V(:,1)'];
ax2=[mean(xyz);mean(xyz)+ellAxLen(2)*V(:,2)'];
ax3=[mean(xyz);mean(xyz)+ellAxLen(3)*V(:,3)'];

Add the axes to the plot

plot3(ax1(:,1),ax1(:,2),ax1(:,3),'-r','LineWidth',2);
plot3(ax2(:,1),ax2(:,2),ax2(:,3),'-g','LineWidth',2);
plot3(ax3(:,1),ax3(:,2),ax3(:,3),'-b','LineWidth',2);

When you do it, you will be able to use the 3d rotate tool to spin the plot around and see the axes better.

Use the eigenvector directions to rotate the data. Scale the eigenvectors to make a rotation matrix.

R1=zeros(3,3);
for i=1:3
    R1(:,i)=V(:,i)/norm(V(:,i));
end

Rotate the data

xyzRot=(inv(R1)*(xyz-mean(xyz))')'+mean(xyz);

Plot the rotated points:

figure(2); plot3(xyzRot(:,1),xyzRot(:,2),xyzRot(:,3),'k.'); 
title('Rotated Point Cloud');
xlabel('X'); ylabel('Y'); zlabel('Z'); grid on;

After viewing the point cloud above from different angles, I think it is over-rotated about the X axis (middle panel below). See screenshot below of 3 views. The rotations about Y and Z look reasonable.

[William Rose]

[edit: correct typo and add description of alternate model]

A more complicated way would be to make a model I-beam with 6 adjustable parameters: L,W,D and orientation, which could be represented as 3 Euler angles, etc. (Assume the beam is made of infinitely thin plates.) Compute the sum squared distances from each point to the nearest part of the model I-beam. This would be tricky since you'd have to check distance from each point to the top plate, the bottom plate, and the middle plate, and take whichever distance is smallest. You may need 3 additional paramers for the center locaiton, if the points are not randomly sampled relative to the center.

Or the model beam is made of 3 infinite plates: 2 parallel, and one perpendicular to the other 2. The beam model has 4 parameters, not including orientation: its center location (x0,y0,z0) and the paprallel plate separation L. Infinite plates will make it easier to compute the distance to each of the three plates.

Then you use a fitting function like fmincon() to find the unknown parameters, including the 3 angles that minimize the sum squared distance. Then you would rotate by the inverse of the rotation specified by the 3 Euler angles.

Better yet, just eyeball it.