


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.

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