Code covered by the BSD License  

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

5.0

5.0 | 1 rating Rate this file 26 Downloads (last 30 days) File Size: 6.11 KB File ID: #43381
image thumbnail

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

by

 

05 Sep 2013 (Updated )

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

| Watch this File

File Information
Description

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 FEX-submission 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.

Acknowledgements

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 (12)
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,
Johannes

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.

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)

17 Jun 2014 yusuf

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.

30 Dec 2013 yusuf

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.

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.

25 Sep 2013 Johannes Korsawe

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

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

24 Sep 2013 Sven

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)
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.

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.

23 Sep 2013 Sven

@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 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 =
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 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!

18 Sep 2013 Johannes Korsawe

@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.

06 Sep 2013 Sven

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.

Updates
18 Sep 2013

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

24 Sep 2013

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

25 Sep 2013

Numerous improvements in speed using Matlab Profiler.

06 Jan 2014

result=1 for testpoints==vertices

Contact us