Problem of rotation of surface on xy plane

8 次查看（过去 30 天）
Elisa
Elisa 2024-7-18
评论： Elisa 2024-7-22
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);

请先登录，再回答此问题。

回答（2 个）

Hassaan
Hassaan 2024-7-18
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 个评论
Elisa
Elisa 2024-7-18
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!!!

请先登录，再进行评论。

Ruchika Parag
Ruchika Parag 2024-7-19
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!
  2 个评论
Elisa
Elisa 2024-7-19
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.
Elisa
Elisa 2024-7-22
Dear @Hassaan, @Star Strider, @Ruchika Parag,
I tried to fix the issue starting from the original point cloud and extracting a series of planes parallel to the maximum dip of the point cloud (i manually calculated the equation of the line of maximum dip and then i extracted 5 planes parallel each others). On those plane I tried to fix a distance and select all the points within it, with the final aim to determine the z values of those planes. However, I get NaN values, also changing the distance. I upload here the code. Can you help me identifying the problem? Thank you so much
clear all
close all
clc
load('XYZ');
X = XYZ(:,1);
Y = XYZ(:,2);
Z = XYZ(:,3);
xyz0=mean(XYZ,1);
A=XYZ-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 original point cloud
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
% line equation parameters
slope = -69 / 41;
intercept = 45999 / 41;
% Define the coordinate ranges
x_min = min(A_rot(1,:));
x_max = max(A_rot(1,:));
z_min = 0;
z_max = max(A_rot(3,:));
% Number of planes
num_planes = 5;
% Extend the length of the line by increasing the range of x values
extended_x_min = x_min - 50;
extended_x_max = x_max + 50;
% Generate points on the initial line for plotting
x_line = linspace(extended_x_min, extended_x_max, 100);
y_line = slope * x_line + intercept;
z_line = zeros(size(x_line)); % Line in the XY plane at z = 0
% Define the spacing between planes
intercept_spacing = linspace(intercept, intercept + (num_planes - 1) * 100, num_planes);
% Generate the grid for the planes
[x_plane, z_plane] = meshgrid(linspace(extended_x_min, extended_x_max, 100), linspace(-100, 100, 100));
% Plot the line
hold on;
plot3(x_line, y_line, z_line, 'r-', 'LineWidth', 2);
min_z_values = zeros(num_planes, 1);
max_z_values = zeros(num_planes, 1);
% Plot the vertical planes and calculate min/max z values
for i = 1:num_planes
current_intercept = intercept_spacing(i);
y_plane = slope * x_plane + current_intercept;
surf(x_plane, y_plane, z_plane, 'FaceAlpha', 0.5, 'EdgeColor', 'none');
% Calculate the distance of each point to the current plane
y_plane_point = slope * A_rot(1,:) + current_intercept;
distances = abs(A_rot(2,:) - y_plane_point);
% Select points within a mm range of the plane
within_range = distances <= 50;
% Extract the z values of the points within the range
z_values_within_range = A_rot(3, within_range);
% Calculate min and max z values for the selected points
if ~isempty(z_values_within_range)
min_z_values(i) = min(z_values_within_range);
max_z_values(i) = max(z_values_within_range);
else
min_z_values(i) = NaN; % No points within range
max_z_values(i) = NaN; % No points within range
end
end
% Add labels and title
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Parallel Vertical Planes and Point Cloud Analysis');
% Add grid and axis equal
grid on;
axis equal;
% Add legend
legend('Line', 'Vertical Planes');
hold off;
% Display the min and max z values
disp('Min Z values for each plane:');
disp(min_z_values);
disp('Max Z values for each plane:');
disp(max_z_values);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Point Cloud Processing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by