www.pudn.com > NumericalComputing.rar > inregion.m, change:2003-09-20,size:4117b

```function [in, on] = inregion(x,y,xv,yv)
%INREGION True for points inside or on a polygonal region.
%   IN = INREGION(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 either strictly inside or
%   on the edge of the polygonal region whose vertices are specified by the
%   vectors XV and YV; otherwise IN(p,q) = 0.
%
%   [IN ON] = INREGION returns a second matrix, ON, which is the size of X
%   and Y.  ON(p,q) = 1 if the point (X(p,q), Y(p,q)) is on the edge of the
%   polygonal region; otherwise ON(p,q) = 0.
%
%   INREGION is a modification of INPOLYGON that uses a roundoff error
%   compensating tolerance in the cross product sign test.
%
%   Example:
%     xv = [-3 -3 1 1 3 1 -1 -1 -3];
%     yv = [-3 -1 3 1 1 -1 -1 -3 -3];
%     [x,y] = meshgrid(-3:1/2:3);
%     [in,on] = inregion(x,y,xv,yv);
%     p = find(in-on);
%     q = find(on);
%     plot(xv,yv,'-',x(p),y(p),'ko',x(q),y(q),'ro')
%     axis([-5 5 -4 4])

% 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

inputSize = size(x);

x = x(:).';
y = y(:).';

mask = (x >= min(xv)) & (x <= max(xv)) & (y>=min(yv)) & (y<=max(yv));
in = zeros(inputSize);
on = in;
return
end

% Choose block_length to keep memory usage of vec_inpolygon around
% 10 Megabytes.
block_length = 1e5;

M = prod(size(x));

if M*Nv < block_length
if nargout > 1
[in on] = vec_inpolygon(Nv,x,y,xv,yv);
else
in = vec_inpolygon(Nv,x,y,xv,yv);
end
else
% Process at most N elements at a time
N = ceil(block_length/Nv);
in = false(1,M);
if nargout > 1
on = false(1,M);
end
n1 = 0;  n2 = 0;
while n2 < M,
n1 = n2+1;
n2 = n1+N;
if n2 > M,
n2 = M;
end
if nargout > 1
[in(n1:n2) on(n1:n2)] = vec_inpolygon(Nv,x(n1:n2),y(n1:n2),xv,yv);
else
in(n1:n2) = vec_inpolygon(Nv,x(n1:n2),y(n1:n2),xv,yv);
end
end
end

if nargout > 1
end

% Reshape output matrix.

%----------------------------------------------
function [in, on] = vec_inpolygon(Nv,x,y,xv,yv)
% vectorize the computation.

% Translate the vertices so that the test points are
% at the origin.

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
% Modified 09/17/03 to use a tolerance in the cross product sign test.
m = 1:Nv-1;
mp1 = 2:Nv;
crossProduct = xv(m,:) .* yv(mp1,:) - xv(mp1,:) .* yv(m,:);
tol = 10*Nv*(max(abs(xv(:)))+max(abs(yv(:))))*eps;
crossProduct(abs(crossProduct)<tol) = 0;
signCrossProduct = sign(crossProduct);
dotProduct = xv(m,:) .* xv(mp1,:) + yv(m,:) .* yv(mp1,:);

% Compute the vertex quadrant changes for each test point.

% 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.