version 3.3.0.0 (22.4 KB) by
Sven

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

**Editor's Note:** This file was selected as MATLAB Central Pick of the Week

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

Sven (2020). inpolyhedron - are points inside a triangulated volume? (https://www.mathworks.com/matlabcentral/fileexchange/37856-inpolyhedron-are-points-inside-a-triangulated-volume), MATLAB Central File Exchange. Retrieved .

Created with
R2012a

Compatible with any release

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Hansong Jiwonderful, thank you Sven.

lunalaraThis is a great function. Thank you immensely for sharing this!

Paul ZhangWorks great! 16x faster than in_polyhedron.

Ilhan BokThis works great. Not sure if posted already, but results from the STL file reader (https://www.mathworks.com/matlabcentral/fileexchange/22409-stl-file-reader) can be fed into this script, rather than using the PDE toolbox.

shmulik eIs it possible to make a modification so it will run on the GPU faster ?

Hamidreza NourzadehHi Sven, Thank you for sharing this code. I was wondering if you could try the following failing examples (some mat files) at https://www.dropbox.com/s/xzb1kcd8cl7z7rp/failedTests.rar?dl=0.

%Test code + visualization

clear variables

clf

load ex5Failed % load on of the example ex2Failed.mat to ex5Failed.mat

TR = triangulation(faces,vertices);

in2= find(inpolyhedron(faces,vertices,unique(points(:,1)),unique(points(:,2)),unique(points(:,3))));

hf= trisurf(TR, 'FaceColor',[0 0 0.8],'FaceAlpha',0.5,'EdgeColor','none');

axis equal, drawnow;

lighting phong;

camlight infinite;

camproj('perspective');

hold on

hScatter2= scatter3(points( in2,1), points( in2,2), points( in2,3),8, 'b', 'fill');

Runze HanHi Sven, thanking so much for the code and it worked perfectly. Just want to know if there is any reference paper this is based on, and how I can cite you?

Daniel LamMaraAlbrecht HaaseGreat tool, perfectly fitting my needs.

Martin TschierskyLine 191: replacing 'clear' with 'clearvars' yields a considerable performance boost

LMSI am trying to use this with a shape model of comet 67P:

>> load('cg-spc-shap5-v1.5-cheops.mat')

>> shp = patch('Faces',redf,'Vertices',redv,'FaceVertexC',[1 1 1],'FaceC','f','LineS','n');

and a set of point-coordinates

>> test = [1.8955 -0.2307 1.2290];

But I get this error:

>> IN = inpolyhedron(shp,test)

Function 'subsindex' is not defined for values of class

'matlab.graphics.primitive.Patch'.

Error in inpolyhedron>parseInputs (line 425)

facets = permute(reshape(facets(:,faces'), 3, 3, []),[2 1 3]);

Error in inpolyhedron (line 105)

[facets, qPts, options] = parseInputs(varargin{:});

....what am I doing wrong, please?

Pradeep ReddyThank you for the code. It works perfectly. I need the principle behind the code. Please can you give me the link/literature/paper based on which it is done? I will be greatfull if you provide me.

Nathaniel H WernerHi again Sven,

This might help. I've tried using

[faces,vertices] = isosurface(X, Y, Z, LEV,-1);

X, Y, Z, and LEV are data that are loaded in

in_2 = inpolyhedron(faces,vertices,TestPoints,'FlipNormals',true);

Points_inside = [TestPoints(in_2,1),TestPoints(in_2,2),TestPoints(in_2,3)];

in_3 = inpolyhedron(faces,vertices,TestPoints);

Points_outside = [TestPoints(in_3,1),TestPoints(in_3,2),TestPoints(in_3,3)];

[intr,i_IN,i_OUT] = intersect(Points_inside,Points_outside,'rows');

Unfortunately, the "extra" points found and saved as in_2, have no intersection with those in in_3. The "extra" points are simply a few scattered lines outside the isosurface.

Sunita SahaError in inpolyhedron (line 185)

cells{yi,xi} = cat(1,tmpInds{xyMinMask & yi >= unqLHgrids(:,2) & yi <= unqLHgrids(:,4)});

Please help with this error

Jay Willisderu jianproczellThis is a great function which saved me a lot of time. Well documented as well!

Ranjan MelpalHi Sven,

Thanks for the upload! Are there any papers/references on which the toolbox is based? Could you share them?

JiayiAmazing function! Just save my ass! Thanks a lot!

Weizhou LiI am using point as center of element. When there is thin walled structure, this code will determine the element is outside. I tried to use tolerance to make it inside. However, the tolerance is only working on z direction. Could changes be made to let user define tolerance in x, y direction?

king june@James Farrell The reason for the error was the 'verLessThan' funtion ,you should modify your code like this：

function options = parseOptions(varargin)

IP = inputParser;

if verLessThan('matlab', '8.2.0.701') % change 'R2013b' to '8.2.0.701'

fcn = 'addParamValue';

else

fcn = 'addParameter';

end

James FarrellHi, I copied the example into a seperate script but when I run the script I get the following error:

No appropriate method, property, or field addParameter for class inputParser.

Error in inpolyhedron>parseOptions (line 453)

IP.(fcn)('gridsize',[], @(x)isscalar(x) && isnumeric(x))

Error in inpolyhedron>parseInputs (line 439)

options = parseOptions(varargin{:});

Error in inpolyhedron (line 105)

[facets, qPts, options] = parseInputs(varargin{:});

I am using MATLAB 2013a

James FarrellNo appropriate method, property, or field addParameter for class inputParser.

isfzadeNot working for cone type. If it is needed I can send you vertex and faces to examine. irashidzade@gmail.com contact to me

Deepa IyerFritzFabian WolfHello Sven,

thank you so much!

That is exactly what I needed.

Sven@Fabian

Try this, adjusting TOL as needed. It detects point ON a 3d triangulation (within a tolerance).

F = [1 2 4;2 3 4;2 1 6;1 5 6;1 4 8;5 1 8;3 2 7;2 6 7;4 3 7;8 4 7;6 5 7;5 8 7];

V = [0 0 0;10 0 0;10 0 10;0 0 10;0 10 0;10 10 0;10 10 10;0 10 10];

points = rand(200,3) * 12 - 1;

TOL = 0.5;

% Make the original triangulation and show it

T1 = triangulation(F, V);

figure, hold on

patch('Faces',F,'Vertices',T1.Points,'FaceColor','g','FaceAlpha',0.2);

% Offset vertices by tiny amounts each way, make new triangulations

T_pos = triangulation(F, T1.Points + T1.vertexNormal*TOL);

T_neg = triangulation(F, T1.Points - T1.vertexNormal*TOL);

patch('Faces',F,'Vertices',T_pos.Points,'FaceColor','r','FaceAlpha',0.2);

patch('Faces',F,'Vertices',T_neg.Points,'FaceColor','b','FaceAlpha',0.2);

% Use inpolyhedron on the offset triangulations

in_pos = inpolyhedron(F,T_pos.Points, points);

in_neg = inpolyhedron(F,T_neg.Points, points);

% Some logic to get inside, outside, and on the edge

on_vol = in_pos & ~in_neg;

in_vol = in_neg & ~on_vol;

out_vol = ~in_pos;

% Display

doplot3 = @(x,varargin)plot3(x(:,1),x(:,2),x(:,3),varargin{:});

inH = doplot3(points(in_vol,:),'ko','MarkerFaceColor','b');

outH = doplot3(points(out_vol,:),'ko','MarkerFaceColor','r');

onH = doplot3(points(on_vol,:),'ko','MarkerFaceColor','g');

legend([inH outH onH],{'IN','OUT','ON'})

Fabian WolfI solved the issue with the nodes above the polyhedron being recognized (user error...)

My other problem changed, as well. The nodes on the edges of the polyhedron aren't recognized that well. I attached the figure to show you (I hope it's allowed to post links here):

https://dl.dropboxusercontent.com/u/12977637/in_polyhedron.jpg

'TOL' didn't help.

Fabian WolfHello everybody!

First of all, thank you Sven for putting so much work in this tool and your work that you put into MatLab. It’s much appreciated.

I have a few difficulties getting inpolyhedron to work. But first, let me explain my situation.

I have a test mesh in the form of a cuboid consisting of around 20000 nodes. This mesh is made of hexa-elements of 8 nodes each.

The first step that I want to do, won’t make a lot of sense but it’s just for testing purpose.

I take the nodes of 1 element and create a polyhedron by using delaunayn and freeboundary, giving me a vertices and faces matrix.

Then, I take the 20000 nodes and check if the nodes are inside that 1 element. Only 8 nodes should be found, the ones creating the 1 element, right?

However, I get a lot more. I get all the nodes that are in the same x- and y- locations as the 8 nodes of the element but which have different z-locations. Can you tell me what I am doing wrong?

My next question is, how can I ignore the elements that are on the surface of that 1 element or rather in the same locations as the 8 nodes. I couldn’t get it to work with vertexnormal() as you suggested to Gabrielle 2 years ago…

I hope you guys can help me! I am rather new to this subject so I feel a little bit lost.

Thanks in advance!

- Fabian

BinuFab BHi Sven,

I can now confirm that the qPtsXY(':',':') bugfix works for me on R2012b and R2015b.

Thanks a lot,

Fab

SvenAh, yes, Tim and Fab - you are right. I had both versions on my machine and didn't realise the link was to 2015a so that's what I was actually running :)

I'll upload a new version with quotes around the semicolons to avoid that error.

Tim LuethqPtsXY(':',':') solved at least my problem on R2015b

Tim LuethHi Sven, we also use R2015b (academic) and have the same problem on PC and MAC. R2015a was working as expected. Tim

Fab BHi Sven,

I'm think I'm using the full release of R2015b. The exact version is 8.6.0.267246 (R2015b).

I've already tried the bugfix you propose ("qPtsXY(':',':')") before and it seems that it did not work properly: in some cases it did in others not (unfortunately I don't have access to another version such that I could compare properly).

Thanks,

Fab

SvenHi Fab, I think there are 2 solutions:

Firstly, are you using the 2015b "Pre-release" or the actual 2015b full release?

I suspect that this issue is with the 2015b pre-release (I saw it too) but I don't actually see it now in the official 2015b so I think that this is where the problem lies.

If you need to continue using the release you have, then I *think* that you can just change instances of "qPtsXY(:,:)" to "qPtsXY(':',':')" and it will work (sorry, I can't verify this for you since I uninstalled the pre-release and can't replicate the problem at the moment).

Thanks,

Sven.

Fab BHey Sven,

Thanks a lot for this great function, I've been using it for a while!

However, after having changed to Matlab R2015b the function does fail with the following error:

"Input arguments to function include colon operator. To input the colon

character, use ':' instead.

Error in inpolyhedron (line 233)

qPtGridXY = floor(bsxfun(@rdivide, bsxfun(@minus, qPtsXY(:,:),

scalingOffsetsXY),..."

Any suggestions?

Thanks,

Fab

JiangpingJiangpingAriel Luthanks Sven, it's perfect to work.

Sven@James: Sure, there should be a "contact" button if you click on my name>profile. A short description (sounds like you're testing inpolyhedron against a dense point cloud) and your email and I'll see what can be done.

James JohnsonHey Sven I have some issues using your program and I was wondering if you can assist me. I tried to email you but it couldn't go through.

Carla Kohoyda-InglisFritzwariasrKHALDON HMESHEHThanks Sven. I will try what you mentioned. I am waiting for a new solution.

Khaldon

SvenHi Khaldon, you've hit an issue that I'm looking to fix. It occurs as you said when query XY points are *exactly* the same as vertex XY points, and those vertices connect faces which "overhang" at a high angle. In those cases machine precision issues (from calculating the distance to an angled surface) can mean that the wrong face is detected as "nearest" and points above/below that face are wrongly set in or out.

I have a slow solution (involving re-testing points detected as coincident with vertices) and I'm looking for a fast(er) solution. In the meantime, I think the best solution is to simply offset your query points by a very tiny amount. I'll email you when I get something.

Thanks,

Sven.

KHALDON HMESHEHHello Sven,

very nice job. I tried your code with an STL file, where the 3D object vertices are also the query points. The result was that some of the points are regarded as OUT. How to resolve this bug or is there something to modify.

Thanks again

Khaldon

Marc LalancetteAmazing 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...

Peter MorovicThanks for sharing, Sven - works great indeed.

yusufThank you Sven.

This code works perfectly

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

RheedaThis works brilliantly! Thank you!!!

GabrieleJohn, 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!

SvenHi 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

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

John D'ErricoBy 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.

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

John D'ErricoGabriele - 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.

GabrieleHi 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!

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.

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

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.

Yan OuI 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;

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.

Dun KirkMy 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.

Benjamin FriedrichExactly what I needed. Worked out of the box.

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

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

APIt worked. Super fast. Great job. Thanks.

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.

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.

API 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(:)];

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.

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.

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

Shawn WalkerAwesome 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.

Shawn Walker