Code covered by the BSD License  

Highlights from
Inhull

4.625

4.6 | 25 ratings Rate this file 80 Downloads (last 30 days) File Size: 3.06 KB File ID: #10226
image thumbnail

Inhull

by

 

04 Mar 2006 (Updated )

Efficient test for points inside a convex hull in n dimensions

| Watch this File

File Information
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 3 D Mesh Transform Using Sparse Control Points.

MATLAB release MATLAB 7.0.1 (R14SP1)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (47)
04 Aug 2014 Jan Froehlich

Very fast, thank you for submitting this.

26 Jul 2014 John D'Errico

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

25 Jul 2014 Andy Horn

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

18 Dec 2013 Francesco

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

05 Sep 2013 QIANG LIU  
29 Jul 2013 John D'Errico

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?

28 Jul 2013 Daniela

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.

23 Jul 2013 John D'Errico

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

23 Jul 2013 Daniela

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

22 Jul 2013 John D'Errico

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

22 Jul 2013 Daniela

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?

24 Jan 2013 Swagatika  
09 Oct 2012 Vincent Jaouen

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

13 Sep 2012 Thomas  
28 Aug 2012 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).

28 Aug 2012 Sven  
24 May 2012 Richard Crozier

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

12 Jun 2011 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.

08 Apr 2011 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

07 Apr 2011 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)

09 Feb 2011 vicky

Thank you John for you valuable suggestions

08 Feb 2011 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

08 Feb 2011 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.

02 Jun 2010 Antonio Pena

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

01 Jun 2010 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.

01 Jun 2010 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]);

01 Jun 2010 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.

01 Jun 2010 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?

01 Jun 2010 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

01 Jun 2010 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.

01 Jun 2010 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...

06 Feb 2010 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.

28 May 2009 David Gingras

Another very usefull code provided by John! Thanks.

02 Mar 2009 Navid Samavati

Excellent Code!

19 Feb 2009 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!!!!!

05 Dec 2008 Nagarjuna

very nice work

02 Dec 2008 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!

02 Dec 2008 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!

23 Jun 2008 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).

22 Feb 2008 Jayveer Thakoor

huge time saver! thank you :)

08 Nov 2007 Lee Shunn

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

31 Oct 2007 Yonas T. Weldeselassie

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

28 Mar 2007 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!

06 Mar 2007 John Galt

Ever heard of commenting your code?

23 Oct 2006 Quan Wen

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

01 May 2006 Mia Ginsburg

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

26 Apr 2006 Liketo Stayanonymous

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

Updates
13 Mar 2006

Added a tolerance on the tests

17 Mar 2006

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

05 Apr 2006

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

08 Apr 2011

Repaired example in the help

06 Sep 2012

minor changes for a tiny speed boost

Contact us