You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
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?
6 views (last 30 days)
Show older comments
Steve
on 9 Dec 2019
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!
Accepted Answer
Matt J
on 9 Dec 2019
Edited: Matt J
on 9 Dec 2019
Using one of my recent FEX postings,
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
24 Comments
Steve
on 9 Dec 2019
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
on 9 Dec 2019
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
on 9 Dec 2019
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
on 9 Dec 2019
Edited: Matt J
on 9 Dec 2019
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
on 10 Dec 2019
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
on 10 Dec 2019
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
on 10 Dec 2019
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
on 10 Dec 2019
Edited: Matt J
on 10 Dec 2019
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
on 11 Dec 2019
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
on 11 Dec 2019
Edited: Matt J
on 11 Dec 2019
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
on 11 Dec 2019
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
on 11 Dec 2019
Edited: Matt J
on 11 Dec 2019
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
on 11 Dec 2019
Edited: Matt J
on 11 Dec 2019
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
Matt J
on 11 Dec 2019
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.
Matt J
on 11 Dec 2019
Wow, really? You managed to compute the areas and perimeters from the [C1,V1,V2,C2] already? Good job!
Steve
on 11 Dec 2019
The number of sides, perimeters, and areas are already listed in the variable "output" (see screenshot). Aren't those correct?
Matt J
on 11 Dec 2019
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
on 11 Dec 2019
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
on 12 Dec 2019
That is correct. Where it says
%%%%CALCULATE STUFF%%%%
is where you would compute the bulgy areas that we talked about above.
Steve
on 12 Dec 2019
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?
Steve
on 12 Dec 2019
Edited: Steve
on 13 Dec 2019
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.
More Answers (1)
Catalytic
on 9 Dec 2019
Edited: Matt J
on 9 Dec 2019
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)
See Also
Categories
Find more on Sparse Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)