In another thread Stuart illustrated an obscure feature of PB’s

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.

**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.

Code:

' 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