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

10 Ansichten (letzte 30 Tage)
Steve am 1 Nov. 2019
Kommentiert: Steve am 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 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Akzeptierte Antwort

Matt J am 1 Nov. 2019
Bearbeitet: Matt J am 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);
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 Kommentare33 ältere Kommentare anzeigen33 ältere Kommentare ausblenden
Steve am 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.
Main code below:
%% Remove border triangles from image
close all;
clearvars;
% 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 = ~TRIA;
% p_tripods = bwskel(TRIA);
% imshow(p_tripods);
% imwrite(p_tripods,'plateau_border_tripods.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
% 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
% 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;
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;
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 am 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.

Melden Sie sich an, um zu kommentieren.

### Weitere Antworten (2)

Matt J am 5 Nov. 2019
Bearbeitet: Matt J am 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 KommentareKeine anzeigenKeine ausblenden
Steve am 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 am 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.

Melden Sie sich an, um zu kommentieren.

Sherif Omran am 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 Kommentar-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden
Steve am 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.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Geometric Transformation and Image Registration finden Sie in Help Center und File Exchange

R2019a

### Community Treasure Hunt

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

Start Hunting!

Translated by