Inhull

Efficient test for points inside a convex hull in n dimensions

Description

Testing if a point is inside a convex hull can be done in a variety of ways. Inhull converts the problem into a dot product. If not supplied, it also computes the convex hull too. Inhull also attempts to avoid memory problems, doing the computation in smaller blocks when appropriate.

Here is a comparison of inhull to tsearchn:
n = 500;
m = 100;
p = 5;
xyz = rand(m,p);
testpts = rand(n,p)-.1;

tic
tess = delaunayn(xyz);
in0 = ~isnan(tsearchn(xyz,tess,testpts));
toc
tic
in1 = inhull(testpts,xyz);
toc

tsearchn: Elapsed time is 0.813646 seconds.
inhull: Elapsed time is 0.242993 seconds.

Acknowledgements

This file inspired Inpolyhedron Are Points Inside A Volume?.

• MATLAB release: MATLAB 7.0.1 (R14SP1)

[Swagatika]

 

[Vincent Jaouen]

This code is much faster than the one I used to use. Thank you very much.

[Thomas]

 

[Sven]

Great submission - very fast. Note that it could be even faster: the calls to repmat(false, X,Y) on lines 109,113,155 could be replaced with false(X,Y).

[Sven]

 

[Richard Crozier]

Another great tool, I'm using it to select arbitrary regions of a finite element mesh.

[Duy]

This is a great tool, but the advantage of tsearchn is that tsearchn also returns the index of the tetrahedron (for 3d case) consisting of the test point. So, if i want that index, i need to run a loop over all tetrahedron and check everytime if the testpoint is inside or not, which would raise the runtime of inhull in an unacceptable way (number of tetrahdrons * inhull). The speed of a single call of inhull with the result of tsearchn would be awesome.

[John D'Errico]

Sorry. I'll fix it. the example should have read

xy = randn(20,2)
tess = convhulln(xy);
testpoints = [ 0 0; 10 10];
in = inhull(testpoints,xy,tess)

in =
1
0

[Roman Voronov]

there is a problem with the example provided in the code:

% xy = randn(20,2)
% tess = convhull(xy(:,1),xy(:,2));
% testpoints = [ 0 0; 10 10];
% in = inhull(testpts,xyz,tess)

[vicky]

Thank you John for you valuable suggestions

[John D'Errico]

Hi Vicky,

So cross gives you a facet normal, but you are unsure if it is outward or inward facing. First of all, you can get all of the normals in one call to cross. Get used to doing these computations in one operation, in a vectorized way. I'll build a test example.

X = rand(10,3);
K = convhulln( X , { 'Qt' } );

So K is the list of triangular facets. Build the normal vectors with cross. All we need do is take the cross product of vectors along a pair of edges. I'll pick two edges, in an arbitrary orientation.

normals = cross(X(K(:,1),:) - X(K(:,2),:),X(K(:,1),:) - X(K(:,3),:));

We could normalize them to have unit length if we wish. I tend to do so for normal vectors.

normals = bsxfun(@times,normals,1./sqrt(sum(normals.^2,2)));

Are they inward or outward pointing? A simple dot product will tell you the answer.

C = mean(X,1);
D = dot(normals,bsxfun(@minus,X(K(:,1),:),C),2)
D =
-0.319410966973992
-0.289171257799806
-0.192437832564587
-0.220416111145392
-0.322367706246851
-0.279897929334218
-0.23737559253969
-0.230111277638159
-0.319713264137273
-0.285713570370542
-0.241935614441904
-0.293453534156934
-0.239664243215607
-0.269863191655601

(I could have done the dot product in other ways too, but this way is clear about what I am doing.)

See that all of these dot products are negative. What does it tell you, if the projection of one vector on another is negative? The answer is the two vectors point in opposite directions. See that I've chosen C in a way such that it MUST be inside the convex hull., then I subtract it from some point on the same facet as each normal.

The negative values indicate every normal vector in this set is inward pointing. If you wanted outward pointing normals, simply negate the normals. If not every normal was consistently pointed in the correct way, then you could as easily correct that.

Regards,
John

[vicky]

Hello John D'Errico,
I read your Inhull matlab file on File Exchange of MathWorks.com. I need a small help from you.

K = convhulln( X , { 'Qt' } );

each row of K contains the indices to the points defined in X, that makeup the triangular facet of the convex hull.
I want to find the outward normals to each of triangular facets in a direction pointing outwards from the convex hull. But I am stuck on how to determine the direction of the normal, whether it is outward or inward.

I think the indices in K (to the points of X) are stored in some meaningful manner(probably in counter clockwise). Here this is how am evaluating the outward normals for each triangular facets(not sure if this is correct way).

%outward normal to facet 1:
temp1 = X( K(1, 2) , : ) - X( K(1, 1) , : );
temp2 = X( K(1, 3) , : ) - X( K(1, 1) , : );
normal = cross( temp1, temp2 );

I need your valuable suggestions for this. I felt in your work, you have encountered this situation. Please guide me. I will be very thankfull for your help.

[Antonio Pena]

Sorry, thanks anyway for your help, i'll try in the forum

[John D'Errico]

What is your question here? Regardless, it is inappropriate for the comments on this file to become an extended set of questions by you on a variety of problems that do not really relate to this code anyway. That is not what these comments are for. I won't do consulting through the comments on this file.

[Antonio Pena]

But i did this for 2D with inpolygon, what I want to get with this are the points inside so as then I can refine the mesh in the internal structure:
x_min=0; x_max = 0.1;
y_min=0; y_max = 0.1;
Nx = 20;
Ny = 20;
x=linspace(x_min,x_max,Nx);
dx = (x_max-x_min)/Nx;
y=linspace(y_min,y_max,Ny);

points = [1 1 x_min 0.158; 2 1 0.053 0.211; 3 1 x_min 0.316; 4 1 0.211 0.526;5 1 0.579 1;6 1 1 0.895;7 1 0.842 0.474;
8 1 0.368 0.158;9 1 0.158 y_min;10 1 0.105 0.053;11 1 0.053 0.105;12 1 x_min 0.158; 13 2 0.368 0.474; 14 2 0.421 0.579;...
15 2 0.474 0.684; 16 2 0.579 0.737; 17 2 0.684 0.632;18 2 0.684 0.526; 19 2 0.579 0.421; 20 2 0.526 0.316;...
21 2 0.474 0.263; 22 2 0.421 0.263; 23 2 0.368 0.368; 24 2 0.368 0.474];

points(:,3) = points(:,3)*x_max;
points(:,4) = points(:,4)*y_max;
xv_medio1=[];
yv_medio1=[];
xv_medio2=[];
yv_medio2=[];

for i= 1:length(points(:,2))
if points(i,2)==1 %MEDIO 1
xv_medio1 = [xv_medio1 points(i,3)];
yv_medio1 = [yv_medio1 points(i,4)];
elseif puntos(i,2)==2 %MEDIO 2
xv_medio2 = [xv_medio2 points(i,3)];
yv_medio2 = [yv_medio2 points(i,4)];
end
end

[X,Y] = meshgrid(x,y);
in_medio1 = inpolygon(X,Y,xv_medio1,yv_medio1);
in_medio2 = inpolygon(X,Y,xv_medio2,yv_medio2);
figure
plot(xv_medio1,yv_medio1,xv_medio2,yv_medio2,X(in_medio1&~in_medio2),Y(in_medio1&~in_medio2),'.r',X(~in_medio1),Y(~in_medio1),'.b',X(in_medio2),Y(in_medio2),'.g')
axis([x_min x_max y_min y_max]);

[John D'Errico]

inhull will NOT solve the question of whether a cell of a rectangular grid is entirely inside a hull. For that matter, neither would inpolygon do that in only two dimensions. It is trivial to show a counterexample.

As far as plotting goes, learn to use the convhull routines in matlab. Then call trisurf to plot a triangulated surface.

[Antonio Pena]

By the way, in this function all i need for define the polygon is the set of vertices, but what if i want to plot it ? how can i do it?

[Antonio Pena]

Thanks for your explanation. What i'm trying to do is find a version of inpolygon in 3D as i need it just to find out if each cell of a regular mesh is inside or outside of this polygon

[John D'Errico]

First of all, I think you misunderstand that this code works with a triangulation, or in higher numbers of dimensions, a tessellation of a domain. Further, you seem to misunderstand how to build those things to define a triangulation/tessellation.

The points in your example define a simple cube in three dimensions. The convex hull of those points is a set of 12 triangles. Each triangle will be given as a set of three INTEGERS.

Before we go any further, let me suggest that you learn to use arrays in matlab, rather than defining separate, numbered variables for every single data point. Thus, define the array variable V, as

V = [0.1 0.1 0.1;
0.1 1 0.1;
0.1 0.1 1;
0.1 1 1;
1 0.1 0.1;
1 1 0.1;
1 0.1 1;
1 1 1];

If you do not learn to do this in matlab, you will find your further work is much more difficult (and VERY SLOW) once you start working with larger sets of data.

Each row of V is ONE vertex of the cube.

tri = convhulln(V)
tri =
6 2 1
5 6 1
7 5 1
3 7 1
7 6 5
6 7 8
2 4 1
4 3 1
6 4 2
4 6 8
4 7 3
7 4 8

See that each row of this array has numbers that run from 1 to 8. Integers, that refer to points or vertices, that are themselves references to rows of V. So the first TRIANGLE in the hull is [6 2 1], defined by the points stored as V([6 2 1],:). There will be 12 triangles, since there are 6 square faces of a cube.

If you will choose to define a triangulation, then you can use inhull. Or, you can let inhull create that triangulation itself. But you cannot create a polyhedron as you have tried to do and try to pass in facets of the cube as 4 sided things, defined in your own arbitrary notation, and hope that inhull will figure out what you mean to do. So, the proper way to call inhull, as already shown in the help, is to do ONE of these things...

% Here, inhull creates the surface triangulation itself.
inhull(testpoints,V)
ans =
1
0

% Here I have used the output of convhulln above.
inhull(testpoints,V,tri)
ans =
1
0

See that either call tells you that the first point is inside the convex hull, and the second is outside the hull.

[Antonio Pena]

could anyone give an example of using this function with a 3d polygon? I've tried with this:
%Test in 3D
clear all,close all,clc;
v1 = [0.1 0.1 0.1];
v2 = [0.1 1 0.1];
v3 = [0.1 0.1 1];
v4 = [0.1 1 1];
v5 = [1 0.1 0.1];
v6 = [1 1 0.1];
v7 = [1 0.1 1];
v8 = [1 1 1];
face1 = [v1 ;v5 ;v6 ;v2 ;v1];
face2 = [v3; v7; v8; v4; v3];
face3 = [v5; v6; v8; v7; v5];
face4 = [v6; v2; v4; v8; v6];
face5 = [v2; v1; v3; v4; v2];
face6 = [v1; v3; v7; v5; v1];
polygon = [face1;face6; face5; face4; face3;face2];
testpoints = [ 0.5 0.2 0.2; 2 2.3 0.4];
x = testpoints(:,1);
y = testpoints(:,2);
z = testpoints(:,3);
in = inhull(testpoints,polygon);
figure
plot3(polygon(:,1),polygon(:,2),polygon(:,3),x(in),y(in),z(in),'.g',x(~in),y(~in),z(~in),'.r')
but the first point in testpoints should be inside the cube but it isn't...

[Pierre Bellec]

Great code. I've used it for head segmentation in MRI. A simple intensity thresholding left lots of holes in the head, and that function filled them easily. It was a problem in low dimension (3) but quite high resolution (256*256*128). In any case, it saved me a lot of time and I really appreciate it.

[David Gingras]

Another very usefull code provided by John! Thanks.

[Navid Samavati]

Excellent Code!

[Florian]

Thank you very much for this code. Particularly thanks to John for explaining me a problem I had with the code:

I have tried the following:
tms=rand(1000,5);
in = inhull(tms,tms);
disp(length(find(in==0)))

which gave roughly ~319

as John points out:
"those 319 points were almost certainly on the surface of the convex hull.

Try this slight modification of your code instead.

tms=rand(1000,5);
in = inhull(tms,tms,[],1.e-13);
disp(length(find(in==0)))
0

See that I've used a tolerance this time. The linear algebra used in the test is done in floating point arithmetic in matlab, as it must be. Without a tolerance, any test done in floating point arithmetic is suspect. By default that tolerance is zero, but this is also why I allowed the user a tolerance at all. "

Thanks John!!!!!

[Nagarjuna]

very nice work

[Luigi Giaccari]

Very good code, readable, uderstandable and even fast.
The n-dimensions also make-it flexible.

The thing I like the most is a wonderfull example of dot product vectorization with memory management. It helps a lot to understdand how to make m-code more efficient.

The only thing i disagree is the comparison with tsearchn.
I thing its purpose is to find where a point is located inside the hull and not whether he is inside or not. The nan output values are just a secondary objective which occurs in the worst case when a point lies outside.

Anyway a five stars, very good work!

[Luigi Giaccari]

Very good code, readable, uderstandable and even fast.
The n-dimensions also make-it flexible.

The thing I like the most is a wonderfull example of dot product vectorization with memory management. It helps a lot to understdand how to make m-code more efficient.

The only thing i disagree is the comparison with tsearchn, I thing its purpose is to find where a point is located inside the hull and whether he is inside or not. The nan output values are just a secondary objective which occurs in the worst case when a point lies outside.

Anyway a five stars, very good work!

[Titus Edelhofer]

Very nice. Just for ease of use: the example given in the help uses different variable names then the call to the function (and should call convhulln instead of convhull).

[Jayveer Thakoor]

huge time saver! thank you :)

[Lee Shunn]

Fast and easy to use. Thanks for the great code.

[Yonas T. Weldeselassie]

Excellent! Great job. Well documented, easy to understand, very similar to inpolygon and what more? Thanks a lot.

[Jon KoS]

Great Job! Works really well and does what it says. I'll be using this function a lot. :) Not sure what the other reviewer was looking at but the code is well commented. Thanks!

[John Galt]

Ever heard of commenting your code?

[Quan Wen]

I agree with what Mia Ginsburg said. This code is very helpful. Thanks a lot.

[Mia Ginsburg]

Just what I wanted. Much faster than tsearchn if you just want to test if a point is inside a convex hull.

[Liketo Stayanonymous]

No, sorry! I think you lost me with this one. What (again) is this software about??
Clear example will be much appreciated!