MATLAB Answers

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

17 views (last 30 days)
Steve
Steve on 1 Nov 2019
Commented: Steve on 8 Sep 2020
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 Comments

Sign in to comment.

Accepted Answer

Matt J
Matt J on 1 Nov 2019
Edited: Matt J on 2 Nov 2019
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 Comments

Steve
Steve on 14 Aug 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
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:
  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
Steve on 8 Sep 2020
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.

Sign in to comment.

More Answers (2)

Matt J
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
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
Matt J on 5 Nov 2019
You're very kind. Just so you know, if you like a particular answer, you can award it extra points by up-voting it.
untitled.png

Sign in to comment.


Sherif Omran
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
Steve on 1 Nov 2019
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.

Sign in to comment.

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!