Code covered by the BSD License  

Highlights from
inpolyhedron - are points inside a triangulated volume?

4.66667

4.7 | 12 ratings Rate this file 66 Downloads (last 30 days) File Size: 22.2 KB File ID: #37856
image thumbnail

inpolyhedron - are points inside a triangulated volume?

by

 

20 Aug 2012 (Updated )

Test if 3d points are inside a mesh. Or, voxelise a mask from a surface. Mesh can be non-convex too!

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

| Watch this File

File Information
Description

INPOLYHEDRON Tests if points are inside a 3D triangulated (faces/vertices) surface

User's note:
inpolyhedron adopts the widely used convention that surface face normals point OUT from the object. If your faces point in, simply call inpolyhedron(...,'flipNormals',true).
(see discussion at http://blogs.mathworks.com/pick/2013/09/06/inpolyhedron/)

  IN = INPOLYHEDRON(FV,QPTS) tests if the query points (QPTS) are inside the
  patch/surface/polyhedron defined by FV (a structure with fields 'vertices' and
  'faces'). QPTS is an N-by-3 set of XYZ coordinates. IN is an N-by-1 logical
  vector which will be TRUE for each query point inside the surface.

  INPOLYHEDRON(FACES,VERTICES,...) takes faces/vertices separately, rather than in
  an FV structure.

  IN = INPOLYHEDRON(..., X, Y, Z) voxelises a mask of 3D gridded query points
  rather than an N-by-3 array of points. X, Y, and Z coordinates of the grid
  supplied in XVEC, YVEC, and ZVEC respectively. IN will return as a 3D logical
  volume with SIZE(IN) = [LENGTH(YVEC) LENGTH(XVEC) LENGTH(ZVEC)], equivalent to
  syntax used by MESHGRID. INPOLYHEDRON handles this input faster and with a lower
  memory footprint than using MESHGRID to make full X, Y, Z query points matrices.

  INPOLYHEDRON(...,'PropertyName',VALUE,'PropertyName',VALUE,...) tests query
  points using the following optional property values:

  TOL - Tolerance on the tests for "inside" the surface. You can think of
  tol as the distance a point may possibly lie above/below the surface, and still
  be perceived as on the surface. Due to numerical rounding nothing can ever be
  done exactly here. Defaults to ZERO. Note that in the current implementation TOL
  only affects points lying above/below a surface triangle (in the Z-direction).
  Points coincident with a vertex in the XY plane are considered INside the surface.
  More formal rules can be implemented with input/feedback from users.

  GRIDSIZE - Internally, INPOLYHEDRON uses a divide-and-conquer algorithm to
  split all faces into a chessboard-like grid of GRIDSIZE-by-GRIDSIZE regions.
  Performance will be a tradeoff between a small GRIDSIZE (few iterations, more
  data per iteration) and a large GRIDSIZE (many iterations of small data
  calculations). The sweet-spot has been experimentally determined (on a win64
  system) to be correlated with the number of faces/vertices. You can overwrite
  this automatically computed choice by specifying a GRIDSIZE parameter.

  FACENORMALS - By default, the normals to the FACE triangles are computed as the
  cross-product of the first two triangle edges. You may optionally specify face
  normals here if they have been pre-computed.

  FLIPNORMALS - (Defaults FALSE). To match a wider convention, triangle
  face normals are presumed to point OUT from the object's surface. If
  your surface normals are defined pointing IN, then you should set the
  FLIPNORMALS option to TRUE to use the reverse of this convention.

  Example:
      tmpvol = zeros(20,20,20); % Empty voxel volume
      tmpvol(5:15,8:12,8:12) = 1; % Turn some voxels on
      tmpvol(8:12,5:15,8:12) = 1;
      tmpvol(8:12,8:12,5:15) = 1;
      fv = isosurface(tmpvol, 0.99); % Create the patch object
      fv.faces = fliplr(fv.faces); % Ensure normals point OUT
      % Test SCATTERED query points
      pts = rand(200,3)*12 + 4; % Make some query points
      in = inpolyhedron(fv, pts); % Test which are inside the patch
      figure, hold on, view(3) % Display the result
      patch(fv,'FaceColor','g','FaceAlpha',0.2)
      plot3(pts(in,1),pts(in,2),pts(in,3),'bo','MarkerFaceColor','b')
      plot3(pts(~in,1),pts(~in,2),pts(~in,3),'ro'), axis image
      % Test STRUCTURED GRID of query points
      gridLocs = 3:2.1:19;
      [x,y,z] = meshgrid(gridLocs,gridLocs,gridLocs);
      in = inpolyhedron(fv, gridLocs,gridLocs,gridLocs);
      figure, hold on, view(3) % Display the result
      patch(fv,'FaceColor','g','FaceAlpha',0.2)
      plot3(x(in), y(in), z(in),'bo','MarkerFaceColor','b')
      plot3(x(~in),y(~in),z(~in),'ro'), axis image

Acknowledgements

This file inspired In Polyhedron.

MATLAB release MATLAB 7.14 (R2012a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (32)
02 Dec 2014 Marc Lalancette

Amazing performance, great to have the grid input option. Sad about meshgrid (Y,X,Z) instead of ndgrid (X,Y,Z) type output though, can't simply convert to list doing IN(:).
Sept 2013 update indicates normals should point in, but help still says out...

31 Oct 2014 Peter Morovic

Thanks for sharing, Sven - works great indeed.

20 Jun 2014 yusuf

Thank you Sven.

This code works perfectly

02 Jun 2014 Vladimir

Hi Sven,
I am working on a project for my research where I would need to write a C++ program that calculates whether-or-not a point lies inside an STL object.

To do this, I would like to port you code over to C++, but I didn't want to do this without your permission since the BSD license says that I should not modify the source code.

If you would be alright with me porting it to C++ and using it for my own research, please let me know by emailing me at vladimir@ufl.edu .

Thanks a lot for making this program available!

Oh, and if it's not too much trouble, could you also include a .txt with any copyright and credits that you'd like me to include in my program.

Thanks again! :)
Vladimir H.

24 Feb 2014 Rheeda

This works brilliantly! Thank you!!!

10 Dec 2013 Gabriele

John, Sven, thanks for the reply.
The reason why I put "exactly" within quotation marks is because I intended implicitly a check considering a tolerance. Since I saw it is possible to specify a TOL parameter to inpolyhedron, I thought it could have been used also for checking whether a point is on the surface within a tolerance, hence the use of "exactly" ;-)
Anyway, again, compliments for the code!

10 Dec 2013 Sven

Hi Bart, I don't think this would be an easy extension since inpolyhedron uses Z-ray (vertical vector) intersections whereas surface distances would be calculated via arbitrary 3D face normal vectors.

Perhaps inpolyhedron would be a useful first filter if you considered all points inside a surface as having zero distance, but it wouldn't save computation for actually computing that distance.

It sounds like an interesting problem with similar performance issues to inpolyhedron. A related question with a proposed algorithm is here: http://stackoverflow.com/questions/18230259/computing-distance-from-a-point-to-a-triangulation-in-3d-with-matlab

10 Dec 2013 Bart Bolsterlee

Works perfectly, thanks! Could this easily be extended to also compute the distance of the point to the triangulated surface?

08 Dec 2013 John D'Errico

By the way, the code looks nice. Superb help. I love to see copious internal comments, enough that it is easy to follow what was done and how the code works.

08 Dec 2013 Sven

Hi Gabriele: John speaks the truth - there's no such thing as exactly on a surface when floating point numbers are involved and the best you can do is check within a tolerance.

To practically get to your query though I'm afraid that even that check will be a little tricky. inpolyhedron uses a separate check on the XY projection of a point (via dot products which range from -1 to 1) and then Z location of a point (using the actual coordinates), so two separate calls to inpolyhedron() with different values for the TOL parameter won't actually do quite what you're looking for.

Instead I suggest this:
1. Run inpolyhedron() on your points.
2. Use the vertexNormal() method of a triangulation object (available from 2012b, I think) to shift every vertex of your surface out/in by a tiny distance
3. Run inpolyhedron() on your points with this shifted surface.

Any points that differed in IN/OUT state between steps 1 and 3 can be considered ON the surface, at least within the tiny distance you specified.

08 Dec 2013 John D'Errico

Gabriele - you will never be able to identify points as EXACTLY on the surface, but only on the surface to within a requested tolerance by some metric.

08 Dec 2013 Gabriele

Hi Sven,
very nice and well working.
I'm wondering if you plan providing a specific output to identify those query points which are "exactly" on the surface.

Compliments!

10 Aug 2013 Sven

@Scott, @Yan,
I've just finished a tool that will help both of your input situations. Please use this: http://www.mathworks.com/matlabcentral/fileexchange/43013-unifymeshnormals

It will take your faces/vertices input and align all of the face normals so that inpolyhedron can do its job.

05 Aug 2013 Cranfield  
01 Jul 2013 Scott  
01 Jul 2013 Scott

Is there a way to find the coordinates of face normals that are pointed in different directions? Majority of my faces are pointed in the right direction but a few are in the opposite direction. I can use 'flipnormals' but that flips every face, not just the few I need flipped. Any help is appreciated.

14 Jun 2013 Sven

@Yan:
The problem here is just that your face normals are inconsistent. Only the bottom triangle (3rd face) is pointed IN, while the other 3 triangles are pointed OUT. If you replace:
face = [1,2,3;3,2,4;1,2,4;1,3,4];
with
face = [1,2,3;3,2,4;2,1,4;1,3,4];

and use (...,'flipnormals',true), then your input will work just fine.

23 May 2013 Yan Ou

I have tried to following case and it does not work
vertex = [0,0,0;2,0,0;1,0,1.75;1,2,0];
face = [1,2,3;3,2,4;1,2,4;1,3,4];
fv.vertices = vertex;
fv.faces = face;

25 Apr 2013 Sven

@Dun Kirk: Please send me your input points and triangles by email. This sounds like there could be something wrong and I would like to see if I can troubleshoot it.

25 Apr 2013 Dun Kirk

My FV has only 63 points and 122 triangles. But it took me 6.9 seconds to analyze 36144 query points and produced the wrong output.

09 Apr 2013 Benjamin Friedrich

Exactly what I needed. Worked out of the box.

07 Nov 2012 Sven

@ Ahmad: The latest version I just uploaded handles your input well now.

It was modified to detect the in/out of query points that are coincident with a mesh node.

I think there's probably a small (I haven't tested exactly how much) speed hit to make an extra comparison, but the new logic actually removes a loop so it will still be quite fast.

I know that there's one more potential issue: query nodes coincident with an edge that is acute (one connecting triangle pointing up, one pointing down)... but that's another issue for another time :)

28 Oct 2012 AP

For my case, inpolyhedron takes less than a second to analyze 60,000 points in the STL I have. STL has 7500 faces and 3800 points. Before inpolyhedron, I used Triangle/Ray Intersection http://www.mathworks.com/matlabcentral/fileexchange/33073-triangleray-intersection which took me about 80 seconds to analyze all 60,000 query points.

So, good job Sven.

28 Oct 2012 AP

It worked. Super fast. Great job. Thanks.

28 Oct 2012 Sven

@Ahmad:

Ok, I've found the issue. For the short term, try this for a solution:

in = inpolyhedron(myStl,X + 0.0001);

The reason you were having issues is that the faces/vertices you have are exactly aligned with the grid of your query points. This means that every query point lies exactly ON an edge. At the moment in my code, ON is considered OUT of the surface (only IN is considered IN). I had only anticipated this would be an issue for the outer surface points, but the same issue arises for inner points too like in your case.

In the help file I mention this as an issue awaiting user feedback. It looks like you're giving that feedback - thanks.

For the moment, just use the little work-around I gave. I think from your example I will try to change the definition so that I add a small tolerance to regard ON query points as actually IN to address this issue in a next version of inpolyhedron.

28 Oct 2012 Binu

@Ahmad:

Ok, I've found the issue. For the short term, try this for a solution:

in = inpolyhedron(myStl,X + 0.0001);

The reason you were having issues is that the faces/vertices you have are exactly aligned with the grid of your query points. This means that every query point lies exactly ON an edge. At the moment in my code, ON is considered OUT of the surface (only IN is considered IN). I had only anticipated this would be an issue for the outer surface points, but the same issue arises for inner points too like in your case.

In the help file I mention this as an issue awaiting user feedback. It looks like you're giving that feedback - thanks.

For the moment, just use the little work-around I gave. I think from your example I will try to change the definition so that I add a small tolerance to regard ON query points as actually IN to address this issue in a next version of inpolyhedron.

27 Oct 2012 AP

I have found the bug of the code. The function does not run the following loop in the code for some query points:

for ptNo = find(qPtHitsTriangles(:))'

which I think is because of the part in the code that gathers the unique XY query locations. In my case, since I have generated a grid of 3D points by meshgrid this would eliminate some of the points and will prevent that aformentioned for-loop. In the example, shown in the code, query points are random.

Now in the example code, instead of generating random points, create the query point according to the following. You will get all the points outside the domain.

[xq, yq, zq] = meshgrid(4:16,4:16,4:16);
pts = [xq(:), yq(:), zq(:)];

27 Oct 2012 AP

@Sven: When I run patch and look at the query points and the triangulated surface, there are many points that are located inside the surface. It is interesting that if the point is inside the surface, calling the function with one single query point will return the correct result.

27 Oct 2012 Sven

@A. Falahatpisheh: Send me a small sample of mesh/query points that don't work for you - I'm sure this can be trouble-shooted. Keep in mind that you should be able to call the following with your input:

figure, hold on, view(3), patch(fv,'FaceColor','g','FaceAlpha',0.2), plot3(pts(:,1),pts(:,2),pts(:,3),'bo','MarkerFaceColor','b')

If you make that call, and you see some of your points inside your volume, inpolyhedron() should be able to tell you which ones.

27 Oct 2012 AP

I flipped the face normals but it did not help. The output of the function is all zero. Any suggestion?

20 Sep 2012 Shawn Walker

Awesome utility. I needed this!

Note to users: if it doesn't look like it is working, you may have to flip the normals of the faces. See help file.

20 Sep 2012 Shawn Walker  
Updates
28 Aug 2012

A speedup tweak to distribute facets to grid locations - way to go, accumarray()

06 Nov 2012

Added logic to consider query points coincident with a vertex as IN (previously considered OUT)

22 Jul 2013

Much improved handling of gridded input. Now no longer unpacks the 3d grid. Saves big memory footprint!

26 Jul 2013

The icon was lost last upload...

11 Sep 2013

NEW CONVENTION ADOPTED to expect face normals pointing IN
Vertically oriented faces are now ignored. Results in speed-up and bug-fix.

23 Sep 2013

v3.1 - Dropped nested unique call made redundant via v2.1 gridded point handling. Refreshed grid size selection via optimisation.

24 Sep 2013

Put an upper limit on gridsize to prevent memory issues for huge trianulations.

27 Feb 2014

Fixed indeterminate behaviour for when query
points lie *exactly* in line with an "overhanging" vertex.

Contact us