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

[Star Strider]

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.

.

[Elisa]

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]

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]

@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]

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]

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
end

distances = 1x3

1.0e+03 * 1.0219 1.1225 1.1183

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

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

z_values_within_range =

  1x0 empty double row vector

distances = 1x3

1.0e+03 * 1.1219 1.0225 1.2183

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

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

z_values_within_range =

  1x0 empty double row vector

distances = 1x3

1.0e+03 * 1.2219 0.9225 1.3183

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

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

z_values_within_range =

  1x0 empty double row vector

distances = 1x3

1.0e+03 * 1.3219 0.8225 1.4183

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

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

z_values_within_range =

  1x0 empty double row vector

distances = 1x3

1.0e+03 * 1.4219 0.7225 1.5183

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

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

z_values_within_range =

  1x0 empty double row vector

% 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:');

Min Z values for each plane:

disp(min_z_values);

   NaN
   NaN
   NaN
   NaN
   NaN

disp('Max Z values for each plane:');

Max Z values for each plane:

disp(max_z_values);

   NaN
   NaN
   NaN
   NaN
   NaN

.

Answer 3

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