Efficient test for points inside a convex hull in n dimensions
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.
1.2  minor changes for a tiny speed boost 

1.1  Repaired example in the help 

Fix a problem in 2d or 3d, catching degenerate facets 

Version 2.0: Repaired a problem when the hull has degenerate facets. 

Added a tolerance on the tests 
Inspired: 3D mesh transform using sparse control points, The Barycentric FixedMass method for estimating fractal dimensions, STORMbased Relative Localization Analysis
Matt J (view profile)
His Teknoloji (view profile)
Instead of using writing own algorithm, one can use inpolygon builtin script with convhull combined. For example:
k = convhull(x,y);
% Here k stores proper indices for convex hull vertices.
% Then we can draw polygon and test if any point inside or
% convex hull
j = 1
for i=1:max(size(k))
xx(j) = x(k(i))
yy(j) = y(k(i))
j = j+1
end
[in on] = inpolygon(xtest,ytest,xx,yy);
numel(xx(in)) % or numel(yy(in))
dhwani contractor (view profile)
I am using the same code but getting this error:
Undefined function or variable 'inhull'.
P.S I am using MATLAB 2016.
Any help on this is appreciated!
Thank u in advance! :)
nilay kant (view profile)
Hii, I am using the inhull function. Is there a way to find if the points lies explicitly inside the convex hull? In other words i want to classify if all the points which lies inside but not on the surface of the convex hull? Currently this code outputs 1 if the points lies on the surface.
Thank you
Bob Roginski (view profile)
Perhaps something I'm not understanding (noob to computational geometry) . . . using the example in the description, if I calculate:
in0 = ~isnan(tsearchn(xyz, tess, xyz));
I obtain all "true" values for each row
If I calculate
in1 = inhull(xyz, xyz);
a substantial portion of the values are "false".
Consider me puzzled . . .
Seth (view profile)
Has anyone been able to create a MEX file or use Matlab Coder with this function? I think it's close to OK, but getting a little hung up on some build errors. If you know how, could you please post a link or instructions?
Thanks!
Seth (view profile)
John D'Errico (view profile)
Oh, gosh, why do you want to work in 6 dimensions at all? :) Just kidding.
Note that things get large quickly enough in higher dimensions. 6 dimensions is starting to get moderately large when you talk about convex hull and triangulations, still very doable though.
X = randn(10,6);
tri = convhulln(X);size(tri)
ans =
40 6
X = randn(40,6);
tri = convhulln(X);size(tri)
ans =
866 6
X = randn(100,6);
tri = convhulln(X);size(tri)
ans =
2871 6
So it is not obscene in only 6 dimensions. And inhull works with no problem.
inhull([.6 .6 .6 .6 .6 .6],X,tri)
ans =
1
sravan venkata (view profile)
Hey John, Can this be used to test for a convex hull of 6 Dimension?
Can you also suggest a function used to compute the distance between origin and facets of Convex Hull.
M. F. (view profile)
Very fast, works like a charm!
Jan Froehlich (view profile)
Very fast, thank you for submitting this.
John D'Errico (view profile)
Andy  there are a few ways to solve your problem.
Once an object is no longer convex, the answer is not as easy. Inhull uses a simple and efficient trick, converting all facets of a convex surface into linear inequalities, so the test becomes not much more than a simple (large & sparse) matrix multiply. However, the linear inequality approach fails when the object is no longer convex, and there is no such thing as only slightly nonconvex. So, some simple alternatives are:
 Form a tessellation of the nonconvex domain, converting it entirely into a list of simplexes that fill the region. A point is inside the surface if it is inside ANY of the component simplexes. This is doable, although a bit slower than inhull would be. The biggest issue is how one converts a surface into a full tessellation.
 Another scheme is to use a ray. A point is inside the domain if any random ray extending from that point to infinity passes through an odd number of facets of the surface. Clearly you need to be careful if the ray happens to be parallel (and coincident with) a facet or edge, or it passes through a vertex. As well, one needs to write careful code to determine what happened if the point fell exactly on a facet, edge, or vertex.
Andy Horn (view profile)
Hello everybody!
I love inhull, it really helps a lot! However, I also have some objects that are not necessarily convex and would need to have a generalized version of this, which could e.g. read the outputs of a isosurface function (triangles and points).. Is anyone aware of a similar function to inhull which does exactly this? (testing if points reside within a surface described by points and triangles)?
Thanks a lot for your help!
Yours, Andy
Francesco (view profile)
Very fast and easy to use! I like how it does the tesselation of surfaces for you #lazy
QIANG LIU (view profile)
John D'Errico (view profile)
Daniella  Sorry, but I really try not to do consulting in the comments, so I won't go any further with this. But ask yourself what you would expect to happen if the three points fell exactly in a straight line?
Daniela (view profile)
Hi
First of all thanks for your reply, it was very helpful. Till now the program was running smoothly but now I tried to find the normal of the plane containing the points(rows of X) X=[134 294 59; 136 295 59; 122 288 59]. The program is calculating the normal but it is giving out a 3x2 matrix rather than a 3x1 matrix causing the program to generate an error later on when projecting the points in Y to the 2D subplot (line of code: YPProj = YP(inplane,:)*B ).
Can you please guide me in what is causing this to happen?
Thanks for your great help once again.
John D'Errico (view profile)
I KNEW you were going to ask this question, and I'm quite sure I'm not even psychic. :)
Consider three points in R^3, here as rows of X. Clearly they form a triangle. Do the points as rows of Y lie inside that triangle?
X = [1 0 0;0 1 0;0 0 1];
Y = [1 1 1;.25 .25 1;.2 .3 .5;4/7 2/7 1/7]
Are the points in Y in the plane of the triangle at all? A plane is defined by a point on the plane, and the normal vector to that plane. We can also define the plane in terms of a point on the plane, and a pair of basis vectors spanning the plane. In fact, we will use both of these forms.
First, choose a point (P) on the plane. Here P=X(1,:) will obviously suffice. Translate the problem so that the point on the plane is the origin. Essentially this means subtracting P from the rows of both X and of Y.
P = X(1,:);
XP = bsxfun(@minus,X,P)
YP = bsxfun(@minus,Y,P)
XP =
0 0 0
1 1 0
1 0 1
YP =
0 1 1
0.75 0.25 1
0.8 0.3 0.5
0.42857 0.28571 0.14286
Now, compute the normal vector to the plane of the triangle. Since this is a three dimensional problem, I COULD do it using the cross product tool found in MATLAB (as cross), but perhaps a more general solution comes from the function null.
N = cross(XP(2,:),XP(3,:))
N =
1 1 1
N = null(XP)
N =
0.5774
0.5774
0.5774
In fact, both vectors are essentially the same, although one is transposed compared to the other, and the latter is normalized to have unit norm. Both aspects are irrelevant here. I'll use the result from null to do my tests.
The test for being inplane is now trivial. If the result of YP*N is zero (to within floating point trash) then the corresponding point is in the plane of the triangle. BE VERY CAREFUL HERE, as unless you take care in the test for zero, it will fail. You need to use an intelligent tolerance.
YP*N
ans =
1.1547
0
0
4.1633e17
OK, so points 2 and 3 are in the plane of the triangle. In fact, point 4 is also in the plane of the triangle, but I chose it so that you would see that a tolerance is needed here. A reasonable tolerance might be:
tol = 10*eps(sqrt(sum(YP.*YP,2)))
tol =
2.2204e15
2.2204e15
1.1102e15
1.1102e15
That tolerance was chosen because N is a UNIT normal vector, so we need only consider the sizes of the numbers in YP.
inplane = find(abs(YP*N) < tol)
inplane =
2
3
4
So we can see that points [2,3,4] are in the plane of the triangle.
Are those points actually inside the triangle? For this, we need to project the problem into the plane of the triangle. Basis vectors that span the plane of the triangle are given in B:
B = null(N')
B =
0.57735 0.57735
0.78868 0.21132
0.21132 0.78868
XPProj = XP*B
XPProj =
0 0
1.366 0.36603
0.36603 1.366
YPProj = YP(inplane,:)*B
YPProj =
0.024519 1.2745
0.59282 0.79282
0.44258 0.29973
Here I've dropped out point 1, because it was of no interest as it lies out of plane. The rows of XPProj define a triangle in the subspace of our plane. Do the points as rows of YPProj lie in that triangle?
inplane(inhull(YPProj,XPProj))
ans =
3
4
So inhull tells us that points 3 and 4 did lie inside the triangle. In fact, inpolygon tells us the same thing.
inplane(inpolygon(YPProj(:,1),YPProj(:,2),XPProj(:,1),XPProj(:,2)))
ans =
3
4
Fairly long winded, but it comes down to only a few short lines of code.
Daniela (view profile)
Thanks for your fast reply. I got the idea of determining whether the point lies on the plane of the triangle or not but I am a bit confused about converting the problem to 2D. Can you guide me to what is meant by this or else indicate me with some useful website? Thanks for your great help
John D'Errico (view profile)
Daniella  I'm sorry, but 3 points forms a TRIANGLE. A triangle is not a closed surface in 3 dimensions.
If you have a point in 3d and wish to test to see if it actually lies inside a triangle, then you would first need to test to see if it even lies inside the plane of the triangle. If it does lie in that plane, then you would need to convert the problem to a 2d problem by projecting it all into the 2d subspace defined by the plane of the triangle. Once you have done so, then inpolygon would do as well as inhull for a triangle, as long as you have done the projection into 2d.
Daniela (view profile)
I am currently using the Inhull code for my project. I am passing to Inhull my set of testing points together with a 3*3 array. The 3*3 array contains the coordinates of the 3 vertices which are to form the convex hull. Upon running the code I am getting an error saying that I did not provide enough points.
This is the actual error I am getting:
Error using cgprechecks (line 41)
Not enough unique points specified.
Error in convhulln (line 42)
cgprechecks(x, nargin, cg_opt);
Error in inhull (line 81)
tess = convhulln(xyz,[]);
What can I do?
Swagatika (view profile)
Vincent Jaouen (view profile)
This code is much faster than the one I used to use. Thank you very much.
Thomas (view profile)
Sven (view profile)
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 (view profile)
Richard Crozier (view profile)
Another great tool, I'm using it to select arbitrary regions of a finite element mesh.
Duy (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
Thank you John for you valuable suggestions
John D'Errico (view profile)
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 (view profile)
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 (view profile)
Sorry, thanks anyway for your help, i'll try in the forum
John D'Errico (view profile)
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 (view profile)
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_maxx_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 (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
Another very usefull code provided by John! Thanks.
Navid Samavati (view profile)
Excellent Code!
Florian (view profile)
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.e13);
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 (view profile)
very nice work
Luigi Giaccari (view profile)
Very good code, readable, uderstandable and even fast.
The ndimensions also makeit 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 mcode 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 (view profile)
Very good code, readable, uderstandable and even fast.
The ndimensions also makeit 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 mcode 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!
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).
huge time saver! thank you :)
Fast and easy to use. Thanks for the great code.
Excellent! Great job. Well documented, easy to understand, very similar to inpolygon and what more? Thanks a lot.
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!
Ever heard of commenting your code?
I agree with what Mia Ginsburg said. This code is very helpful. Thanks a lot.
Just what I wanted. Much faster than tsearchn if you just want to test if a point is inside a convex hull.
No, sorry! I think you lost me with this one. What (again) is this software about??
Clear example will be much appreciated!