Find cycles in an undirected graph
Show older comments
Hi, I need to find cycles in a graph, exactly as it was asked here (and apparently without fully clear/working solutions!):
- find cycle in array https://ch.mathworks.com/matlabcentral/answers/425321-find-cycle-in-array
- find a cycles in undirected graph https://ch.mathworks.com/matlabcentral/answers/421435-find-a-cycles-in-undirected-graph
Here my example/code:
clear all; clc; close all;
figure('Color',[1 1 1]);
s = [1 1 1 2 3 3 4 4 5 6 7 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26];
G = graph(s,t);
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0];
h = plot(G,'XData',x,'YData',y,'linewidth',2,'MarkerSize',7);
nl = h.NodeLabel;
h.NodeLabel = '';
xd = get(h, 'XData');
yd = get(h, 'YData');
text(xd, yd, nl, 'FontSize',17, 'FontWeight','bold', ...
'HorizontalAlignment','left', 'VerticalAlignment','top')
set(gca,'Fontsize',15,'FontWeight','Bold','LineWidth',2, 'box','on')
% Remove "branches"
xy = [x' y'];
while ~isempty(find(degree(G)==1))
degreeone = find(degree(G)==1);
G = rmnode(G,degreeone);
xy(degreeone,:) = [];
end
Here the corresponding Figure (after removal of "branches"):
My goal would be to find the following 5 cycles as output (i.e. lists of nodes composing each cycle):
- 1-2-3-4-5-6-7-8-1
- 6-7-8-9-10-11-6
- 1-8-9-10-12-14-18-1
- 1-18-19-20-1
- 12-13-15-16-17-18-14-12
Note 1:
This is the Sergii Iglin 's idea, that I found in https://ch.mathworks.com/matlabcentral/fileexchange/4266-grtheory-graph-theory-toolbox
This method is partially working for my purposes.
Unfortunately, the 2nd and 4th cycles are not what I needed/expected.
% Sergii Iglin
% https://iglin.org/All/GrMatlab/grCycleBasis.html
E = table2array(G.Edges);
Output_SI = grCycleBasis(E);
% [my part] From the Sergii Iglin's output to cycles nodes
for i = 1 : size(Output_SI,2)
w = [];
u = E(find(Output_SI(:,i)),:); % edges list
w(1) = u(1,1);
w(2) = u(1,2);
u(1,:) = [];
j = 2;
while ~isempty(u)
[ind,~] = find(u==w(j));
[~,ind2] = ismember(u, u(ind,:), 'rows');
g = u( ind2==1 ,:) ~= w(j);
w(j+1) = u( ind2==1 , g);
u( ind2==1 ,:) = [];
j = j + 1;
end
cycles_SI{i} = w;
end
% Sergii Iglin's results
>> cycles_SI{:}
1 2 3 4 5 6 7 8 1
1 2 3 4 5 6 11 10 9 8 1
1 8 9 10 12 14 18 1
1 8 9 10 12 13 15 16 17 18 1
1 18 19 20 1
Note 2:
This is the Christine Tobler 's idea that I found in https://ch.mathworks.com/matlabcentral/answers/353565-are-there-matlab-codes-to-compute-cycle-spaces-of-graphs
This method is partially working for my purposes.
Unfortunately, the 2nd and 4th cycles are not what I needed/expected.
% Christine Tobler
% https://ch.mathworks.com/matlabcentral/answers/353565-are-there-matlab-codes-to-compute-cycle-spaces-of-graphs
t = minspantree(G, 'Type', 'forest');
% highlight(h,t)
nonTreeEdges = setdiff(G.Edges.EndNodes, t.Edges.EndNodes, 'rows');
cycles_CT = cell(size(nonTreeEdges, 1), 1);
for i = 1 : length(cycles_CT)
src = nonTreeEdges(i, 1);
tgt = nonTreeEdges(i, 2);
cycles_CT{i} = [tgt shortestpath(t, src, tgt)];
end
% Christine Tobler's results
>> cycles_CT{:}
8 7 6 5 4 3 2 1 8
11 10 9 8 1 2 3 4 5 6 11
18 14 12 10 9 8 1 18
18 17 16 15 13 12 10 9 8 1 18
20 19 18 1 20
Note 3:
Methods from Sergii Iglin and Christine Tobler give the same result!
Note 4:
The ideas / FileExchange submissions
- Count all cycles in simple undirected graph version 1.2.0.0 (5.43 KB) by Jeff Howbert
- Count Loops in a Graph version 1.1.0.0 (167 KB) by Joseph Kirk
kindly suggested here
are not working for my case...
Any further idea / suggestion ?
Thanks a lot!
6 Comments
Steven Lord
on 8 Oct 2019
Why isn't 1-2-3-4-5-6-11-10-12-13-15-16-17-18-19-20-1 on your list of cycles?
Sim
on 9 Oct 2019
I found also this idea/method "Polygons From Set of Line segments", explained in
and based on Ferreira A., Fonseca M.J., Jorge J.A. (2003) Polygon detection from a set of lines
The basic idea can be summarised as: "Finding small polygons is the same as searching for a Minimum Cycle Basis"
% Pseudocode of the algorithm
% Minimum basis cycle [3]
% Initialize empty sets F,P
% P=All-Pairs-Shortest-Paths(G)
% for each v ∈ V
% for each (x,y) ∈ E
% if Px,v ∩ Pv,y = {v}
% C=Px,v ∪ Pv,y ∪ (x,y)
% add C to F
% Order-By-Length
% return Select-Cycles(F)
% Polygons from cycles [4]
% Initialize empty set P
% for each cycle C in F
% p=new polygon
% for each vertex v ∈ V
% add vertex v to p
% add polygon p to set P
% return P
Sim
on 9 Oct 2019
Can Chen
on 5 Jun 2020
Hi Sim, I work at MathWorks on graphs. If you have a few minutes, I would very much appreciate hearing more about your workflow using cycles. Would you please contact me directly?
Sim
on 27 Dec 2020
Accepted Answer
Because this sounds like a generally useful thing, I cooked up the attached polyregions class to do the partitioning that you described. It uses graph theoretic functions only.
Here is its application to the data example that you provided. Each partitioned polygon is contained in the polyshape array, pgon.
s = [1 1 1 2 3 3 4 4 5 6 7 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26];
G = graph(s,t);
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0];
obj=polyregions(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off
43 Comments
Matt J
on 8 Oct 2019
I've updated polyregions.polyshape to allow you to return a list of node indices for each of the polygons:
>> [pgon,nodeIndices]=polyshape(obj);
>> nodeIndices{:}
ans =
1 8 7 6 5 4 3 2
ans =
9 10 11 6 7 8
ans =
12 14 18 1 8 9 10
ans =
13 15 16 17 18 14 12
ans =
19 20 1 18
darova
on 9 Oct 2019
Does it mean that my solution is worse?
% Matt J time
Elapsed time is 0.907870 seconds.
% Darova time
Elapsed time is 21.533292 seconds.
darova
on 9 Oct 2019
Thanks you bro =)
Feel better now
Sim
on 10 Oct 2019
Hi Matt, a question about the removal of "branches" in your code:
function obj = prune(obj)
%Remove dangling branches
....In case I have a small loop inside one of the dangling "branches" that we need to remove, like the loop 25-26-27, would you have some suggestion to workaround it ?
The code corresponding to that Figure, with the looped-branch is the following:
s = [1 1 1 2 3 3 4 4 5 6 7 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25 26 27 27];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26 27 25 28];
G = graph(s,t);
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5 5 5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0 0 0.5];
xy = [x' y']
h = plot(G,'XData',xy(:,1),'YData',xy(:,2),'linewidth',2,'MarkerSize',7);
axis([-3 7 -1.5 2.5]);
nl = h.NodeLabel;
h.NodeLabel = '';
xd = get(h, 'XData');
yd = get(h, 'YData');
text(xd, yd, nl, 'FontSize',17, 'FontWeight','bold', ...
'HorizontalAlignment','left', 'VerticalAlignment','top')
set(gca,'Fontsize',15,'FontWeight','Bold','LineWidth',2, 'box','on')
If I apply your code...
% Matt J Polygons/Cycles
obj=polyregions(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off
[pgon,nodeIndices]=polyshape(obj);
nodeIndices{:}
Inside the function
function [P,E]=polyshortest(obj,arg1,arg2)
I get this error:
nodeA =
17
nodeB =
21
P =
[]
Index exceeds array bounds.
Error in polyregions/polyshortest (line 90)
E=G.findedge(P,[P(2:end),P(1)])
Error in polyregions/polyshape (line 118)
[P,E]=polyshortest(obj,E0); % F4
Error in GraphCycles (line 22)
pgon=polyshape(obj);
Thanks a lot!
Matt J
on 10 Oct 2019
I think this revised prune() method should do it.
function obj = prune(obj)
%Remove dangling branches
Gp=obj.G;
xp=obj.x;
yp=obj.y;
tips=find(degree(Gp)<=1);
while ~isempty(tips)
Gp=rmnode(Gp,tips);
xp(tips)=[];
yp(tips)=[];
tips=find(degree(Gp)<=1);
end
bins=biconncomp(Gp);
[~,cmaj]=max(histcounts(bins,1:max(bins)+1));
ekeep=find(bins==cmaj);
[sk,tk]=findedge(Gp,ekeep);
uk=unique([sk,tk]);
xp=xp(uk); yp=yp(uk);
obj=polyregions(graph(sk,tk),xp,yp);
end
I think the attached version (rename to jigsaw.m) may fix it, but it would help to have your G,x,y attached in a .mat file for testing purposes. In the meantime, the test below verifies that scrambling the order of the node labeling doesn't interfere with things.
L=randperm(28);
Li(L)=1:28;
s = [1 1 1 2 3 3 4 4 5 6 7 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25 26 27 27];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26 27 25 28];
G = graph(L(s),L(t));
x = ([0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5 5 5]);
y = ([0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0 0 0.5]);
x=x(Li); y=y(Li);
obj=jigsaw(G,x,y);
plot(obj)
[pgon,cycles]=polyshape(obj);
hold on;
plot(pgon);
hold off
>> cycles{:}
ans =
3 16 12 19 20 17 11
ans =
4 26 5 14 28 15 10 18
ans =
5 14 28 22 8 23
ans =
8 22 28 15 20 19 12
ans =
13 25 20 15
Sim
on 11 Oct 2019
In general, is MATLAB "preserving" confidential data/material if we exchange data using this official channel (with private message/email) ?!
No, everything is totally public. If you can make some modification of the data such that it is releasable to the public, but still produces the same error, that would be optimal. In the meantime, you can try this replacement
pgon=[pgon(:).',hgon(:).'];
did you use any method from the Literature?
I didn't refer to any literature, but that doesn't mean I didn't re-invent the wheel inadvertently :-).
In general, I would just cite the URL if I want to give credit to a File Exchange submission, but because you are still making tweaks to the code that I don't understand, I am not sure if jigsaw.m is perfected enough to be posted on the FEX.
In addition, I replaced this nodeIndices{end+1}=lp(P.'); with nodeIndices{end+1}=P.';
If you do that, then the example below doesn't work properly. Notice how cycle{5} incorrectly shows 6-5-8-7.
s = [1 1 1 2 2 2 3 3 3 4 4 4 6 7 8 9 4];
t = [2 7 4 8 3 1 9 2 4 1 6 3 7 8 9 6 5];
G = simplify(graph(s,t));
x = [0 3 3 0 0.5 1 1 2 2];
y = [3 3 0 0 1 1 2 2 1];
obj=jigsaw(G,x,y);
[pgon,cycles]=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off
>> cycles{:}
ans =
1 7 8 2
ans =
1 7 6 4
ans =
2 8 9 3
ans =
3 9 6 4
ans =
6 5 8 7
But I am just wondering why the coordinates of some "vertices" are not perfectly overlapping during the operation of "ismember", since all of those coordinates come from the same initial dataset...
I think I understand why. The union() operation
U=union(pgon);
causes a recomputation of the vertices and that introduces floating point errors. In theory, we know the union of these polygons should have vertices from the original data set, but that won't be true for a general union of polygons and so polyshape.union() doesn't recognize that.
I think ismembertol is a good solution for now, but you should restore this line to the way it was,
nodeIndices{end+1}=lp(P.');
I attach here the jigsaw.m file with the latest small adjustment...
Unfortunately, this version still has a lot of problems, see attached example. I'll continue to think about solutions as well.
load DTexample
plot(obj);
tic;
[pgon,cycles]=polyshape(obj);
toc
hold on
plot(pgon);
hold off
This version seems to be doing well. I've tested it on a 1000-node graph. Unfortunately, it now requires the Statistics Toolbox (to make available the commands knnsearch and pdist2).
N=1000;
x = randn(N,1)*100;
y = randn(N,1)*100;
DT=delaunayTriangulation(x,y);
EdgeTable=table(DT.edges,'VariableNames',{'EndNodes'});
G=graph(EdgeTable);
obj=jigsaw(G,x,y);
obj.tol
tic;
[pgon,cycles]=polyshape(obj);
toc %Elapsed time is 16.075788 seconds.
hold on
plot(pgon);
hold off
my network is squeezed to the right side (right half) of the Figure window (with the previous version of jigsaw.m it did not happen.. previously, it was well centered..)
Maybe there is an isolated node somewhere that you can't see?
I get these warnings (repeated many times in the command window.. I dont know why..)...
I don't get those. If you use,
>>dbstop if warn
you can pause execution where the warning occurs. Then, save the vertex data as fed to polyshape for further study.
Sim
on 14 Oct 2019
Many thanks Matt :) I used what you kindly suggested (something similar)
dbstop in jigsaw if warning
in my list of commands... i.e. in this way...
dbstop in jigsaw if warning
obj=jigsaw(G,x,y);
obj.tol
[pgon,cycles]=polyshape(obj);
plot(obj)
hold on
plot(pgon);
hold off
....and I found out something, that I explaine here below...
Let's start from this warning
Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected
results. Input data has been modified to create a well-defined polyshape.
which disappeares if, inside the
function [pgon,labelIndices]=polyshape(obj)
I replace this line
pgon(end+1)=polyshape(Vall(P,:),'Simplify',true);
by this one
pgon(end+1)=polyshape(Vall(P,:),'Simplify',false);
In your opinion, is this replacement a problem/issue for the calculation of polygons ?
Note: in general, I can see that the warnings disappear if I use "Simplify, false"
polyshape(..,'Simplify',false);
Let's go now to the "network squeezed towards the right side of the Figure window" and to the other warnings, which occur all together at once, as follows
Warning: Output contains an empty polyshape due to duplicate vertices, intersections, or other inconsistencies.
Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected
results. Input data has been modified to create a well-defined polyshape.
Warning: Boundaries with less than 3 points were removed.
....All these warnings disappear and the picture is restored to the center as previously, if I replace the last loop of the code
for i=1:numel(hgon)
P= knnsearch(Vall,hgon(i).Vertices);
pgon(end+1)=polyshape(Vall(P,:),'Simplify',true);
labelIndices{end+1}=lp(P.'); %#ok<*AGROW>
end
with that one we used in a previous version of jigsaw.m
pgon=[pgon(:).',hgon(:).'];
for i=1:numel(hgon)
[~,P]=ismembertol(hgon(i).Vertices,[x(:),y(:)],0.02,'ByRows',true);
labelIndices{end+1}=lp(P.');
end
By using the jigsaw.m code with the above mentioned modifications (see attached file too) I checked the "visual results" for both your network and mine.... At a first glance, it looks like that both networks give correct visual results, i.e. all polygons/holes are detected... However, I am aware that "ismembertol" is not the best, while the usage of "knnsearch" is a more robust solution.. but employing "knnsearch", I got many warnings and a squeezed picture which I would not know how to fix for my network/case...
Just for information, please see below the code I used - that you kindly provided - and the corresponding Figure I got with the attached (slightly) modified jigsaw.m
N=1000;
x = randn(N,1)*100;
y = randn(N,1)*100;
DT=delaunayTriangulation(x,y);
EdgeTable=table(DT.edges,'VariableNames',{'EndNodes'});
G=graph(EdgeTable);
% dbstop in jigsaw if warning
obj=jigsaw(G,x,y);
obj.tol
tic;
[pgon,cycles]=polyshape(obj);
plot(obj)
toc %Elapsed time is 16.075788 seconds.
hold on
plot(pgon);
hold off
Let's start from this warning
Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpectedresults. Input data has been modified to create a well-defined polyshape.
which disappeares if, inside the ... I replace this line
It disappears, but that doesn't tell us why it was there, and that means there could be lingering hidden problems... Why would you have duplicate vertices, intersections, or other inconsistencies? I encourage you to post an example set of vertices that was triggering the warning so that it can be studied.
Sim
on 15 Oct 2019
The way to extract a relevant polygon is using "dbstop in jigsaw if warning". It should have stopped execution at this line,
pgon(end+1)=polyshape(Vall(P,:),'Simplify',true);
Now you just need to save and post the particular example input Vall(P,:) that is triggering the warning:
K>> vertices=Vall(P,:); nodenums= lp(P);
K>> save Debug_Data.mat vertices nodenums
Using our jigsaw class, we can plot the sub-section of the graph that is triggering the warning,
load Debug_Data
N=numel(nodenums);
s=1:N;
t=circshift(s,-1);
obj=jigsaw(graph(s,t), vertices(:,1), vertices(:,2),nodenums);
obj.plot
Somewhere in your graph, you have the following, which does indeed show an intersection of two node connections.
The algorithm I use to divide the graph into sub-polygons assumes no intersections. In general, I wouldn't know how to predict behavior if you allow them.
Assuming cross-overs/intersections are supposed to be forbidden, then it is best to use,
polyshape(______,'Simplify',true);
in the code, because that will alert you that a cross-over is present when it shouldn't be.
Sim
on 22 Oct 2019
Matt J
on 9 Dec 2019
I finally got around to posting it on the FEX
though I renamed the class yet again to spatialgraph2D. At some point, I'll get around to expanding its capabilities.
Sim
on 9 Dec 2019
Tarik Chowdhury
on 24 Dec 2019
Hi, I have also used the same code for my data.
s = [1 1 2 2 2 3 4 4 4 5 6 6 6 7 7 9 9 10 12 13];
t = [2 5 3 4 5 4 5 7 9 6 11 12 13 8 9 10 14 11 13 14];
G = graph(s,t);
x = [15 13 14.5 11.5 11.5 8 9.5 20 7.0 5.5 5.9 5.5 4.5 2.7];
y = [03 04 5.50 06.5 02.5 2 9.5 20 6.5 5.0 3.5 1.0 3.0 6.0];
obj=polyregions(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%
I am getting the loops as following:
1 5 2
2 4 3
4 2 5
4 8 7
5 4 8 9 10 6
6 12 11
8 9 10 6 12 13
at the 4th loop why node * is coming also the node 8 is creating problem with some other loops too.
Matt J
on 26 Dec 2019
@Tarik,
It looks like you are using an early, non-perfected version of the code. Here is what I get using spatialgraph2D from the link mentioned above,
s = [1 1 2 2 2 3 4 4 4 5 6 6 6 7 7 9 9 10 12 13];
t = [2 5 3 4 5 4 5 7 9 6 11 12 13 8 9 10 14 11 13 14];
G = graph(s,t);
x = [15 13 14.5 11.5 11.5 8 9.5 20 7.0 5.5 5.9 5.5 4.5 2.7];
y = [03 04 5.50 06.5 02.5 2 9.5 20 6.5 5.0 3.5 1.0 3.0 6.0];
obj=spatialgraph2D(G,x,y);
[pgon,loops]=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off
>> loops{:}
ans =
1 5 2
ans =
2 4 3
ans =
4 2 5
ans =
4 9 7
ans =
5 4 9 10 11 6
ans =
6 13 12
ans =
9 10 11 6 13 14
Matt J
on 13 Apr 2020
Am I doing something wrong or using the mismatch version of Matt's code?
Yes, you should be using the version posted here,
@Cheng,
The reason for the excluded polygons is that they dangle off the main mosaic by a trivial one-vertex long branch. As mentioned in the Overview of the FEX submission:
"Note that currently, graphs containing multiple connected or disconnected mosaics are outside the scope of what the class code will analyze."
This is the effect you're seeing and is the behavior that the OP wanted (though possibly Sim would not have regarded a 1-vertex long branch as a branch). In any case, you can try replacing the class' prune method with the following experimental alternative. It gives the behavior you're looking for on your two examples, but it needs to go through more extensive beta-testing before I officially release it.
function obj = prune(obj)
%Remove dangling branches
Gp=obj.G;
xp=obj.x;
yp=obj.y;
lp=obj.labels;
tips=find(degree(Gp)<=1);
while ~isempty(tips)
Gp=rmnode(Gp,tips);
xp(tips)=[];
yp(tips)=[];
lp(tips)=[];
tips=find(degree(Gp)<=1);
end
obj=spatialgraph2D(Gp,xp,yp,lp);
end
Matt J
on 16 Apr 2020
Well, the second error is easy. You have 420 nodes,
>> G.numnodes
ans =
420
but your x,y vectors only provide 408 values. The other error, I will have to take some time to examine.
OK, here's a new beta-version to try. It's worked on everything I've tried so far. I've also made it so that it will now throw more meaningful error messages in cases of improper input like your 420 node example above.
Matt J
on 17 Apr 2020
Never mind the last version. Try the one attached here. It should be faster.
Matt J
on 5 Jun 2020
By the way, there is a modest suggestion for someone who want to use "spatialgraph2D.m" for other purposes: any node in your graph needs to be linked by at least one edge!!
No, that is not correct. Below is an example to show that spatialgraph2D will work for graphs with isolated nodes, like in the attached example.mat.
load example
hg=plot(G);
obj=spatialgraph2D(G,hg.XData,hg.YData);
mosaic(obj)
@Matt J
Hi Matt,
Thanks for your awesome work. Could you please do a 3D version ? Here're a testing sample.
vertices = [6805.95000000000,-942.260000000000,497.960000000000;6805.95000000000,942.260000000000,497.960000000000;6862.09000000000,-946.640000000000,548.720000000000;6862.09000000000,946.640000000000,548.720000000000;6875.25000000000,-732.770000000000,503.040000000000;6875.25000000000,732.770000000000,503.040000000000;6903.34000000000,978.720000000000,617.260000000000;6917.19000000000,-1008.46000000000,847.040000000000;6917.19000000000,1008.46000000000,847.040000000000;6930.70000000000,-1007.19000000000,687.570000000000;6930.70000000000,1007.19000000000,687.570000000000;6938.81000000000,-726.260000000000,555.870000000000;6938.81000000000,726.260000000000,555.870000000000;6943.92000000000,-600.580000000000,534.310000000000;6943.92000000000,600.580000000000,534.310000000000;6953.83000000000,-827.580000000000,619.900000000000;6953.83000000000,827.580000000000,619.900000000000;6962.81000000000,-581.470000000000,570.950000000000;6962.81000000000,581.470000000000,570.950000000000;6968.29000000000,-868.050000000000,681.270000000000;6968.54000000000,-618.520000000000,597.530000000000;6968.54000000000,618.520000000000,597.530000000000;6968.55000000000,-508.330000000000,584.870000000000;6968.55000000000,508.330000000000,584.870000000000;6979.30000000000,-535.330000000000,643.970000000000;6979.30000000000,535.330000000000,643.970000000000;6980.01000000000,-650.780000000000,625.820000000000;6980.01000000000,650.780000000000,625.820000000000;6980.04000000000,-706.420000000000,620.380000000000;6980.04000000000,706.420000000000,620.380000000000;6980.28000000000,-761.010000000000,636.010000000000;6980.28000000000,761.010000000000,636.010000000000;6990.21000000000,-631.490000000000,664.060000000000;6990.21000000000,631.490000000000,664.060000000000;6990.67000000000,-776.270000000000,678.530000000000;6990.67000000000,776.270000000000,678.530000000000;7019.02000000000,-544.690000000000,770.930000000000;7019.02000000000,544.690000000000,770.930000000000];
edge_list = [9,11;11,26;26,38;9,38;37,38;25,26;25,37;8,10;8,37;10,25;3,12;3,10;10,20;12,20;1,5;1,3;5,12;31,35;16,31;16,20;20,35;29,31;12,29;12,16;27,31;33,35;27,33;27,29;21,27;12,21;5,6;12,13;6,13;18,21;14,18;12,14;25,33;21,25;18,23;23,25;18,19;15,19;14,15;24,26;23,24;19,22;19,24;22,26;28,34;22,28;26,34;13,22;13,15;28,30;13,30;28,32;34,36;32,36;30,32;17,32;13,17;7,32;11,36;7,11;2,4;2,6;4,13;7,17;4,7];
I tried on spatialgraph2D, but it's missing two pgons, 1-12-35-33-25 and 12-16-20.
Hi JUN WANG,
In 3D, what would be the criterion be for deciding which vertices form a tile? In 2D, a polygon is a tile if it doesn't overlap with any other polygons in the graph. How would that generalize to 3D?
Also, in your example, there are only 23 vertices, so it is not clear how you could have a polygon with vertex indices 1-12-35-33-25.
More Answers (2)
darova
on 7 Oct 2019
Just use for loop and cells since you already know indices of each polygon
s = [1 1 1 2 3 3 4 4 5 6 6 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26];
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0];
ind = {{1 2 3 4 5 6 7 8 1}
{6 7 8 9 10 11 6}
{1 8 9 10 12 14 18 1}
{1 18 19 20 1}
{12 13 15 16 17 18 14 12}};
cla
% plot([x(s);x(t)],[y(s);y(t)],'.b')
hold on
for i = 1:length(ind)
ix = cell2mat(ind{i});
plot(x(ix),y(ix),'color',rand(1,3))
end
hold off
axis equal
5 Comments
Sim
on 7 Oct 2019
darova
on 7 Oct 2019
Here is an idea:
Create an empty image something like
% make origin (0,0)
x = x - min(x);
y = y - min(y);
% size of image
x = 10*x;
y = 10*y;
width = max(x);
height= max(y);
I = zeros(height,width); % pixels
Create more points on each edge (more pixels)
dilate an image using imdilate() to get solid lines. Use bwlabel to separate each region:
Loop through all regions. imdilate some region a bit to identificate which points belong to area
Branches remain. Maybe loop through all points and see some point has only one neighbour?
darova
on 8 Oct 2019
Here is an example. Any ideas how to remove branches?
See the attached script
... I have a quite silly question..... I was trying to use your function in a loop for, in order to get automatically the number of polygons / cycles.... However, sometimes, it happens that there are not polygons / cycles (i.e. there are only tree-like graphs) to detect, as in this example, which returns an error...
% (1) use the function "spatialgraph2D"
>> obj = spatialgraph2D(SG, SG.Nodes.X, SG.Nodes.Y)
obj =
spatialgraph2D with properties:
G: [1×1 graph]
x: [12×1 double]
y: [12×1 double]
labels: [1 2 3 4 5 6 7 8 9 10 11 12]
pruneType: 'basic'
% (2) plot the graph
>> plot(obj)
% (3) calculate "polyshape" of "obj"
>> pgon = polyshape(obj)
Index in position 1 exceeds array bounds.
Error in spatialgraph2D (line 70)
obj.tol=min(D(2,:))/1000;
Error in spatialgraph2D/pruneBranches (line 253)
obj=spatialgraph2D(Gp,xp,yp,lp);
Error in spatialgraph2D/polyshape (line 98)
obj=pruneBranches(obj);
Is there a way to workaround this error in spatialgraph2D ?
Maybe just giving an empty "pgon" or something, but not an error (otherwise the loop for breaks..) ?
All the best,
Sim
P.S.: Just for a sake of completness.... here following the edge list of my example:
>> obj.G.Edges(:,1)
ans =
11×1 table
EndNodes
________
1 6
2 3
2 9
3 8
4 7
5 6
5 7
7 11
9 10
10 11
11 12
7 Comments
Great!
Your updated function seems to be working fine!
So, if I have understood correctly, you have set "pgon "as zero.. Thank you very much!!
...However, I got another error :-(
Error using spatialgraph2D
Input vectors x and y must have lengths equal to the number of graph nodes.
Error in spatialgraph2D/pruneBranches (line 258)
obj=spatialgraph2D(Gp,xp,yp,lp);
Error in spatialgraph2D/polyshape (line 98)
obj=pruneBranches(obj);
Error in Zofingen_temporal_path_activation_connected_components (line 289)
pgon=polyshape(obj);
I got it for this case....
>> plot(obj)
And I have also got an error just typing "plot(obj)":
Error using matlab.graphics.chart.primitive.GraphPlot
Expected XData to be an array with number of elements equal to 160.
Error in matlab.graphics.chart.primitive.GraphPlot/set.XData
Error in matlab.graphics.chart.primitive.GraphPlot
Error in graph/plot (line 112)
hObj = matlab.graphics.chart.primitive.GraphPlot('BasicGraph', ...
Error in spatialgraph2D/plot (line 182)
h = plot(obj.G,'XData',obj.x,'YData',obj.y,'linewidth',2,'MarkerSize',7);
Do you know what I can do to fix it ? :-)
Could be possible that spatialgraph2D does not work when we have isolated nodes ? (just my hypothesis..)
Just for a sake of completeness, here you are (1) the details of "obj" and (2) the corresponding "list of edges" :-)
obj =
spatialgraph2D with properties:
G: [1×1 graph]
x: [161×1 double]
y: [161×1 double]
labels: [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 … ]
pruneType: 'basic'
>> obj.G
ans =
graph with properties:
Edges: [164×2 table]
Nodes: [160×0 table]
>> obj.G.Edges(:,1)
ans =
164×1 table
EndNodes
__________
1 36
1 90
2 40
2 80
2 160
:
137 138
142 151
143 148
151 156
153 159
Display all 164 rows.
ans =
164×1 table
EndNodes
__________
1 36
1 90
2 40
2 80
2 160
3 16
3 41
4 27
4 115
5 31
5 94
6 146
6 157
7 26
7 102
7 145
8 45
8 109
8 112
9 42
9 128
10 125
10 150
11 59
11 60
11 122
12 23
14 64
14 100
15 20
16 28
16 84
17 44
17 121
18 139
19 75
20 48
20 53
21 61
21 135
21 158
22 124
22 134
23 73
23 104
24 143
24 159
25 98
26 69
27 38
27 123
28 125
29 60
29 147
30 78
30 93
31 77
31 139
32 85
32 103
32 106
33 34
33 115
33 143
34 35
34 74
36 83
37 47
37 51
38 100
38 148
39 48
39 53
41 89
42 49
42 113
43 131
43 133
44 79
44 156
45 67
45 107
48 124
49 55
51 97
52 112
52 125
53 61
53 81
55 68
55 70
56 79
56 116
57 120
57 127
57 157
58 146
59 119
60 109
60 133
61 149
64 130
65 96
65 119
65 138
66 71
66 82
66 144
69 83
70 90
70 114
71 136
72 95
72 136
74 115
74 155
75 141
76 93
77 139
77 157
78 159
79 91
80 86
81 141
82 134
84 110
86 121
88 98
89 142
90 114
91 96
91 156
94 117
95 104
95 127
96 138
97 151
98 119
100 105
101 155
102 137
102 145
103 105
103 116
104 136
105 130
108 118
109 122
111 129
111 144
114 147
117 146
118 149
123 155
123 160
129 131
129 134
131 147
133 140
137 138
142 151
143 148
151 156
153 159
>>
Matt J
on 16 Jun 2022
P.S. it looks like matlab does not allow me to edit my post in the Discussion section of the File Exchange, therefore, if something coems to my mind, I write here... sorry..
Yes, I'm not sure why that hasn't been fixed. What I do as a workaround is to copy the post to the clipboard, then delete and recompose it.
Sim
on 16 Jun 2022
Categories
Find more on Convert Image Type in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
