Problem of rotation of surface on xy plane (2024)

58vues (au cours des 30derniers jours)

Afficher commentaires plus anciens

Elisa le 18 Juil 2024 à 10:53

Lien

Utiliser le lien direct vers cette question

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane

⋮

Lien

Utiliser le lien direct vers cette question

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane

Commenté: Star Strider le 22 Juil 2024 à 11:23

xyz_c1p1.mat

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

0commentaires
Afficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

Connectez-vous pour répondre à cette question.

Réponses (3)

Hassaan le 18 Juil 2024 à 11:14

Lien

Utiliser le lien direct vers cette réponse

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#answer_1487491

⋮

Lien

Utiliser le lien direct vers cette réponse

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#answer_1487491

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

1commentaire
Afficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens

Elisa le 18 Juil 2024 à 13:05

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3214801

⋮

Lien

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3214801

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

Connectez-vous pour commenter.

Star Strider le 18 Juil 2024 à 14:33

Lien

Utiliser le lien direct vers cette réponse

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#answer_1487591

⋮

Lien

Utiliser le lien direct vers cette réponse

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#answer_1487591

Modifié(e): Star Strider le 18 Juil 2024 à 16:39

Ouvrir dans MATLAB Online

Ran in:

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

5commentaires
Afficher 3 commentaires plus anciensMasquer 3 commentaires plus anciens

Elisa le 18 Juil 2024 à 15:31

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3214961

⋮

Lien

Utiliser le lien direct vers ce commentaire

https://fr.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 le 18 Juil 2024 à 16:52

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215061

⋮

Lien

Utiliser le lien direct vers ce commentaire

https://fr.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 le 19 Juil 2024 à 7:22

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215551

⋮

Lien

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215551

Modifié(e): Elisa le 19 Juil 2024 à 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 le 19 Juil 2024 à 10:47

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215681

⋮

Lien

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3215681

Modifié(e): Star Strider le 19 Juil 2024 à 12:08

Ouvrir dans MATLAB Online

Ran in:

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.

Star Strider le 22 Juil 2024 à 11:23

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3217366

⋮

Lien

Utiliser le lien direct vers ce commentaire

https://fr.mathworks.com/matlabcentral/answers/2138421-problem-of-rotation-of-surface-on-xy-plane#comment_3217366

Ouvrir dans MATLAB Online

Ran in:

XYZ.mat

Running yur code, checking ‘distances’ reveals that they are all on the order of so none of them are . The following if block evaluates correctly and fills ‘min_z_values’ with NaN. (I also checked to be certain that there are no NaN values or 0 values in your data, since that can cause NaN values to appear as the results of calculatitons. There aren’t so that’s not the problem.)

I don’t understand how your code works, so that’s the best I can do with respect to ttroubleshooting it.

My analysis —

clear all

close all

clc

load('XYZ');

whos('-file','XYZ')

Name Size Bytes Class Attributes XYZ 325728x3 7817472 double

% XYZ

% Lm = any(isnan(XYZ))

% Lm0 = any(XYZ == 0)

% return

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

distances = 1x3

1.0e+03 * 1.0219 1.1225 1.1183

<mw-icon class=""></mw-icon>