version 1.2.0.0 (3.06 KB) by
John D'Errico

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.

John D'Errico (2020). Inhull (https://www.mathworks.com/matlabcentral/fileexchange/10226-inhull), MATLAB Central File Exchange. Retrieved .

Created with
R14SP1

Compatible with any release

**Inspired:**
3D mesh transform using sparse control points, Agglomorative Clustering for Fault Network Reconstruction, The Barycentric Fixed-Mass method for estimating fractal dimensions, inpolyhedron - are points inside a triangulated volume?, wrench2d.zip, STORM-based Relative Localization Analysis

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.

Jonathan McGuckinI am trying to find out if query points are located inside a computed convex hull. I am having trouble finding the vertices needed to use this function. it says as used by convhulln, but I cannot find where vertices are used in this function. Could you provide some indication to how I calculate the vertices?

Thanks

Rafael GaticaJimWorks well most of the time but I am not getting the error "Not enough unique points specified". I am testing points in convex hulls with 20 normally distributed points in 10 dimension. Might I need to change TOL or does the error mean some other problem?

LaurentYannisRobust and fast :-))

LukasSimon MuellerSunita SahaError using qhullmx

The data is degenerate in at least one dimension - ND set of points lying in (N+1)D space.

Error in delaunayn (line 101)

t = qhullmx(x', 'd ', opt);

Error in teest (line 32)

tess = delaunayn(M);

why I am getting this error

Sunita SahaError in convhulln (line 64)

[k,vv] = qhullmx(x', opt);

Please help me solving this error

Stefan SpelitzWorks great! Thanks a lot!

@His Teknoloji: Matlab's inpolygon only works in 2D space.

Matt JHis TeknolojiInstead of using writing own algorithm, one can use inpolygon built-in 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 contractorI 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 kantHii, 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 RoginskiPerhaps 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 . . .

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

SethJohn D'ErricoOh, 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 venkataHey 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.Very fast, works like a charm!

Jan FroehlichVery fast, thank you for submitting this.

John D'ErricoAndy - 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 non-convex. So, some simple alternatives are:

- Form a tessellation of the non-convex 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 HornHello 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

FrancescoVery fast and easy to use! I like how it does the tesselation of surfaces for you #lazy

QIANG LIUJohn D'ErricoDaniella - 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?

DanielaHi

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'ErricoI 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 in-plane 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.1633e-17

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.2204e-15

2.2204e-15

1.1102e-15

1.1102e-15

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.

DanielaThanks 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'ErricoDaniella - 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 3-d 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 2-d problem by projecting it all into the 2-d 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 2-d.

DanielaI 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?

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

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

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

DuyThis 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'ErricoSorry. 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 Voronovthere 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)

vickyThank you John for you valuable suggestions

John D'ErricoHi 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

vickyHello 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 PenaSorry, thanks anyway for your help, i'll try in the forum

John D'ErricoWhat 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 PenaBut 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'Erricoinhull 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 PenaBy 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 PenaThanks 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'ErricoFirst 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 Penacould 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 BellecGreat 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 GingrasAnother very usefull code provided by John! Thanks.

Navid SamavatiExcellent Code!

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

Nagarjunavery nice work

Luigi GiaccariVery 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 GiaccariVery 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 EdelhoferVery 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 Thakoorhuge time saver! thank you :)

Lee ShunnFast and easy to use. Thanks for the great code.

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

Jon KoSGreat 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 GaltEver heard of commenting your code?

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

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

Liketo StayanonymousNo, sorry! I think you lost me with this one. What (again) is this software about??

Clear example will be much appreciated!