How do we have the code eliminate edges whose angles are too dissimilar? Also, why are we getting dangling edges in our output graph?

Does anyone know a good way to compare similar angles (theta1 and theta2) in this code, where we could set a tolerance on the delta, above which, the edges belonging to the associated angles will be discarded?

In the code, obj.ArcPoints contains all four points needed to define the angles; however, even before the implementation of this criterion on dissimilar angles, some of the edges are being deleted when they shouldn't be. We still do not understand why this is happening.

Also, we seem to be getting dangling, colored, edges (even in MattJ's version too)--along the borders of the MATLAB output graph (see attached). Firstly, we're are not sure why these edges are not black, like the other edges. Secondly, we're not sure why we are even getting this dangling edges and what this means. Lastly, We are not sure if these dangling boundary edges are being counted in the output data. If so, we really need to fix that. To do this, we're thinking that the code needs to identify and/or define:

1. All boundary edges.
2. That any remaining vertex that has two remaining edges, is also a boundary edge.
3. Those vertices which are connected to the proper amount of edges (i.e., three), but where one of those edges is already defined as a boundary edge.

Any help you could offer would be greatly appreciated.

Aforementioned code below:

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=50;
            
            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
            
            U=(C2-C1)./L;
            S1=C1+CTDist*U;   
            S2=C2-CTDist*U; %Synthetic endpoints inserted along DT edges
            
            
            %%Pass 1 - throw away DT edges whos 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
            
            E(junk,:)=[];S1(junk,:)=[];S2(junk,:)=[]; %corresponding discards in other arrays
            L(junk)=[];
 

            
            
            %%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,:); 
            
            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,:);
          
            
            
        end
        
        function varargout=plotlines(obj)
            
            C=obj.C; T=obj.T; %#ok<*PROP>
            G=obj.Cgraph;
            N=numnodes(G); %#ok<*NASGU>
            h=plot(G,'XData',C(:,1),'YData',C(:,2),...
                'Marker','.','MarkerSize',.1,'NodeColor','k');
%             h=plot(G,'XData',C(:,1),'YData',C(:,2),...
%                 'Marker','.','MarkerSize',5,'NodeColor','k','NodeLabel',(1:N));
            hold on
            plotpts_2D([],reshape(T(:,:).',2,[]))
            
            hold off
            if nargout, varargout={h}; end
            
        end
        
%         function plotspline(obj)
%             
%            
%             h=obj.plotlines;
%             h.LineStyle='none';
%             
%             AP=obj.ArcPoints;
%             N=size(AP,3);
%             
%             sq=linspace(0,1,100);
%             X=nan(100,N);
%             Y=nan(100,N);
%             
%             for i=1:N
%                 
%                 C1=AP(:,1,i); C2=AP(:,4,i);
%                 V1=AP(:,2,i); V2=AP(:,3,i);
%                 
%                 L=norm(C2-C1);
%                 U=(C2-C1)/L;
%                 
%                 s=[0,  dot(V1-C1,U)/L ,  dot(V2-C1,U)/L  , 1];
%                 
%                 APi=interp1(s.',AP(:,:,i).',sq.','spline');
%                 X(:,i)=APi(:,1);
%                 Y(:,i)=APi(:,2);
%                 
%                 
%             end
%             
%             hold on
%             plot(X,Y,'-');
%             hold off
%             
%         end
%         
%         function plotquad(obj)
%             
%            
%             h=obj.plotlines;
%             h.LineStyle='none';
%             
%             AP=obj.ArcPoints;
%             N=size(AP,3);
%             
%             sq=linspace(0,1,100)-0.5;
%             X=nan(100,N);
%             Y=nan(100,N);
%             
%             for i=1:N
%                 
%                 C1=AP(:,1,i); C2=AP(:,4,i);
%                 V1=AP(:,2,i); V2=AP(:,3,i);
%                 
%                 L=norm(C2-C1);
%                 U=(C2-C1)/L;
%                 dV1=(V1-C1)/norm(V1-C1);
%                 dV2=(V2-C2)/norm(V2-C2);
%                 theta1=acosd(dot( dV1, U) );
%                 theta2=acosd(dot( dV2, -U) );
%                 theta=(theta1+theta2)/2;
% 
%                 
%                 
%                 W=cross( cross([U;0],[dV1;0]), [U;0]); 
%                 W=W(1:2)/norm(W);
%                 D=L/2;
%                 a=-tand(theta)/2/D;
%                 
%                 t=2*D*sq;
%                 s=polyval([a,0,-a*D^2],t);
%                 
%                XY=(C1+C2)/2+U*t+W(1:2)*s;
%                X(:,i)=XY(1,:);
%                Y(:,i)=XY(2,:);
%                 
%                 
%             end
%             
%             hold on
%             plot(X,Y,'-');
%             hold off
%             
%         end    
        
        function plotcirc(obj)
            
           
            h=obj.plotlines;
            h.LineStyle='none';
            
            AP=obj.ArcPoints;
            N=size(AP,3);
            
            sq=linspace(0,1,100)-0.5;
            X=nan(100,N);
            Y=nan(100,N);
            
            for i=1:N
                
                C1=AP(:,1,i); C2=AP(:,4,i);
                V1=AP(:,2,i); V2=AP(:,3,i);
                
                L=norm(C2-C1);
                U=(C2-C1)/L;
                dV1=(V1-C1)/norm(V1-C1);
                dV2=(V2-C2)/norm(V2-C2);
                theta1=acosd(dot( dV1, U) );
                theta2=acosd(dot( dV2, -U) );
                theta=(theta1+theta2)/2;
                arcradius=norm(C1-C2)/(2*sind(theta));
                             
                
                W=cross( cross([U;0],[dV1;0]), [U;0]); 
                W=W(1:2)/norm(W);
                D=L/2;

                t=2*D*sq;
                Dcot=D*cotd(theta);
                s=sqrt(Dcot.^2-(t.^2-D.^2))-Dcot;
                
               XY=(C1+C2)/2+U*t+W(1:2)*s;
               X(:,i)=XY(1,:);
               Y(:,i)=XY(2,:);
                
                
            end
            
            hold on
            plot(X,Y,'-');
            hold off
            
        end   
    end
end

%% Remove border triangles from image

close all;
clearvars;

% TBimage = imread('new_vertex_triangles2_low.jpg');
% TBW = imbinarize(TBimage);
% TBWI = ~TBW;
% TBWI2 = imclearborder(TBWI);
% imshow(TBWI2);
% title('Objects touching image borders removed');
% imwrite (TBWI2,'border_triangles_removed_by_matlab77.tif');
% 
% %% Find tripods from triangles in cropped foam image
% 
% TRIA = imread('border_triangles_removed_by_matlab77.tif');
% TRIA = ~TRIA;
% p_tripods = bwskel(TRIA);
% imshow(p_tripods);
% imwrite(p_tripods,'plateau_border_tripods.tif');
T = imread('Mask_of_img019_tria_inv_cropped2.tif');
TBW = imbinarize(T);
% TBW = T;
% T = ~T;
% tripods77 = bwskel(TBW);
tripods77 = bwmorph(TBW,'thin',inf);
% tripods77 = bwskel(TBW,'MinBranchLength',1);
% imshow(tripods77);
imwrite(tripods77,'tri_pods.tif');

%% Find endpoints from tripods

IMM = imread('tri_pods.tif');
% TRIP = imbinarize(IMM);
TRIP = IMM;
% TRIP = ~TRIP;
TBW2 = bwmorph(TRIP,'endpoints');
%TBW3 = bwmorph(TRIP,'branchpoints');
imwrite(TBW2,'plateau_border_endpoints.png');
% imwrite(TBW3,'tripodbranchpoints2a.png');

%% Find Triplets from endpoints

IME = imread('plateau_border_endpoints.png');
% B = imbinarize(IM);
% BW = ~B;
thresh=150; %distance threshold
[I,J]=find(TBW2);
X=[I,J];
[D,K]=pdist2(X,X,'euclidean','Smallest',3);  X=X.';
 K=unique(  sort(K).',  'rows').';    K=K(:);
 assert(numel(unique(K))==numel(K) ,'Points not well-separated enough for threshold'  );
Triplets=num2cell(  reshape( X(:,K), 2,3,[] )  ,[1,2]);
Triplets=Triplets(:).';
Cm=cell2mat(cellfun(@(z)mean(z,2),Triplets,'uni',0)   );
% imshow(IM); hold on; plot_pts2D(flipud(Cm)); hold off;

%% find fpep from Triplets

N = length(Triplets);
xA = zeros(1,N);  
xB = zeros(1,N);
xC = zeros(1,N);
yA = zeros(1,N);
yB = zeros(1,N);
yC = zeros(1,N);
xF = zeros(1,N);
yF = zeros(1,N);
vABx = zeros(1,3,N);
vABy = zeros(1,3,N);
vCBx = zeros(1,3,N); 
vCBy = zeros(1,3,N);
A = zeros(1,N);
B = zeros(1,N);
C = zeros(1,N);
xyF = zeros(2,N);
xye = zeros(2,3,N);

for i = 1 : N
  
        xC(i) = Triplets{i}(2,1);
        xB(i) = Triplets{i}(2,2);
        xA(i) = Triplets{i}(2,3);
        yC(i) = Triplets{i}(1,1);
        yB(i) = Triplets{i}(1,2);
        yA(i) = Triplets{i}(1,3);
        
        vABx(i) = xA(i)-xB(i);
        vABy(i) = yA(i)-yB(i);
        vCBx(i) = xC(i)-xB(i);
        vCBy(i) = yC(i)-yB(i);
        
        A(i) = dot(cross([vABx(i) vABy(i) 0],[vCBx(i) vCBy(i) 0]), [0 0 1]); %% normal vector dotted with z hat to check sign
            
    if   A(i) > 0
               
            xA(i) = Triplets{i}(2,1);
            xB(i) = Triplets{i}(2,2);
            xC(i) = Triplets{i}(2,3);
            yA(i) = Triplets{i}(1,1);
            yB(i) = Triplets{i}(1,2);
            yC(i) = Triplets{i}(1,3);
            
          
        
                    xF(i) = xA(i) + ((-2*xA(i) + xB(i) + xC(i) + sqrt(3)*(-yB(i) + yC(i)))*(sqrt(3)*xA(i)^2 + (yA(i) - yB(i))*(xC(i) + sqrt(3)*(yA(i) - yC(i))) + xB(i)*(sqrt(3)*xC(i) - yA(i) + yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - yB(i) + yC(i))))/(2*(sqrt(3)*xA(i)^2 + sqrt(3)*xB(i)^2 + sqrt(3)*xC(i)^2 + 3*xC(i)*(yA(i) - yB(i)) - xB(i)*(sqrt(3)*xC(i) + 3*yA(i) - 3*yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - 3*yB(i) + 3*yC(i)) + sqrt(3)*(yA(i)^2 - yA(i)*yB(i) + yB(i)^2 - (yA(i) + yB(i))*yC(i) + yC(i)^2)));
                    yF(i) = yA(i) + ((sqrt(3)*xB(i) - sqrt(3)*xC(i) - 2*yA(i) + yB(i) + yC(i))*(sqrt(3)*xA(i)^2 + (yA(i) - yB(i))*(xC(i) + sqrt(3)*(yA(i) - yC(i))) + xB(i)*(sqrt(3)*xC(i) - yA(i) + yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - yB(i) + yC(i))))/(2*(sqrt(3)*xA(i)^2 + sqrt(3)*xB(i)^2 + sqrt(3)*xC(i)^2 + 3*xC(i)*(yA(i) - yB(i)) - xB(i)*(sqrt(3)*xC(i) + 3*yA(i) - 3*yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - 3*yB(i) + 3*yC(i)) + sqrt(3)*(yA(i)^2 - yA(i)*yB(i) + yB(i)^2 - (yA(i) + yB(i))*yC(i) + yC(i)^2)));
                    xyF(1,i) = yF(i);
                    xyF(2,i) = xF(i);
                    xye(1,:,i) = [yA(i) yB(i) yC(i)];
                    xye(2,:,i) = [xA(i) xB(i) xC(i)];
                    
                    B(i) = dot(cross([vABx(i) vABy(i) 0],[vCBx(i) vCBy(i) 0]), [0 0 1]); %% checking sign again
            
        
    else
      
                    xF(i) = xA(i) + ((-2*xA(i) + xB(i) + xC(i) + sqrt(3)*(-yB(i) + yC(i)))*(sqrt(3)*xA(i)^2 + (yA(i) - yB(i))*(xC(i) + sqrt(3)*(yA(i) - yC(i))) + xB(i)*(sqrt(3)*xC(i) - yA(i) + yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - yB(i) + yC(i))))/(2*(sqrt(3)*xA(i)^2 + sqrt(3)*xB(i)^2 + sqrt(3)*xC(i)^2 + 3*xC(i)*(yA(i) - yB(i)) - xB(i)*(sqrt(3)*xC(i) + 3*yA(i) - 3*yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - 3*yB(i) + 3*yC(i)) + sqrt(3)*(yA(i)^2 - yA(i)*yB(i) + yB(i)^2 - (yA(i) + yB(i))*yC(i) + yC(i)^2)));
                    yF(i) = yA(i) + ((sqrt(3)*xB(i) - sqrt(3)*xC(i) - 2*yA(i) + yB(i) + yC(i))*(sqrt(3)*xA(i)^2 + (yA(i) - yB(i))*(xC(i) + sqrt(3)*(yA(i) - yC(i))) + xB(i)*(sqrt(3)*xC(i) - yA(i) + yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - yB(i) + yC(i))))/(2*(sqrt(3)*xA(i)^2 + sqrt(3)*xB(i)^2 + sqrt(3)*xC(i)^2 + 3*xC(i)*(yA(i) - yB(i)) - xB(i)*(sqrt(3)*xC(i) + 3*yA(i) - 3*yC(i)) - xA(i)*(sqrt(3)*xB(i) + sqrt(3)*xC(i) - 3*yB(i) + 3*yC(i)) + sqrt(3)*(yA(i)^2 - yA(i)*yB(i) + yB(i)^2 - (yA(i) + yB(i))*yC(i) + yC(i)^2)));
                    xyF(1,i) = yF(i);
                    xyF(2,i) = xF(i);
                    xye(1,:,i) = [yA(i) yB(i) yC(i)];
                    xye(2,:,i) = [xA(i) xB(i) xC(i)];
                    
                    C(i) = dot(cross([vABx(i) vABy(i) 0],[vCBx(i) vCBy(i) 0]), [0 0 1]); %% checking sign again
                    
%                   pA(i) = [yA(i) xA(i)];
%                   pB(i) = [yB(i);xB(i)];
%                   fC(i) = [yC(i);xC(i)];
%                   fF(i) = [yF(i);xF(i)];
                    
    end
    
                    axy = [yA;xA]'; %% endpoint A coords
                    bxy = [yB;xB]'; %% endpoint B coords
                    cxy = [yC;xC]'; %% endpoint C coords
                    fxy = [yF;xF]'; %% Fermat point that goes with above endpoints A, B, & C
                    fpep = [axy bxy cxy fxy]; %% All endpoint coords shown with their corresponding Fermat points
                   
end
  

% plot(xA,yA, 'ro','markersize',5,'MarkerFaceColor','r');
% hold on
% plot(xB,yB, 'ro','markersize',5,'MarkerFaceColor','r');
% hold on
% plot(xC,yC, 'ro','markersize',5,'MarkerFaceColor','r');
% hold on
% plot(xF,yF, 'ro','markersize',5,'MarkerFaceColor','b');
% axis equal
% grid on

%% Find bubble info from fpep

% create the object
obj=TripletGraph2(fpep);
% draw the figure
% obj.plotspline;
obj.plotcirc;
% hold the figure and draw upon the previous figure
axis equal;
hold on;
% get polygons from the object
[polygons, vertices_index] = getPolyshapes(obj);
% get the arc point information
arcpoint = obj.ArcPoints;
% interpolation length
sq=linspace(0,1,100)-0.5;
[polygons,vertexLists]=getPolyshapes(obj);
N=numel(polygons);
output(N).numsides=[];   %Return output as a structure array
output(N).perimeter=[];
output(N).area=[];
% for each polygon, calculate its area and draw it
for i= 1:1:length(polygons)
    p=polygons(i);
    v=vertexLists{i};
% get the vertices of current polygon
vertices = polygons(i).Vertices;
% calculate its approximation area
area(i) = polyarea(vertices (:,1), vertices (:,2));
% plot the current polygon
plot(polygons(i));
axis equal;
hold on;
% for each edge of this polygon, build the circular segment from the arc
for edge = 1:1:size(vertices,1)

    % use start and end points C1, C2 of this edge as a key, find the two
    % bridge points V1 V2 between them
    start_point = vertices(edge, :);
    if edge + 1 > size(vertices,1)
        end_point = vertices(mod(edge+1, size(vertices,1)), :);
    else
        end_point = vertices(edge+1, :);
    end

    temp = arcpoint(:,:,sum(arcpoint(:,1,:)==start_point')==2);
    res = temp(:,:,sum(end_point'==temp(:,4,:))==2);
    if isempty(res)
        temp = arcpoint(:,:,sum(arcpoint(:,1,:)==end_point')==2);
        res = temp(:,:,sum(start_point'==temp(:,4,:))==2);
    end
    if isempty(res)
        error('can not find V1 and V2!!')
    end

    % C1 start point, C2 end point, V1 V2 bridge points
    C1 = res(:,1,:);
    C2 = res(:,4,:);
    V1 = res(:,2,:);
    V2 = res(:,3,:);

    % calcuate by polyarea
    L=norm(C2-C1);
    U=(C2-C1)/L;
    dV1=(V1-C1)/norm(V1-C1);
    dV2=(V2-C2)/norm(V2-C2);
    theta1=acosd(dot( dV1, U) );
    theta2=acosd(dot( dV2, -U) );
    theta=(theta1+theta2)/2;
    arcradius=norm(C1-C2)/(2*sind(theta));
    W=cross( cross([U;0],[dV1;0]), [U;0]);
    W=W(1:2)/norm(W);
    D=L/2;
    t=2*D*sq;
    Dcot=D*cotd(theta);
    s=sqrt(Dcot.^2-(t.^2-D.^2))-Dcot;
    XY=(C1+C2)/2+U*t+W(1:2)*s;
    X = XY(1,:);
    Y = XY(2,:);
    % estimate arc segment as polygon
    circle_segement = polyshape(X, Y);
    plot(circle_segement);
    area_by_polygon = polyarea(X, Y);
    % end here

    % calculate arc segment w/ algebra
    ca = 2* theta;
    r = arcradius;
    circle = pi * r ^ 2;
    triangle = r ^ 2 * sind(ca) / 2;
    area_by_algebra =  circle * ca / 360 - triangle;
    % algebra ends here

    % if V1 V2 live inside the polygonal area, we should substract the circle_segement
    if inpolygon([V1(1),V2(1)], [V1(2),V2(2)], vertices (:,1), vertices (:,2))
        area(i) = area(i) - area_by_algebra;
    % if V1 V2 live outside the polygonal area, we should add the circle_segement
    else
        area(i) = area(i) + area_by_algebra;
    end

end
    output(i).numsides=numsides(p);
    output(i).perimeter=perimeter(p);
    output(i).area=area(i); 
end

% display the area values
disp(area);

[Steve]

Can anyone tell us why this graph failed to close all the shapes in the interior of the graph? Also, can someone explain why there are dangling edges along the border? Thanks!

[Steve]

Is anyone capable of helping us solve these issues? If not, does anyone know how to have this code pause and ask the user which vertices should be chosen to create an edge that it is unable to determine on its own?