How can I use graph object info to find: areas, number of sides, and perimeter lengths of a each connected polygon within a network of polygons?

I was told by an MVP in the MathWorks community, "You can save the graph object, G, which contains all the connection information. You should read up on the capabilities of graph objects, as they allow you to do many things!

https://www.mathworks.com/help/matlab/ref/graph.html#d117e525223". However, I was never told how to do this and the link does not provide information on how to do this either. Can someone explain how I could use the attached graph object info to find: areas, number of sides, and perimeter lengths of a each connected polygon within a network of polygons?

Thanks in advance for your help!

Answer 1

[Matt J]

Using one of my recent FEX postings,

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

you should be able to convert a TripletGraph object tg to a polyshape array.

function [pshape,pgonNodes] = getPolyshapes(tg)
%
%IN:
%
%  tg: TripletGraph object
%
%OUT:
%
% pshape: polyshape array. Each pshape(i) is a constituent polygon.
% pgonNodes: A cell array such that pgonNodes{i} contain the node indices of the
%            i-th polygon.

  G=tg.CGraph;
  x=tg.C(:,1);
  y=tg.C(:,2);

  [pshape,pgonNodes]=polyshape( spatialgraph2D(G,x,y) );

end

[Steve]

Thanks Matt. Would this be implemented in your existing code? If so, where would it be placed and how would it be used to extract the needed information (i.e., polygon areas, number of sides, and perimeter lengths)? Attached is the latest version of your code (TripletGraph2.m) that I have been using.

[Matt J]

No, you are free to put this function wherever you wish. Once you get the array of pshape objects from the function, you can use polyshape methods, as demonstrated for you by Catalytic, to get areas, perimeters, etc...

Be mindful, however, that this analysis will treat your network as consisting of actual polygons. In other words, the area and perimeter calculations will be those of true polygons (i.e., with straight sides) and not the shapes you get by connecting the points with circular arcs, as you have done in other threads.

[Steve]

Is it possible to get the true values (i.e., for the polygons with circular arc edges) using this or any method that you know of?

[Matt J]

Well, first of all you need to come up with better terminology for it than "polygons", because polygons by definition have straight sides. You will invite much confusion otherwise. I will call them piecewise circular cells (PCCs).

For the PCC case, getPolyshapes will still give you useful information. It will output for you all the node numbers that participate in each PCC. Then, it's just a geometry game. I believe you already have the radius and central angle (in radians) of each of the circular arcs, so each arc length will be,

ArcLength = *Radius.

Then you just add up the arc lengths of all the sides of the PCC to get the perimeter.

For area calculation, it will be a bit more complicated, but not too bad. You can get the bulk area of the polygon (with straight sides) from the pshape output of getPolyshapes and then add/subtract the area of the bulgy orange part of any given circular arc like in the figure below. The area will need to be added or subtracted, depending on whether the arc bulges out of the polyshape or inward. You can check which way the arc bulges by testing whether one of the end points V1,V2 lies inside the polyshape or not. This can be done with the isinterior() function for polyshapes.

[Steve]

Matt,

Here's one way we could get the required info--please let me know what you think and how we could implement this in Matlab, in conjunction with your existing scripts:

We have a list of all edges (IDs 1...2750 or so) with the corresponding

IDs of two central points (obj.C) belonging to each. The goal is to construct

n-tuples of edge IDs that describe PCCs (e.g. maybe the 5 edge IDs

10, 35, 47, 98, 110 are the edges of one of the closed PCCs). One algorithm

that does that is:

• Make a second list of the three (at most) central points connected to each other central point.

• Start a for loop with the first edge ID, say E1, choose one end's central point (C ID, say C1); make that point the head of a vector pointing there from the other central point (chord vector).

• Find the other two central points connected to this head point (from the other list) and the corresponding vectors that point from the head to them (two edges, two chords).

• Take cross products of the first chord with each of the other two.

• Only one of these will be a positive cross product (counterclockwise turning); retain the ID of that edge and concatenate it to the edge ID we started with; we now have a 2-tuple [E1 E2]; also retain the ID of its central point to make [C1 C2].

• Keep doing this iteratively (in a while loop), with E2 taking the role of E1, etc. After every step [E1 ... Ek], check if Ek=E1. If yes, we have gone too far, and BE=[E1 ... E(k-1)] is the polygon edge set, i.e., n=k-1; likewise, you get BV=[C1 ... Cn] as the vertex set. Store BE,BV as new variables and exit the while loop.

• Close the for loop by taking the next edge ID, and do it for all of them.

• This will generate a lot of duplicate polygon variables BE(), BV(), but that's OK. Discard duplicates.

• Once we have the BE, BV vectors of edge IDs, the length of the vectors is the topology of the PCC.

• The area is computed in two steps: the polygonal area (straight edges) is obtained by the standard formula using the (x,y) coordinates of the central points belonging to the IDs of BV, or by the method you had described above.

• The circle segment areas within the arcs are either added or subtracted from that. Here we may have to go back to the way the theta is determined: we computed inner products of chords and end point directors (right?), but in the end we must have discarded their sign. We should keep the sign, because it tells us whether the arc curves counterclockwise or clockwise with respect to our edges as we go around the PCC. If we consistently concatenate our edges (as described above) by turning ccw at each central point, a ccw turning arc will subtract, a cw arc will add.

If this makes sense, how can we implement this into the existing code? Thanks again for all your help!

[Matt J]

I'm not sure if you saw my last comment. Using getPolyshapes in conjunction with things you've already computed, you pretty much have everything you need.

[Steve]

I like your approach better. I'm unsure of how to write that up in Matlab code though. Mainly because I am unsure of how to call on and/or pass objects from your code as opposed to passing standard, traditional, variables.

[Matt J]

Mainly because I am unsure of how to call on and/or pass objects from your code as opposed to passing standard, traditional, variables.

Well, let's take stock of what you have. For a given pair of node numbers, s and t, you have the means to extract an edge ID for the arc connecting them:

edgeID=tg.Cgraph.findedge(s,t);

Here, tg is once again a TripletGraph object from our earlier endeavors.

Once you've obtained the edge ID, you have the means to extract the coordinates C1,V1,V2,C2 corresponding to the center points and end points in my diagram above

C1=tg.ArcPoints(:,1,edgeID);
V1=tg.ArcPoints(:,2,edgeID);
V2=tg.ArcPoints(:,3,edgeID);
C2=tg.ArcPoints(:,4,edgeID);

Now, notice that C1,V1,V2,C2 are all just standard numeric Matlab variables, so these should be very friendly and familiar to you. But once you have these 4 points, you have everything you need to analyze the geometry of one of the arcs in a PCC. The rest of it should really just be a matter of repeating these operations in appropriate loops.

[Steve]

Matt, as usual...looks great. you're amazing. I can see we are almost there. I keep getting an "undefined variable" error pertaining to getPolyshapes. Can you explain what getPolyshapes is? Is it a built-in Matlab command? I'm also getting errors with tg. Was tg defined somewhere else?

[Matt J]

tg is the TripletGraph object that you created from your fpep data, e.g., the one in the .mat file you attached to your post.

getPolyshapes is not a native Matlab function. It is the function that I wrote for you in my original answer,

function [pshape,pgonNodes] = getPolyshapes(tg)
%
%IN:
%
%  tg: TripletGraph object
%
%OUT:
%
% pshape: polyshape array. Each pshape(i) is a constituent polygon.
% pgonNodes: A cell array such that pgonNodes{i} contain the node indices of the
%            i-th polygon.
  G=tg.CGraph;
  x=tg.C(:,1);
  y=tg.C(:,2);
  [pshape,pgonNodes]=polyshape( spatialgraph2D(G,x,y) );
end

[Steve]

Thanks for your response and your continued effort Matt. I still cannot get anything to run though, mainly (i think) because I am not proficient enough to piece all the snippets of code together into their proper places. I downloaded your spatialgraph2D file, but I am unsure of where you got the tg file as I don't believe I ever uploaded a file with that name--did you rename the graph_object_G.mat file (that I uploaded in the beginning of this thread) to tg? I would really love to see your code work over here on my end. Are you able to attach all the required code as a single file so I could run it? Thanks again for all your help!

[Matt J]

Nothing we've been discussing involves any code other than TripletGraph.m (orTripletGraph2.m, but you really should rename it), spatialgraph2D.m, and the above short getPolyshape() function.

tg is not the name of a file. It is just an arbitrary hypothetical name I chose for a variable generated using TripletGraph.m, for example,

load fpep.mat
tg=TripletGraph(fpep);

You can of course choose any names you wish for your variables.

Since I have told you several times that tg is a "TripletGraph object", I am going to assume that you did not understand that terminology. It is important that you do, however. Every variable in Matlab has a type or "class". These include traditional kinds of Matlab variables like doubles, chars, cells, and structs. It is also possible to define our own classes, however. TripletGraph is not just the name of the file that created, tg. It is also the name of its class - a class defined by us using the classdef command. You can see the classes of different variables using either the whos or class commands, for example,

>> a=1:3; b='ABC';

>> class(a)

ans =

    'double'

>> class(b)

ans =

    'char'

>> class(tg)

ans =

    'TripletGraph'

>> whos a b tg
  Name      Size             Bytes  Class           Attributes

  a         1x3                 24  double                    
  b         1x3                  6  char                      
  tg        1x1             191213  TripletGraph              

When I told you that tg is a TripletGraph object, this simply means that it is a variable of type 'TripletGraph'.

[Matt J]

So, here is how I envision the code looking at a high level,

load fpep

tg=TripletGraph(fpep);

[polygons,vertexLists]=getPolyshapes(tg);

N=numel(polygons);

output(N).numsides=[];   %Return output as a structure array
output(N).perimeter=[];
output(N).area=[];

for i=1:N  %For each of the PCCs
    
    p=polygons(i);
    v=vertexLists{i};  
           v=[v,v(1)]; %make list wrap back to the beginning
    
    for j=1:numel(v)-1  %For each edge of the PCC
    
     s=v(j); 
     t=v(j+1);
     edgeID=tg.Cgraph.findedge(s,t);
     
        
     C1=tg.ArcPoints(:,1,edgeID);
     V1=tg.ArcPoints(:,2,edgeID);
     V2=tg.ArcPoints(:,3,edgeID);
     C2=tg.ArcPoints(:,4,edgeID);
    
     
       %%%%CALCULATE STUFF%%%%
     
    end
    
    output(i).numsides=numsides(p);
    output(i).perimeter=_________;
    output(i).area=________;   
    
end

[Steve]

I'm still getting an error--see attached screenshot:

Any idea what's going on?

[Matt J]

Yes, you must put getPolyshapes somewhere that Matlab can see it... As I mentioned earlier, you are free to put this function anywhere that you would put one of your own functions.

[Steve]

Matt, It worked. You, my friend, are the best! Thank you so much!

[Matt J]

Wow, really? You managed to compute the areas and perimeters from the [C1,V1,V2,C2] already? Good job!

[Steve]

The number of sides, perimeters, and areas are already listed in the variable "output" (see screenshot). Aren't those correct?

[Matt J]

Depends how you got them. I left the area and perimeter calculations as blanks for you to fill in, so I don't know what you put there.

[Steve]

Ahhh, ok. I just put in perimeter(p) and area(p) as I thought those dots represented ditto marks. I guess this would give the values for only straight edges?

[Matt J]

That is correct. Where it says

 %%%%CALCULATE STUFF%%%%

is where you would compute the bulgy areas that we talked about above.

[Steve]

How close are we to being able to calculate the actual areas? Do you have a preferred approach for how to tackle this with respect to the code?

[Matt J]

Like I said here, I don't think it's too bad.

[Steve]

Well, we know the chord-lengths (c) and arc-lengths (s) of each edge. We also know the radii (r), and thus the central angles for the same. So, we can definitely calculate the segment areas via the following equation: circular segment area = (1/2)*(r^2)*((180/Pi)*alpha - sin(alpha)). I am just not sure how to have Matlab to do it (in terms of syntax, etc.). Also, I'm not sure how to set up isinterior() to determine concavity (for +/- areas). Have any suggestions on how to set this up? Thanks again.

Answer 2

[Catalytic]

Your attachment doesn't contain a graph object. However, perhaps it will help to mention that if you have the vertices of a polygon, you can use it to create a polyshape object.

With the polygon represented in polyshape form,

>> p=polyshape([0,0; 0 1;1 2; 1 0])

p = 

  polyshape with properties:

      Vertices: [4×2 double]
    NumRegions: 1
      NumHoles: 0
      

you can query all the kinds of things that you mentioned. For instance,

      
>> area(p)

ans =

    1.5000

>> perimeter(p)

ans =

    5.4142

>> numsides(p)

ans =

     4
     
>> plot(p)

[Steve]

Thanks for your response Matt, but I don't know how to input the vertices from your previous code's output files into the polyshape code. Can you explain how to do this?

[Steve]

Thanks also to Catalytic.