Two answers in one.
First answer :
I would like here to consider the question under a different angle : how can "Sage" be used for the computation of the winding number of any point on a grid with respect to a self-crossing polygon.
This can be done using the so-called "Sunday's algorithm" well described here , being understood that a representation such as Fig. 1 can be easily converted into a colored 'pixelized" representation.
Fig. 1. The principle of Sunday's algorithm can be explained on the example of the red point situated in $(5,5.5)$ : its winding index is $2$ because the half-line issued from it to the right crosses 4 arrows : $3$ of them are upwards, $1$ is downwards ; total $3 \times (+1) + 1 \times (-1) = 2$.
Sage program for the generation of Fig. 1 :
L = 10 # plotting region is [0,L]^2 n = 12 # number of vertices p = [vector([L*random(),L*random()]) for _ in range(n)] # polygon's vertices G = Graphics() for i in range(n): G+=arrow(p[i], p[(i+1)%n]) r = srange(0, L, 0.25) for x0 in r: for y0 in r: wn = 0 # winding number for i in range(n): x1, y1 = p[i] x2, y2 = p[(i+1)%n] if (y0-y1)*(y0-y2)<0: xi = x1 + (y0-y1)*(x2-x1)/(y2-y1) # intersection abscissa if xi > x0: wn += sign(y2-y1) G += text(str(wn),(x0,y0),fontsize=8) G.show(aspect_ratio=1)
Remark : the central part of the previous algorithm could be replaced by
if (y0-y1)*(y0-y2)<0: a1=matrix([[x1,x2,L],[y1,y2,y0],[1,1,1]]).determinant() a2=matrix([[x1,x2,x0],[y1,y2,y0],[1,1,1]]).determinant() if a1*a2<0: wn+=sign(y2-y1)
Indeed, Let us call $P_k=(x_k,y_k), \ k=0,1,2$ and $P_3=(x_0,L)$ (the point at the extremity of the horizontal half-line $(\ell)$ issued on the right of $P_0$), a way to express that line $P_1P_2$ intersects $(\ell)$ is to consider the third barycentric coordinate of $P_0$ with respect to triangle $P_1P_2P_3$. This barycentric coordinate whose value is equal to the ratio $\tfrac{\text{area}(P_1P_2P_0)}{\text{area}(P_1P_2P_3)}=\tfrac{a_2}{a_1}$ must be negative. Equivalently, the product $a_1a_2$ must be negative. Using this modification of the program is important because it avoids in particular the division by $y_2-y_1$ which can be $0$.
Second answer
(Closer to your objective but still not systematic)
Fig. 2.
The winding number of a curve $\gamma(t)$ around a point $z_0=x_0+iy_0$, in complex function theory, is obtained as the integral
$$W_{z_0}=\frac{1}{2i\pi}\int_0^{2 \pi} \frac{\gamma(t)'}{\gamma(t) - z_0}dt\tag{1}$$
(where the prime notation means "derivative").
[Integral (1) comes from the contour integral of $\frac{1}{z-z_0}dz$ along circuit $\gamma$ whose primitive function is a branch of $\log(z-z_0)$, explaining the connection with "angle sweeping" ; this remark assumes a certain acquaintance with complex function theory, theorem of residues, etc.]
The following "Sage" program is based on (1). It assumes that the different regions delimited by the curve have been previously identified "by hand".
(questions about this program are welcome).
var('x y t') # 2 equivalent ways to define the lemniscate : fx(t)=cos(t);fy(t)=sin(2*t); # (parametric way) f(x,y)=4x^4-4x^2+y^2 # equal to 0 (implicit way) # winding number definition around z0=x0+Iy0: def wn(x0,y0): gamma(t) = (fx(t)-x0)+I*(fy(t)-y0) g(t) = (diff(gamma(t), t)/gamma(t)).imag() R = numerical_integral(lambda t: g(t).n(),0,2*pi)[0] R=(R/(2*pi)).n() return round(R,0) G=Graphics() for k in range(10): a=2*random()-1;b=2*random()-1; w=wn(a,b) if w==1: G+=region_plot([f(x,y)<0,x>0],(x,-1,1),(y,-1,1),incol="lightgreen",plot_points=[100,100]) if w==-1: G+=region_plot([f(x,y)<0,x<0],(x,-1,1),(y,-1,1),incol="lightblue",plot_points=[100,100]) G+=parametric_plot((fx(t),fy(t)),(t,0,2*pi),color="black") show(G)
