This post explains the key findings and algorithm of our 2023 SIGGRAPH paper “Winding Numbers on Discrete Surfaces” .

What problem does this paper solve?

For a given point on a surface, we would like to know whether it is inside or outside a given curve.

For example, coloring points on this cow according to inside vs. outside would look like this:

Of course, before asking how to determine inside/outside, we should ask whether it’s always possible to define inside/outside. At least in the Euclidean space $\mathbb{R}^2$, the Jordan curve theorem guarantees that every simple planar closed curve separates $\mathbb{R}^2$ into an inside region and an outside region. And a similar statement is also true for “simple” surfaces like the above cow.

But on some surfaces, not all closed curves bound regions:

On the left, a fence on Planet Cow forms a closed curve, and the cows are confined to one region. But the fence on the right also forms a closed curve with no gaps, and yet the cows are free to roam wherever they’d like.

We’ll call curves like the one on the left bounding , because they bound a region, and curves like the one on the right nonbounding because they don’t bound regions. So far, we’ve established that while inside/outside is not well-defined for nonbounding curves, inside/outside is at least well-defined for bounding curves.

However, an algorithm for computing inside/outside is not so simple as classifying loops as “bounding” vs. “nonbounding”. In particular, nonbounding loops can combine to bound valid regions (figure reproduced from Riso et al. 2022, Fig. 4 ):

And more generally, there are numerous ways inside/outside can fail to be well-defined if we’re willing to consider all types of surface domains:

Computing inside/outside gets even harder if some of the curves have been corrupted somehow and are now “broken”, for example the “dashed” curves in the rightmost example above, which represent partial observations of an originally whole curve. Ideally, we would like to answer “inside or outside?” for the point $p$ in all of the examples shown above.

Our starting point to reasoning about all of the above is the winding number , a quantity which appears widely across math and physics. Intuitively, the winding number at a point $p$ measures the infinitesimal angle subtended by the curve at $p$, and sums up all these angles over the curve. For a closed curve in the plane, this gives us a piecewise constant function whose value at each point tells us how many times the curve “winds” around it. Hence winding number is fundamentally a topological quantity (though it is often defined through a geometric formula). One can also consider a generalization of winding number called solid angle , which considers open curves in addition to closed curves.

Researchers in computer graphics and geometry processing have in fact already used the properties of solid angle to great success. Perhaps the most well known examples are the equivalent* methods of Poisson Surface Reconstruction [Kazhdan et al. 2006] and generalized winding number [Jacobson et al. 2013] (*see Chen et al. 2024 for a precise explanation of the equivalence). However, all past methods are only applicable to Euclidean domains like $\mathbb{R}^n$, not general surface domains. For deeper discussion of why usual interpretations of solid angle don’t apply in surface domains, see our supplemental, “Perspectives on Winding Numbers” .

At the time of our work, there had been relatively little mathematical or algorithmic treatment of winding number on surfaces. Riso et al. 2022 had recently developed the algorithm “BoolSurf” for computing booleans on surfaces — but assumed that the input curves are closed, and already partitioned into loops. It remained an open question of how one could optimally infer the curve’s topology from unstructured input. More generally, we aimed for a principled handling of surfaces with boundary, and other types of surfaces.

Concretely, our algorithm does the following: Given an unorganized collection of curves with orientations – meaning each curve has a specified direction from one endpoint to the other, but otherwise with no additional structure – we determine inside vs. outside on general surfaces. Importantly, we do so robustly , meaning the algorithm reliably outputs correct answers even in the presence of noise and errors. In the process, we determine the bounding vs. nonbounding components of the curve, and produce a closed, completed version of the input curve.

In short, we obtain an algorithm that answers the fundamental question of inside/outside. In turn, the algorithm can form a basic building block of many geometry processing tasks ranging from region selection to curve completion to boolean operations on shapes.

An example of using our algorithm to compute boolean operations on regions defined by imperfect, broken curves on surfaces.

How does the paper solve the problem?

In a nutshell : the key idea of our method is to turn the harder problem of determining inside/outside of arbitrary curves , into an easier, but equivalent problem about vector fields .

As mentioned earlier, we start from the established concept of winding number . Ordinary winding number is often presented as an integral: for example, in differential geometry, the winding number $w_\Gamma(x)$ of a curve $\Gamma$ at a point $x\in\mathbb{R}^2$ is usually defined as $w_\Gamma(x)=\frac{1}{2\pi}\int_\Gamma \frac{\langle n_\Gamma(y), y-x\rangle}{|x-y|^2} {\rm d}y$, or equivalently in complex analysis as the complex contour integral $w_\Gamma^\mathbb{C}(p) = \frac{1}{2\pi i}\int_\Gamma \frac{1}{z-p} {\rm d}z$ for $p\in\mathbb{C}$. You could, for example, extend these integral definitions to surfaces that admit planar or spherical parameterizations. But in general, these contour integrals can’t easily be computed on surfaces. There are several other characterizations of winding numbers (see “Perspectives on Winding Numbers” ), but they also fail to generalize on surfaces.

One insight is that the winding number is a special type of harmonic function . A harmonic function is one whose Laplacian vanishes; you can just think of it as a special type of smooth function. In particular, winding number is a jump harmonic function (for more detail, see the section below, “The algorithm in more detail”).

Jump harmonic functions can be computed using the jump Laplace equation , which is the partial differential equation (PDE) version of solid angle , a generalization of winding number to open curves $\Gamma$; from this point on, we’ll use the terms “winding number” and “solid angle” interchangeably. Simply solving the jump Laplace equation, however, doesn’t yield a winding number function on surfaces, because of the influence of nonbounding loops. In particular, we get spurious discontinuities that do not correspond to anything in the input:

Formally, (non)bounding loops are loops that are (non-)nullhomologous , meaning congruent to zero in the first homology group $H_1(M) = \rm{ker}(\partial_1)\setminus\rm{im}(\partial_2)$ of the surface $M$. In words, nonbounding loops are closed curves which don’t bound regions.

Great, we’ve identified the core challenge of our problem! We’ve been able to formalize the meaning of “nonbounding” through the lens of homology: in particular, we just need to determine the homology class of the curve (meaning the nonbounding components of the curve which belong to the first homology group of $M$). But homology only directly applies to closed curves. So we use de Rham cohomology to instead process functions dual to curves.

At a high level, we start by taking an appropriate derivative of our jump harmonic function. Once in the “gradient domain”, our problem of decomposing a curve into its bounding and nonbounding components — meaning components that can and cannot be explained as region boundaries — becomes the more well-studied problem of decomposing a vector field into components that can and cannot be explained by a potential, which we can solve using the classical tool of Hodge decomposition .

This is a key insight of our method: Given some unstructured collection of broken curves, it’s really hard to reason about them directly. So we instead transform the curves into vector fields , which are objects we do know how to reason about – and in particular, know how to decompose using Hodge decomposition. Once we have the decomposition of the vector field, we transform it into a decomposition of the curve.

Further, because we work with smooth functions defined globally on the domain, our algorithm is much more robust than if we had tried to work directly with sparse, singular curves. In summary, the final algorithm is able to robustly compute inside/outside by determining the curve’s nonbounding components, in particular by reasoning about vector fields associated with the curve rather than the curve itself. The algorithm is guaranteed to work if the input is “perfect” (no noise, gaps, etc.), and otherwise degrades gracefully in the presence of imperfections (Figures 12–15, 20–26 in the paper.)

The algorithm in more detail

We use $M$ to denote a surface domain, and $\Gamma$ some collection of oriented curves on $M$. As input, the algorithm takes in $M$ and $\Gamma$. As output, the algorithm outputs

• an integer labeling of regions bounded by $\Gamma$
• a decomposition of $\Gamma$ into bounding components that induce valid regions vs. nonbounding components
• a closed, completed version of the input curve $\Gamma$.

As mentioned above, winding number is a jump harmonic function , meaning a function that is harmonic everywhere — including at points on the curve — modulo the jumps. (The function only fails to be harmonic at the curve’s endpoints where there is a branch point singularity, shown below.)

Hence our starting point to winding numbers on surfaces is the jump Laplace equation defined on a domain $M$, $$\begin{align*} \Delta u(x) &= 0 &x\in M\setminus\Gamma \newline u^+(x) - u^-(x) &= 1^* &x\in \Gamma \newline \frac{\partial u^+}{\partial n}(x) &= \frac{\partial u^-}{\partial n}(x) &x\in \Gamma \end{align*}$$ *More generally, $u^+(x) - u^-(x) =$ number of times $\Gamma$ covers $x$.

Basically, the jump Laplace equation corresponds to the standard Laplace equation $\Delta u=0$ (which defines a standard harmonic function), plus a condition that stipulates the solution must jump by a specified integer across the curve $\Gamma$. Finally, the third equation is a compatibility condition that ensures the solution is indeed harmonic “up to jumps” everywhere on $M\setminus\partial\Gamma$. Intuitively, for all $x\notin\Gamma$ the function $u(x)$ is harmonic in the usual sense, and for $x\in\Gamma$ (where the jump occurs) it’s still possible to locally shift $u$ such that it once again becomes harmonic in the usual sense (above, right).

Assuming $M$ is simply-connected , meaning there cannot exist nonbounding loops on $M$, then we can simply solve the jump Laplace equation above for the function $u(x)$. If $\Gamma$ is closed, this case is analogous to ordinary winding number in the plane, which yields a piecewise constant, integer-valued function; every component of $\Gamma$ is a bounding component.

If $\Gamma$ is not closed (but $M$ is still simply-connected), then $u(x)$ yields a “soft” real-valued (rather than integer-valued) indicator function, analogous to solid angle.

Solving the jump Laplace equation above is Step $\raisebox{.5pt}{\textcircled{\raisebox{-.9pt}{1}}}$. If $M$ is not simply-connected, we need to take additional steps. The key insight here is that jump harmonic functions form the “bridge” through which we can transform curves into functions, and vice versa. In particular, jump harmonic functions can both encode curves (where the function jumps in value) as well as a curve’s homology class (through its derivative). Our algorithm will amount to a round-trip around the following diagram, where we first translate the input curve $\Gamma$ into a jump harmonic function, then into a differential 1-form , which you can think of as a vector field. We’ll then translate the resulting 1-form back into a jump harmonic function and curve, which will correspond to our final winding number function and curve decomposition, respectively. Performing these translations amounts to using appropriate notions of differentiation and integration.

At a high level, these notions of differentiation and integration simply correspond to their usual notions from ordinary calculus, except we have to take extra care when dealing with the discontinuities of jump harmonic functions. For simplicity, we’ll first consider a periodic 1-dimensional function $f(x)$ defined on the unit interval $[0,1]$:

In general, the gradient of such a function $f(x)$ with respect to $x$ yields a continuous part $\omega:=\mathcal{D}f$, obtained from considering each continuous interval of $f(x)$, as well as a “jump” part $\mathcal{J}f$ encoding the location and magnitude of the jumps in $f(x)$, expressed here as a sum of Dirac-deltas.

For a jump harmonic function in particular, the continuous part $\omega:=\mathcal{D}f$, which we term the Darboux derivative , “forgets” jumps:

(If you’re familiar with differential geometry, you can think of $\mathcal{D}f$ as “$df$ modulo jumps”, where $d$ denotes the standard exterior derivative. Just as the standard exterior derivative $df$ measures the failure of $f$ to be constant, $\mathcal{D}f$ measures the failure of $f$ to be piecewise constant.) In essence, even though jump harmonic functions have discontinuities, the condition of harmonicity forces the Darboux derivative of a jump harmonic function to be globally continuous because the derivatives must “match up” across the discontinuity.

Because the Darboux derivative of a jump harmonic function “forgets” the jumps, there is no unique inverse to Darboux differentiation. In particular, we can only integrate “up to jumps”: to continue the previous example, both the original jump harmonic function $f(x)$ and the function $g(x)$ shown below are valid possibilities of integrating $\omega$:

The moral of this story is that to map from derivatives back to curves, we can integrate $\omega$ and choose the jumps. (This choice of jumps is what will allow us to decompose curves into bounding vs. nonbounding components.)

Our understanding of the 1D case more or less extends directly to the 2D case. Instead of a 1D function on a periodic interval, we can, for example, consider a solution $u(x)$ of the jump Laplace equation on a 2D periodic domain — namely, a torus:

We arrive at Step $\raisebox{.5pt}{\textcircled{\raisebox{-.9pt}{2}}}$: taking the Darboux derivative of the jump harmonic function obtained in Step $\raisebox{.5pt}{\textcircled{\raisebox{-.9pt}{1}}}$. The Darboux derivative $\omega$ of $u$ is now a 1-form, which can be thought of as a vector field (shown above). If $\omega=0$, then $u(x)$ is piecewise constant, meaning that $u(x)$ is already a valid (piecewise constant) region labeling. Conversely, if $w\neq 0$, then there are nonbounding components of $\Gamma$. In particular, nonbounding components of $\Gamma$ are encoded by the harmonic component $\gamma$ of the 1-form $\omega$. More formally, (non)bounding components of $\Gamma$ correspond to 1-forms (non)congruent to zero in the first cohomology group $H^1(M)=\rm{ker}(d_1)\setminus\rm{im}(d_0)$.

Luckily for us, extracting the harmonic component of 1-form can be done via Hodge decomposition , which states that our 1-form $\omega$ can be decomposed as $\omega=d\alpha+\delta\beta+\gamma$, where $\gamma$ is the harmonic 1-form we seek to isolate (and the differential forms $\alpha$ and $\beta$ are some scalar and vector potentials). So Step $\raisebox{.5pt}{\textcircled{\raisebox{-.9pt}{3}}}$ amounts to doing Hodge decomposition to extract the harmonic 1-form $\gamma$.

In summary, we have transformed our original problem of decomposing the input curve into components that can and cannot be explained as region boundaries, into the well-studied problem of decomposing a vector field into components that can and cannot be explained by a potential.

To recap, here is where we are in our round-trip journey so far, using as a visual example a collection of (broken) bounding and nonbounding loops $\Gamma$ on a torus:

We’ve used differentiation to obtain a 1-form associated with the original curve $\Gamma$, then use Hodge decomposition to process this 1-form into a harmonic 1-form $\gamma$ encoding the nonbounding component of $\Gamma$. We arrive at Step $\raisebox{.5pt}{\textcircled{\raisebox{-.9pt}{4}}}$: To determine our curve decomposition, we just have to use integration to translate $\gamma$ back into a jump harmonic function, and then a curve that will represent the nonbounding component of $\Gamma$.

More specifically, we search for a scalar potential $v$ that could have generated $\gamma$, meaning $\mathcal{D}v=\gamma$. Since $\gamma$ is harmonic, $v$ must jump somewhere; but remember, we get to choose these jumps — and we want the jump derivative $g:=\mathcal{J}v$ to represent the nonbounding component of $\Gamma$. We formulate the following optimization that does just this:

In words, the first constraint $\mathcal{D}v=\gamma$ stipulates that the solution $v$ should integrate $\gamma$ in the “Darboux” sense, while the second term in the objective function encourages the jumps $g=\mathcal{J}v$ in $v$ to concentrate on a subset of the original curve $\Gamma$. In the case that $\Gamma$ is a broken curve, the first term in the objective function also encourages $g=\mathcal{J}v$ (which necessarily must form closed loops, since $\gamma$ is harmonic and thus not globally integrable) to be a shortest closed-curve completion of the nonbounding components of $\Gamma$. Finally, the second constraint discourages spurious extraneous loops in the case that $\Gamma$ contains multiple nonbounding components of the same homology class.

The solution of this optimization yields a so-called “residual function” $v$ such that the solution $w$ to a jump Laplace equation using the curve $\Gamma - \mathcal{J}v$ is our final winding number function. In other words, we solve for the function $w$ for which we have removed the influence of $\Gamma$’s nonbounding components.

Step $\raisebox{.5pt}{\textcircled{\raisebox{-.9pt}{5}}}$ is simply reading off the curve decomposition from the residual function $v$. The locus of points $\mathcal{J}v$ where $v$ jumps represents a completion $G$ of the nonbounding components of the input curve $\Gamma$; the bounding components of $\Gamma$ are simply the complement of the nonbounding components, taken as a subset of $\Gamma$. Additionally, the function $w$ can be rounded to yield integer region labels $W$. The entire algorithm is summarized in the following diagram:

With our algorithm described, all that remains is to discretize the algorithm.

In the discrete setting, we assume the surface domain $M$ is specified as a triangle mesh. We represent curves and regions as chains on $M$, meaning as signed integer-value functions on mesh edges and mesh faces, respectively. The jump Laplace equation can be discretized on the mesh as a sparse, positive definite linear system. Performing Hodge decomposition on the gradient of the solution $u$ amounts to solving another Laplace-type sparse linear system. Finally, performing the integration procedure necessary to map the decomposition of our 1-form to a decomposition of our curve amounts to solving sparse linear program. (The last two steps are only necessary if $M$ is topologically nontrivial.) For more details, check out our paper , which also contains detailed pseudocode, or our code repository .

Takeaway

In summary, in this paper we addressed the problem of computing inside/outside of curves on surfaces, including the cases where curves fundamentally might not even have a well-defined inside or outside. We were able to do this by using a duality between curves and vector fields, using insights from differential geometry and algebraic topology.