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.