isinpoly

PURPOSE ^

ISIN - indicate if points are inside or outside a polygon

SYNOPSIS ^

function isin = isinpoly(x,y,xp,yp)

DESCRIPTION ^

ISIN - indicate if points are inside or outside a polygon

 ISIN = ISINPOLY(X,Y,XP,YP)   Finds whether points with
           coordinates X and Y are inside or outside of
           a polygon with vertices XP, YP.
           Returns matrix ISIN of the same size as X and Y
           with 0 for points outside a polygon, 1 for
           inside points and  0.5 for points belonging
           to a polygon XP, YP itself.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function  isin = isinpoly(x,y,xp,yp)
0002 %ISIN - indicate if points are inside or outside a polygon
0003 %
0004 % ISIN = ISINPOLY(X,Y,XP,YP)   Finds whether points with
0005 %           coordinates X and Y are inside or outside of
0006 %           a polygon with vertices XP, YP.
0007 %           Returns matrix ISIN of the same size as X and Y
0008 %           with 0 for points outside a polygon, 1 for
0009 %           inside points and  0.5 for points belonging
0010 %           to a polygon XP, YP itself.
0011 
0012 %  Copyright (c) 1995  by Kirill K. Pankratov
0013 %    kirill@plume.mit.edu
0014 %    4/10/94, 8/26/94.
0015 
0016  % Handle input ::::::::::::::::::::::::::::::::
0017 if nargin<4
0018   fprintf('\n  Error: not enough input arguments.\n\n')
0019   return
0020 end
0021 
0022  % Make the contour closed and get the sizes
0023 xp = [xp(:); xp(1)];
0024 yp = [yp(:); yp(1)];
0025 sz = size(x);
0026 x = x(:); y = y(:);
0027 
0028 lp = length(xp); l = length(x);
0029 ep = ones(1,lp); e = ones(1,l);
0030 
0031  % Calculate cumulative change in azimuth from points x,y
0032  % to all vertices
0033 A = diff(atan2(yp(:,e)-y(:,ep)',xp(:,e)-x(:,ep)'))/pi;
0034 A = A+2*((A<-1)-(A>1));
0035 isin = any(A==1)-any(A==-1);
0036 isin = (abs(sum(A))-isin)/2;
0037 
0038  % Check for boundary points
0039 A = (yp(:,e)==y(:,ep)')&(xp(:,e)==x(:,ep)');
0040 fnd = find(any(A));
0041 isin(fnd) = .5*ones(size(fnd));
0042 isin = round(isin*2)/2;
0043 
0044  % Reshape output to the input size
0045 isin = reshape(isin,sz(1),sz(2));
0046

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