It's often useful to know if a point is inside of a polygon (two dimensions) and sometimes if a point is inside of a polyhedron (three dimensions). For two dimensions, MATLAB ships with inpolygon, a nice function to handle this. However, this same operation in three dimensions becomes more complicated. It is simplified if the object is convex. In this case, the vertices can be represented as a set of constraints and we can apply these constraints to the points to test whether they are inside or not. Matt J does this in his vert2lcon function (The runner-up, who will also get some MathWorks swag, for this week's post!).
But what about when the polydehron is not necessarily convex? Now we no longer have this luxury and have to resort to more complex algorithms.
Fortunately, Sven has done the work for us!
Let's use Paul's three dimensional L-Shaped Membrane.
% Load shellVertices, shellfaces from Paul's Blog: load PaulsMembrane.mat % Draw the membrane: patch('Vertices',shellVertices,... 'Faces',shellFaces,... 'FaceColor','r'); axis tight view(-51,24)
Add 10000 random points to it
pts = rand([10000,3]); %Generate a random 10000pts hold on; plot3(pts(:,1),pts(:,2),pts(:,3),'b*'); %add points to plot
Find the points in the polyhedron
in = inpolyhedron(shellFaces,shellVertices,pts,'FlipNormals',true); delete(hLine); %clean up plot3(pts(in,1),pts(in,2),pts(in,3),'b*'); %add points to plot set(hPatch,'FaceAlpha',0.3,... 'EdgeColor','none'); %make the points inside visible
Get the MATLAB code
Published with MATLAB® R2013a
5 CommentsOldest to Newest
Hei ho, thanks Sean!
I’d like to ask the opinion of anyone here who has an opinion on the topic: “what IS the expected face orientation of a closed polyhedron?”
I have a feeling (since many other softwares use this) that face normals should point OUT of a solid object.
In inpolyhedron, I went with the opposite convention, but only because a call to isosurface() with a binary mask produces faces pointing IN to the solid object.
It turns out that isosurface() actually has the (quite reasonable) convention of producing faces pointing towards the higher region in its original mask (ie, for a binary mask it points from zeros outside the object to ones inside the object).
So… I’m beginning to think that I should update inpolyhedron to fit the wider convention of expecting face normals pointing OUT.
So two questions:
1. What do YOU expect as face/normal convention. Faces IN or faces OUT?
2. If it’s faces OUT, should I update inpolyhedron to match that convention?
Looks like a great function to extend inpolygon. I’ll certainly check it out.
On the topic of inpolygon, though, it works in cartesian coordinates and I’m often interested in knowing whether points on the surface of the earth are inside a polygon. Using latitudes and longitudes is there a good alternative to inpolygon? The option I see first is to simply interpolate the polygon edges using either great circle or rhumb lines as appropriate for the case and then use the inpolygon cartesian results as a close approximation. Any better methods?
@Sven, You’re welcome!
As far as the normals, I don’t really have an opinion. I was surprised when the first run through yesterday yielded all of the points outside of the membrane. However, the FlipNormals option painlessly took care of that.
Interesting question and it lead to a good whiteboard discussion here with the neighbors. I think your idea is the most straightforward approach. Use INTERPM from the Mapping Toolbox to interpolate between lat/lon vertices (probably using the ‘gc’ or great circle interpolation option). You can refine your polygons to the desired tolerance so that you can then use INPOLYGON.
Sven, I’d probably vote for face normals pointing out. Cheers!
Ok, I’ve just submitted a new version that has a small speedup/bugfix plus the FACE NORMAL CONVENTION CHANGE.
From now, inpolyhedron() expects face normals to point out. This seems to be the wider convention, and will avoid the “surprise” that Sean encountered when first running.
Alert to current users: this new version is better (with the bugfix) but be aware of this change and how it affects your code. For example, Sean’s post above no longer needs (‘FlipNormals’,true) option to be added.