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?

[Sim]

Hi Can, sorry I have just noticed your post! I am available for answering your questions. How can I contact you 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

[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

[Matt J]

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

[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

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

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

[Jun W]

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

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.

[NA]

What is the time complexity of this algorithm?

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!

Answer 3

[Sim]

Hey @Matt J! hope this message finds you well! I am back to my post after a while :-)

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

[Matt J]

Thanks @Sim, I've uploaded a new version which hopefully resolves the bug.

[Sim]

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]

I think I know where it's coming from, but please attach your test case data (G,x,y) in a .mat file. Also, please use the Discussion section of the FIle Exchange page for spatialgraph2D to report problems, rather than the Answers forum.

[Sim]

Ok, thanks a lot!

Please see attached the following file for you :-)

obj.mat

I have just added this thread to the Discussion section of the File Exchange as well ! :-)

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]

P.S. it looks like matlab does not allow me to edit my post in the Discussion section of the File Exchange, therefore, if something coems to my mind, I write here... sorry..

Yes, I'm not sure why that hasn't been fixed. What I do as a workaround is to copy the post to the clipboard, then delete and recompose it.

[Sim]

Yes, I do the same, copy to the clipboard, delete from the discussion panel and paste it there with my modifications... but I think that the matlab staff already knows it, and I am confident they will fix it soon :-)

[Sim]

@Matt J, any idea on how to solve this small issue ? :)