Polygon winding number

Compute the winding number of a polygon around a reference point or multiple reference points.

% Syntax:

% w = windnum(poly);

% w = windnum(poly,v);

% w = windnum(poly,v,tol);

%

% Inputs:

%

% poly: 2-column matrix of polygon vertices ([x,y] coordinate pairs),

% size(poly,1)>=3.

%

% v: 2-column matrix of reference points ([x,y] coordinate pairs),

% OPTIONAL, default = [0,0].

%

% tol: scalar, non-negative, proximity tolerance for point equality

% testing, OPTIONAL, default = 0. Points [x1,y1] and [x2,y2] are

% considered to be proximate if norm([x1-x2,y1-y2])<=tol. (CAUTION: The

% result might be in error due to numerical precision limits if tol = 0

% and a reference point is on or very close to the polygon boundary.)

%

% Output:

%

% w: column vector of winding numbers (length(w)==size(v,1)), integer or

% nan. For a simple closed, non-self-intersecting polygon and vector v

% (size-[1,2]), w is one of the following:

% * w = 0 if v is outside of the polygon with proximity clearance tol.

% * w = +/1 if v is inside of the polygon with proximity clearance tol,

% +1 if the polygon is CCW-oriented and -1 if it is CW-oriented.

% * nan if v is on the polygon boundary within tolerance tol.

% More generally, w is the angle swept by the polygon at point v, divided

% by 2*pi, or nan is point v is proximate to any polygon point. If v has

% multiple rows, then w(i) = windnum(poly,v(i,:),tol).