Problem of rotation of surface on xy plane

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);

Answer 1

[Hassaan]

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)]);

[Elisa]

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!!!

Answer 2

[Ruchika Parag]

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!

[Elisa]

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]

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);