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.)
See
help intriangulation
for an example and details.
In contrast to the FEXsubmission inpolyhedron, you do not need to compute and/or preprocess any face normals in intriangulation.
The mesh may not have internal faces!
ACKNOWLEDGMENT:
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.
Johannes Korsawe (2020). intriangulation(vertices,faces,testp,heavytest) (https://www.mathworks.com/matlabcentral/fileexchange/43381intriangulationverticesfacestestpheavytest), MATLAB Central File Exchange. Retrieved .
1.5.0.0  Added function drehmatrix, which is needed for optional parameter heavytest.


1.4.0.0  result=1 for testpoints==vertices 

1.3.0.0  Numerous improvements in speed using Matlab Profiler. 

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

1.1.0.0  Following the idea of Sven's hint. Testing successive directions only if former tests left undeterminable results. 
Inspired by: Mesh voxelisation
Create scripts with code, output, and formatted text in a single executable document.
Jurjen Leer (view profile)
Ioannis Nemparis (view profile)
Hello Johannes,
Is there also a way to get in which triangle the actual testpoint is on, instead of the simple the boolean yes or no value?
Dampf Hans (view profile)
Rik (view profile)
I had loads of points with unexpected results (about 5%), until I thought of using heavytest. This reduced the number of unexpected outcomes to nearly 0.
Your code was a great help in my research, thanks a lot!
Sascha Duczek (view profile)
Sascha Duczek (view profile)
Yanrong Chen (view profile)
Johannes Korsawe (view profile)
Fabian,
intriangulation includes checks for rays crossing triangles. The big question is: If a ray crosses the edge of a triangle, should this be considered a "hit" for this triangle or not?
I can construct examples and counter examples for both assumed cases, where intriangulation will fail or succeed in determining inside/outside.
And your example just resembles one of these "fail"examples.
If you still want to use this submission, you could alter the way, the optional input heavytest works. By now, i rotate the problem heavytesttimes and assume a testpoint to be inside just if it is inside everytime. (see lines 117 to 123). One could alter this to:
 considered inside, if inside at least once or
 considered inside, if inside the most times.
I think, you will be able to apply this change.
Regards,
Johannes
Fabian Wolf (view profile)
Hello Johannes,
first of all, thank you for providing us with this tool!
I have a problem, which I hope you can help me with.
I have a cube consisting of the nodes:
x = [0 10 10 0 0 10 10 0 5];
y = [0 0 0 0 10 10 10 10 5];
z = [0 0 10 10 0 0 10 10 5];
vertices(:,:)=[x;y;z]';
I create a polyhedron with it by:
tetra = delaunayn(vertices);
faces = freeBoundary(TriRep(tetra,vertices));
The points that I want to check are basically smaller cubes inside the bigger cube:
points=[0 0 0; 10 0 0; 10 0 10; 0 0 10; 0 10 0; 10 10 0; 10 10 10; 0 10 10;...
5 0 0; 10 0 5; 5 0 10; 0 0 5; 5 0 5; 10 5 0; 10 10 5; 10 5 10;...
10 5 5; 5 10 10; 0 10 5; 5 10 0; 5 10 5; 0 5 10; 0 5 0; 0 5 5;...
5 5 5; 5 5 0; 5 5 10];
The problem is that inpolyhedron.m doesn't reliably recognize the points on the surface of the big cube.
Do you know a way to make it work?
Your help is much appreciated!
Regards,
Fabian
Johannes Korsawe (view profile)
Hi Jiahuan,
i have no reference for this. Please refer to the ideas in “Mesh Voxelisation”. But there is also no reference, so you will have to take it from the code.
Regards, Johannes
jiahuan cui (view profile)
Hi Johannes,
Thanks for sharing the code.
I need to use this function but in Fortran language. I am thinking to implement this method myself.
Do you have any reference for this code please?
Best wishes
Jiahuan
James Johnson (view profile)
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.
James
Daniel Leib (view profile)
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!
John (view profile)
Johannes Korsawe (view profile)
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,
Johannes
Rene Labounek (view profile)
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.
Best,
Rene
[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];
h=trisurf(CYL.faces,CYL.vertices(:,1),CYL.vertices(:,2),CYL.vertices(:,3));
set(h,'FaceColor','black','FaceAlpha',1/4,'EdgeColor','none');
hold on
plot3(testp(:,1),testp(:,2),testp(:,3),'b.');
hold off
in = intriangulation(CYL.vertices,CYL.faces,testp)
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.
Johannes Korsawe (view profile)
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.
yusuf (view profile)
Hi,
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.
Johannes Korsawe (view profile)
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 (onetime) effort to do at least with my realworld examples.
Johannes Korsawe (view profile)
Sven,
your numbers are right. I did not have the latest version of your submission.
So i still have to accelerate intriangulation.... :)
Sven (view profile)
Hi Johannes, are you using latest versions of each? After adding voxeliseability 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)
end
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.
Johannes Korsawe (view profile)
@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.
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)
end
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 divideandconquer 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 =
0
>> inpolyhedron(fv.faces,fv.vertices,[3 3 3])
ans =
1
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 trianglecrossings approach, but changed to the current triangledistance 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 trianglecrossing 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!
Johannes Korsawe (view profile)
@Sven:
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.
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.