Computing the winding number of a point vis-à-vis a polygon

[Mark Hunter]

In another thread Stuart illustrated an obscure feature of PB’s GRAPHIC POLYGON statement. It can color the polygon in two colors: regions with an odd winding number get one color, regions with an even winding number the other color.

Why odd and even? A more general feature would allow you to specify an array of colors for winding numbers from minus whatever to plus whatever, the whatevers given by the upper and lower bounds of the array.

The problem I consider here is this: Given a polygon with a specific ordering of its vertices, and any point in the plane, what is the winding number at that point?

My example polygon is Stuart’s septigram except that I reversed the order of the vertices so the tendency of the path is counter-clockwise. That way the winding number of the central region (and in this case all the other regions) is positive rather than negative.

This method uses vector algebra (trigonometry and determinants). I imagine a simpler method exists but anyway this one works.

One intriguing fact you’ll discover is that if the point happens to be lie right on an edge, between, say, regions having a winding number of 2 and 3 respectively, the winding number is exactly 2.5 — a curious solution to Leibnitz’s doctrine of cause-effect continuity. On that subject, see also this old post quoting a program I wrote.

' Compute winding numbers
'
' For PBCC4 or 5, and probably 6.
' To convert to PBWin, replace Print with MsgBox and delete WaitKey$.
'
' Given a polygon whose vertices are in a specific order,
' and a point in the plane,
' what is that point's winding number with respect to the polygon?
'
' That is, if a man stood at the point and watched someone walk
' all the way around the polygon in the order of the vertices, how
‘ many times would be turn completely around counter-clockwise,
‘ understanding that turning clockwise undoes turning counter-clockwise.
'------------------------------------------

Global PI, PH, PD As Single
'------------------------------------------

' The inverse of Cos().
' -1 <= x <= 1,  PI <= ArcCos(x) <= 0
Function ArcCos(ByVal x As Single) As Single
 If     x >=  1 Then
  Function = 0
 ElseIf x <= -1 Then
  Function = PI
 Else
  Function = PH - Atn(x / Sqr(1 - x * x))
 End If
End Function
'------------------------------------------

' n  is the number of vertices x(), y().
' a, b is a point in the plane.
' Return the winding number of the point.
Function WindingNumber(ByVal n As Long, x() As Single, y() As Single, ByVal a As Single, ByVal b As Single) As Single
 Local i As Long                         'index into vertices, starts at 0
 Local vx, vy As Single                  'vector from a, b to x, y
 Local norm As Single                    'norm of vx, vy
 Dim ux(n) As Single, uy(n) As Single    'unit vector from a, b to x, y
 Local theta As Single                   'signed angle subtended by a side of the polygon
 Local det As Single                     'determinant
 Local Angle As Single                   'signed total angle
 Local fracky As Single                  'fractional part
 Local WN As Single

 ' Add an extra vertex  x(n), y(n)  duplicating  x(0), y(0).
 x(n) = x(0) : y(n) = y(0)

 ' Create unit vectors pointing from the given point (a, b)
 ' to each of the polygon's vertices.
 For i = 0 To n
  vx = x(i) - a : vy = y(i) - b
  norm = Sqr(vx *vx + vy * vy)
  ux(i) = vx / norm : uy(i) = vy / norm
 Next

 ' The dot product of two vectors v and w  is  vx * wx + vy * wy
 ' and that equals the product of their norms times the cosine
 ' of the angle between them.  Thus the angle between two unit
 ' vectors is the ArcCos of their dot product.  That, however,
 ' is the unsigned angle.  To get the sign, compute the determinant
 ' vx * wy - vy * wx  which is the signed area spanned by v and w,
 ' and use its sign.
 Angle = 0
 For i = 0 To n - 1
  theta = ArcCos( ux(i) * ux(i + 1) + uy(i) * uy(i + 1) )
  det = ux(i) * uy(i + 1) - uy(i) * ux(i + 1)
  Angle = Angle + theta * Sgn(det)
 Next

 ' Angle is in radians and will be a multiple of PD.
 ' Check that it is.
 
 WN = Angle / PD
 fracky = Frac(WN)
 If Abs(fracky - .5) < .001 Then
  Print "This point is on an edge (WN half-integral):  ";
 ElseIf .001 < fracky And fracky < .999 Then
  Print "Back to the drawing board:  ";
 End If

 ' Avoid printing out minus zero
 If Abs(WN) < .001 Then WN = 0
 Function = WN

End Function
'------------------------------------------

%n = 7                                  'total number of vertices

Function PBMain
 Dim x(%n) As Single, y(%n) As Single   'the polygon, vertices 0 to %n - 1, then %n being the first again
 Local a, b As Single                   'a point anywhere in the plane
 Local i As Long                        'count
 Local WN As Single                     'the winding number
 Local fracky As Single                 'fractional part of WN

 PI = 4 * Atn(1)                        'global, for subroutines
 PH = PI / 2                            '
 PD = PI * 2                            '

 Array Assign x() =  58, 180,  0, 200,  22, 144, 100     'index 0 to %n - 1
 Array Assign y() =   0, 162, 72,  72, 162,   0, 200      '

 a = 101       'in the central region, WN should be 3 for this polygon
 b =  95
 Print "(";Format$(a);", ";Format$(b);")", :GoSub PrintWN

 a = 300       'way outside the polygon, WN should be 0
 b = 300

 Print "(";Format$(a);", ";Format$(b);")", :GoSub PrintWN

 a = 110       'in a region near central
 b = 120
 Print "(";Format$(a);", ";Format$(b);")", :GoSub PrintWN

 a =  80       'in a region near outside
 b = 120
 Print "(";Format$(a);", ";Format$(b);")", :GoSub PrintWN

 Print "-----------------"
 ' A few random points.
 For i = 1 To 16
  a = Rnd(50, 100) : b = Rnd(50, 100)
  Print "(";Format$(a);", ";Format$(b);")", :GoSub PrintWN
 Next

 WaitKey$
 Exit Function
 
PrintWN:
 WN = WindingNumber(%n, x(), y(), a, b)
 fracky = Frac(WN)
 If Abs(fracky - .5) < .01 Then
  Print Format$(WN, "#.00")
 Else
  Print Format$(WN,"#")
 End If
 Return
End Function