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

4 次查看（过去 30 天）
Marina Baltz
Marina Baltz 2021-10-15
回答： Vidhi Agarwal 2024-8-28
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.

请先登录，再回答此问题。

回答（1 个）

Vidhi Agarwal
Vidhi Agarwal 2024-8-28
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:
I hope that was helpful

类别

在 Help Center 和 File Exchange 中查找有关 Point Cloud Processing 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by