function in = inpolygon(x,y,xv,yv)
%INPOLYGON True for points inside polygonal region.
%   IN = INPOLYGON(X, Y, XV, YV) returns a matrix IN the size of
%   X and Y.  IN(p,q) = 1 if the point (X(p,q), Y(p,q)) is
%   strictly inside the polygonal region whose vertices are
%   specified by the vectors XV and YV;  IN(p,q) is 0.5 if
%   the point is on the polygon; otherwise IN(p,q) = 0.
%
%   Example:
%       xv = rand(6,1); yv = rand(6,1);
%       xv = [xv ; xv(1)]; yv = [yv ; yv(1)];
%       x = rand(1000,1); y = rand(1000,1);
%       in = inpolygon(x,y,xv,yv);
%       plot(xv,yv,x(in),y(in),'.r',x(~in),y(~in),'.b')

%   Steven L. Eddins, July 1994
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.5 $  $Date: 1997/11/21 23:40:39 $

% If (xv,yv) is not closed, close it.
xv = xv(:);
yv = yv(:);
Nv = length(xv);
if ((xv(1) ~= xv(Nv)) | (yv(1) ~= yv(Nv)))
  xv = [xv ; xv(1)];
  yv = [yv ; yv(1)];
  Nv = Nv + 1;
end

% Make matrices to vectorize the computation.
% Translate the vertices so that the test points are
% at the origin.

inputSize = size(x);
x = x(:).';
y = y(:).';
Np = length(x);
x = x(ones(Nv,1),:);
y = y(ones(Nv,1),:);
xv = xv(:,ones(1,Np)) - x;
yv = yv(:,ones(1,Np)) - y;

% Compute the quadrant number for the vertices relative
% to the test points.
posX = xv > 0;
posY = yv > 0;
negX = ~posX;
negY = ~posY;
quad = (negX & posY) + 2*(negX & negY) + ...
  3*(posX & negY);

% Compute the sign() of the cross product and dot product
% of adjacent vertices.
m = 1:Nv-1;
mp1 = 2:Nv;
signCrossProduct = sign(xv(m,:) .* yv(mp1,:) ...
        - xv(mp1,:) .* yv(m,:));
dotProduct = xv(m,:) .* xv(mp1,:) + yv(m,:) .* yv(mp1,:);

% Compute the vertex quadrant changes for each test point.
diffQuad = diff(quad);

% Fix up the quadrant differences.  Replace 3 by -1 and -3 by 1.
% Any quadrant difference with an absolute value of 2 should have
% the same sign as the cross product.
idx = (abs(diffQuad) == 3);
diffQuad(idx) = -diffQuad(idx)/3;
idx = (abs(diffQuad) == 2);
diffQuad(idx) = 2*signCrossProduct(idx);

% Find the inside points.
in = (sum(diffQuad) ~= 0);

% Find the points on the polygon.  If the cross product is 0 and
% the dot product is nonpositive anywhere, then the corresponding
% point must be on the contour.
on = find(any((signCrossProduct == 0) & (dotProduct <= 0)));

% Indicate points on the polygon with a value of 0.5.
in(on) = 0.5*ones(size(on));

% Reshape output matrix.
in = reshape(in, inputSize);
