Can someone propose some code that will "connect the dots" to produce the correct geometric shapes (i.e., hexagons, pentagons, rectangles) from these points?

Hi,

I'm wondering if someone could propose some clever code that will automatically connect the dots (as shown in the fabricated images below) to make the correct geometric shapes below from the given points in the file (fpep.mat)? Note: the edges are actually arcs that emanate from the center points (coordinates for these center points are contained in the last two colomns of fpep.mat). Important note: these arcs must also pass through the endpoints of the triplets (triplets = the 3 points surrounding each center point).The coords of each endpoint are contained in the first 6 columns of fpep.mat). I'm excited to see what you can come up with. Thanks in advance for your help!

Center points (shown with red crosshairs in image above) are connected to each neighboring center point to produce the geometric shapes as shown in the image below:

Answer 1

[Matt J]

This seems to do it. For the display part, it uses plotpts2d, which you've seen before.

maxDistLine=10;     %User tolerance settings
maxLineLength=1500;

P=reshape(fpep(:,7:8).',2,1,[])/1000; P(3,:)=1;
V=reshape(fpep(:,1:6).',2,3,[])/1000; V(3,:)=1;
 
DT=delaunayTriangulation(P(1:2,:).');
E=edges(DT);

P1=P(:,:,E(:,1));
P2=P(:,:,E(:,2));
V1=V(:,:,E(:,1));
V2=V(:,:,E(:,2));

L=cross(P1,P2,1);  
L=L./vecnorm(L(1:2,:,:),2,1);

sdot=@(a,b) squeeze(sum(a.*b,1));
dL=abs( sdot(L,[V1,V2]) ) ;
dL1=dL(1:3,:); dL2=dL(4:6,:);

D=1000*vecnorm(P1(1:2,:)-P2(1:2,:),2,1).';


d=maxDistLine/1000;
Q1=sdot(P2-P1,V1-P1);
Q2=sdot(P1-P2,V2-P2);
C1=dL1<d & Q1>0;
C2=dL2<d & Q2>0;


keep=sum(C1)==1 & sum(C2)==1 & (D<maxLineLength).';

E=E(keep,:);
P1=P1(:,:,keep); P2=P2(:,:,keep); 
V1=V1(:,:,keep); V2=V2(:,:,keep); 
dL1=dL1(:,keep); dL2=dL2(:,keep);
Q1=Q1(:,keep); Q2=Q2(:,keep);
D=D(keep,:);

EdgeTable=table(E,D,'VariableNames',{'EndNodes','Weights'});

G=graph(EdgeTable);


J=find(sum(adjacency(G),2)>3).';

cut=cell(1,numel(J));
for k=1:numel(J)
    
    j=J(k);
    ej=find(any(E==j,2));
    
    Ej=E(ej,:);
    
    dj=dL1(:,ej).*(Ej(:,1)==j).' + dL2(:,ej).*(Ej(:,2)==j).';
    qj=Q1(:,ej).*(Ej(:,1)==j).' + Q2(:,ej).*(Ej(:,2)==j).';
    dj(qj<=0)=inf;
    
    [val,keep]=min(dj,[],2);
    keep=keep(isfinite(val));
    ej(keep)=[];
    
    cut{j}=ej(:).';
    
end

G=G.rmedge([cut{:}]);



    N=size(P,3);
    h=plot(G,'XData',P(1,:)*1000,'YData',P(2,:)*1000,...
           'Marker','o','MarkerSize',5,'NodeColor','k','NodeLabel',(1:N));
    hold on
    plotpts2d([],reshape(V(1:2,:)*1000,2,[]))
    hold off

[Steve]

Thanks so much Matt; it looks great!

I named the file arc_plot and ran it, but I can't seem to get it to work on my end...I'm getting an error that says:

Undefined function or variable 'P'.

Error in arc_plot (line 3)

P1=P(:,:,E(:,1));

Any idea what I'm doing wrong?

[Matt J]

My copy/paste missed the first 2 lines. Try now.

[Steve]

I re-ran it with the 2 new lines, but now I'm getting an error that says:

Index in position 2 exceeds array bounds (must not exceed 1).

Error in arc_plot (line 6)

P2=P(:,:,E(:,2));

[Matt J]

Try now.

[Steve]

Thanks again Matt. It worked! I'm just getting a few anomalies (see attached screenshot). The endpoint "triplets" are not showing in the plot (even when I zoom all the way in); only the center points show, and there are some lines that were created that should not be there (marked by stars). Any idea how I could fix this?

[Steve]

Just a follow up for clarification: No nodes can have a degree other than degree 3 (because of the triplets).

[Steve]

Correction: The nodes along the outer borders could have degrees less than 3 but all of the internal nodes should all have degrees of precisely 3. No nodes (anywhere) should ever have a degree higher than 3. I'm not sure how to have this controlled in the code. Any further help you could offer would be greatly appreciated.

[Matt J]

The endpoint "triplets" are not showing in the plot (even when I zoom all the way in);

I'm guessing you do not have plotpts2d anymore and an error was thrown when you tried to invoke it. You can get it again from the link I gave you in my initial answer.

there are some lines that were created that should not be there (marked by stars).

The udpated version will remove the anomolous edges [562,595] and [609,636] . The edge [633,925] I argue is not an anomaly under the connection criteria that you've stipulated so far. That edge line does pass through a vertex at each of these triplets and therefore should be included.

[Steve]

I still have and have been running your plotpts2d file with no error thrown pertaining to that. However, now, when I run the code, for some reason, it's throwing an error that says:

Undefined function or variable 'iseven'.

Error in plotpts2D (line 26)

if ~iseven(nargin)

Error in arc_plot2 (line 53)

plotpts2D([],reshape(V(1:2,:)*1000,2,[]))

With regard to the long connectors/lines connecting 633 to 925, and the other border edge connectors, Could we set a distance/length limit so those are disqualified?

Thanks for your continued effort--it is truly appreciated.

[Matt J]

Could we set a distance/length limit so those are disqualified?

See the updated version.

[Steve]

Thanks again Matt! You are awesome! Your modification worked perfectly to remove the erroneous edges/lines.

By the way, after much effort, I finally figured out how to fix the error I was getting (pertaining to the plotpts2d file) as described in my previous comment. It seems that the command "iseven" used in plotpts2d (i.e., if ~iseven(nargin)) was not a recognized function in my version/configuration of Matlab, and since it was also negated (~), I just replaced it with "mod(nargin,2)==1" and everything worked just fine.

I just have a few other questions for you:

1. Is this purely graphical or can I store the arrangement of edges/lines as vectors in a file?
2. Is it possible to have the node numbers correspond to the indices in the fpep.mat file?
3. In the original question, I had specified that the edges must pass through the triplet end points (and hence must actually be arcs), is there any way to do this with this method you have shown?

In any event, thank you again so much for all your help! You are (by far) the most talented Matlab MVP I have ever worked with.

[Matt J]

Is this purely graphical or can I store the arrangement of edges/lines as vectors in a file?

You can save the graph object, G, which contains all the connection information. You should read up on the capabilities of graph objects, as they allow you to do many things!

https://www.mathworks.com/help/matlab/ref/graph.html#d117e525223

Is it possible to have the node numbers correspond to the indices in the fpep.mat file?

That is already the case.

In the original question, I had specified that the edges must pass through the triplet end points (and hence must actually be arcs), is there any way to do this with this method you have shown?

Not sure. It isn't even clear what the connection rules would be in that case. Once you allow the connections to be arbitrary curves, there are infinite ways to connect the points while passing through the vertices.

[Steve]

Thanks again for your response Matt! Regarding your first answer: That is awesome, I will definitely read up on that. By the way, I just noticed that one of the vertices (22) failed to connect to its paired vertex (58)--not sure what happened there (see last image below).

Regarding the second answer, I just realized that the plot (of the original fpep points) was transposed.

Regarding your third answer:

"Not sure. It isn't even clear what the connection rules would be in that case. Once you allow the connections to be arbitrary curves, there are infinite ways to connect the points while passing through the vertices."

Isn't it true that a circle is fully defined by three points on its circumference? Therefore, the arc of a circle passing through 4 points is not arbitrary. I have double-checked on this here:

Below are some images for clarification--the first image shows the arc of a unique circle passing through all 4 points; the second image shows an arc through two chosen vertices and their corresponding end points, and also includes a superimposed screenshot of the "double-check" for these particular points:

Please let me know what you think. Thanks again for all your help!

[Matt J]

By the way, I just noticed that one of the vertices (22) failed to connect to its paired vertex (58)--not sure what happened there (see last image below).

To me, they don't look like they are supposed to be connected under our current rules. The rule followed by the current code is that two triplets A and B are connected (by a line segment) if and only if the line segment connects the centers of the triplets and 2 of their vertices in the order CenterA-VertexA-VertexB-CenterB. The tolerance parameter maxDistLine determines how close a vertex needs to be to the line segment in order to be considered colinear with CenterA and CenterB. Pair (22,58) does not satisfy this connection rule, or at least not within the tolerance parameter value currently chosen. You would have to increase maxDistLine rather generously to have that pair joined and that would probably create many unwanted connections elsewhere (but you could always try).

Isn't it true that a circle is fully defined by three points on its circumference? Therefore, the arc of a circle passing through 4 points is not arbitrary.

There are a few problems with this.

First, there are more kinds of curves in the world than just circles (parabolas, ellipses, splines, etc...) and you have not mentioned up until now that the connection curves should be circular arcs.

Second, while 3 points uniquely define a circle, a circle joining 4 arbitrary points generally does not exist. For example, 4 colinear points will never all lie on a common circle. Since most of the points you are connecting are colinear, or very nearly so, circular arcs would be an odd choice for joining them.

Third, this choice does not make it unambiguously clear what the connection rule should be. There may be multiple ways of choosing which combinations of triplets should be joined. There also may be multiple choices of the circle that connects the same 2 triplets. In the image below, I show a second circle that might be used to connect the 2 triplets in your example above. Keep in mind that while the amber arc doesn't contain all 4 points exactly, it does to within pretty good tolerances, and we must always apply tolerances whether we are connecting the triplets with lines, circles, or whatever.

[Steve]

Wow. Thanks again for all your input Matt. I hear what you are saying; specifically speaking, with respect to your second response, the circle could come from the back side. However, this would not produce the geometric shapes (polygons) we desire. Therefore, the chords generated by your code, should necessarily make the smallest angle with the circular arc that passes through the four points. One way to look at it is that the arc passing from one vertex to the next has the initial condition (initial slope) defined by the angle (alpha) that exists between the chord and the "endpoint vector" (i.e., line from the vertex to the corresponding endpoint) as shown in the image below:

[Matt J]

Just to be clear, what you call endpoints, I've been calling vertices. They are vertices of the triangle formed by the triplets...

Regarding your diagram, though, what happens when the alpha at the two triplets are different? And if one is positive and one is negative, i.e., the endpoints are on different sides of the line connecting the centerslike in this section of your drawing?

How do you join them in that case? And what are the rules supposed to be on how big alpha1,alpha2, and alpha1-alpha2 can be before being disqualified as a connection?

And after you've established all these rules, what form is final output of the code supposed to take? An array listing the connected centers and endpoints? We have that already. What has changed now that we're contemplating arc connections rather than lines?

[Image Analyst]

If you have 4 points, you can fit a cubic. I bet it would look pretty reasonable.

[Steve]

Matt,

Well, in theory, that should never happen. If you have found an instance where that is the case, then it is the result of a an error in processing the original image.

By the way, I have an important note to make that I failed to mention before:

The length of the "endpoint vector" is not a factor here, rather it is the angle that the "endpoint vector" makes with the chord that matters. In processing the original image, each line that was used to generate each of the endpoints (see image of "tripods" below) should have been trimmed back so that all the endpoints in each triplet are at the same distance from the common center points. So, fitting an arc through the current end points would be pointless (n.p.i.). I believe this is the source of the error that would be causing anomalies like in the 22/58 case, because the lower endpoint is abnormally far from its center point (see second image below).

[Steve]

Hi Image Analyst,

They must be arcs of circles (in theory). Why can't a circular arc be fit through four points? Think of the arc being fit through any three of the four points, then it will automatically fit through any other three points of the four. This is true because the points are not arbitrary, and did originally belong to an arc of a circle.

Best,

Steve

[Matt J]

Think of the arc being fit through any three of the four points, then it will automatically fit through any other three points of the four.

That is only apparent now that we know the center-to-endpoint distance can be assumed the same for all triplets. Were this not the case, it would not be true in general. Four arbitrary points will not generally be fit by a circle.

There is still the question though of what output you are still looking for from the code. If it is the listing of connected center points, you already have that. That does not change, regardless of whether we are connecting things with lines or circular arcs.

[Steve]

There is still the question though of what output you expect from the code. If it is the listing of connected center points, you already have that. That does not change, regardless of whether we are connecting things with lines or circular arcs.

How would I actually go about connecting the dots with circular arcs instead of straight lines? In addition, I would really love to see comments in your code so I could follow along. I am still amazed at how well it works and how quickly you whipped it up. Thanks again Matt!

[Steve]

Hi Matt,

If it's not too much trouble could you possibly apply comments to your code (above)? I am having trouble figuring out what was done and I would really like to learn.

Thanks,

Steve

[Steve]

Hi People,

Could anyone help solve the issue we are having with the code proposed (above) by MattJ? Seems we keep getting a lot of discontinuous/incomplete shapes in the output graph--as seen in the attached picture. Red dots were added to show where the lines failed to connect. Thanks for any help you could offer.

Best,

Steve

[Catalytic]

Did you play with the tolerance parameters? I see a few in the first lines of the code

maxDistLine=10;     %User tolerance settings
maxLineLength=1500;

[Steve]

Yes, I did try playing with those two parameters but to no avail. So far, I tried:

maxDistLine=15;

maxLineLength=1800

then,

maxDistLine=20;

maxLineLength=1800

then,

maxDistLine=15;

maxLineLength=2000

Any idea how we could fix this?

[Catalytic]

And you get the exact same result? What if you increase maxDistLine more significantly, like to 100?

Anyway, I don't understand why some of the regions are considered to have missing connections. Why would the green line be considered a missing connection in this blue section, for example? There doesn't appear to be any vertex to connect to on the south-east side of the polygon.

And similarly, where would the missing connections be drawn in this yellow region? Would you again try to symmetrically cut it with a line running north-west to south-east? But if so, I again don't see any vertex on the north-west side for the cutting line to join to.

And while we're on the subject, what criterion do you use to decide when two vertices should be joined or not? From your question post, the lines are supposed to join vertices whose 'endpoints' are approximately colinear with the joining line, it appears. But since there will never be exact collinearity, what kind of tolerance do you use to decide whether they are "colinear enough"?

[Steve]

And you get the exact same result? What if you increase maxDistLine more significantly, like to 100?

Anyway, I don't understand why some of the regions are considered to have missing connections. Why would the green line be considered a missing connection in this blue section, for example? There doesn't appear to be any vertex to connect to on the south-east side of the polygon.

I got slightly different results (see attached). When I go out to say 100, other trhings get much worse (see attached).

You raised another good point...those blue-region "twin" shapes should not be there at all. They should not be able to form any complete shapes in that region because that area was where the original border was cut off.

And similarly, where would the missing connections be drawn in this yellow region? Would you again try to symmetrically cut it with a line running north-west to south-east? But if so, I again don't see any vertex on the north-west side for the cutting line to join to.

Here, there was supposed to be a tripet set that divides the area--for some reason it did not survive the transfer to code.

And while we're on the subject, what criterion do you use to decide when two vertices should be joined or not?

This was determined by MattJ's initial version of the code.

From your question post, the lines are supposed to join vertices whose 'endpoints' are approximately colinear with the joining line, it appears. But since there will never be exact collinearity, what kind of tolerance do you use to decide whether they are "colinear enough"?

The tolerance is set by MattJ's next version of the code, which compares the angle made between the line from the centerpoint to an endpoint and the line connecting the two neighboring centerpoints.

There are also other errors that I noticed along the borders--As alluded to before, there should not always be complete shapes there, but also, there should not be any stray lines emanating from the complete shapes.

I'm not sure if MattJ had any other knobs that he was planning to build into the code that would prevent the types of errors that are showing up. If he didn't, can you think of a way to prevent these errors? Thanks for any help you could offer.

[Catalytic]

This was determined by MattJ's initial version of the code.

I wasn't asking what MattJ's idea of the correct criterion was. I was asking yours. If you cannot articulate a criterion of your own, then as a naive onlooker, I have no basis for seeing MattJ's omissions as "errors". Have you zoomed in to look at the endpoints that you think should be joined by lines? Maybe they are not very colinear at all...

[Steve]

Well, after addressing some anomalies found in the original input image, we were able to fix all but one missing edge--this was found to be the result of that particular edge not being able to meet the prescribed tolerances as currently set up in the code. Also, upon inspection, we realized that the code is not using the right criteria (closeness between the starting and ending angles, i.e., theta1 and theta2, of each potential arc) to choose and/or check the correct centerpoint/node neighbor to create an edge with. We always felt that this was a necessary criterion and thought we had mentioned this early on. In any event, would anyone be able to suggest the best way to get this feature into the code? Thanks in advance for your help!

[Catalytic]

Well, this section of the code seems to be the part that governs which potential arcs will be kept and which will be rejected.

d=maxDistLine/1000;

Q1=sdot(P2-P1,V1-P1);
Q2=sdot(P1-P2,V2-P2);
C1=dL1<d & Q1>0;
C2=dL2<d & Q2>0;

keep=sum(C1)==1 & sum(C2)==1 & (D<maxLineLength).';

It appears that each potential connection examines points P1 and P2 and associated endpoints V1 and V2. The criterion dL1<d requires that the distance of V1 to the line [P1,P2] must be less than d, in other words that V1 and V2 be nearly colinear with P1 and P2. The quantity acos(Q1) seems to be your theta1 angle and the requirement Q1>0 ensures that this angle is acute. So, if you have additional criteria on theta1=acos(Q1) and theta2=acos(Q2) you can impose them in the expression for keep.

It doesn't appear to me, however, that you can just impose a similarity criterion between theta1 and theta2 and throw the other criteria away. It looks like you need a colinearity criterion as well, because there are bound to be some potential arcs with similar, but very large theta that cut the hexagons in half. I can't imagine you would want to keep those if the general honeycomb pattern is what you're aiming to get.

[Steve]

Thank you for your input Catalytic. I believe the missing edge was the result of the primary straight edge being discarded via the tolerance settings, which happens way before the introduction of the arcs in the code. I have included a screenshot of the area (arrow is pointing to the missing black edge), with all of the initial Delaunay Triangulation edges (DT edges, i.e., the straight lines) plotted for clarity. Please let me know what you think. By the way, I agree with what you said, "It doesn't appear to me, however, that you can just impose a similarity criterion between theta1 and theta2 and throw the other criteria away. It looks like you need a colinearity criterion as well". I just think that inclusion of a similar-angles test would make the code a bit more robust. However, this would have to be implemented earlier in the code than where we currently define theta1 and theta2 though; I'm not sure what the best way to do this is. Would you happen to have any suggestions? In any event, thanks again!

[Catalytic]

I don't see anything in your comment that affects my earlier suggestion that you modify these lines. They are the ones that control which edges from the Delaunay triangulation are kept, and which are thrown away.

d=maxDistLine/1000;

Q1=sdot(P2-P1,V1-P1);
Q2=sdot(P1-P2,V2-P2);
C1=dL1<d & Q1>0;
C2=dL2<d & Q2>0;

keep=sum(C1)==1 & sum(C2)==1 & (D<maxLineLength).';

[Steve]

Thanks again for your response Catalytic. First, I have to apoloogize. I failed to mention that we had some other issues early on with the code originally proposed by MattJ; after not being able to reach Matt again for assistance, we were offered a different approach in the code by another amazing MathWorks Technical Support Department staff member, Zhikang Zhang (please see his version below). Now, we are unable to reach Zhikang for further help. So, any input on the issues we are having with this code would be very much appreciated.

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
        
            
            

            
        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

[Steve]

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.

Main code below:

%% 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]

By the way, one very important concept that needs to be in the code is that there should never be allowed, a polygon that has edges that turn outward (see attached images) at any point. In other words, all of the polygons must never have edges that (depending upon the direction one traverses along the edges of the polygon) make an angle greater than 180 degrees with any of its conjoined edges (i.e., having a common vertex). If anyone could provide any input on how to do this, that would be great.

Answer 2

[Matt J]

Here is a refinement of my earlier answer which I think performs better. It uses the attached classdef file to create objects representing these special kinds of triplet graphs. You will find the code commenting more detailed than the earlier version.

For display purposes, the class contains methods that will plot the graph joining the center points with either lines or arcs. The usage is simple,

obj=TripletGraph(fpep);

figure(1);
obj.plotarcs;

figure(2);
obj.plotlines;

The plotarcs method joins the points using cubic spline fitting. I think you will find this gives essentially the same as what you would get with circular arc fitting. The 4th and higher order Taylor expansion terms of the circular arcs should be negligibly small for such large radii of curvature as we see here.

[Steve]

Looks amazing Matt! Truly, a big thanks to you for all your help and effort. I wish there was a way (on this site) to commend you properly--not only for your stellar MATLAB skills, but also for your willingness to go above and beyond when helping others. You are (by far) the most knowledgeable and helpful MATLAB expert I have ever come across, here, or anywhere. If you ever need a job recommendation, I would be happy to give you one. Please feel free to contact me anytime.

All the best,

Steve D.

[Matt J]

You're very kind. Just so you know, if you like a particular answer, you can award it extra points by up-voting it.

Answer 3

[Sherif Omran]

Yes, i propose using the least square root used for shortest distance, see the Dijkstra algorithm to connect your dots

Implementation is not that difficult, i am sure you can do it.

[Steve]

Thank you Sherif. I don't think that this would work though because these connector lines are not necessarily the shortest distances. Also, they are actually arcs, not straight lines. However, perhaps you can run some code on this to show otherwise. I would be curious to see what comes of it.