Points and Polygons

March 31st, 2011

Ever wanted to know if some points of interest were inside a region? You can answer that if the region is specified as a polygon in MATLAB. The key is the inpolygon function.

Polygons
Some Points of Interest
Test the Points with the Polygon
Multiply-connected Polygon
Do You Work with Points and Polygons?

Polygons

Polygons in two dimensions are generally represented in MATLAB with two arrays, locations for the X vertices and Y vertices. There is no need to have the final points in these match the initial points; that is, when arrays as described are used in situations where they are interpreted as polygon vertices, the polygon is automatically closed.

Let's define a very simple polygon to start.

X = [0 0.5 1]';
Y = [0 0.5 0]';

And look at it.

patch(X,Y, [0.2, 0.7, 0.8], 'edgecolor','r',...
    'facealpha', 0.2, 'linewidth',2);
axis([0 1 -0.2 1])

Some Points of Interest

Let's create some points of interest now. Some inside, some outside, and some on the boundary.

xin = [0.3 0.75 0.82]';
yin = [0.25 0.1, 0.05]';
xout = [0.3 0.75 0.82]';
yout = [-0.15 0.6, 0.3]';
xedge = [0.3 0.75 0.82]';
yedge = [0.3 0.25, 0.18]';
hold all
plot(xin, yin, 'g*', 'markersize',5)
plot(xout, yout, 'mx', 'markersize',9)
plot(xedge, yedge, 'bd', 'markersize',7)
hold off

Test the Points with the Polygon

[inIN onIN] = inpolygon(xin,yin, X, Y)

[inOUT onOUT] = inpolygon(xout,yout, X, Y)

inOUT =
     0
     0
     0
onOUT =
     0
     0
     0

[inEDGE onEDGE] = inpolygon(xedge,yedge, X, Y)

inEDGE =
     1
     1
     1
onEDGE =
     1
     1
     1

Multiply-connected Polygon

Here's an example from the help for inpolygon, with a square containing a square hole. The outer loop is counterclockwise, the inner clockwise.

xv = [0 3 3 0 0 NaN 1 1 2 2 1];
yv = [0 0 3 3 0 NaN 1 2 2 1 1];
x = rand(1000,1)*3; y = rand(1000,1)*3;
in = inpolygon(x,y,xv,yv);
plot(xv,yv,x(in),y(in),'.r',x(~in),y(~in),'.b')

Do You Work with Points and Polygons?

Have you used inpolygon? I'd love to hear the contexts of the problems you solve with it. Let me know here.

Get the MATLAB code

Published with MATLAB® 7.11

By Loren Shure

You can follow any responses to this entry through the RSS 2.0 feed. Both comments and pings are currently closed.

12 Responses to “Points and Polygons”

Rodney Thomson replied on April 1st, 2011 at 01:15 UTC :

I have used inpolygon for data trimming purposes. Typically using ButtonDownFcn’s to draw and select a bounded region on a plot, then using inpolygon to identify and modify any points within said reason.

JuanPi replied on April 3rd, 2011 at 20:22 UTC :

Hi!
Thanks again for all your super interesting post.

I use inpolygon and related functions quite often (I find geom2d and 3d from David Legland and excellent extension to MATLAB core functions http://www.mathworks.com/matlabcentral/fileexchange/7844-geom2d).

What I am really looking forward to is a MATLAB version of alpha shapes http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Alpha_shapes_2/Chapter_main.html
I searched Central, but so far nothing really satisfying.
Do you know any source?

Thanks again

Damian Sheehy replied on April 4th, 2011 at 21:20 UTC :

Hi Juan,

I am not aware of a fully functional Alpha-shape algorithm for MATLAB, but this is a feature that is occasionally requested by our users. I would be happy to file an enhancement request for this feature on your behalf. To help us assess real need from a casual “it would be nice to have”, it is always helpful for us to know how our customers would use the feature. If you can share details about your use case we would be happy to hear from you. We would also like to understand your preference for alpha shapes as opposed to a curve reconstruction algorithm. You may respond to me directly or file the request with customer support.

Thanks for your feedback!

Damian

Chao-Wei Chen replied on April 7th, 2011 at 01:16 UTC :

I didn’t know inpolygon before, but love it!
Alternatively I would use roipoly to create a mask M.
To know the number of points inside M, maybe

ROI= and(IMG, M);
Num = sum(ROI(:));

Matt Kindig replied on April 7th, 2011 at 15:22 UTC :

@JuanPi

There is an alpha shape function on the FEX that I’ve used. I didn’t develop it, but have used it and found it very useful. Check http://www.mathworks.com/matlabcentral/fileexchange/6760-ashape-a-pedestrian-alpha-shape-extractor

@Loren

I have found that the built-in inpolygon is fairly slow. I have found the FEX submission “inpoly” to be much faster. It is here: http://www.mathworks.com/matlabcentral/fileexchange/10391-fast-points-in-polygon-test

-Matt

JuanPi replied on April 8th, 2011 at 10:10 UTC :

Hi all,
Thank you very much for reading my post an your comments.

@Damian Sheehy: I would gladly follow your suggestion if I knew your e-mail or had a link to exactly where to post the request. I am sorry for my ignorance, I tried searching and could not identify the right place. Thank you very much.
I have already discuss it on this forum, but basically I am studying the “topology” (lets call it geometry for the moment) of solutions of nonlinear hybrid dynamical systems (nonlinear sets of ODE’s connected with “if” statements). The “properties” I am interested in are represented as scattered points in R² or R³. Of course, the regions where the different properties appear are connected (not simply) and I am interested on their boundaries. I have no prior knowledge of the geometry of that border and therefor the alpha shapes provide a way to approximate it without mayor complications. The more simulations I made, the better define the boundary.
I hope you got an idea form that compressed description.
In other words you could think that you are given a set scattered points ad somebody ask you “I know this points belongs to some unknown regions(they are grouped). Could you find the geometrical boundary of those regions for me?”. I think this is almost the motivation of the alpha shape concept.
I have a post in my blog explaining some more details http://lavidasegunjuanpi.blogspot.com/2010/06/cgal-python-octave.html
I manged to use CGAL library with python and from there do a “system call” to use it in Matlab or Octave (I use Matlab only at the University).

@Matt Kindig: Thank you very much for the link, I am testing it, as I did with the previous ones. I will post feedback in the File Exchange post.

Loren replied on April 8th, 2011 at 16:39 UTC :

JuanPi,

On the right side of this blog, you will see a link for technical support. Follow that and you get the chance to enter a bug report or enhancement suggestion.

–Loren

Damian Sheehy replied on April 8th, 2011 at 17:23 UTC :

Juan,

Thanks for a providing an example that outlines your use of Alpha shapes. You can follow the steps Loren suggested to file an enhancement request or you reach me at dsheehy#mathworks*com (make the usual substitutions # = @ and * = . ). Feedback like this helps us to gauge real user need. Example applications such as research papers and reports help prioritize aspects of the feature that make it useful.

Damian

Noosh replied on November 16th, 2011 at 22:42 UTC :

I would like to use inpolygon() for a set of data points that do not form a well defined geometrical shape. Since inpolygon requires the vertices of the shape, I have to first find a way to detect the approximate vertices of my dataset. I have tried:
(1) convhull(): due to the asymmetrical shape of my data it generates large gaps. Here’s a simplified example:

SIGMA1 =[0.05 0;0 0.05];
MU1 =[1 2];
X = [mvnrnd(MU1,SIGMA1,2000)];
X=[X;2,1; 1.7 1.6; 1.69 1.5; 1.75 1.6; 1.5 1.5;1.6 1.5 ];
scatter(X(:,1),X(:,2)); ylim([- 3 3]);
hold on;
hullIndexes = convhull(X(:,1),X(:,2));
hullx = X(hullIndexes,1);
hully = X(hullIndexes,2);
plot(hullx,hully);
text(1.4,1.2,'THIS GAP IS BIG!','Color','r');

(2): alpha shape solution that Juan mentioned: but my system goes to not responding mode when I run ASHAPE() on my large dataset (>150K data points)
(3) I also tried converting the data points to image and then use BWBOUNDARIES() to find the boundary of my shape, but the image looks like sparse pixels (each corresponding to one data point) and BWBOUNDARIES treats each pixel (data point) as a separate object, and I don’t know how to resolve this issue.

So my question is can you help me find a way I can find vertices in my dataset that approximately follow the shape of my input data, so I can use the inpolygon() function?

I really appreciate your help.
Noosh

Damian Sheehy replied on November 17th, 2011 at 19:15 UTC :

Hi Noosh,

I would recommended Alpha Shapes for your particular problem. Since the application on the File Exchange didn’t work out for you I would suggest the following alternative.
If you triangulate the set of points you will notice the triangles lying near the convex hull are large and will have large circumcircles. If we were to identify and remove these triangles, we would get a better approximation to your region of interest.

We can use DelaunayTri to triangulate the data and compute the radius of the circumcircles for all triangles.
RCC is of size 1:NumOfTriangles
Here are the steps:

dt = DelaunayTri(X)
[CC, RCC] = dt.circumcenters();

I did some checks and chose a threshold radius of 0.2.
So, let’s find the indices of the triangles whose radius is greater than that.

idx = find(RCC > 0.2)

Next , we want to get hold of the triangulation array and delete the triangles that exceed the threshold.

tri = dt.Triangulation;
tri(idx,:) = [];

OK, create a TriRep to represent this triangulation and if we plot the TriRep we can see it represents a domain that is closer to the one you want.

tr = TriRep(tri,X)
triplot(tr)
hold off

The TriRep can give you the boundary of this triangulation using the freeBoundary method and you could use that to represent your polygon. However, it will only chain the segments together if you have a single region with no holes. Alternatively, we can impose these boundaries on the DelaunayTri we just created and then use it to perform an in-polygon test.
Here’s how:

constrainededges = tr.freeBoundary();
dt.Constraints = constrainededges;

These constraints separate the domain into interior and exterior triangles.
Get in/out status of each triangle as follows:

inout = dt.inOutStatus();

OK, let’s see if this works. Generate a grid of points and perform the in/out test and plot the inside and outside points against the triangulation.

[Xq, Yq] = ndgrid(0:0.2:2, 1:0.2:3);
Xq = Xq(:);
Yq = Yq(:);

Find the indices of the triangles that contain these points. If the index is NaN it means it’s outside the convex hull, so get rid of those points.

tids = dt.pointLocation(Xq,Yq);
nanidx = isnan(tids);
Xq(nanidx) = [];
Yq(nanidx) = [];
tids(nanidx) = [];

Find the indices of the remaining points that lie outside our region of interest

outidx = inout(tids)==0;
Xq(outidx) = [];
Yq(outidx) = [];

figure
plot(hullx,hully);
text(1.4,1.2,'THIS GAP IS BIG!','Color','r');
hold on
triplot(tr)
plot(Xq,Yq,'or')
hold off

Catherine Robin replied on December 8th, 2011 at 10:01 UTC :

Hi,

Damian, thank you for that, it has been very useful.

One question – this works well for a set of very evenly spaced grid points. The difference between large and small triangles is clear. What about the case where some holes we would like to identify are smaller than some of the grid separations? I am looking at FEM ocean models, with large elements offshore and small elements with islands near shore. Any suggestions?

Thanks
Catherine

Damian Sheehy replied on December 8th, 2011 at 13:21 UTC :

Hi Catherine,

Your FE model defines the domain of your polygon. Why throw away this information and then attempt to recreate it? You can import the nodes and elements into MATLAB. Then create a TriRep from this data, the TriRep/freeBoundary will give you the outline of the polygon. If you wish to do a point-in-polygon test on this region, you could create a constrained DelaunayTri from the boundary points and segments. Triangulate the points and constrain the segments. You do not need to triangulate interior points. Is your FE model is composed of quadrilateral elements, split them into triangles. For example, a quad with nodes (5, 27, 16, 35) can be split into two triangles, (5, 27 16) and (16, 35, 5). In general, if your quads are in an N-by-4 array format, your triangles will be:

tri = [quads(:,1:3); quads(:, (3, 4, 1))

Damian

Loren on the Art of MATLAB

March 31st, 2011