Code covered by the BSD License  

Highlights from
intriangulation - which points are inside a 3d watertight triangulation?

4.7 | 3 ratings Rate this file 25 Downloads (last 30 days) File Size: 6.11 KB File ID: #43381 Version: 1.4
image thumbnail

intriangulation - which points are inside a 3d watertight triangulation?



05 Sep 2013 (Updated )

Are 3D-testpoints located inside or outside an arbitrary watertight mesh with vertices and faces?

| Watch this File

File Information

How often does this happen: You have a nice mesh in 3d, which is described by an np by 3 array of vertices and an nt by 3 array of indices into this array, which describe the faces. lets further say, the mesh is closed, i.e. it divides the total 3d space into a bounded and an unbounded domain in a way that it is impossible to connect any point inside with any point outside the domain without crossing the triangulation.

Now you have another fine set of points in 3d and you need to know which of these points are inside the triangulated volume and which aren't.

Intriangulation is the solution.

With the command

in = intriangulation(vertices,faces,testpoints)

you get an array of the same length as the testpoints with an entry 1 for testpoints inside the triangulated volume and an entry 0 for those points that aren't. (OK, you get -1 for points where the algorithm failed.)


help intriangulation

for an example and details.

In contrast to the FEX-submission inpolyhedron, you do not need to compute and/or preprocess any face normals in intriangulation.

The mesh may not have internal faces!


I have to announce that the core algorithm is 99% identical to the algorithm from Adam Aitkenhead in his submission "Mesh voxelisation". I use it here with his permission. Please go and rate "mesh voxelisation" also if you like intriangulation.

Also thanks to Sven's comments on algorithm and special cases.


Mesh Voxelisation inspired this file.

Required Products MATLAB
MATLAB release MATLAB 7.11 (R2010b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (15)
04 Jun 2015 James Johnson

Hey Johannes

I ran into some issues using your code. There are some points inside the volume that are recognized as out when they are in. The characteristic of the problem is that if I draw a volume and saturate it with points there is crevasse of points that is recognized as outside of the volume when it is in. My polyhedron has some very thin and long triangles which may be a factor, but I would like your advice.


Comment only
22 Apr 2015 Daniel Leib

Hi Johannes,

I hope you're still monitoring this submission! I'm using it in what might be a bit of a fringe case and am having a problem with memory. The overall use is that I have created two distinct surface meshes and also a tetrahedral mesh that represents the total volume of those two surface meshes. I want to be able to tell if the individual tetrahedra fall between those surface meshes or just within one, so I'm running a check for the centroid of each tetrahedron. If I run 10 checks, no problem. 100, no problem. When I have solid meshes up around 100k elements, though, I run into memory allocation errors, despite the fact that each call to intriangulation only involves a single surface mesh (which can be rather large) and one 3D point. Is there possibly a memory leak that might be happening somewhere? Or is there a fringe case that might cause one of the algorithms to freak out that I might be encountering by chance?

Thanks for any help you might be able to provide!

28 Nov 2014 John

John (view profile)

01 Sep 2014 Johannes Korsawe

The mesh generated by cylinder is not closed (or watertight). Close it by adding middlepoints on top and bottom and generate additional appropriate triangles to form the caps on the two ends.

Best regards,

Comment only
28 Aug 2014 Rene Labounek

Hi everyone,
I am trying to separate points outside and inside the cylinder, but it is working as I wish.

Please, do you have any idea what I should change?

I suppose that cylinder function generates only surface. I was trying to generate my own faces as all possible combinations between vertices but it did not help. The output of intriangulation function was same as for faces generated by surf2patch function.


[CYL1 CYL2 CYL3]=cylinder(5,200);
CYL = surf2patch(5+CYL1, 5+CYL2, 150*CYL3,'triangles');

testp=[5 5 50;2 1 20; 0.5 0.5 35; 1 1 100;0 5 0;10 5 150; 3 3 80];

hold on
hold off

in = intriangulation(CYL.vertices,CYL.faces,testp)

Comment only
17 Jun 2014 yusuf

yusuf (view profile)

Hi Johannes,

Thank you for your interest.

Could you rearrange the code to give the points exactly on the facets and on the vertex of the facet.

06 Jan 2014 Johannes Korsawe

Hi yusuf,
you are right, you can easily extend the program or make a check prior to intriangulation to give results for this case as a quick help.
I will fix this as soon as i come to time.

Comment only
30 Dec 2013 yusuf

yusuf (view profile)


I faced with some problem about program. When the point coincides the vertex point of the geometry the program can not determine whether the point is inside or outside. Moreover when I change the number of the points (increased) for the same geometry the program gives same error.

Comment only
31 Oct 2013 Johannes Korsawe

Acceleration is done but without domain decomposition methods as of Sven's submission, which i currently have no time to do.

The main advantage of intriangulation is the fact i do not need any orientation of the triangles which can be a significant (one-time) effort to do at least with my realworld examples.

Comment only
25 Sep 2013 Johannes Korsawe

your numbers are right. I did not have the latest version of your submission.

So i still have to accelerate intriangulation.... :-)

Comment only
24 Sep 2013 Sven

Sven (view profile)

Hi Johannes, are you using latest versions of each? After adding voxelise-ability to inpolyhedron, there was a version that essentially had duplicate calls to unique() (which is always the slowest part). I've got the latest versions of both on my computer now and I get the following.

load tetmesh, verts = X; faces = freeBoundary(TriRep(tet,verts));
for n = [1 10 50 100 500 1000 5000 10000 50000 100000 500000]
fprintf('For %6d test points: ',n); pts=bsxfun(@plus,45*rand(n,3),[0 0 40]);
tic, int=intriangulation(verts,faces,pts); fprintf('intri: %fs. ', toc)
tic, inp=inpolyhedron(faces,verts,pts); fprintf('inpol: %fs.\n',toc)

For 1 test points: intri: 0.021168s. inpol: 0.034748s.
For 10 test points: intri: 0.009391s. inpol: 0.017144s.
For 50 test points: intri: 0.007063s. inpol: 0.012926s.
For 100 test points: intri: 0.010297s. inpol: 0.015172s.
For 500 test points: intri: 0.045396s. inpol: 0.035927s.
For 1000 test points: intri: 0.095363s. inpol: 0.060970s.
For 5000 test points: intri: 0.468793s. inpol: 0.236190s.
For 10000 test points: intri: 0.926118s. inpol: 0.437674s.
For 50000 test points: intri: 4.608397s. inpol: 1.976781s.
For 100000 test points: intri: 9.166739s. inpol: 3.824124s.
For 500000 test points: intri: 46.20851s. inpol: 18.702466s.

I would be interested if you get particularly different results on a different machine - I've only optimised the grid selection based on my own computer but I think it should be quite general.

Comment only
24 Sep 2013 Johannes Korsawe

@Sven. I fail to recover your statement that for this very example inpolyhedron is faster. With my PC, it is exactly the other way round. Here are my results of your script:

For 1 test points: intri: 0.114747s. inpol: 0.124156s.
For 10 test points: intri: 0.013370s. inpol: 0.021715s.
For 100 test points: intri: 0.010008s. inpol: 0.028562s.
For 1000 test points: intri: 0.071619s. inpol: 0.127528s.
For 5000 test points: intri: 0.351135s. inpol: 0.912710s.
For 10000 test points: intri: 0.704125s. inpol: 2.767123s.
For 100000 test points: intri: 7.039906s. inpol: 209.625841s.

Could it be that there is some cuda acceleration which could result in such differences in speed judgement?

You are right, that intriangulation gives an uncorrect result for your second example. I will think about that. For the meantime i hope that for real meshes, this exact case will not occur.

Comment only
23 Sep 2013 Sven

Sven (view profile)

@Johannes: glad you got the speed increase from that last tweak.

As for speed comparison, you're right that it depends on problem type but I'm afraid inpolyhedron will still be faster for all but the simplest of cases:

load tetmesh, vertices = X; faces = freeBoundary(TriRep(tet,vertices));
for n = [1 10 100 1000 5000 10000 100000]
fprintf('\nFor %d test points:\n',n);
testp = bsxfun(@plus, 45*rand(n,3), [ 0 0 40]);
tic, int = intriangulation(vertices,faces,testp); fprintf('intri: %fs.\n',toc)
tic, inp = inpolyhedron(faces,vertices,testp); fprintf('inpol: %fs.\n',toc)

intriangulation is indeed faster when the total computation time is tiny, but inpolyhedron takes over as more points are queried or the number of faces increases. For 2000 faces and 5000 points inpolyhedron is around 2x faster and this scales with complexity. The reason for this is that inpolyhedron runs with a divide-and-conquer algorithm that takes some small cost to set up, but recoups that cost by reducing the number of triangles to compare to each query point.

And I'm afraid that there's a class of problem where the crossings method doesn't give the right answer. Consider the simple test:

>> BW=true(5,5,5);BW(2:4,2:4,2:4)=false;fv=isosurface(BW)
>> intriangulation(fv.vertices,fv.faces,[3 3 3])
ans =
>> inpolyhedron(fv.faces,fv.vertices,[3 3 3])
ans =

intriangulation fails here because the query point aligns with the facet edge. In this case, counting the crossings doesn't tell you if a point is inside because it's possible for ANY number of triangles to meet at a given point. Perhaps you can start taking the unique crossings, but that tends to slow things down quite a bit. I actually originally wrote inpolyhedron using this triangle-crossings approach, but changed to the current triangle-distance method for this exact reason.

I think you're on the right track and probably the fastest implementation would have a combination of my grid division and your triangle-crossing once this issue is ironed out, so I still encourage you to dig into intriangulation to see if there's a fast solution to be found!

Comment only
18 Sep 2013 Johannes Korsawe

Thank you for your suggestion, i updated the file based on that idea.
The nice thing about intriangulation is, that you need not care about face normals or their in/outward direction. And i think intriangulation is - at least after the last improvement - now faster than inpolyhedron, but this may depend on the type of problem.

Comment only
06 Sep 2013 Sven

Sven (view profile)

Hi Johannes, nice submission. It looks like your work looks at arbitrarily placed points Adam's was restricted to grids.

In terms of speed, it seems the biggest culprit is the call to unique() on line 220. You might be able to get a little boost in speed if you can somehow reduce the number of points being sent to each successive ray direction (x, y, z). I think (although I'm not sure) that after checking the x direction, there will be *some* points which you are certain aren't in the volume, so you can avoid asking about them in the y and z directions.

Anyway, also check out inpolyhedron. It performs a similar task with a different approach, and I've tried lots of different ways to shave off calculation time and memory footprint... you might be able to adopt some of them.

Comment only
18 Sep 2013 1.1

Following the idea of Sven's hint. Testing successive directions only if former tests left undeterminable results.

24 Sep 2013 1.2

Algorithm update: No more fails if testpoint is projected on edge of mesh.

25 Sep 2013 1.3

Numerous improvements in speed using Matlab Profiler.

06 Jan 2014 1.4

result=1 for testpoints==vertices

Contact us