inhull

PURPOSE ^

INHULL - tests if a set of points are inside a convecx hull

SYNOPSIS ^

function in = inhull(testpts,xyz,tess,tol)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003