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

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');