


INHULL - tests if a set of points are inside a convecx hull
inhull: tests if a set of points are inside a convex hull
usage: in = inhull(testpts,xyz)
usage: in = inhull(testpts,xyz,tess)
usage: in = inhull(testpts,xyz,tess,tol)
arguments: (input)
testpts - nxp array to test, n data points, in p dimensions
If you have many points to test, it is most efficient to
call this function once with the entire set.
xyz - mxp array of vertices of the convex hull, as used by
convhulln.
tess - tessellation (or triangulation) generated by convhulln
If tess is left empty or not supplied, then it will be
generated.
tol - (OPTIONAL) tolerance on the tests for inclusion in the
convex hull. You can think of tol as the distance a point
may possibly lie outside the hull, and still be perceived
as on the surface of the hull. Because of numerical slop
nothing can ever be done exactly here. I might guess a
semi-intelligent value of tol to be
tol = 1.e-13*mean(abs(xyz(:)))
In higher dimensions, the numerical issues of floating
point arithmetic will probably suggest a larger value
of tol.
DEFAULT: tol = 0
arguments: (output)
in - nx1 logical vector
in(i) == 1 --> the i'th point was inside the convex hull.
Example usage: The first point should be inside, the second out
xy = randn(20,2)
tess = convhull(xy(:,1),xy(:,2));
testpoints = [ 0 0; 10 10];
in = inhull(testpts,xyz,tess)
in =
1
0
A non-zero count of the number of degenerate simplexes in the hull
will generate a warning (in 4 or more dimensions.) This warning
may be disabled off with the command:
warning('off','inhull:degeneracy')
See also: convhull, convhulln, delaunay, delaunayn, tsearch, tsearchn
Author: John D'Errico
e-mail: woodchips@rochester.rr.com
Release: 3.0
Release date: 10/26/06

0001 function in = inhull(testpts,xyz,tess,tol) 0002 %INHULL - tests if a set of points are inside a convecx hull 0003 % 0004 % inhull: tests if a set of points are inside a convex hull 0005 % usage: in = inhull(testpts,xyz) 0006 % usage: in = inhull(testpts,xyz,tess) 0007 % usage: in = inhull(testpts,xyz,tess,tol) 0008 % 0009 % arguments: (input) 0010 % testpts - nxp array to test, n data points, in p dimensions 0011 % If you have many points to test, it is most efficient to 0012 % call this function once with the entire set. 0013 % 0014 % xyz - mxp array of vertices of the convex hull, as used by 0015 % convhulln. 0016 % 0017 % tess - tessellation (or triangulation) generated by convhulln 0018 % If tess is left empty or not supplied, then it will be 0019 % generated. 0020 % 0021 % tol - (OPTIONAL) tolerance on the tests for inclusion in the 0022 % convex hull. You can think of tol as the distance a point 0023 % may possibly lie outside the hull, and still be perceived 0024 % as on the surface of the hull. Because of numerical slop 0025 % nothing can ever be done exactly here. I might guess a 0026 % semi-intelligent value of tol to be 0027 % 0028 % tol = 1.e-13*mean(abs(xyz(:))) 0029 % 0030 % In higher dimensions, the numerical issues of floating 0031 % point arithmetic will probably suggest a larger value 0032 % of tol. 0033 % 0034 % DEFAULT: tol = 0 0035 % 0036 % arguments: (output) 0037 % in - nx1 logical vector 0038 % in(i) == 1 --> the i'th point was inside the convex hull. 0039 % 0040 % Example usage: The first point should be inside, the second out 0041 % 0042 % xy = randn(20,2) 0043 % tess = convhull(xy(:,1),xy(:,2)); 0044 % testpoints = [ 0 0; 10 10]; 0045 % in = inhull(testpts,xyz,tess) 0046 % 0047 % in = 0048 % 1 0049 % 0 0050 % 0051 % A non-zero count of the number of degenerate simplexes in the hull 0052 % will generate a warning (in 4 or more dimensions.) This warning 0053 % may be disabled off with the command: 0054 % 0055 % warning('off','inhull:degeneracy') 0056 % 0057 % See also: convhull, convhulln, delaunay, delaunayn, tsearch, tsearchn 0058 % 0059 % Author: John D'Errico 0060 % e-mail: woodchips@rochester.rr.com 0061 % Release: 3.0 0062 % Release date: 10/26/06 0063 0064 % get array sizes 0065 % m points, p dimensions 0066 p = size(xyz,2); 0067 [n,c] = size(testpts); 0068 if p ~= c 0069 error 'testpts and xyz must have the same number of columns' 0070 end 0071 if p < 2 0072 error 'Points must lie in at least a 2-d space.' 0073 end 0074 0075 % was the convex hull supplied? 0076 if (nargin<3) || isempty(tess) 0077 tess = convhulln(xyz); 0078 end 0079 [nt,c] = size(tess); 0080 if c ~= p 0081 error 'tess array is incompatible with a dimension p space' 0082 end 0083 0084 % was tol supplied? 0085 if (nargin<4) || isempty(tol) 0086 tol = 0; 0087 end 0088 0089 % build normal vectors 0090 switch p 0091 case 2 0092 % really simple for 2-d 0093 nrmls = (xyz(tess(:,1),:) - xyz(tess(:,2),:)) * [0 1;-1 0]; 0094 0095 % Any degenerate edges? 0096 del = sqrt(sum(nrmls.^2,2)); 0097 degenflag = (del<(max(del)*10*eps)); 0098 if sum(degenflag)>0 0099 warning('inhull:degeneracy',[num2str(sum(degenflag)), ... 0100 ' degenerate edges identified in the convex hull']) 0101 0102 % we need to delete those degenerate normal vectors 0103 nrmls(degenflag,:) = []; 0104 nt = size(nrmls,1); 0105 end 0106 case 3 0107 % use vectorized cross product for 3-d 0108 ab = xyz(tess(:,1),:) - xyz(tess(:,2),:); 0109 ac = xyz(tess(:,1),:) - xyz(tess(:,3),:); 0110 nrmls = cross(ab,ac,2); 0111 degenflag = repmat(false,nt,1); 0112 otherwise 0113 % slightly more work in higher dimensions, 0114 nrmls = zeros(nt,p); 0115 degenflag = repmat(false,nt,1); 0116 for i = 1:nt 0117 % just in case of a degeneracy 0118 nullsp = null(xyz(tess(i,2:end),:) - repmat(xyz(tess(i,1),:),p-1,1))'; 0119 if size(nullsp,1)>1 0120 degenflag(i) = true; 0121 nrmls(i,:) = NaN; 0122 else 0123 nrmls(i,:) = nullsp; 0124 end 0125 end 0126 if sum(degenflag)>0 0127 warning('inhull:degeneracy',[num2str(sum(degenflag)), ... 0128 ' degenerate simplexes identified in the convex hull']) 0129 0130 % we need to delete those degenerate normal vectors 0131 nrmls(degenflag,:) = []; 0132 nt = size(nrmls,1); 0133 end 0134 end 0135 0136 % scale normal vectors to unit length 0137 nrmllen = sqrt(sum(nrmls.^2,2)); 0138 nrmls = nrmls.*repmat(1./nrmllen,1,p); 0139 0140 % center point in the hull 0141 center = mean(xyz,1); 0142 0143 % any point in the plane of each simplex in the convex hull 0144 a = xyz(tess(~degenflag,1),:); 0145 0146 % ensure the normals are pointing inwards 0147 dp = sum((repmat(center,nt,1) - a).*nrmls,2); 0148 k = dp<0; 0149 nrmls(k,:) = -nrmls(k,:); 0150 0151 % We want to test if: dot((x - a),N) >= 0 0152 % If so for all faces of the hull, then x is inside 0153 % the hull. Change this to dot(x,N) >= dot(a,N) 0154 aN = sum(nrmls.*a,2); 0155 0156 % test, be careful in case there are many points 0157 in = repmat(false,n,1); 0158 0159 % if n is too large, we need to worry about the 0160 % dot product grabbing huge chunks of memory. 0161 memblock = 1e6; 0162 blocks = max(1,floor(n/(memblock/nt))); 0163 aNr = repmat(aN,1,length(1:blocks:n)); 0164 for i = 1:blocks 0165 j = i:blocks:n; 0166 if size(aNr,2) ~= length(j), 0167 aNr = repmat(aN,1,length(j)); 0168 end 0169 in(j) = all((nrmls*testpts(j,:)' - aNr) >= -tol,1)'; 0170 end 0171 0172 0173