Find cycles in an undirected graph

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

https://ch.mathworks.com/matlabcentral/answers/425321-find-cycle-in-array

are not working for my case...

Any further idea / suggestion ?

Thanks a lot!

[Steven Lord]

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?

[Matt J]

Because the ulterior motive is to partition the graph into non-overlapping polygons.

[Sim]

I found also this idea/method "Polygons From Set of Line segments", explained in

http://ac.aua.am/sergey_yeranosyan/Web/algorithm%20project/

and based on Ferreira A., Fonseca M.J., Jorge J.A. (2003) Polygon detection from a set of lines

http://web.ist.utl.pt/alfredo.ferreira/publications/12EPCG-PolygonDetection.pdf

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]

@Steven Lord : Kindly Matt J gave the answer... In this moment, I would just need to know the "cycles" corresponding to the polygons which compose my network...

[Can Chen]

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?

Answer 1

[Matt J]

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

[Matt J]

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]

Hi Matt J, thanks a lot for you answer! It works!

I accepted your answer since it gives exactly what I was looking for in a short time..

..Thanks a lot Matt J and Darova for your efforts and ideas/solutions!

[darova]

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.

[Sim]

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!

[darova]

Thanks you bro =)

Feel better now

[Sim]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

Undefined function or variable 'Vall'.

When execution pauses, you need to be sure you first change to the workspace of jigsaw.polyshape.

[Sim]

Ops...thanks... Please see attached the file

Debug_Data.mat        

(and sorry for my late reply)

[Matt J]

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]

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!

[Matt J]

I finally got around to posting it on the FEX

https://www.mathworks.com/matlabcentral/fileexchange/73630-spatialgraph2d

though I renamed the class yet again to spatialgraph2D. At some point, I'll get around to expanding its capabilities.

[Sim]

Cool!! I am sure it will be very useful to many people :) ..Hopefully I will cite you soon :)

[Tarik Chowdhury]

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]

@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

[Cheng Ting Tsai]

Hi, Sim

I have the same error about this when I apply Matt's code to my real network (about 400 nodes) :

Index exceeds the number of array elements (415).

Error in jigsaw (line 16)
              w=sqrt((x(s)-x(t)).^2 + (y(s)-y(t)).^2);

Error in loadingFiles_CY_senpai (line 610)
obj=jigsaw(G,x,y);

I have alreadly read the discussions between you and Matt, and I tried the solutions in your discussions. However, I still can't figure out where the problem is. Could you tell me how you solve this problem in your previous discussions with Matt J on 11 Oct 2019? Thanks a lot!

(The below is my method to set s, t, x and y because of 400 nodes in my case...)

In my real network, I use "for" loop and "adjacency matrix" to put elements into s, t, x and y :

s = [];
t = [];
x = [];
y = []; 

n = 1;
p = 1;
for row=1:AllParticles
    for col=1:AllParticles
        if col>row && AdjMat{i}(row,col)==1  % I know all elements in AdjMat(adjacency matrix)
            s(1,n) = row;    % s = [ ... ];
            t(1,n) = col;    % t = [ ... ];
            n = n + 1;
        end
    end
end

G = graph(s,t);              % G = graph(s,t);

for j=1:AllParticles
    if AdjMat{i}(j,j)~=0  
       x(1,p) = Plot_Position(j,1);    % x = [ ... ];
       y(1,p) = Plot_Position(j,2);    % y = [ ... ];
       p = p + 1;
    end 
end

% Matt J Polygons/Cycles
obj=jigsaw(G,x,y);
plot(obj)
[pgon,cycles]=polyshape(obj);
hold on
plot(pgon);
hold off

I try to figure out my error :

After running the for loop from (s & t) to (x & y), I got the number of elements in s = 1* 576 double, t = 1* 576 double, x = 1* 404 double, y = 1* 404 double. I check the network and data, and they are correct with my network figure.

[these are not code]
cell in worksapce :
s = 1x576 double
t = 1x576 double
x = 1x404 double
y = 1x404 double

However, when I keep running the "Matt J Polygons/Cycles" code, the number of elements change into s = 1* 613 double, t = 1* 613 double, x = 1* 415 double, y = 1* 415 double.

[these are not code]
cell in worksapce :
s = 1x613 double
t = 1x613 double
x = 1x415 double
y = 1x415 double

Then, I got the error below because the number 415(fake) exceeds 404(real) ?

Index exceeds the number of array elements (415).

Error in jigsaw (line 16)
              w=sqrt((x(s)-x(t)).^2 + (y(s)-y(t)).^2);

Error in loadingFiles_CY_senpai (line 610)
obj=jigsaw(G,x,y);

Am I doing something wrong or using the mismatch version of Matt's code?

Thanks a lot!

[Matt J]

Am I doing something wrong or using the mismatch version of Matt's code?

Yes, you should be using the version posted here,

https://www.mathworks.com/matlabcentral/fileexchange/73630-spatialgraph2d

[Cheng Ting Tsai]

Hi, Matt J

Really thanks for your kindly help!! Your "spatialgrapg2D.m" is really a great work!!

It does well on the examples of previous discussions you gave.

...However, when I try to apply it on my real case (about 400 nodes), I still got the same error.

By using the attached file from your link :

s = [];
t = [];
x = [];
y = []; 

n = 1;
p = 1;
for row=1:AllParticles
    for col=1:AllParticles
        if col>row && AdjMat{i}(row,col)==1 
            s(1,n) = row;     % s = [ ... ];
            t(1,n) = col;     % t = [ ... ];
            n = n + 1;
        end
    end
end

G = graph(s,t);               % G = graph(s,t);

for j=1:AllParticles
    if AdjMat{i}(j,j)~=0
       x(1,p) = Plot_Position(j,1);    % x = [ ... ];
       y(1,p) = Plot_Position(j,2);    % y = [ ... ];
       p = p + 1;
    end 
end

% spatialgraph2D.m
obj=spatialgraph2D(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

I got the same error in the command window

Index exceeds the number of array elements (418).

Error in spatialgraph2D (line 18)
              w=sqrt((x(s)-x(t)).^2 + (y(s)-y(t)).^2);

Error in loadingFiles_CY_senpai (line 708)
obj=spatialgraph2D(G,x,y);

[ sorry for my repetition ]

I try to figure out my error :

After running the for loop from (s & t) to (x & y), I got the number of elements in s = 1* 595 double, t = 1* 595 double, x = 1* 408 double, y = 1* 408 double. I check the network and data, and they are correct with my network figure.

cell in worksapce :

s = 1x595 double

t = 1x595 double

x = 1x408 double

y = 1x408 double

However, when I keep running the "spatialgraph2D" section, the number of elements change into s = 1* 632 double, t = 1* 632 double, x = 1* 418 double, y = 1* 418 double.

cell in worksapce :

s = 1x632 double

t = 1x632 double

x = 1x418 double

y = 1x418 double

Then, I got the error below because the number 418(fake) exceeds 404(real)?

Index exceeds the number of array elements (418).

I will be really appreciated for your guidance and kind help. Thank you!

[Matt J]

I suggest you attach your s,t,x,y data in a .mat file, so that they can be examined.

[Cheng Ting Tsai]

Dear Matt, thanks for your suggestion!

[Cheng Ting Tsai]

@Matt J

Hi, Matt

This is my matlab code of figure 1 (using your spatialgraph2D.m ) :

% Example. 1
s = [1 2 3 4 5 6 7 8 8 9 9 9 9 10 11 12 13 13 13 14 15 15 16 17 18 18 19 20 20 21 21 22 23 24 24 25 26 27 28 29 30 32 32 34 35 36 36 38];
t = [3 3 6 5 6 8 8 11 12 10 13 33 35 11 14 15 14 17 20 16 16 23 20 18 19 31 20 28 29 22 23 24 24 26 25 26 27 28 29 30 31 33 34 35 36 37 38 39];
G = graph(s,t);
x = [1 1.5 2 3 2.5 3 3.5 4 4 4 4.5 5 5 5 6 6.5 5.5 6 6.5 7 7 8 7.5 8.5 10 9 8.5 8 7.5 7.5 7 4 4 3 3 2.5 2 2 1.5];
y = [0 1 0 2 1 0 1 0 3 2 1 0 3 2 0 1.5 3.5 4 3.5 2.5 1 1 0 0 0 1 2 2.5 3.5 4.5 5 5 4 4.5 3.5 3 3.5 2.5 2];
obj=spatialgraph2D(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

Then, I got the 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.

Results:

I have some questions when i use the version from your link:

Why cycle 9-33-32-34-35-9 in figure 1 can't be recognized as a polygon?

[Figure. 1]

However, cycle 9-33-32-34-35-9 in figure 2 can be recognized as a polygon?

[Figure. 2]

According to above questions, I test a more basic case to examine (sorry but i need to say it) your code:

% Example. 2
s = [1 1 1 1 2 4];
t = [2 3 4 5 3 5];
G = graph(s,t);
x = [3 4 2 4 6];
y = [2 0 0 4 4];
obj=spatialgraph2D(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

...and I got the answer below : cycle 1-4-5-1 in figure 3 can't be recognized as a polygon.

[Figure. 3]

It is really closed to my real case. Only one step away!

You help me a lot. Thank you!

[Matt J]

@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

[Cheng Ting Tsai]

@Matt J

Hi, Matt

I am really appreciated for your kindly help! It works well after using your alternative class' prune method.

However, I always get an error when I use the code on my real network case:

[Case 1 ] : more complex than previous examples (64 nodes)

s = [1 2 3 4 5 6 7 8 8 9 9 9 9 10 11 12 13 13 13 14 15 15 16 17 18 18 19 20 20 21 21 22 23 24 24 25 25 26 27 27 28 29 30 31 32 32 34 35 36 36 38 40 41 41 42 43 44 44 45 46 47 47 48 49 50 50 51 52 52 53 54 56 57 59 60 61 61 62 63];
t = [3 3 6 5 6 8 8 11 12 10 13 33 35 11 14 15 14 17 20 16 16 23 20 18 19 31 20 28 29 22 23 24 24 26 25 26 62 27 28 55 29 30 31 40 33 34 35 36 37 38 39 41 42 44 43 44 46 47 46 50 48 50 49 52 51 55 52 53 54 54 58 57 59 62 62 63 64 63 64];
G = graph(s,t);
x = [1 1.5 2 3 2.5 3 3.5 4 4 4 4.5 5 5 5 6 6.5 5.5 6 6.5 7 7 8 7.5 8.5 10 9 8.5 8 7.5 7.5 7 4 4 3 3 2.5 2 2 1.5 7 8 8.5 9 8.5 7.5 8 9 9.5 11 8.5 9.5 11 12 11.5 8.5 10.5 11 12 11.5 10.5 12.5 11 12 13];
y = [0 1 0 2 1 0 1 0 3 2 1 0 3 2 0 1.5 3.5 4 3.5 2.5 1 1 0 0 0 1 2 2.5 3.5 4.5 5 5 4 4.5 3.5 3 3.5 2.5 2 6 6.5 8 7 6 5.5 5 5 5.5 5 4 4.5 4 4 3 3 3 2 2 1 1 1 0 0 0];
obj=spatialgraph2D(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

and the following error:

Index exceeds the number of array elements (0).

Error in spatialgraph2D/polyshortest (line 103)
            E=G.findedge(P,[P(2:end),P(1)]);

Error in spatialgraph2D/polyshape (line 132)
                [P,E]=polyshortest(obj,E0);

Error in loadingFiles_CY_senpai (line 730)
pgon=polyshape(obj);

[ Case 2 ] : my real network (about 400 nodes)

s = [1,1,2,2,3,3,4,4,4,5,5,5,7,7,7,8,8,8,9,9,9,10,10,12,12,12,12,12,13,13,13,14,14,14,14,16,17,17,17,18,18,18,19,19,19,19,20,20,20,21,21,21,21,22,22,22,22,23,23,24,24,24,24,24,25,25,26,26,27,27,27,27,28,28,29,29,30,30,30,30,31,31,31,32,32,33,33,33,34,34,35,35,36,36,37,37,38,38,38,38,39,39,39,39,40,40,40,41,42,42,42,43,43,44,44,44,45,45,46,46,46,46,47,48,48,48,49,49,50,50,51,51,51,52,53,54,54,54,54,55,55,55,55,56,56,57,57,58,59,59,60,60,60,61,62,62,62,63,63,63,63,64,64,65,65,65,65,66,66,67,67,67,67,68,68,69,69,69,70,70,71,71,72,72,72,72,73,74,74,74,75,75,75,75,76,76,77,77,77,78,79,79,79,80,81,81,81,82,82,82,82,83,84,84,84,84,85,85,86,86,86,86,87,87,87,88,89,89,90,90,91,92,92,92,93,93,93,94,94,95,95,95,96,97,97,98,98,98,99,99,100,100,101,102,102,102,103,103,104,104,105,105,106,107,107,107,107,109,109,110,111,112,112,112,113,114,114,114,115,115,116,116,117,117,118,119,119,119,119,120,120,121,121,121,122,122,123,123,124,125,126,127,128,128,129,130,130,130,131,132,132,132,133,133,133,134,135,135,135,135,136,136,136,137,137,138,138,138,138,140,140,143,143,145,146,146,146,147,147,148,148,149,150,151,151,152,152,153,153,153,154,154,154,155,155,155,156,156,156,157,157,157,157,157,158,158,159,159,160,161,161,162,162,163,163,163,164,164,165,165,165,166,166,167,168,168,169,169,170,170,170,171,172,173,174,175,175,176,177,177,178,178,179,179,179,180,180,182,182,182,182,183,184,185,185,186,187,187,188,188,188,189,189,189,190,191,191,192,192,193,193,194,194,195,195,196,196,197,197,198,199,199,201,201,201,202,202,203,205,205,207,207,207,208,209,209,212,213,214,214,215,216,217,217,218,219,220,220,220,220,221,222,222,222,223,223,224,224,225,225,226,226,227,227,228,228,228,229,229,230,233,234,235,235,236,238,238,239,242,242,244,244,245,245,248,248,249,250,250,250,250,252,254,254,255,256,257,259,260,260,261,261,264,269,269,270,270,271,272,273,273,276,276,277,278,282,282,283,284,285,288,293,293,295,298,298,300,301,302,308,308,309,311,312,312,314,318,327,327,328,330,331,332,338,340,348,350,352,353,353,357,370,372,381,386,386,401,402,403,404,405,406,407,408,409,410,411,412,413,414,415,416,417,418,419];
t = [16,59,35,64,389,416,279,303,315,66,102,312,227,395,396,32,90,362,225,296,411,165,181,50,52,76,128,297,211,266,304,36,116,167,358,328,190,322,363,48,289,393,94,142,175,342,33,106,221,26,28,47,103,120,258,266,303,121,183,36,40,108,289,356,132,238,186,262,144,233,286,318,51,73,218,376,171,291,355,397,45,173,277,348,373,57,226,291,124,159,131,174,289,356,169,271,47,101,118,300,104,159,369,382,80,108,393,209,92,185,317,223,365,93,258,326,323,404,52,148,241,354,247,163,170,199,180,254,297,354,111,247,343,53,332,97,127,313,397,173,215,236,323,271,378,160,270,310,235,364,218,376,394,119,122,232,361,156,202,299,400,161,174,106,221,269,379,156,302,125,333,341,342,108,203,101,118,355,112,384,97,299,144,152,172,347,224,285,317,365,83,97,230,268,174,251,114,217,331,110,239,291,397,88,172,269,309,111,290,325,343,280,105,214,240,352,162,345,194,219,356,388,147,200,329,391,148,234,216,384,266,153,358,365,205,385,390,341,368,215,301,322,148,230,400,112,320,377,120,258,166,214,270,104,130,134,117,278,345,375,274,352,270,183,261,315,346,196,281,403,118,129,277,377,253,126,329,389,261,303,195,204,192,200,247,206,264,287,388,279,303,274,275,346,172,286,300,374,274,311,152,239,251,284,173,134,333,368,251,232,257,285,185,281,317,162,267,282,311,316,175,308,342,345,375,230,268,304,357,214,299,244,377,370,220,233,286,360,392,175,212,155,404,319,373,309,347,195,263,358,306,325,357,273,290,327,282,312,316,184,191,259,328,378,381,408,352,382,355,265,307,353,370,167,272,330,181,264,292,297,307,274,346,289,260,311,273,385,187,336,348,313,379,276,297,212,332,329,254,310,344,420,338,344,419,186,399,189,268,280,304,197,290,263,283,200,301,322,216,301,348,211,243,321,336,224,343,242,329,253,359,219,264,198,234,226,257,315,349,219,231,391,248,349,366,240,302,388,245,327,239,267,281,417,262,399,227,296,243,280,253,337,414,415,251,241,232,314,339,362,257,258,350,385,314,339,262,328,298,347,291,387,283,395,320,409,410,243,261,313,335,283,288,399,301,281,317,334,278,374,362,384,266,398,315,340,360,306,313,325,355,356,310,344,331,302,361,324,281,396,321,346,354,351,374,300,351,324,358,378,390,337,384,406,367,334,400,396,394,339,310,338,344,402,412,413,367,337,369,311,395,374,316,316,333,373,320,357,390,364,348,413,394,386,366,371,385,369,368,376,398,376,411,409,392,417,402,403,404,405,406,407,408,409,410,411,412,413,414,415,416,417,418,419,420];
G = graph(s,t);
x = [23.8262490000000,2.02198900000000,17.6007040000000,22.1728820000000,11.8965620000000,8.77028500000000,6.58055500000000,12.7853790000000,0.127039000000000,3.39282500000000,19.4155370000000,4.50498400000000,23.2761410000000,0.540076000000000,2.81948100000000,7.83719800000000,13.0894750000000,18.8797340000000,21.6978350000000,19.3008590000000,2.20991000000000,10.2695750000000,19.9011580000000,10.3218370000000,18.8931840000000,7.02514300000000,14.5889410000000,5.07614600000000,5.63567200000000,13.5088600000000,15.7719520000000,3.10665800000000,3.34152200000000,23.3685740000000,16.8145410000000,13.3881850000000,1.18769500000000,21.8117840000000,7.45503600000000,7.55301100000000,21.9379060000000,4.39822900000000,4.68505100000000,17.8417940000000,3.29164600000000,22.1191020000000,2.56654800000000,18.7887790000000,4.37248700000000,5.48611100000000,14.7835530000000,3.88648200000000,22.7905010000000,14.3781330000000,22.6718300000000,22.9295220000000,6.70256400000000,0.721472000000000,11.1034610000000,14.0853860000000,1.59887100000000,13.1593430000000,12.8854460000000,10.0016590000000,0.512816000000000,16.1988150000000,6.97700000000000,15.3555850000000,12.9257420000000,19.9961390000000,8.32222000000000,16.3027430000000,3.76925200000000,15.9737090000000,2.63972400000000,13.3530500000000,0.491128000000000,13.5071600000000,18.8557900000000,16.3549460000000,15.9122610000000,10.1781340000000,3.00094900000000,18.1579760000000,0.217156000000000,6.16577000000000,6.07346000000000,20.0792770000000,6.58842100000000,21.6659620000000,8.71617900000000,2.72904500000000,5.55486400000000,15.4007190000000,8.72196200000000,23.5306500000000,17.1877480000000,15.7935140000000,11.5027000000000,18.0158900000000,12.0120310000000,16.2699530000000,14.0172160000000,20.3379350000000,1.04836000000000,11.1909620000000,2.08783500000000,18.1904910000000,7.35183200000000,0.865135000000000,16.4088210000000,20.7469850000000,5.05494600000000,18.0786340000000,17.1712470000000,1.27971700000000,22.7851390000000,18.2613720000000,11.5317210000000,15.6192440000000,16.7308930000000,10.5014310000000,15.2036840000000,13.6563040000000,4.62452000000000,6.11339400000000,10.4793880000000,4.36705100000000,10.1460550000000,9.25804300000000,10.5650120000000,11.9320510000000,8.21742200000000,11.7582940000000,17.7647450000000,15.8885390000000,7.60380900000000,8.26044900000000,11.5763410000000,8.77338800000000,9.71889100000000,18.9129600000000,5.87945900000000,19.1859170000000,3.20001100000000,6.41687800000000,14.1912770000000,6.09623000000000,18.2394390000000,20.2990830000000,12.7137420000000,21.8021880000000,8.00047600000000,14.6513600000000,14.9958150000000,0.443226000000000,9.32962600000000,4.09895900000000,1.79216400000000,0.894920000000000,18.0989240000000,3.90086800000000,10.0856400000000,22.1582080000000,3.49342200000000,15.3576900000000,12.5136630000000,5.07052700000000,2.49471400000000,6.99430500000000,15.8905120000000,23.7175480000000,23.6005880000000,22.0115040000000,21.5723960000000,1.35434800000000,17.6302270000000,20.1516940000000,20.8828640000000,8.16051700000000,20.4563290000000,3.02612700000000,4.47119200000000,18.1532190000000,1.23722000000000,20.6682690000000,17.0827830000000,0.786685000000000,3.09581800000000,6.08702800000000,11.6909580000000,21.2300270000000,5.29950900000000,2.17876100000000,19.2016430000000,22.3981580000000,14.1551320000000,1.02881900000000,4.17430100000000,20.6446030000000,0.548409000000000,11.2370250000000,18.8002530000000,21.2017430000000,19.1174820000000,7.13742100000000,11.8645570000000,16.7517880000000,2.76032400000000,5.30933500000000,16.4000210000000,5.90604000000000,4.16725700000000,8.84954600000000,12.1923150000000,23.2737980000000,7.89966900000000,21.3169330000000,13.5599660000000,12.6378570000000,7.55060900000000,10.4004970000000,18.3682550000000,16.6253170000000,1.01248000000000,10.0131560000000,9.37017200000000,6.65655700000000,23.0642720000000,4.01434900000000,9.68787000000000,12.5848630000000,14.8621990000000,4.02262400000000,16.2126720000000,17.7763060000000,7.51519800000000,20.2359330000000,17.9442460000000,22.2540500000000,20.5073770000000,16.7838310000000,4.88057800000000,1.84869900000000,1.67440000000000,22.6742820000000,14.6334910000000,12.7129340000000,11.2758170000000,22.4833240000000,22.7421430000000,9.94355400000000,19.5479790000000,20.5183240000000,7.07582700000000,2.57989000000000,0.829763000000000,20.4855390000000,11.6609180000000,17.2483810000000,14.1088810000000,15.0190290000000,23.9686850000000,4.92073800000000,21.3656520000000,17.2629020000000,17.6588680000000,5.64387700000000,6.34787900000000,17.1727090000000,23.2191010000000,16.9765030000000,10.4004500000000,13.2395820000000,7.65393000000000,5.70430700000000,8.96569900000000,10.5411090000000,0.0520960000000000,23.7979340000000,2.85491800000000,20.0827040000000,13.5799320000000,1.28271900000000,21.5899690000000,0.863213000000000,12.1744180000000,2.03592600000000,14.0001460000000,15.2064700000000,15.8406140000000,3.97489400000000,13.5319730000000,21.6877370000000,18.3987610000000,17.2647990000000,0.135615000000000,9.44502100000000,14.3678120000000,23.2253050000000,10.5715470000000,11.6871710000000,16.0180230000000,7.84432600000000,21.3619400000000,11.5176910000000,8.51347400000000,9.95105400000000,6.02767200000000,9.81953400000000,18.8785260000000,1.82446100000000,3.77215500000000,23.8144170000000,18.0347320000000,21.2344670000000,19.7837510000000,22.3385650000000,16.9498250000000,4.27363700000000,15.0927970000000,6.47773500000000,10.8434470000000,12.8163630000000,8.49995300000000,2.28178100000000,4.98506500000000,20.9715070000000,9.02079800000000,23.2726740000000,9.61444100000000,8.91996700000000,19.5660080000000,22.7151120000000,11.1404890000000,19.2033750000000,13.5836820000000,4.59413000000000,21.5113390000000,23.8463520000000,14.7265670000000,15.1886840000000,8.50657400000000,3.51056200000000,15.9031810000000,3.02819600000000,18.6613020000000,5.54265000000000,0.316063000000000,19.7036550000000,11.1437910000000,7.79524400000000,0.104565000000000,22.4756360000000,7.20263500000000,23.4627340000000,16.7886400000000,9.48998600000000,14.1379390000000,8.18261900000000,5.29880300000000,11.6000170000000,6.85149800000000,15.0723400000000,12.6889330000000,7.55882400000000,8.02458300000000,21.7717560000000,12.6436230000000,9.20094500000000,13.7782180000000,6.75831800000000,22.7143350000000,19.9611680000000,12.2863790000000,1.93241600000000,16.9918400000000,20.9935110000000,1.24426700000000,19.3997310000000,1.85927000000000,6.52322400000000,8.51172700000000,8.82919900000000,14.2390700000000,19.4174150000000,21.9803610000000,14.4166060000000,0.200000000000000,1.40000000000000,2.60000000000000,3.80000000000000,5,6.20000000000000,7.40000000000000,8.60000000000000,9.80000000000000,11,12.2000000000000,13.4000000000000,14.6000000000000,15.8000000000000,17,18.2000000000000,19.4000000000000,20.6000000000000,21.8000000000000,23];
y = [8.27149400000000,22.8299740000000,1.03487900000000,20.0546930000000,19.2103280000000,14.7750580000000,6.42103200000000,1.01359700000000,18.7852280000000,19.6803720000000,16.8936150000000,11.9554910000000,9.15535400000000,6.19994900000000,10.0948380000000,18.7491890000000,9.16107300000000,7.38501000000000,17.8476540000000,22.6412240000000,12.3600810000000,9.88090500000000,6.88591700000000,3.91872600000000,8.52067300000000,22.2936940000000,11.9940200000000,1.87666900000000,7.38949400000000,10.2569080000000,21.9511130000000,22.1620010000000,12.2618520000000,13.5183470000000,8.60393700000000,21.3791740000000,11.5179340000000,7.79904200000000,11.8556190000000,9.47037400000000,15.6470170000000,0.916959000000000,17.4079240000000,7.90957600000000,9.10067300000000,4.09875200000000,18.5498430000000,9.70204200000000,18.6542240000000,19.3013690000000,14.0847500000000,3.12925400000000,12.3485010000000,9.51594100000000,5.12104400000000,7.44411500000000,20.7767640000000,17.6588510000000,6.76581700000000,17.1098090000000,21.7733890000000,7.58029600000000,18.4341100000000,18.1369460000000,13.7431250000000,10.6283210000000,3.41029600000000,16.4713750000000,4.00411500000000,8.96659200000000,10.4209570000000,15.9857440000000,20.9736380000000,1.83937600000000,2.18828200000000,12.4550630000000,10.5284330000000,5.59251400000000,11.7612640000000,17.1066630000000,19.8589320000000,22.8493930000000,14.5493850000000,4.11371400000000,9.52281400000000,15.5542000000000,5.25591300000000,18.5741980000000,11.2464830000000,14.6825930000000,19.5953720000000,5.15596800000000,18.1284070000000,15.2400480000000,2.27055700000000,17.5536720000000,20.1710720000000,9.50016200000000,20.2952010000000,6.77860400000000,21.3916510000000,20.9916980000000,8.40915500000000,20.8021630000000,12.7667770000000,11.2672370000000,1.11633400000000,10.7665430000000,2.26990300000000,4.45679000000000,2.88218700000000,19.7495050000000,13.0457320000000,5.70845300000000,9.91958500000000,16.5441940000000,18.5985060000000,21.9599140000000,5.72859500000000,6.82081200000000,22.3661370000000,17.2040730000000,3.20289700000000,13.9622900000000,20.0879010000000,2.66171200000000,19.7985820000000,22.0569680000000,8.50171500000000,12.4189420000000,21.1123850000000,15.6487700000000,17.1688440000000,23.0228760000000,15.2819700000000,18.3658750000000,20.0213030000000,4.32993400000000,4.38824000000000,22.7975760000000,5.92083500000000,3.05175500000000,16.9077970000000,12.9425160000000,0.898377000000000,8.92553000000000,3.84580400000000,12.2765950000000,13.6623360000000,13.0885100000000,17.1159570000000,10.5419810000000,1.10297500000000,21.7801020000000,10.7109910000000,21.5359460000000,21.7657770000000,9.80340400000000,17.6345570000000,19.7282250000000,20.8679400000000,10.8716320000000,14.6527040000000,13.7027500000000,7.98967800000000,12.9576420000000,5.20071800000000,3.04819500000000,20.9740130000000,17.7157650000000,4.18315400000000,3.20230200000000,1.11494700000000,1.12401000000000,5.05104000000000,18.5685210000000,17.2492620000000,21.8496220000000,11.1509070000000,12.7255590000000,5.73256800000000,6.74688600000000,6.84393000000000,18.1455340000000,7.25905200000000,10.0468860000000,5.04492600000000,3.03704200000000,15.7570710000000,13.5343030000000,9.98380300000000,21.6740430000000,14.5361940000000,9.08261300000000,5.05246900000000,22.0856990000000,18.4197540000000,14.6207530000000,13.8134370000000,15.1776320000000,15.5368460000000,13.2223930000000,1.01477100000000,6.85284800000000,17.8710000000000,16.4300980000000,3.06437100000000,19.0622360000000,3.83369400000000,6.11431300000000,0.914799000000000,21.6127090000000,15.1804060000000,6.72756800000000,8.32849100000000,15.6845990000000,8.34705400000000,8.88637800000000,1.86362500000000,10.7701620000000,15.2177700000000,0.974548000000000,19.8720860000000,14.9900490000000,8.54527100000000,7.21067400000000,4.78900400000000,14.5054530000000,6.16154800000000,4.43375400000000,11.1648080000000,13.4441280000000,19.3645260000000,16.4345940000000,5.84313800000000,19.0472560000000,5.16680800000000,16.2196120000000,8.99087700000000,21.0262640000000,4.51775900000000,12.8871490000000,21.1295470000000,13.7122250000000,3.57139700000000,3.25048700000000,2.17155600000000,20.0652100000000,8.83866900000000,16.7091020000000,11.1636910000000,13.3718540000000,19.9093040000000,7.86912400000000,12.8382670000000,16.8184990000000,22.7136530000000,17.3175830000000,14.4195520000000,16.2739180000000,6.62138400000000,8.56556600000000,12.5276820000000,10.4664050000000,12.7862860000000,21.4513520000000,22.9492050000000,4.05134900000000,1.30549300000000,6.28068400000000,19.7557150000000,18.0517350000000,12.2486520000000,15.8767700000000,13.8868380000000,20.4239300000000,9.23163100000000,5.13908100000000,16.4981460000000,5.21262000000000,11.2883360000000,11.9686850000000,11.3268900000000,20.7029110000000,2.39564400000000,1.06976700000000,1.92465300000000,19.7709740000000,0.899167000000000,17.5235850000000,7.90861400000000,5.73165100000000,19.3781240000000,19.1092960000000,16.3571880000000,14.0119510000000,20.4233510000000,16.7082000000000,4.94581600000000,4.19384600000000,15.9577040000000,17.9667250000000,13.9796180000000,7.24474900000000,20.6118670000000,16.7956300000000,11.5967540000000,2.81918100000000,10.2054930000000,1.83115100000000,18.9082890000000,6.13346100000000,1.82835600000000,11.4253990000000,12.5816050000000,16.3855440000000,14.1608070000000,9.55124700000000,3.89530000000000,8.77404600000000,1.11328400000000,18.7917340000000,18.7455560000000,14.7004220000000,5.48735100000000,7.88760400000000,5.05202600000000,1.50810900000000,7.93643500000000,21.0473280000000,19.0777760000000,18.1786190000000,10.7223780000000,2.08765400000000,22.1667310000000,21.1513720000000,2.95926700000000,7.85898100000000,22.7172550000000,14.6453530000000,7.61366400000000,20.7518110000000,20.8092420000000,17.6495750000000,11.7747030000000,13.2957060000000,14.6103290000000,11.2736330000000,2.06767800000000,3.88528300000000,7.83178800000000,6.21878100000000,7.26437400000000,8.54074500000000,10.4555830000000,22.0444230000000,7.28181600000000,20.3301190000000,20.3122140000000,22.0214830000000,8.53602500000000,0.923407000000000,7.70504700000000,5.83058800000000,22.4114980000000,21.2778760000000,3.26649500000000,11.6957580000000,6.37504400000000,0.982365000000000,22.5633050000000,4.41817000000000,14.6907990000000,1.01447600000000,11.7674930000000,15.3408130000000,1.94439600000000,13.8845050000000,9.75343400000000,1.98109000000000,10.5761400000000,19.8022300000000,15.9582030000000,13.6056120000000,13.0858430000000,15.4579860000000,6.05202000000000,15.9511590000000,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
obj=spatialgraph2D(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

and the following error:

Index exceeds the number of array elements (408).

Error in spatialgraph2D (line 18)
              w=sqrt((x(s)-x(t)).^2 + (y(s)-y(t)).^2);

Error in loadingFiles_CY_senpai (line 755)
obj=spatialgraph2D(G,x,y);

My view on these 2 error could be:

1. Too many nodes?

2. There are some imperfect polygons in my undirected network?

Hope these will not cause a misunderstanding when I recognize that what I am confused are the basic knowledge I should already know.

I'm really appreciated for your guidance and kind help. Thank you a lot!

[Matt J]

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.

[Matt J]

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]

Never mind the last version. Try the one attached here. It should be faster.

[Cheng Ting Tsai]

Hi, Matt~

Thanks for your attention on my case and questions.

I'll try your new version code as soon as possible.

[Cheng Ting Tsai]

Hi Matt~

Finally, I complete my real case by using your newest version "spatialgraph2D.m" on MathWorks file exchange! Really really thanks a lot. I will cite you in my research. Cheers!

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

Take a simple case (Figure 1.) as an example:

there are 10 nodes (blue) in my graph, but there are only 7 nodes which are linked by edges. In other words, I have 3 nodes without linking as "marginal men". Now, code "spatialgraph2D.m" would not recognize these 3 marginal men as "nodes". However, because of the ID number (blue) of the last nodes is "10", G = graph(s,t) would mistankly believed that the total number of nodes is 10 ( is 7 actually). Therefore, the lable (ID number) of my nodes in x=[] and y=[] will get some trouble:

Error using spatialgraph2D (line 45)
Input vectors x and y must have lengths equal to the
number of graph nodes.

Error in loadingFiles_CY_senpai (line 845)
obj=spatialgraph2D(G,x,y);

In other words, I will get an error if my graph have nodes without linking. Therefore, renewing ID numbers (red) in my case is necessary. Hope this will help someone who is confused like me.Thanks!

[Figure 1.]

[Matt J]

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)

Answer 2

[darova]

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

[Sim]

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]

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]

Thanks Darova!

[darova]

Here is an example. Any ideas how to remove branches?

See the attached script

[Sim]

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!