inside

 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 near (but outside) the polygon  : IND(i)=2

         Points classified as "near" are actually outside
         the polygon, but less than approx TOL away (This is sometimes
         useful for identifying "peninsulas" when gridding data).

         XV,YV is an ordered list of polygon vertices. The last point
         in the list (identical to the first point) may be omitted.
         There is no restriction on the shape of the polygon.

         IND=INSIDE(...,TOL) sets a tolerance for identifying "nearby"
         points. (the default is 1e-6).

0001 function ind=inside(x,y,xv,yv,tol);
0002 % INSIDE  determines if points are inside/outside a polygon.
0003 %         IND=INSIDE(X,Y,XV,YV) is a vector of flags corresponding to
0004 %         the vector of points X,Y. If
0005 %
0006 %       -  X(i),Y(i) is outside the polygon XV,YV       : IND(i)=0
0007 %       -  X(i),Y(i) is inside   "   "       "  "       : IND(i)=1
0008 %       -  X(i),Y(i) is near (but outside) the polygon  : IND(i)=2
0009 %
0010 %         Points classified as "near" are actually outside
0011 %         the polygon, but less than approx TOL away (This is sometimes
0012 %         useful for identifying "peninsulas" when gridding data).
0013 %
0014 %         XV,YV is an ordered list of polygon vertices. The last point
0015 %         in the list (identical to the first point) may be omitted.
0016 %         There is no restriction on the shape of the polygon.
0017 %
0018 %         IND=INSIDE(...,TOL) sets a tolerance for identifying "nearby"
0019 %         points. (the default is 1e-6).
0020 
0021 %RP (WHOI) 19/Dec/91
0022 %   - this routine uses the fact that the sum of angles between
0023 %     lines from a point to successive vertices of a polygon is 0
0024 %     for a point outside the polygon, and +/- 2*pi for a point
0025 %     inside. A point on an edge will have one such succesive angle
0026 %     equal to +/- pi.
0027 %RP (WHOI) 5/Dec/92
0028 %   - some bulletproofing for really wierd checkerboard polygons and
0029 %     single point "polygons".
0030 %RP (WHOI) 6/Dec/92
0031 %   - Redid edge stuff a little better to allow  "rounded gridding"
0032 %     i.e. getting grid points within some TOL of the polygon.
0033 %
0034 %    updated to version 5.1  AN
0035 
0036 if (nargin<5),
0037    tol=1e-6;
0038 end;
0039 
0040 % Put things into rows or columns...
0041 x=x(:)';
0042 y=y(:)';
0043 xv=xv(:);
0044 yv=yv(:);
0045 
0046 NP=length(x);
0047 NV=length(xv);
0048 
0049 if (NP~=length(y)), error('Point x,y vectors of unequal length!'); end;
0050 if (NV~=length(yv)), error('Polygon vertex x,y vectors of unequal length!');end;
0051 
0052    % If last point is missing append it
0053 
0054 if ( (xv(NV) ~= xv(1)) | (yv(NV) ~= yv(1))),
0055    xv=[xv;xv(1)];
0056    yv=[yv;yv(1)];
0057    NV=NV+1;
0058 end;                      
0059 
0060     % distances from points to vertices
0061 
0062 dx=ones(NV,1)*x-xv*ones(1,NP);  
0063 dy=ones(NV,1)*y-yv*ones(1,NP);
0064 
0065     % If we have only 1 point in the polygon, we can only search
0066     % for matches to the "vertex". Otherwise we do all the tests.
0067 
0068 if (NV==1),
0069    vertexpt=( abs(dx)<tol & abs(dy)<tol );
0070    ind=2*vertexpt;
0071 else
0072 
0073         % Test for points at vertices
0074 
0075    vertexpt=any( abs(dx)<tol & abs(dy)<tol );
0076 
0077         % Compute angular intervals between vertex points. We shift
0078         % to put things in the range -pi<= ang <pi, with some slop for
0079         % possible edge points
0080 
0081    angs=rem( diff(atan2(dy,dx)) + 3*pi+eps , 2*pi ) - pi-eps;
0082 
0083         % Look for points "near" the edge - acts to smear out the polygon
0084         % edge by tol - useful for(maybe) neglecting shallow water?
0085  
0086         % If a point is closest to a given segment at an endpoint, it
0087         % may be a vertex - since we have found these points above, we can
0088         % ignore them. If the closest point is in the middle of a segment,
0089         % classify it as an edge if the dx or dy displacements are < tol
0090         % (strictly speaking I suppose we should test sqrt(dx.^2+dy.^2)<tol
0091         % but this is faster).
0092 
0093    dxv=(xv(2:NV)-xv(1:NV-1))*ones(1,NP);
0094    dyv=(yv(2:NV)-yv(1:NV-1))*ones(1,NP);
0095    dxl=dx(1:NV-1,:);
0096    dyl=dy(1:NV-1,:);
0097    
0098         % Projection onto segment. if 0<lgt<1 closest point is in the
0099         % middle of the segment!
0100 
0101    lgt=(dxl.*dxv + dyl.*dyv)./(dxv.*dxv+dyv.*dyv);
0102    distx=2*tol*ones(NV-1,NP); disty=distx;
0103 
0104    ii=(lgt>0 & lgt < 1);
0105    distx(ii)=dxl(ii)-lgt(ii).*dxv(ii);
0106    disty(ii)=dyl(ii)-lgt(ii).*dyv(ii);
0107    
0108    edgept= any( abs(distx)<tol & abs(disty)<tol );
0109 
0110         % If the angles sum to 0, 4*pi, etc., we are outside - if they
0111         % sum to  2*pi, 6*pi, etc. we are inside. This will result in a
0112         % "checkerboard" inside/outside classification when polygon
0113         % segments intersect.
0114 
0115    interiorpt=rem( abs(round( sum(angs)/(2*pi) )) , 2);
0116 
0117         % Classify points by priority
0118 
0119    ind=2*edgept;
0120    ind(logical(vertexpt))=2*ones(sum(vertexpt),1);
0121    ind(logical(interiorpt))=ones(sum(interiorpt),1);
0122 end;