How to calculate angle difference between 3D vectors and conditionally plot differences

I have a 3D image of an object with holes in it. The holes appear more as indents when plotting the image - this is acceptable for me. When looking at the surfnorm and quiver3 plots of the image, I noticed the majority of the surface vectors appear to be at the same angle, while the areas with indents, are at varying angles. I have not been able to verify if this is true or not, however I am working on the assumption that it is true.

I have calculated the normal vectors [Nx, Ny, Nz] using surfnorm, and I have attempted to calculate the angles between each and every vector. The results are not really what I was going for though - the plots look similar to a quiver3 plot of the gradients which makes sense in retrospect.

Below is what I have come up with so far.

% Importing data from .xyz
Data = Drill2Edited;
x = Data(:,1);
y = Data(:,2);
z = Data(:,3);

% Defining step size
xv = linspace(min(x), max(x), 100);
yv = linspace(min(y), max(y), 100);
[X,Y] = meshgrid(xv, yv);
Z = griddata(x,y,z,X,Y);

%% Stem Plots
figure()
subplot(2, 1, 1)
stem3(X, Y, Z)
subplot(2, 1, 2)
stem3(x, y, z)
grid on
xlabel('x')
ylabel('y')
zlabel('z')

%% Surface Plot
figure()
surf(X, Y, Z);
grid on
xlabel('x')
ylabel('y')
zlabel('z')

%% Surface Normal Calculations and Vector Plots
figure()
subplot(2,1,1)
surfnorm(X,Y,Z)
[Nx,Ny,Nz] = surfnorm(X,Y,Z);
subplot(2,1,2)
quiver3(X,Y,Z,Nx,Ny,Nz);
xlabel('x')
ylabel('y')
zlabel('z')

%% Angle Calculations
[theta] = zeros(100);

% All three matrices have the same NaN/non-NaN positions, therefore I am assuming the
% positions are transferrable from matrix to matrix and the positions only
% need to be found in one matrix
% Note: I am only referring to the cell positions, not the values stored
% within those cells/positions

% Creates a for loop iterating through each position in the Nx matrix
for k = 1:numel(Nx)
    pos = k;
    [r3]=[Nx(pos),Ny(pos),Nz(pos)];
    
    % Uses position loop to find neighbors of every position
    for j = pos
        szx = size(Nx);
        [rx,cx] = ind2sub(szx,pos);
        neighx(1:8,1:2) = [rx+[-1;0;1;-1;1;-1;0;1] cx+[-1;-1;-1;0;0;1;1;1] ];
        neighx = neighx(all(neighx,2) & neighx(:,1) <= szx(1) & neighx(:,2) <= szx(2),:);
        idx = (neighx(:,2)-1)*szx(1) + neighx(:,1); %positions of neighbor        
    end
    
    % Values of selected positions' neighbors
    for i =1:1:length(idx)
            [r4]=[Nx(idx(i)),Ny(idx(i)),Nz(idx(i))];
            
            % Calculating Theta
            ps = idx(i);
            c = cross(r3,r4);
            c1 = abs(c);
            res = real(acos((dot(r3,r4))/sqrt(c1(1)^2+c1(2)^2+c1(3)^2)));

            if (~isnan(res) && res>0)
                if (res>theta(ps))
                    theta(ps) = res;
                end
            end             
    end
end

figure()
plot3(theta, Y, Z);
plot3(X, Y, theta); % no indents appear which leads me to believe XY isn't the correct plane to plot on
return

My intended result is a plot with the indented areas and the edges (as I expect the edges to also have changes in the vector angle) highlighted in some way. Eventually, I would like to add a threshold for the vector angle difference but that may be a problem for another time.

I'm using MatLab 2021a

The coordinate points are stored in an xyz file ('Drill2Edited') but it is too large to attach to this question. Instead I have attached 'Drill.mat' with the coordinate data. If another format is preferred let me know.

Answer 1

[Vidhi Agarwal]

Hi Marina,

I understand that you are trying to highlight the indented areas and edges based on the angle differences of the surface normal.

To highlight the indent areas, you can follow the below given steps:

1. Use “surfnorm” to compute the normal vectors at each point on the surface.
2. Calculate the angle between each normal vector and its neighbours using the dot product. The angle can be computed using the arccosine of the dot product of two normalized vectors.
3. Highlight areas where the angle difference exceeds a certain threshold.

The updated code, incorporating the points described above, is provided below:

% Load your data
load('Drill.mat');  % Ensure you have the correct filename

% Extract coordinates from 'Drill2Edited'
x = Drill2Edited(:, 1);
y = Drill2Edited(:, 2);
z = Drill2Edited(:, 3);

% Combine x, y, z into a single matrix
data = [x, y, z];

% Use unique to find unique x, y pairs and their indices
[uniqueXY, ~, idx] = unique(data(:, 1:2), 'rows');

% Preallocate a vector for averaged z values
averageZ = zeros(size(uniqueXY, 1), 1);

% Average z values for each unique x, y pair
for i = 1:size(uniqueXY, 1)
    % Find all z values corresponding to the current x, y pair
    zValues = data(idx == i, 3);
    % Compute the average
    averageZ(i) = mean(zValues);
end

% Use the unique x, y pairs and their averaged z values
xUnique = uniqueXY(:, 1);
yUnique = uniqueXY(:, 2);
zUnique = averageZ;

% Define step size for grid
xv = linspace(min(xUnique), max(xUnique), 100);
yv = linspace(min(yUnique), max(yUnique), 100);
[X, Y] = meshgrid(xv, yv);
Z = griddata(xUnique, yUnique, zUnique, X, Y);

% Surface Normal Calculations
[Nx, Ny, Nz] = surfnorm(X, Y, Z);

% Initialize angles matrix
theta = zeros(size(Nx));

% Compute angles between each vector and its neighbors
for i = 2:size(Nx, 1) - 1
    for j = 2:size(Nx, 2) - 1
        % Current normal vector
        r3 = [Nx(i, j), Ny(i, j), Nz(i, j)];
        
        % Neighboring normal vectors
        neighbors = [
            Nx(i-1, j-1), Ny(i-1, j-1), Nz(i-1, j-1);
            Nx(i-1, j), Ny(i-1, j), Nz(i-1, j);
            Nx(i-1, j+1), Ny(i-1, j+1), Nz(i-1, j+1);
            Nx(i, j-1), Ny(i, j-1), Nz(i, j-1);
            Nx(i, j+1), Ny(i, j+1), Nz(i, j+1);
            Nx(i+1, j-1), Ny(i+1, j-1), Nz(i+1, j-1);
            Nx(i+1, j), Ny(i+1, j), Nz(i+1, j);
            Nx(i+1, j+1), Ny(i+1, j+1), Nz(i+1, j+1);
        ];
        
        % Calculate angles with neighbors
        for k = 1:size(neighbors, 1)
            r4 = neighbors(k, :);
            cosTheta = dot(r3, r4) / (norm(r3) * norm(r4));
            angle = acosd(min(1, max(-1, cosTheta))); % Clamp value to avoid NaN due to floating-point errors
            theta(i, j) = max(theta(i, j), angle);
        end
    end
end

% Plotting the areas with significant angle changes
threshold = 15; % Define a threshold for angle differences
figure;
surf(X, Y, Z, theta, 'EdgeColor', 'none');
colormap(jet);
colorbar;
caxis([0, threshold]);
title('Highlighted Indented Areas Based on Normal Vector Angles');
xlabel('x');
ylabel('y');
zlabel('z');
view(3);
grid on;

For detailed understanding of “surfnorm” and “meshgrid” refer to the following documentation:

• https://www.mathworks.com/help/releases/R2021a/matlab/ref/surfnorm.html
• https://www.mathworks.com/help/releases/R2021a/matlab/ref/meshgrid.html

I hope that was helpful