Can anyone tell us why this graph failed to close all the shapes in the interior of the graph? Also, can someone explain why there are dangling edges along the border? Thanks!
How do we have the code eliminate edges whose angles are too dissimilar? Also, why are we getting dangling edges in our output graph?
2 次查看（过去 30 天）
显示 更早的评论
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.
Aforementioned code below:
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
VE1=ArcPoints(:,2,:);
VE2=ArcPoints(:,3,:);
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
%% 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);
回答（0 个）
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Computational Geometry 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
您也可以从以下列表中选择网站：
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)