No announcement yet.

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

  • Filter
  • Time
  • Show
Clear All
new posts

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

    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.
    Click image for larger version  Name:	Clipboard01.gif Views:	1 Size:	10.9 KB ID:	816403

    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
      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
     ' 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)
     ' 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
     Exit Function
     WN = WindingNumber(%n, x(), y(), a, b)
     fracky = Frac(WN)
     If Abs(fracky - .5) < .01 Then
      Print Format$(WN, "#.00")
      Print Format$(WN,"#")
     End If
    End Function
    Algorithms - interesting mathematical techniques with program code included.