**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

5 views (last 30 days)

Show older comments

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:

### Accepted Answer

Matt J
on 1 Nov 2019

Edited: Matt J
on 2 Nov 2019

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

##### 35 Comments

Steve
on 1 Nov 2019

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?

Steve
on 1 Nov 2019

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

Steve
on 1 Nov 2019

Steve
on 1 Nov 2019

Steve
on 2 Nov 2019

Matt J
on 2 Nov 2019

Edited: Matt J
on 2 Nov 2019

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
on 2 Nov 2019

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
on 2 Nov 2019

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

See the updated version.

Steve
on 2 Nov 2019

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:

- Is this purely graphical or can I store the arrangement of edges/lines as vectors in a file?
- Is it possible to have the node numbers correspond to the indices in the fpep.mat file?
- 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
on 2 Nov 2019

Edited: Matt J
on 2 Nov 2019

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!

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
on 2 Nov 2019

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
on 2 Nov 2019

Edited: Matt J
on 2 Nov 2019

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
on 2 Nov 2019

Matt J
on 2 Nov 2019

Edited: Matt J
on 3 Nov 2019

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
on 3 Nov 2019

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

Steve
on 3 Nov 2019

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
on 3 Nov 2019

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
on 3 Nov 2019

Edited: Matt J
on 3 Nov 2019

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
on 3 Nov 2019

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
on 4 Nov 2019

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
on 30 Jul 2020

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
on 30 Jul 2020

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
on 30 Jul 2020

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
on 30 Jul 2020

Edited: Catalytic
on 30 Jul 2020

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
on 31 Jul 2020

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
on 31 Jul 2020

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
on 11 Aug 2020

Catalytic
on 12 Aug 2020

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
on 13 Aug 2020

Catalytic
on 14 Aug 2020

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
on 14 Aug 2020

Edited: Steve
on 8 Sep 2020

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
on 21 Aug 2020

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:

- All boundary edges.
- That any remaining vertex that has two remaining edges, is also a boundary edge.
- 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
on 8 Sep 2020

### More Answers (2)

Matt J
on 5 Nov 2019

Edited: Matt J
on 5 Nov 2019

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.

##### 2 Comments

Steve
on 5 Nov 2019

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
on 5 Nov 2019

Sherif Omran
on 1 Nov 2019

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.

##### 1 Comment

Steve
on 1 Nov 2019

### See Also

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)