Problem of rotation of surface on xy plane
56 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (3 件)
Hassaan
2024 年 7 月 18 日 11:14
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)]);
Star Strider
2024 年 7 月 18 日 14:33
編集済み: Star Strider
2024 年 7 月 18 日 16:39
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
Yr = hp31.YData
Zr = hp31.ZData
Zfitr = hp32.ZData
% 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)
[El1,El2] = bounds(Elf)
[R1,R2] = bounds(Rf)
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.
.
5 件のコメント
Star Strider
2024 年 7 月 19 日 10:47
編集済み: Star Strider
2024 年 7 月 19 日 12:08
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)
meanZfitr = mean(Zfitr)
[minZr,maxZr] = bounds(Zr-meanZr)
[minZfitr,maxZfitr] = bounds(Zfitr-meanZfitr)
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)
[El1,El2] = bounds(Elf)
[R1,R2] = bounds(Rf)
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
約12時間 前
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')
% 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
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);
.
Ruchika Parag
2024 年 7 月 19 日 7:41
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!
参考
カテゴリ
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!