# Extracting profiles from a point cloud surface along the mazimum dip of the plane

46 views (last 30 days)
Elisa on 19 Jul 2024 at 9:40
Commented: Elisa on 29 Jul 2024 at 13:59
Dear all,
Based on the result of this code (blue cloud rotated), I need to extract a series of profiles (with a fixed diastance) to get the Z values of all the points constituting the profile. These profiles should be extracted from vertical planes oriented according to the maximum slope of the cloud.
you have other ideas, they are accepted!! I have included a sketch of what I would like.
clear all
close all
clc
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 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
##### 3 CommentsShow 1 older commentHide 1 older comment
Elisa on 22 Jul 2024 at 8:21
dear @Umar thank you so much!! i get an error: "At least one END is missing. The statement beginning here does not have a matching end. "
Umar on 22 Jul 2024 at 11:15
Edited: Umar on 23 Jul 2024 at 14:14
Hi Elisa,
Please see Avni response below. She has already provided the solution you are seeking.

Avni Agrawal on 23 Jul 2024 at 11:50
Hi @Elisa,
I understand that you are trying to extract the profiles oriented according to the maximum slope of the cloud. To achieve this:
• Compute the maximum slope direction and rotate the point cloud:
% code already mentioned in the question
[~,~,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
• Extract Vertical planes: Create vertical planes at fixed intervals and extract points that lie on these planes
• Extract Z values: For each plane, extract the Z values of the points that lie on that
clear all
close all
clc
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-direction
% 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
% Extract profiles
fixed_distance = 1; % Define the fixed distance between profiles
x_min = min(A_rot(:,1));
x_max = max(A_rot(:,1));
profile_planes = x_min:fixed_distance:x_max;
for i = 1:length(profile_planes)
plane_x = profile_planes(i);
% Find points within a small tolerance of the plane
tolerance = 0.05; % Adjust the tolerance as needed
profile_points = A_rot(abs(A_rot(:,1) - plane_x) < tolerance, :);
if ~isempty(profile_points)
% Plot the profile points
scatter3(profile_points(:,1), profile_points(:,2), profile_points(:,3), 10, 'filled');
% Extract Z values
Z_values = profile_points(:,3);
% Display Z values (or process further as needed)
disp(['Profile at X = ', num2str(plane_x), ': Z values = ', num2str(Z_values')]);
end
end
hold off
I hope this helps!
Elisa on 29 Jul 2024 at 13:59
dear @Avni Agrawal thank you so much for the help! I tried my self to implement the code to do that but I am myssing the way to determine the maximum slope and thus the equation of the line (that becomes a plane then) parallel to it. If you see in my code I inserted an equation that I manually calculated based on the coordinates of selected points on the point cloud plot. Can you help me implementing this code to automatically find the planes aligned on the maxium slope of the projected plane? thank you so much
clear all
close all
clc
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(:,2) + current_intercept;
distances = abs(A_rot(:,2) - y_plane_point);
% Select points within a mm range of the plane
within_range = distances <= 150;
% Extract the z values of the points within the range
% pippo=A_rot(:,3)
z_values_within_range = squeeze(A_rot(within_range,3));
% 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); %,'omitnan');
max_z_values(i) = max(z_values_within_range); %,'omitnan');
else
min_z_values(i) = NaN; % No points within range
max_z_values(i) = NaN; % No points within range
end
end
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Parallel Vertical Planes and Point Cloud Analysis');
% Add grid and axis equal
grid on;
axis equal;
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);
Max_z_values_mm = 70.69/1000*max_z_values
Min_z_values_mm = 70.69/1000*min_z_values
Rz = Max_z_values_mm - Min_z_values_mm
E5 = 1.5715*Rz+4.0318
mean_E5 = mean(E5)
E6 = 4.4192*Rz.^0.6482
mean_E6 = mean(E6)