Testing for Points on Border or Interior of a Polygon

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2007

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1];
Needs["Graphics`"]

Introduction

It is often necessary to numerically test whether a point is within a specified polygon or on its border.  Here we will develop some functions for this purpose and later we will use them in the pde relaxation notebook in problems that compute the electrostatic potential outside systems of electrodes in two dimensions.
We define a polygon 𝒫 = {(Subscript[z, 1]), (Subscript[z, 2]), ⋯, 
          (Subscript[z, n]), (Subscript[z, n + 1])} as a list of vertices whose first and last points coincide, such that (Subscript[z, n + 1]) = (Subscript[z, 1]).  Suppose that the point p is an interior point of 𝒫.  Let \(TraditionalForm\`{α\_1, α\_2, ⋯, α\_n}\) represent the n interior angles formed between lines from p to each (Subscript[z, i]), as sketched below.  If we interpret the points as lying in the complex plane, these angles can be computed easily according to

(Subscript[α, i])  =  
      Arg[(((Subscript[z, i + 1]) - p)/((Subscript[z, i]) - p))]

The identification of interior points is based upon the winding number

W = (1/(2  π)) (∑_(i = 1)^n) (Subscript[α, i])

that counts the number of times 𝒫 winds around p in a positive (counter-clockwise) sense; a negative winding number corresponds to a clockwise loop.  A point p is interior to 𝒫 iff W ≠ 0 and is exterior otherwise.

p = {1, 2};
poly = {{3, 4}, {1, 5}, {-2, 3}, {2, -1}, {3, 
            4}};
ListPlot[Append[poly, p], PlotStyle -> AbsolutePointSize[4], 
        Epilog -> {Line[poly], Dashing[{0.02, 0.02}], 
            Table[Line[{p, 
                  poly[[k]]}], {k, 
                Length[poly]}]}];

Implementation

winding number, wn

The following function evaluates the winding number.  The dot products with {1, i} convert the Cartesian coordinates to the complex numbers needed by Arg.  Chop ensures that vanishing winding numbers are returned as 0 instead of a numerically insignificant number.  Note that this function is undefined for points that coincide with any of the vertices.

wn[a_, poly_] := 
      Chop[(1/(2  π)) 
          Sum[Arg [(poly[[k + 1]] - a) . {1. , I}/((poly[[k]] - a)) . {1. , I}], {k, 1, Length[poly] - 1}]]

To test this function, we construct a closed polygon with several loops.  Selected points are numbered in the diagram.  Point 1 is clearly surrounded once in a positive sense, point 2 once in a negative sense, point 3 twice in a negative sense, and point 4 is exterior.  Also observe that point 5 is not enclosed while point 6 is encompassed once in a negative sense.  As shown below, our function successfully computes these winding numbers.

poly = {{0, 4}, {7, 4}, {7, -3}, {-2, -3}, {0, 2}, {5, 
            2}, {3, -2}, {-4, 3}, {-5, 0}, {-3, -3}, {0, 
            4}};
pts = {{-3, -1}, {5, 3}, {2, 1}, {-2, 
            4}, {-1.85, -1.5}, {-0.6, 1.5}};

TextListPlot[pts, PlotStyle -> AbsolutePointSize[4], 
        PlotRange -> {{-7.5, 7.5}, {-5, 5}}, Ticks -> None, 
        Epilog -> {Line[poly], 
            Table[Arrow[
                poly[[
                  k]], ((poly[[k]] + poly[[k + 1]])/2)], {k, 1, Length[poly] - 1}]}];

wn[#, poly] & /@ pts

OnBoundaryQ

Notice that the winding number vanishes for points that are on the boundary but not on any vertex.  For example,

wn[((2  poly[[2]] + 3  poly[[3]])/5), poly]

The following function determines whether there is a solution 0 ≤ t ≤ 1 to the equation

p ==  
      t (Subscript[z, i]) + (1 - t) (Subscript[z, i + 1])

for any vertex 1 ≤ i ≤ n.  The Table produces a list of True|False for each vertex and MemberQ asks whether there is a True in that list.

OnBoundaryQ[pt_, poly_] := 
      Module[{t}, 
        MemberQ[Table[
            If[0 ≤ (t /. 
                    Flatten[
                      Solve[pt == 
                          t*poly[[
                                k]] + (1 - t)*
                              poly[[k + 
                                  1]], t]]) ≤ 
                1, True, False, False], {k, 1, Length[poly] - 1}], 
          True]]

OnBoundaryQ[#, poly] & /@ 
      Join[pts, {poly[[2]], 
          0.27  poly[[2]] + 
            0.73  poly[[3]]}]

Notice that vertices are on the boundary.  The boundary need not be closed.

InPolygonQ

Finally, InPolygonQ combines these functions into a single query that returns True for points enclosed by the polygon or on its boundary and False otherwise.  If given an open curve, it closes the polygon by connecting the first and last points.  The interior and boundary tests are included explicitly so that this function is self-contained.

InPolygonQ[pt_, poly_] := 
      Module[{t, cpoly, n, k}, 
 (* 
          Prepare closed figure. *) 
cpoly = 
          If[Last[poly] ≠ First[poly], Append[poly, First[poly]], 
            poly]; 
n = 
          Length[cpoly] - 1; 
If[
 (* 
            Check whether point is on the boundary. *) 
MemberQ[
Table[
              If[0 ≤ (t /. 
                      Flatten[
                        Solve[pt == 
                            t*cpoly[[
                                  k]] + (1 - t)*
                                cpoly[[
                                  k + 1]], 
                          t]]) ≤ 1, True, False, False], {k, 1, 
                n}], True], 
          True, 
 (* 
            Use winding number to test whether the point is inside. *) 
Chop[(1/(2  π)) 
                Sum[Arg [(cpoly[[k + 1]] - pt) . {1. , I}/((cpoly[[k]] - pt)) . {1. , I}], {k, 1, n}]] ≠ 0]
]

InPolygonQ[#, poly] & /@ 
      Join[pts, {poly[[2]], 
          0.27  poly[[2]] + 
            0.73  poly[[3]]}]