function ind=inside(x,y,xv,yv,tol);
% INSIDE  determines if points are inside/outside a polygon.
%         IND=INSIDE(X,Y,XV,YV) is a vector of flags corresponding to
%         the vector of points X,Y. If
%
%       -  X(i),Y(i) is outside the polygon XV,YV           : IND(i)=0
%       -  X(i),Y(i) is inside   "   "       "  "           : IND(i)=1
%       -  X(i),Y(i) is on an exterior edge of the polygon  : IND(i)=2
%       -  X(i),Y(i) is an exterior vertex of the polygon   : IND(i)=4
%
%         XV,YV is an ordered list of polygon vertices. The last point
%         in the list (identical to the first point) may be omitted.
%
%         IND=INSIDE(...,TOL) sets a tolerance for assigning points
%         to edges or vertices. (the default is 1e-6).

%RP (WHOI) 19/Dec/91
%   - this routine uses the fact that the sum of angles between
%     lines from a point to successive vertices of a polygon is 0
%     for a point outside the polygon, and +/- 2*pi for a point
%     inside. A point on an edge will have one such succesive angle
%     equal to +/- pi. 

if (nargin<5),
   tol=1e-6;
end;

% Put things into rows or columns...
x=x(:)';
y=y(:)';
xv=xv(:);
yv=yv(:);

NP=length(x);
NV=length(xv);

   % If last point is missing append it

if ( (xv(NV) ~= xv(1)) | (yv(NV) ~= yv(1))),
   xv=[xv;xv(1)];
   yv=[yv;yv(1)];
   NV=NV+1;
end;                      

    % distances from points to vertices

dx=ones(NV,1)*x-xv*ones(1,NP);  
dy=ones(NV,1)*y-yv*ones(1,NP);

    % Test for points at vertices

vertexpt=any( abs(dx)<tol & abs(dy)<tol );

    % Compute angular intervals between vertex points. We shift
    % to put things in the range -pi<= ang <pi, with some slop for
    % possible edge points

angs=rem( diff(atan2(dy,dx)) + 3*pi+tol , 2*pi ) - pi-tol;

    % If ang= -pi then we are on an edge

edgept=any( abs(angs+pi)<tol );

    % If the angles sum to zeros, we are outside - if they sum to 2*pi,
    % we are inside

interiorpt=abs( sum(angs) )>1.5;

    % If a point is on an edge or vertex which is inside the polygon
    % we classify this as interior (this handles some wierd special cases).

ind=4*vertexpt;
ind(edgept)=2*ones(sum(edgept),1);
ind(interiorpt)=ones(sum(interiorpt),1);

