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

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.

Thanks in advance!!

you have other ideas, they are accepted!!  I have included a sketch of what I would like.
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 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

[Umar]

Hi Elisa,

To achieve the extraction of profiles from the rotated 3D point cloud, I had to follow a series of steps after rotating the data to align with the maximum slope direction. Here's an updated version of the code that includes the extraction of profiles.

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 the z-direction

% Extract profiles along vertical planes oriented based on the maximum slope

% Maximum slope direction

max_slope_direction = V(:,3);

% Fixed distance for profile extraction

fixed_distance = 0.1;

profiles = cell(1, size(A_rot, 1));

for i = 1:size(A_rot, 1) profile_points = A_rot(i, :); profile_points = repmat(profile_points, size(A_rot, 1), 1); distances = sqrt(sum((A_rot - profile_points).^2, 2)); [~, sorted_indices] = sort(distances); profile = A_rot(sorted_indices(2:end), :); % Exclude the point itself profiles{i} = profile; end

% Plot the extracted profiles

figure;

hold on;

for i = 1:size(profiles, 2) scatter3(profiles{i}(:,1), profiles{i}(:,2), profiles{i}(:,3), 0.1, 'green'); end

xlabel('X-Axis', 'FontSize', 14, 'FontWeight', 'bold');

ylabel('Y-Axis', 'FontSize', 14, 'FontWeight', 'bold');

zlabel('Z-Axis', 'FontSize', 14, 'FontWeight', 'bold');

axis equal;

So, first I needed to determine the maximum slope direction from the rotated data. Then, iterate through each point in the rotated cloud, extract profiles at a fixed distance, and plot these extracted profiles which makes sure that profiles are extracted along vertical planes oriented according to the maximum slope of the cloud.

I just tried to execute the code utilizing your XYZ.mat file but my MATLAB session kept disconnected. Please see attached.

So, I need your help to execute the code and let me know if this helps resolve your problem. Thanks for your cooperation.

[Elisa]

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]

Hi Elisa,

Please see Avni response below. She has already provided the solution you are seeking.

Answer 1

[Avni Agrawal]

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
% Load the data
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-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]

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

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