How should I modify this code so that it will discard the correct Delaunay Triangulation (DT) edges?

17 ビュー (過去 30 日間)
Steven D
Steven D 2022 年 6 月 30 日
編集済み: Steven D 2022 年 7 月 1 日
The attached code (TripletGraph2) was written a while back with the help of a few experts in this community. It uses a routine that discards delaunay triangulated edges, which fail to meet a maximum-distance critereon. I need to change this to a routine based on maximum-angle-difference between the two vectors (as shown in the attached image):
This is really only half of the picture; C2 actually has a vector coming off of it in the opposite direction of v1. So there will actually be two angles (i.e., alpha1 and alpha2), as opposed to just the one (alpha). Also, alpha1 and alpha2 should be very close numerically; the maximum difference between these two angles is the criterion I need to use to discard unwanted edges.
Any help you could offer would be very much appreciated.
Relevant Portion of the Code:
classdef TripletGraph2
properties
C; %triplet center points
T; %triplet endpoints
Cgraph %The undirected graph of the center points
ArcPoints %2x4xN list of arc points CenterA-EndpointA-EndpointB-CenterB
end
methods
function obj=TripletGraph2(fpep)
% maxDistLine=20;
maxLineLength=1800;
CTDist=25;
% maxDistLineDiff=7;
l2norm=@(z,varargin) vecnorm(z,2,varargin{:});
C=fpep(:,7:8); %Central points
T=reshape( fpep(:,1:6) ,[],2,3); %Triplet end-points Nx2x3
T=(T-C)./l2norm(T-C,2)*CTDist + C; %Normalize center to endpoint distances
DT=delaunayTriangulation(C); %Delaunay triangulation of centers
E=edges(DT); %Edge list
C1=C(E(:,1),:); %Center points listed for all DT edges
C2=C(E(:,2),:);
L=l2norm(C1-C2,2); %Center-to_center lengths for DT edges
if L<500
maxDistLine=30;
maxDistLineDiff=7;
elseif (500<L)&(L<1000)
maxDistLine=35;
maxDistLineDiff=7;
else
maxDistLine=20;
maxDistLineDiff=7;
end
U=(C2-C1)./L;
S1=C1+CTDist*U;
S2=C2-CTDist*U; %Synthetic endpoints inserted along DT edges
%%Pass 1 - throw away DT edges whose nodes aren't both sufficiently close to
%%any endpoint
clear D1 D2
for j=3:-1:1
D1(:,:,j)=pdist2(T(:,:,j),S1,'Euclidean','Smallest',1);
D2(:,:,j)=pdist2(T(:,:,j),S2,'Euclidean','Smallest',1);
end
% junk=min(D1,[],3)>maxDistLine & min(D2,[],3)>maxDistLine; %discard these DT edges
junkmaxDistLineDiff=abs(abs(min(D1,[],3)) - abs(min(D2,[],3)))>maxDistLineDiff; %discard these DT edges
% E(junk,:)=[];S1(junk,:)=[];S2(junk,:)=[]; %corresponding discards in other arrays
% L(junk)=[];
E(junkmaxDistLineDiff,:)=[];S1(junkmaxDistLineDiff,:)=[];S2(junkmaxDistLineDiff,:)=[]; %corresponding discards in other arrays
L(junkmaxDistLineDiff)=[];
%%all data used in statistics (i.e., numsides) must only come from valid edges
%%identify all edges whose vertices have less than 3 edges as a boundary edge
%%also Identify edges whose vertices have 3 edges, but are connected to an
%%edge that has already been identified as a boundary edge
%%Pass 2 - keep only DT edges **minimally** distant from an endpoint at both
%%ends
clear EID1 EID2
for j=3:-1:1
[D,I]=pdist2(S1,T(:,:,j),'Euclidean','Smallest',1);
I(D>maxDistLine)=nan;
EID1(:,:,j)=I;
[D,I]=pdist2(S2,T(:,:,j),'Euclidean','Smallest',1);
I(D>maxDistLine)=nan;
EID2(:,:,j)=I;
end
Eset=1:size(E,1);
keep=ismember(Eset,EID1) & ismember(Eset,EID2) & L.'<=maxLineLength; %Keep these DT edges
E=E(keep,:);
L=L(keep);
S1=S1(keep,:); S2=S2(keep,:);
% keep2=abs(ismember(Eset,EID1) - ismember(Eset,EID2))<=maxDistLineDiff & L.'<=maxLineLength; %Keep these DT edges
% E=E(keep2,:);
% L=L(keep2);
% S1=S1(keep2,:); S2=S2(keep2,:);
EdgeTable=table(E,L,'VariableNames',{'EndNodes','Weights'});
obj.Cgraph=graph(EdgeTable);
obj.C=C;
obj.T=T;
%%Pass 3 - All spurious DT edges should be gone. Do one last
%%distance comparison to associate endpoints with edges pass
C1=C(E(:,1),:);
C2=C(E(:,2),:);
Tr=reshape(T(:,:).',2,[]);
ArcPoints=nan(2,4,obj.Cgraph.numedges);
ArcPoints(:,1,:)=C1.';
ArcPoints(:,4,:)=C2.';
[~,I1]=pdist2(Tr.',S1,'Euclidean','Smallest',1);
[~,I2]=pdist2(Tr.',S2,'Euclidean','Smallest',1);
ArcPoints(:,2,:)=Tr(:,I1);
ArcPoints(:,3,:)=Tr(:,I2);
obj.ArcPoints=ArcPoints; %All 4-point data for each edge
VE1=ArcPoints(:,2,:);
VE2=ArcPoints(:,3,:);
% VEW = [VE1;VE2];
% writematrix(VEW,'VEW.xls');

回答 (0 件)

カテゴリ

Help Center および File ExchangeDelaunay Triangulation についてさらに検索

製品


リリース

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by