**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Find cycles in an undirected graph

44 views (last 30 days)

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

### Accepted Answer

Matt J
on 8 Oct 2019

Edited: Matt J
on 8 Oct 2019

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

Sim
on 9 Oct 2019

Hi Darova, absolutely not worse! ..Your idea is promising and cool!

In my case study, I am dealing with one network composed of around 3000 nodes and 6000 edges, which is changing over time, and I would need to find the polygons for several time steps. Therefore I would prefer, when possible, the fastest available solution to save time.... and with this "speed criterion" I decided to accept the Matt's answer... This does not mean that your code is worse or better... Probably, in my question, I should have mentioned something about the real size of my network and my specific application.... I think both solutions are very nice and clever, and if I could I would assign the "accept answer" equally to both of you! Thanks again for your nice code!

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

Sim
on 11 Oct 2019

Hi Matt, your solution perfectly works on this network..! Thanks a lot! :)

...However, when I try to apply this revised prune() method to my real case (a larger network), I get this error:

Index exceeds array bounds.

Error in polyregions (line 15)

w=sqrt((x(s)-x(t)).^2 + (y(s)-y(t)).^2);

Error in polyregions/prune (line 44)

obj=polyregions(graph(sk,tk),xp,yp);

Error in polyregions/polyshape (line 85)

obj=prune(obj);

Error in floods (line 98)

pgon=polyshape(obj);

At a first glance, it looks like that the components of the "w" array (i.e. x(s), x(t), y(s), y(t)) cannot be computed due to

Index exceeds array bounds.

Trying to understand the reasons of this error, I have noticed that the x and y arrays, which contain the nodes coordinates, reduce their sizes when the function prune() is called. This occurs via the xp and yp arrays, that are initialised as x and y arrays, at the beginnning of the prune() function. Indeed, when we remove the "dangling branches", the xp and yp arrays are reduced:

xp(tips)=[];

yp(tips)=[];

The xp and yp arrays are then "exported" to the new "obj", at the end of the prune() function

obj=polyregions(graph(sk,tk),xp,yp);

Therefore, in the following/remaining part of the code, the updated graph (without dangling branches) contains the reduced x and y arrays.

Now, if we change the last line of the prune() function as

obj=polyregions(graph(sk,tk),obj.x,obj.y);

where obj.x and obj.y have the original x and y arrays, the code runs further, and almost at the very end, where pgon is calculated

U=union(pgon) %#ok<*LTARG>

if ~U.NumHoles, return; end

hgon=holes(U)

pgon=[pgon,hgon]

for i=1:numel(hgon)

[~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows')

nodeIndices{end+1}=P.' %#ok<*AGROW>

end

....and exactly there I get this error:

U =

polyshape with properties:

Vertices: [349×2 double]

NumRegions: 1

NumHoles: 15

hgon =

15×1 polyshape array with properties:

Vertices

NumRegions

NumHoles

Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in polyregions/polyshape (line 117)

pgon=[pgon,hgon]

Error in floods (line 98)

pgon=polyshape(obj);

Am I doing something wrong?

Thanks a lot!

Matt J
on 11 Oct 2019

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

Hi Matt, I really really appreciate your huge help!!

...I sent you a private message/email since my data are confidential.. If you have not receive it, I can re-send it... :)

In general, is MATLAB "preserving" confidential data/material if we exchange data using this official channel (with private message/email) ?!

....I think we are getting very close to the perfect code for my case / real network... I mean, your code "somehow" works for my real network only if I comment "pgon" at the end of your code, in order to avoid the previous error

U=union(pgon); %#ok<*LTARG>

if ~U.NumHoles, return; end

hgon=holes(U);

% pgon=[pgon,hgon]

for i=1:numel(hgon)

[~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows');

nodeIndices{end+1}=P.' %#ok<*AGROW>

end

...this is the previous error....

Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in polyregions/polyshape (line 117)

pgon=[pgon,hgon]

.....but doing this, some polygons are not "detected".............. Btw, very thankful for your huge support!! Really grateful!

Matt J
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(:).'];

Sim
on 11 Oct 2019

Perfect! Thank you very much!

I replaced

pgon=[pgon,hgon]

with

pgon=[pgon(:).',hgon(:).'];

In addition, I replaced this

nodeIndices{end+1}=lp(P.'); %#ok<*AGROW>

with

nodeIndices{end+1}=P.'; %#ok<*AGROW>

And it is now working for my real network!

Please see following the exact code I used, and that you kindly provided...

s = [ ... ];

t = [ ... ];

G = graph(s,t);

x = [ ... ];

y = [ ... ];

% Matt J Polygons/Cycles

obj=jigsaw(G,x,y);

plot(obj)

[pgon,cycles]=polyshape(obj)

hold on;

plot(pgon);

hold off

....together with the dimensions of my graph G, the final result and the elapsed time...

G =

graph with properties:

Edges: [6038×1 table]

Nodes: [3111×0 table]

pgon =

1×2495 polyshape array with properties:

Vertices

NumRegions

NumHoles

cycles =

1×2495 cell array

Columns 1 through 7

{1×3 double} {1×23 double} etc.......

Elapsed time is 18.913817 seconds.

If possible, I would suggest to create a "MATLAB File Exchange"....maybe with a Reference about the method you used (did you use any method from the Literature?) I partially understood what you did, but not completely... I would be curious to know if this is somehow related/connected to the "Minimum Cycle Basis" of Graph Theory... How to cite your code in research papers ? Thanks a lot for your precious help!

Matt J
on 12 Oct 2019

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

Sim
on 12 Oct 2019

Thanks... By using this line of code

nodeIndices{end+1}=lp(P.');

I was getting the following error

Array indices must be positive integers or logical values.

Error in jigsaw/polyshape (line 137)

nodeIndices{end+1}=lp(P.');

Error in program (line 94)

[pgon,cycles]=polyshape(obj)

....Then, I checked the previous line

[~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows');

and the three quantites inside it, i.e.

P

hgon(i).Vertices

[x(:),y(:)]

As example, I got something like this:

P =

188

1873

191

0

0

1902

245

281

0

0

0

177

178

hgon(1).Vertices =

719206.761 262076.829

720011.091 260159.216

720306.696 259519.618

720812.14899999 258050.23099999

722903.21299998 256225.431

726046.209 256652.186

727673.923 255337.524

732122.269 253113.636

736764.74299995 256159.99799996

728781.39999995 267902.73300001

725492.548 269719.18599999

722128.621 271490.111

723849.229 268561.544

[x(:),y(:)] =

704726.255 281756

711700.953 265115.835

712542.533 264359.824

713410.689 263107.846

713769.936 262368.76

714125.573 261718.406

715803.539 261306.605

717089.203 260113.503

717635.691 259625.701

694183.174 253518.56

etc... etc...

Then, I thought to replace this line

[~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows');

with (ismember with a tolerance)

[~,P]=ismembertol(hgon(i).Vertices,[x(:),y(:)],0.02,'ByRows',true)

and I did not get any error anymore... 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...

In summary, I modified the very last loop of the code from this

for i=1:numel(hgon)

[~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows');

nodeIndices{end+1}=lp(P.'); %#ok<*AGROW>

end

to this:

for i=1:numel(hgon)

[~,P]=ismembertol(hgon(i).Vertices,[x(:),y(:)],0.02,'ByRows',true)

nodeIndices{end+1}=lp(P.'); %#ok<*AGROW>

end

and the last error disappeared....However, I am not sure if the usage of "ismembertol" is fully correct, instead of the previous "ismember"

Matt J
on 14 Oct 2019

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.');

Sim
on 14 Oct 2019

Thanks a lot.. I am sorry, there was a typo in my previous comment (I correct it immediately in my previous reply).... I used what you rightly proposed:

nodeIndices{end+1}=lp(P.');

and the very last loop can be summarised as

for i=1:numel(hgon)

[~,P]=ismembertol(hgon(i).Vertices,[x(:),y(:)],0.02,'ByRows',true)

nodeIndices{end+1}=lp(P.');

end

...I will try to provide a sample of my network...and to think how to "adjust" the issue about the union() operation... In order to avoid confusion, I attach here the jigsaw.m file with the latest small adjustment...

Matt J
on 14 Oct 2019

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

Matt J
on 14 Oct 2019

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

Sim
on 14 Oct 2019

Thats cool! It looks perfect!

I am trying the newer version of jigsaw.m on my network, by using exactly your commands

% my real network

s = [ ... ];

t = [ ... ];

G = graph(s,t);

x = [ ... ];

y = [ ... ];

% Matt J Polygons/Cycles

obj=jigsaw(G,x,y);

obj.tol

tic;

[pgon,cycles]=polyshape(obj);

toc

hold on

plot(pgon);

hold off

% my elapsed time is

% Elapsed time is 16.687441 seconds.

..It looks like that all the polygons are well detected and the code is faster, but I have just noticed a couple of things..

- 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..)
- I get these warnings (repeated many times in the command window.. I dont know why..)...

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: Output contains an empty polyshape due to duplicate vertices, intersections, or other inconsistencies.

Warning: Boundaries with less than 3 points were removed.

:-)

Matt J
on 14 Oct 2019

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

Matt J
on 15 Oct 2019

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

Thanks for the reply.. Yes, I might know why I have some intersections or general inconsistencies in my case... I am trying to select and "cut" a small slice of my network to bring it here as sample/example.. but, so far, my attempts to cut and isolate a portion (some polygons) from the rest of my network failed.. In this moment, I am trying to remove all the nodes in the network with

rmnode(G,nodeIDs)

except those ones composing the polygons I would like to bring here... Actually, I also tried the other way round, just picking up the polygons I want, but it looks like MATLAB is reassigning the nodes IDs, messing up the edges... I mean, in my network, a polygon can appear as this chain of nodes IDs: 1032, 1033, 560, 1044, 1045, 1046, 799, etc.. (and the nodes IDs are not always sequential), then, when I make the graph of just those nodes IDs, in that order, MATLAB is reassigning those nodes IDs from 1, to the lenght of such a polygon, as, 1, 2, 18, 6, 7, 8, 13, ...but it assigns the edges in this order: 1, 2, 3, 4, 5, 6, ... therefore the edges are in different order than in my polygon... Well, I will find some solution :)

Matt J
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

Sim
on 15 Oct 2019

By using the attached file (the same you sent me around 16 hours ago) and these commands

dbstop in jigsaw if warning

obj=jigsaw(G,x,y);

obj.tol

plot(obj)

[pgon,cycles]=polyshape(obj);

hold on

plot(pgon);

hold off

I got the following output in the command window

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.

K>> vertices=Vall(P,:); nodenums= lp(P);

Undefined function or variable 'Vall'.

K>>

In addition, the debugging indicates the warning line

--> warning(message('MATLAB:polyshape:repairedBySimplify'));

inside the file (line 421)

polyshape.m

Matt J
on 21 Oct 2019

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

Thanks a lot Matt J for this clever analysis!

You are rigth about the usage of

polyshape(______,'Simplify',true);

....now I need to think a bit about how to handle with those intersections... Btw, jigsaw class is an excellent tool, and I am really grateful for your help!

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 14 Apr 2020

@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 17 Apr 2020

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)

Jun W
on 27 Dec 2020

@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.

Matt J
on 27 Dec 2020

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

Hi Darova, what you indicated as

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}};

for me is unknown, and it is exactly what I asked for in my question.. That should be the output and not the input.. I wrote it just to indicate exactly what I would need...

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?

Sim
on 9 Oct 2019

Hi Darova, to remove branches I used this:

% Remove branches

xy = [x' y'];

while ~isempty(find(degree(G)==1))

degreeone = find(degree(G)==1);

G = rmnode(G,degreeone);

xy(degreeone,:) = [];

end

About your idea/proposal, it looks cool and promising! Also, a nice animation!

However, is there any way to extrapolate the nodes composing each "cycle"/"polygon" that you were able to isolate ?

Thanks a lot for your efforts!

Sim
on 15 Jun 2022

Edited: Sim
on 15 Jun 2022

... 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

Sim
on 16 Jun 2022

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

Sim
on 16 Jun 2022

Ok, thanks a lot!

Please see attached the following file for you :-)

obj.mat

Btw, It looks like that the difference is just of one node between obj.G and obj.x and obj.y

>> numnodes(obj.G)

ans =

160

>> size(obj.x)

ans =

161 1

>> size(obj.y)

ans =

161 1

I am confident that keeping 160 nodes for both obj.G and obj.x and obj.y will solve the problem... I have tried to play around obj.x and obj.y in order to identify what/where is this "extra" value, but I did not solve it yet..

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..

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.

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- 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)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)