Code covered by the BSD License

Decimate Polygon

Anton Semechko (view profile)

17 Jan 2012 (Updated )

Simplify a 2D closed, piecewise linear contour by specifying boundary offset tolerance.

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

inpoly(p,node,edge,TOL)
```function [cn,on] = inpoly(p,node,edge,TOL)

%  INPOLY: Point-in-polygon testing.
%
% Determine whether a series of points lie within the bounds of a polygon
% in the 2D plane. General non-convex, multiply-connected polygonal
% regions can be handled.
%
% SHORT SYNTAX:
%
%   in = inpoly(p,node);
%
%   p   : The points to be tested as an Nx2 array [x1 y1; x2 y2; etc].
%   node: The vertices of the polygon as an Mx2 array [X1 Y1; X2 Y2; etc].
%         The standard syntax assumes that the vertices are specified in
%         consecutive order.
%
%   in  : An Nx1 logical array with IN(i) = TRUE if P(i,:) lies within the
%         region.
%
% LONG SYNTAX:
%
%  [in,on] = inpoly(p,node,edge);
%
%  edge: An Mx2 array of polygon edges, specified as connections between
%        the vertices in NODE: [n1 n2; n3 n4; etc]. The vertices in NODE
%        do not need to be specified in connsecutive order when using the
%        extended syntax.
%
%  on  : An Nx1 logical array with ON(i) = TRUE if P(i,:) lies on a
%        polygon edge. (A tolerance is used to deal with numerical
%        precision, so that points within a distance of
%        eps^0.8*norm(node(:),inf) from a polygon edge are considered "on"
%        the edge.
%
% EXAMPLE:
%
%   polydemo;       % Will run a few examples
%
% See also INPOLYGON

% The algorithm is based on the crossing number test, which counts the
% number of times a line that extends from each point past the right-most
% region of the polygon intersects with a polygon edge. Points with odd
% counts are inside. A simple implementation of this method requires each
% wall intersection be checked for each point, resulting in an O(N*M)
% operation count.
%
% This implementation does better in 2 ways:
%
%   1. The test points are sorted by y-value and a binary search is used to
%      find the first point in the list that has a chance of intersecting
%      with a given wall. The sorted list is also used to determine when we
%      have reached the last point in the list that has a chance of
%      intersection. This means that in general only a small portion of
%      points are checked for each wall, rather than the whole set.
%
%   2. The intersection test is simplified by first checking against the
%      bounding box for a given wall segment. Checking against the bbox is
%      an inexpensive alternative to the full intersection test and allows
%      us to take a number of shortcuts, minimising the number of times the
%      full test needs to be done.
%
%   Darren Engwirda: 2005-2007
%   Email          : d_engwirda@hotmail.com
%   Last updated   : 23/11/2007 with MATLAB 7.0
%
% Problems or suggestions? Email me.

%% ERROR CHECKING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<4
TOL = 1.0e-12;
if nargin<3
edge = [];
if nargin<2
error('Insufficient inputs');
end
end
end
nnode = size(node,1);
if isempty(edge)                                                           % Build edge if not passed
edge = [(1:nnode-1)' (2:nnode)'; nnode 1];
end
if size(p,2)~=2
error('P must be an Nx2 array.');
end
if size(node,2)~=2
error('NODE must be an Mx2 array.');
end
if size(edge,2)~=2
error('EDGE must be an Mx2 array.');
end
if max(edge(:))>nnode || any(edge(:)<1)
error('Invalid EDGE.');
end

%% PRE-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n  = size(p,1);
nc = size(edge,1);

% Choose the direction with the biggest range as the "y-coordinate" for the
% test. This should ensure that the sorting is done along the best
% direction for long and skinny problems wrt either the x or y axes.
dxy = max(p,[],1)-min(p,[],1);
if dxy(1)>dxy(2)
% Flip co-ords if x range is bigger
p = p(:,[2,1]);
node = node(:,[2,1]);
end
tol = TOL*min(dxy);

% Sort test points by y-value
[y,i] = sort(p(:,2));
x = p(i,1);

%% MAIN LOOP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cn = false(n,1);     % Because we're dealing with mod(cn,2) we don't have
% to actually increment the crossing number, we can
% just flip a logical at each intersection (faster!)
on = cn;
for k = 1:nc         % Loop through edges

% Nodes in current edge
n1 = edge(k,1);
n2 = edge(k,2);

% Endpoints - sorted so that [x1,y1] & [x2,y2] has y1<=y2
%           - also get xmin = min(x1,x2), xmax = max(x1,x2)
y1 = node(n1,2);
y2 = node(n2,2);
if y1<y2
x1 = node(n1,1);
x2 = node(n2,1);
else
yt = y1;
y1 = y2;
y2 = yt;
x1 = node(n2,1);
x2 = node(n1,1);
end
if x1>x2
xmin = x2;
xmax = x1;
else
xmin = x1;
xmax = x2;
end

% Binary search to find first point with y<=y1 for current edge
if y(1)>=y1
start = 1;
elseif y(n)<y1
start = n+1;
else
lower = 1;
upper = n;
for j = 1:n
start = round(0.5*(lower+upper));
if y(start)<y1
lower = start;
elseif y(start-1)<y1
break;
else
upper = start;
end
end
end

% Loop through points
for j = start:n
% Check the bounding-box for the edge before doing the intersection
% test. Take shortcuts wherever possible!

Y = y(j);   % Do the array look-up once & make a temp scalar
if Y<=y2
X = x(j);   % Do the array look-up once & make a temp scalar
if X>=xmin
if X<=xmax

% Check if we're "on" the edge
on(j) = on(j) || (abs((y2-Y)*(x1-X)-(y1-Y)*(x2-X))<tol);

% Do the actual intersection test
if (Y<y2) && ((y2-y1)*(X-x1)<(Y-y1)*(x2-x1))
cn(j) = ~cn(j);
end

end
elseif Y<y2   % Deal with points exactly at vertices
% Has to cross edge
cn(j) = ~cn(j);
end
else
% Due to the sorting, no points with >y
% value need to be checked
break
end
end

end

% Re-index to undo the sorting
cn(i) = cn|on;
on(i) = on;

end      % inpoly()
```

Contact us