This notebook demonstrates computation of winding numbers in two dimensions based on Jacobson et al.'s Robust Inside-Outside Segmentation using Generalized Winding Numbers as well as the subsequent optimization of Barill et al. in Fast Winding Numbers for Soups and Clouds.
The winding number of a curve relative to a point is a straightforward concept. We trace the curve, add up the angle swept out relative to the point—positive for counterclockwise motion and negative for clockwise—and divide by . For any point not on the curve itself, this procedure counts count how many times the curve encircles the point, for each counterclockwise encircling, for each clockwise encircling, and for points entirely outside the curve.
Stating this procedure mathematically, we write the winding number relative to point as the integral
For a curve comprised of discrete line segments, the integral becomes the summation
Jacobson et al. offer a formula for computing the angle swept out relative to point by the line segment with endpoints and . Defining and , the angle swept out by the segment is
The curve does not need to be closed to use the above equation. The figure below plots the winding number induced by a single segment of length . Color represents this generalized winding number. (The contours are scaled by to better represent the shape of the field.) The winding number is exactly on one side of the segment, on the other, and varying shades between as we move around the segment from one side to the other.
The summation above represents a perfectly adequate means of computing a winding number. We may superimpose the above segment solutions, allowing us to iterate over the segments which make up a curve and directly compute the total winding number relative to a given point.
The plot below illustrates this summation for a circle made up of 32 segments. Note that when the curve is fully closed, the winding number becomes exactly 0 or 1 everywhere.
This procedure works great for toy examples, but if we want to visualize an entire field, the global nature of this algorithm presents a major problem: every point depends on all geometry! Imagine we have geometry with 100,000 line segments. A 512 × 512 image contains about 260,000 pixels. Rasterizing such an image requires 26 billion evaluations of the above equation!
If you have a sharp eye, you might realize that the field in the above example seems to depend only on the endpoints of the curve. More specifically, for any point outside the convex hull of the geometry, knowledge of the endpoints alone is adequate to compute the field exactly. Jacobson et al. use this fact to construct an efficient scheme for exact far-field evaluation, but as I’m particularly interested in a soup of disorganized line segments without clean, well-maintained endpoints, the approximate algorithm of Barill et al. seemed preferable.
The key observation is that a line segment from far away looks like a point. More specifically, the effect of a source or collection of sources approaches in the far field that of a dipole defined solely by a location, an orientation, and a strength. You can think of a dipole like the limiting behavior of two oppositely charged particles brought infinitesimally close together.
Adjust the scale of the geometry below and observe that from far away, the structure disappears and the winding number looks identical to that of the single short line segment shown above.
This suggests a hierarchical approximation: nearby, we iterate over all geometry explicitly; far away, we treat the geometry in aggregate. This method of iterating over a tree of clustered geometry is called the Barnes-Hut approximation. Barnes-Hut can be a somewhat harsh approximation, but with respect to the number of segments , its performance for a single evaluation is such a speedup compared to the naive that we consider ourselves lucky to get so far with relatively little effort.
(The Fast Multipole Method (FMM) performs a more careful and accurate expansion, though it’s a lot more complicated and I’ve never tried it. Barill et al. mention that it didn’t prove advantageous for this particular problem.)
For a deeper dive into the Barnes-Hut approximation, see Jeffrey Heer’s excellent notebook The Barnes-Hut Approximation.
Since we’re dealing with oriented dipole-like segments, our case is a bit different than the more typical gravitational simulation. To that end, there’s one more improvement worth discussing, and at this point we have no choice but to break out the math. As Barill et al. discuss in Fast Winding Numbers for Soups and Clouds, instead of aggregating geometry into single point dipoles, we can perform a Taylor series expansion to account for internal structure of each point cluster. They focus primarily on the 3D case. I found the math a bit tedious and challenging to work through, so I’ll offer some details about the 2D case here.
Finally, the main feature! The simulation below computes the winding number of a curve using a Bounding Volume Hierarchy (BVH) for acceleration. On construction, the segments which make up a leaf node are aggregated into a single dipole. For non-leaf nodes, it combines the two dipoles of the respective child nodes into a single dipole. During evaluation, far-away nodes are accounted for only in the aggregate sense, as a dipole. Nearby nodes are either then recursed into or in the case of a leaf node, their segments treated individually.
In the end, I’m still not completely content with the Taylor series expansion and small improvements. I think this is because I’ve sort of combined the Taylor series expansion for a point-dipole cloud with the Green’s function of a finite-length segment. The resulting improvements sorta help, but not completely. Fortunately, the end result I seek—a thresholded winding number for loosely-but-not-exactly-watertight geometry—appears quite tolerant to error, such that the unimproved solution seems good enough for a pretty good inside-outside query.