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

Question

Steve 2019-11-1

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/488755-can-someone-propose-some-code-that-will-connect-the-dots-to-produce-the-correct-geometric-shapes

评论： Steve 2020-9-8

fpep.mat

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:

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2019-11-1

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/488755-can-someone-propose-some-code-that-will-connect-the-dots-to-produce-the-correct-geometric-shapes#answer_399415

编辑：Matt J 2019-11-2

在 MATLAB Online 中打开

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

35 个评论
显示 33更早的评论隐藏 33更早的评论

Steve 2019-11-2

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.

Steve 2019-11-2

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 2019-11-2

编辑：Matt J 2019-11-2

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.

Catalytic 2020-7-30

编辑：Catalytic 2020-7-30

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

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 2020-8-12

在 MATLAB Online 中打开

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 2020-8-13

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!

Steve 2020-8-14

编辑：Steve 2020-9-8

在 MATLAB Online 中打开

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 2020-8-21

在 MATLAB Online 中打开

missing-edges-matlab.PNG

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 2020-9-8

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 2019-11-5

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/488755-can-someone-propose-some-code-that-will-connect-the-dots-to-produce-the-correct-geometric-shapes#answer_399847

编辑：Matt J 2019-11-5

在 MATLAB Online 中打开

TripletGraph.m

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 个评论
显示无隐藏无

Steve 2019-11-5

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 2019-11-5

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 2019-11-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/488755-can-someone-propose-some-code-that-will-connect-the-dots-to-produce-the-correct-geometric-shapes#answer_399391

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Steve 2019-11-1

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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

35 个评论
显示 33更早的评论隐藏 33更早的评论

更多回答（2 个）

2 个评论
显示无隐藏无

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

35 个评论 显示 33更早的评论隐藏 33更早的评论

更多回答（2 个）

2 个评论 显示 无隐藏 无

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

35 个评论
显示 33更早的评论隐藏 33更早的评论

2 个评论
显示无隐藏无

1 个评论
显示 -1更早的评论隐藏 -1更早的评论