function ind=insidepts(xpoly, ypoly, xpts, ypts);
% INSIDEPTS:  finds all the points strictly inside a polygon.
%         
%   USAGE:  ind = insidepts(xpoly,ypoly,xpts,ypts);
%             xpoly and ypoly define the polygon
%             xpts and ypts define the set of points to test
%
%   See also:  INSIDEP, upon which this is based.
%
%   John Evans (jevans@sakhalin.er.usgs.gov)
%

n = length(xpoly);

for k = 1:length(xpts)
  c = 0;

  x = xpts(k);
  y = ypts(k);
  

  j = n;
  for i = 1:n

    if ( ( (ypoly(i)<=y) & (y<ypoly(j)) ) | ( (ypoly(j)<=y) & (y<ypoly(i)) ) ) 
      if (x < (xpoly(j) - xpoly(i)) * (y - ypoly(i)) / (ypoly(j) - ypoly(i)) + xpoly(i) )
        c = ~c;
      end
    end
  
    j = i;
  end

  if (c)
    ind = [ind; k];
  end
  
end

return;




