50 views (last 30 days)
Show older comments
Elisa on 18 Jul 2024 at 10:53
Edited: Star Strider on 19 Jul 2024 at 12:08
- xyz_c1p1.mat
Open in MATLAB Online
Hi everyone,
I am using this code on this point cloud imported in XYZ format, in which I would like to calculate the difference between the maximum and minimum points in the profile. To do this, the cloud must be rotated on the XY plane (see output figure). However, it seems to have an incorrect rotation. Can someone help me fix the rotation?
Thank you very much
clear all
close all
clc
load('xyz_c1p1');
X = xyz_c1p1(:,1);
Y = xyz_c1p1(:,2);
Z = xyz_c1p1(:,3);
xyz0=mean(xyz_c1p1,1);
A=xyz_c1p1-xyz0; % center the data at zero
% Find the direction of most variance using SVD and rotate the data to make
% that the x-axis
[~,~,V]=svd(A,0);
a=cross(V(:,3),[0;0;1]);
T=makehgtform('axisrotate', a, -atan2(norm(a),V(3,3)));;
R=T(1:3,1:3);
A_rot = A*R;
A_rot = A_rot + [xyz0(1) xyz0(2) 0]; % move so the centers are aligned in z-diection
% Plot the raw data
close all
scatter3(X,Y,Z,0.1,"magenta")
hold on
scatter3(A_rot(:,1),A_rot(:,2),A_rot(:,3),0.1,'blue');
xlabel('X-Axis','FontSize',14,'FontWeight','bold')
ylabel('Y-Axis','FontSize',14,'FontWeight','bold')
zlabel('Z-Axis','FontSize',14,'FontWeight','bold')
axis equal
Zmax = max(A_rot(:,3))
Zmin = min(A_rot(:,3))
Rz = Zmax - Zmin
hold on
% Alpha Shape
shpINT=alphaShape(A_rot,5);
figure(3)
plot(shpINT)
VolAlphaShapeINT=volume(shpINT);
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
Answers (3)
Hassaan on 18 Jul 2024 at 11:14
Open in MATLAB Online
clear all;
close all;
clc;
% Load your point cloud data
load('xyz_c1p1'); % Assuming 'xyz_c1p1' is in the format [X, Y, Z]
X = xyz_c1p1(:,1);
Y = xyz_c1p1(:,2);
Z = xyz_c1p1(:,3);
% Center the data at zero
xyz0 = mean(xyz_c1p1, 1);
A = xyz_c1p1 - xyz0;
% Find the direction of most variance using SVD and align it with the Z-axis
[~, ~, V] = svd(A, 0);
% Rotation to align the principal component with the Z-axis
axis = cross(V(:,3), [0; 0; 1]); % Cross product to find the rotation axis
angle = acos(dot(V(:,3), [0; 0; 1]) / (norm(V(:,3)) * norm([0; 0; 1]))); % Angle calculation
R = axang2rotm([axis' angle]); % Create a rotation matrix from axis-angle
% Apply rotation
A_rot = (R * A')'; % Rotate and transpose back
A_rot = A_rot + repmat(xyz0, size(A_rot, 1), 1);
% Plot the raw data
figure(1);
scatter3(X, Y, Z, 0.1, "magenta");
hold on;
% Plot the rotated data
scatter3(A_rot(:,1), A_rot(:,2), A_rot(:,3), 0.1, 'blue');
xlabel('X-Axis', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Y-Axis', 'FontSize', 14, 'FontWeight', 'bold');
zlabel('Z-Axis', 'FontSize', 14, 'FontWeight', 'bold');
axis equal;
title('Raw and Rotated Point Cloud');
% Calculate the range of Z-values after rotation
Zmax = max(A_rot(:,3));
Zmin = min(A_rot(:,3));
Rz = Zmax - Zmin;
disp(['Range in Z after rotation: ', num2str(Rz)]);
% Alpha Shape to visualize the boundary of the point cloud (optional)
shpINT = alphaShape(A_rot, 5);
figure(2);
plot(shpINT);
title('Alpha Shape of Rotated Point Cloud');
VolAlphaShapeINT = volume(shpINT);
disp(['Volume of Alpha Shape: ', num2str(VolAlphaShapeINT)]);
1 Comment Show -1 older commentsHide -1 older comments
Show -1 older commentsHide -1 older comments
Elisa on 18 Jul 2024 at 13:05
Direct link to this comment
https://nl.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3214801
Open in MATLAB Online
dear @Hassaan thank so much for your kind reply and help!! is not easy to explain what I need. I try uploading an example of what I need:
the original point cloud has specific coordinates x,y,z and I want to rotate my point cloud in the same way you see in the above picture, i.e. on the x-y plane, with z = 0. In that way I can calculate the roughness as
Rz = Zmax - Zmin;
I hope to be more clear now. Thank you so so much!!!
Sign in to comment.
Star Strider on 18 Jul 2024 at 14:33
Edited: Star Strider on 18 Jul 2024 at 16:39
Open in MATLAB Online
- xyz_c1p1.mat
I am not certain what you want. Yopur data do not appear to define a surface, instead they appear almost linear.
An alternative approach culd be something like this —
% clear all
% close all
% clc
load('xyz_c1p1');
X = xyz_c1p1(:,1);
Y = xyz_c1p1(:,2);
Z = xyz_c1p1(:,3);
B = [X Y ones(size(X))] \ Z; % Simple Multitple Linear Regeression
Zfit = [X Y ones(size(X))] * B; % Fit To Data
[Azf,Elf,Rf] = cart2sph(X, Y, Zfit); % Convert Fit
[Azd,Eld,Rd] = cart2sph(X, Y, Z); % Convert Data
figure
scatter3(X, Y, Z, '.', 'DisplayName','Data')
hold on
plot3(X, Y, Zfit, '-r', 'LineWidth',2, 'DisplayName','Regression')
hold off
axis('equal')
xlabel('X')
ylabel('Y')
zlabel('Z')
legend('Location', 'best')
figure
hp31 = plot3(X, Y, Z, '.', 'DisplayName','Data');
hold on
hp32 = plot3(X, Y, Zfit, '-r', 'LineWidth',2, 'DisplayName','Regression');
hold off
grid on
axis('equal')
xlabel('X')
ylabel('Y')
zlabel('Z')
title('Rotated Plot')
legend('Location', 'best')
Azc = mean(rad2deg(Azf));
Elc = mean(rad2deg(Elf));
orign = mean([X Y Z]);
rotate(hp31, [Elc/90 -Azc/90 0], 90, orign)
rotate(hp32, [Elc/90 -Azc/90 0], 90, orign)
Xr = hp31.XData
Xr = 1x3188
534.3975 512.9273 540.2768 537.8913 536.3009 535.5057 534.7105 533.9153 535.0235 534.2283 533.4331 532.9510 534.0592 543.1992 542.4041 540.8137 540.0185 539.2233 538.4281 537.6329 533.6569 532.8618 532.0666 531.2714 512.9820 512.1868 533.9700 533.1748 532.3796 531.5844
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Yr = hp31.YData
Yr = 1x3188
901.4372 885.0655 906.9204 905.1013 903.8886 903.2822 902.6759 902.0695 903.9146 903.3082 902.7018 903.3342 905.1792 906.3863 905.7799 904.5672 903.9608 903.3545 902.7481 902.1418 899.1100 898.5036 897.8972 897.2909 883.3446 882.7383 900.3486 899.7423 899.1359 898.5296
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Zr = hp31.ZData
Zr = 1x3188
676.5572 676.5572 677.1635 677.1635 677.1635 677.1635 677.1635 677.1635 677.7699 677.7699 677.7699 678.3762 678.9826 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 677.3524 677.3524 677.3524 677.3524
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Zfitr = hp32.ZData
Zfitr = 1x3188
676.5572 676.5572 677.1635 677.1635 677.1635 677.1635 677.1635 677.1635 677.7699 677.7699 677.7699 678.3762 678.9826 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 676.7460 677.3524 677.3524 677.3524 677.3524
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
% figure
% scatter3(Xr, Yr, Zfitr, '.', 'DisplayName','Regression')
% hold on
% scatter3(Xr, Yr, Zr, '.', 'DisplayName','Data')
% hold off
% xlabel('Az')
% ylabel('El')
% zlabel('R')
% legend('Location', 'best')
[Az1,Az2] = bounds(Azf)
Az1 = 0.7615
Az2 = 1.2281
[El1,El2] = bounds(Elf)
El1 = 0.3869
El2 = 1.0590
[R1,R2] = bounds(Rf)
R1 = 783.8440
R2 = 1.2100e+03
Rrsd = Zfitr - Zr; % Calculate Residuals
figure
stem3(Xr, Yr, Rrsd, '.', 'LineWidth',0.01)
axlims = axis;
hold on
patch([axlims([1 2]) axlims([2 1])], [axlims([3 3]) axlims([4 4])], zeros(1,4), 'r', 'FaceAlpha',0.25 )
hold off
xlabel('Az')
ylabel('El')
zlabel('R')
title('Residuals')
view(27,30)
return
xyz0=mean(xyz_c1p1,1);
A=xyz_c1p1-xyz0; % center the data at zero
% Find the direction of most variance using SVD and rotate the data to make
% that the x-axis
[~,~,V]=svd(A,0);
a=cross(V(:,3),[0;0;1])
T=makehgtform('axisrotate', a, -atan2(norm(a),V(3,3)));;
R=T(1:3,1:3);
A_rot = A*R;
A_rot = A_rot + [xyz0(1) xyz0(2) 0]; % move so the centers are aligned in z-diection
% Plot the raw data
close all
scatter3(X,Y,Z,0.1,"magenta")
hold on
% scatter3(A_rot(:,1),A_rot(:,2),A_rot(:,3),0.1,'blue');
xlabel('X-Axis','FontSize',14,'FontWeight','bold')
ylabel('Y-Axis','FontSize',14,'FontWeight','bold')
zlabel('Z-Axis','FontSize',14,'FontWeight','bold')
axis equal
Zmax = max(A_rot(:,3))
Zmin = min(A_rot(:,3))
Rz = Zmax - Zmin
hold on
% Alpha Shape
shpINT=alphaShape(A_rot,5);
figure(3)
plot(shpINT)
VolAlphaShapeINT=volume(shpINT);
EDIT — (18 Jul 2024 at 16:39)
Added Rotated Plot and recalculated Residuals from the rotated data.
.
4 Comments Show 2 older commentsHide 2 older comments
Show 2 older commentsHide 2 older comments
Elisa on 18 Jul 2024 at 15:31
Direct link to this comment
https://nl.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3214961
dear @Star Strider
Thank you so much for your efforts to help me! The method of using residuals works well for me, and I can definitely use it.
To answer your question, my task involves segmenting a small section of a point cloud (similar to a profile) and calculating the roughness. I need to rotate the cloud according to the reference image, which means setting z = 0.Basically the cloud should be put on the xy plane (see the example attached before)
Star Strider on 18 Jul 2024 at 16:52
Direct link to this comment
https://nl.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215061
My pleasure!
I reconsidered this problem and then I edited my earlier code, adding:
figure
hp31 = plot3(X, Y, Z, '.', 'DisplayName','Data');
hold on
hp32 = plot3(X, Y, Zfit, '-r', 'LineWidth',2, 'DisplayName','Regression');
hold off
grid on
axis('equal')
xlabel('X')
ylabel('Y')
zlabel('Z')
title('Rotated Plot')
legend('Location', 'best')
Azc = mean(rad2deg(Azf));
Elc = mean(rad2deg(Elf));
orign = mean([X Y Z]);
rotate(hp31, [Elc/90 -Azc/90 0], 90, orign)
rotate(hp32, [Elc/90 -Azc/90 0], 90, orign)
and re-calculated and plotted the residuals based on the values of the rotated data. It is not absolutely flat in the Z-direction, however tthe Z-axis in that plot only goes from 660 to 675, and the resuduals calculated from the rotated values are essentially zero. (I extracted the relevant values from the rotated line objects in the plots. Changing them from scatter3 to plot3 is the only significant change, because rotate does not work with scatter objects. The effect is the same, otherwise.) Everything else is still there (some commented-out). See the rotate function documentation for details on it. Experiment with it, since you may be able to get even better results.
.
Elisa on 19 Jul 2024 at 7:22
Direct link to this comment
https://nl.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215551
Edited: Elisa on 19 Jul 2024 at 7:26
@Star Strider thank you again for your effort! Certainly, I must have explained myself poorly, and I apologize. In the "Rotated plot" figure, I still see the same issue. What I would like is for the graph to be simply "flattened" on the z-axis (z = 0), lying on the xy plane without appearing rotated. Right now, it looks like the graph is rotated by 90 degrees. I insert an example:
Starting from the point cloud profile aligned in this way:
i want it rotated in this way:
Sorry for misunderstanding
Star Strider on 19 Jul 2024 at 10:47
Direct link to this comment
https://nl.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215681
Edited: Star Strider on 19 Jul 2024 at 12:08
Open in MATLAB Online
- xyz_c1p1.mat
The Rotated Plot is not absolutely flat (that may not be possible since it is of coures a point cloud), however it is reasonably close.
This is the best I can do —
% clear all
% close all
% clc
load('xyz_c1p1');
X = xyz_c1p1(:,1);
Y = xyz_c1p1(:,2);
Z = xyz_c1p1(:,3);
B = [X Y ones(size(X))] \ Z; % Simple Multitple Linear Regeression
Zfit = [X Y ones(size(X))] * B; % Fit To Data
[Azf,Elf,Rf] = cart2sph(X, Y, Zfit); % Convert Fit
[Azd,Eld,Rd] = cart2sph(X, Y, Z); % Convert Data
figure
scatter3(X, Y, Z, '.', 'DisplayName','Data')
hold on
plot3(X, Y, Zfit, '-r', 'LineWidth',2, 'DisplayName','Regression')
hold off
axis('equal')
xlabel('X')
ylabel('Y')
zlabel('Z')
legend('Location', 'best')
figure
hp31 = plot3(X, Y, Z, '.', 'DisplayName','Data');
hold on
hp32 = plot3(X, Y, Zfit, '-r', 'LineWidth',2, 'DisplayName','Regression');
hold off
grid on
axis('equal')
xlabel('X')
ylabel('Y')
zlabel('Z')
title('Rotated Plot')
legend('Location', 'best')
Azc = mean(rad2deg(Azf));
Elc = mean(rad2deg(Elf));
orign = mean([X Y Z]);
orign = [0 0 0];
rotate(hp31, [Elc/45 -Azc/45 0], 45, orign)
rotate(hp32, [Elc/45 -Azc/45 0], 45, orign)
Xr = hp31.XData;
Yr = hp31.YData;
Zr = hp31.ZData;
Zfitr = hp32.ZData;
meanZr = mean(Zr)
meanZr = 896.2878
meanZfitr = mean(Zfitr)
meanZfitr = 896.2878
[minZr,maxZr] = bounds(Zr-meanZr)
minZr = -266.3257
maxZr = 277.1837
[minZfitr,maxZfitr] = bounds(Zfitr-meanZfitr)
minZfitr = -256.6166
maxZfitr = 266.5059
figure
scatter3(Xr, Yr, Zfitr-meanZfitr, '.', 'DisplayName','Regression')
hold on
scatter3(Xr, Yr, Zr-meanZr, '.', 'DisplayName','Data')
axlims = axis;
patch([axlims([1 2]) axlims([2 1])], [axlims([3 3]) axlims([4 4])], zeros(1,4), 'r', 'FaceAlpha',0.25, 'DisplayName','Zero Plane')
hold off
% xlabel('Az')
% ylabel('El')
% zlabel('R')
xlabel('X')
ylabel('Y')
zlabel('Z')
legend('Location', 'best')
view(-27,30)
grid on
zlim([minZfitr maxZr])
[Az1,Az2] = bounds(Azf)
Az1 = 0.7615
Az2 = 1.2281
[El1,El2] = bounds(Elf)
El1 = 0.3869
El2 = 1.0590
[R1,R2] = bounds(Rf)
R1 = 783.8440
R2 = 1.2100e+03
Rrsd = Zfitr - Zr; % Calculate Residuals
figure
stem3(Xr, Yr, Rrsd, '.', 'LineWidth',0.01)
axlims = axis;
hold on
patch([axlims([1 2]) axlims([2 1])], [axlims([3 3]) axlims([4 4])], zeros(1,4), 'r', 'FaceAlpha',0.25 )
hold off
xlabel('X')
ylabel('Y')
zlabel('Z')
% xlabel('Az')
% ylabel('El')
% zlabel('R')
title('Residuals')
view(27,30)
Stopping here, unless I can get some further insight into this.
Good luck!
EDIT — (19 Jul 2024 att 12:08)
Note thtat you can interactively rotate your figure in the user interface by clicking on the appropriate icon in the top toolstrip (alternatively by invoking the rotate3d function). The azimutth and elevattion will appear in the lower left corner of the figure UI as you do this. (Record those somewhere.) When you find a rotation that gives you the result you want, you can save that figure (saving it as a .fig file is best, since that preserves all the information) and then use that information in your code. You can also do that in the figure UI by clicking on the ‘File’ option in the top toolstrip.
.
Sign in to comment.
Ruchika Parag on 19 Jul 2024 at 7:41
Open in MATLAB Online
Hi Elisa, to address the rotation issue in your point cloud data, we need to ensure that the rotation aligns the profile correctly along the XY plane. The current code uses Singular Value Decomposition (SVD) to find the principal components, but the rotation might not be applied correctly to achieve the desired orientation.Please modify your code as follows that ensures the rotation aligns the point cloud properly along the XY plane:
clear all
close all
clc
% Load the point cloud data
load('xyz_c1p1');
% Extract X, Y, Z coordinates
X = xyz_c1p1(:,1);
Y = xyz_c1p1(:,2);
Z = xyz_c1p1(:,3);
% Center the data at zero
xyz0 = mean(xyz_c1p1, 1);
A = xyz_c1p1 - xyz0;
% Perform SVD to find the principal axes
[~, ~, V] = svd(A, 'econ');
% Calculate the rotation matrix to align the first principal component with the X-axis
rotationAngle = atan2(V(2,1), V(1,1));
R = [cos(rotationAngle) -sin(rotationAngle) 0;
sin(rotationAngle) cos(rotationAngle) 0;
0 0 1];
% Rotate the data
A_rot = A * R;
% Translate back to the original center
A_rot = A_rot + xyz0;
% Plot the original and rotated data
figure;
scatter3(X, Y, Z, 0.1, "magenta");
hold on;
scatter3(A_rot(:,1), A_rot(:,2), A_rot(:,3), 0.1, 'blue');
xlabel('X-Axis', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Y-Axis', 'FontSize', 14, 'FontWeight', 'bold');
zlabel('Z-Axis', 'FontSize', 14, 'FontWeight', 'bold');
axis equal;
% Calculate the difference between max and min Z values in the rotated data
Zmax = max(A_rot(:,3));
Zmin = min(A_rot(:,3));
Rz = Zmax - Zmin;
% Display the results
disp(['Maximum Z: ', num2str(Zmax)]);
disp(['Minimum Z: ', num2str(Zmin)]);
disp(['Difference (Rz): ', num2str(Rz)]);
% Alpha Shape
shpINT = alphaShape(A_rot(:,1:2), 5); % Use only X and Y for alpha shape
figure;
plot(shpINT);
VolAlphaShapeINT = volume(shpINT);
disp(['Volume of Alpha Shape: ', num2str(VolAlphaShapeINT)]);
This revised code should give you the correct rotation of your point cloud data on the XY plane. Hope this helps!
1 Comment Show -1 older commentsHide -1 older comments
Show -1 older commentsHide -1 older comments
Elisa on 19 Jul 2024 at 8:16
Direct link to this comment
https://nl.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215581
dear @Ruchika Parag tahnk you so much for your help! unfortunately this is still not what I need. I try to insert a skecth below. I would like to rotate the point cloud based on the same origin of the original one, without changing coordinates x and y. At the end it should be on z = 0 (see pictures uploaded in the last answer). Hope to be more clear.
Sign in to comment.
Sign in to answer this question.
See Also
Categories
Image Processing and Computer VisionComputer Vision ToolboxPoint Cloud Processing
Find more on Point Cloud Processing in Help Center and File Exchange
Tags
- rotation
- axis rotation
- xyz file
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
