Technical Report

Winder

A Library for Fully Differentiable Fast Winding Number Fields on GPU

⚠️
Early Alpha • Active Development Only the forward pass is implemented so far. The backward pass is currently under active research and engineering—please do not expect a finalized or production-stable product yet.
By Raphael Braun
🛠️ Point Clouds & Triangle Soups Millisecond BVH8 Generation 🧠 Second-Order Taylor Accuracy 🔄 Full Partial Derivatives 💾 Optimized for Speed and Flexibility

01.Winding Numbers: Introduction & Motivation

🧭 The Inside-Outside Dilemma: A Story of Broken Geometry

Imagine you are slicing a crappy 3D model for a 3D printer, converting a complex video game asset into a signed distance field, or you are training a neural network to understand 3D shapes from real-world 3D scans. All those cases, are so much easier to handle if the 3D representation has a property that every object in the physical world possesses: a closed surface. Even a piece of paper has a front side, a thickness, and a back side.

Most meshes - even hand-crafted ones - are built to look good when rendered into images, not to be topologically correct. Native 3D scans consist only of disconnected point clouds, which are processed into more or less closed surfaces. The simple question: "Is a given point inside the object, or outside in empty air?" can be exceptionally hard to answer reliably.

If your 3D model is a perfect, watertight (manifold) mesh, this is easy. You cast a ray from your query point; if it crosses the mesh boundary an odd number of times, you're inside. But perfect 3D models are rare. Real-world 3D data—whether harvested from raw LiDAR scans, chaotic CAD conversions, or raw triangle meshes—is inherently broken. It is full of cracks, missing patches, self-intersections, and floating geometry.

Classical ray-casting can fail catastrophically with any microscopic topological flaw, completely breaking downstream processing pipelines.

This is where the seminal work of Robust Inside-Outside Segmentation using Generalized Winding Numbers (Jacobson et al., 2013) comes in. Instead of a brittle binary test, it treats your 3D shape as a source for a continuous field, similar to how mass acts as a source for gravitational potential fields in physics. Generalized Winding Numbers aggregate the contribution - or "vote" - of every single surface element into a smooth, continuous field. Instead of a rigid binary answer, we get a continuous scalar factor tracking how many times the 3D shape wraps around the query point.

❌ Classical Boundaries

Asks: "Did I cross a line?" Fragile, binary, and blind to broken geometry. A single missing triangle breaks the global spatial state.

🔮 The Winding Number View

Asks: "How much of the object wraps around me?" Global, continuous, and robust to defects. If 99% of a shape surrounds you, you are 99% inside.

For a perfect, watertight surface, a point inside has a winding number of $1$, indicating it is completely enclosed by geometry, while a point outside has a winding number of $0$. Remarkably, this formulation remains completely stable when confronting holes, intersecting surfaces, or severe topological defects. For instance, due to a missing patch or crack, an interior query point might capture an $80\%$ vote for being inside and a $20\%$ vote for being outside, yielding a graceful winding number of $0.8$. Conversely, an exterior point close to that same hole might register a winding number of $0.2$ by getting positive votes for surface elements visible through the opening. On highly corrupted meshes featuring self-intersections or duplicate faces, the field accommodates this by naturally producing values going below $0$ or above $1$.

Crucially, away from the boundary geometry, the generalized winding number field is strictly harmonic, meaning it satisfies Laplace's equation ($\nabla^2 \mathcal{W} = 0$). In physical terms, this means the field behaves like steady-state heat distribution or a gravitational potential: it is as smooth as mathematically possible, completely devoid of sudden local spikes, maxima, or minima in empty space. It acts as an optimal global interpolator that bridges topological gaps automatically.

This mathematical elegance elevates the winding field into a powerful, implicit 3D shape representation. Because the field equation can be queried at any arbitrary coordinate $x \in \mathbb{R}^3$, it possesses infinite spatial resolution, entirely unfettered by voxel grids or discrete memory arrays until evaluation. The true underlying surface of the shape is elegantly preserved as an iso-surface (a constant-value contour line) slicing through this continuous volume field, allowing downstream reconstruction algorithms to handle broken geometry.

Winder fully unlocks the beauty of this robust 3D representation for large-scale GPU optimizations. By adapting the Fast Winding Number algorithm for the GPU and extending it to partial derivatives, Winder provides fully differentiable, millisecond-level winding fields on raw, broken datasets at scale. Here is how the math actually unrolls.

Figure 1.1: Generalized Winding Numbers Demonstration. The video above shows a 2.1M 3D printing asset through three geometric states: clean, with high-frequency Gaussian noise, and open non-manifold topology (holes). Winder was used to computes the winding number field for all three configurations over a $512^3$ query grid in under 1 second each. The scalar winding number field goes from 0 meaning outside (blue) to 1 meaning inside (green). Colors in between show regions that are not 100% surrounded by geometry. We show a marching cubes reconstruction with varying thresholds on the $512^3$ winding number grid to demonstrate the smooth iso-surfaces inherent to the winding field representation.

The Continuous Formulation

For any continuous, closed manifold surface $\mathcal{M}$, the continuous winding field maps any spatial location to a continuous scalar field $\mathbb{R}^3 \rightarrow \mathbb{R}$, calculating precisely $1$ for any point enclosed by the manifold and $0$ outside it:

$$\mathcal{W}(x) = \frac{1}{4\pi} \iint_{\mathcal{M}} \frac{(p - x) \cdot n(p)}{\|p - x\|^3} \, dA(p)$$

Where $x$ is the 3D query point, $p$ is an integration coordinate lying on the continuous surface geometry $\mathcal{M}$, $n(p)$ denotes the outward-facing unit normal vector at position $p$, and $dA(p)$ represents the infinitesimal area element.

Discretization: Triangle Soups

When dealing with digital geometry, this continuous formulation is discretized. For standard Triangle Soups (any set of triangles, connected or not), the total winding number at query $x$ is evaluated as the sum of signed solid angles $\Omega_t$ subtended by each individual triangle $t$ within the triangle collection $\mathcal{T}$:

$$\begin{align}\mathcal{W}_{\text{mesh}}(x) &= \frac{1}{4\pi} \sum_{t \in \mathcal{T}} \Omega_t(x)\end{align}$$

To guarantee numerical stability near triangle edges, vertices, or coplanar configurations, $\Omega_t(x)$ is computed robustly using the reliable signed solid angle formulation established by Van Oosterom and Strackee (1983). For a triangle with vertex positions defined relative to the query point as $a = v_1 - x$, $b = v_2 - x$, and $c = v_3 - x$, the solid angle is computed via:

$$\tan\left(\frac{\Omega_t(x)}{2}\right) = \frac{a \cdot (b \times c)}{\|a\|\|b\|\|c\| + (a \cdot b)\|c\| + (b \cdot c)\|a\| + (c \cdot a)\|b\|}$$

With this we can compute the winding field for any triangle mesh. As an example we show the winding number field of a single triangle in the following video.

Winding Number Field Representations (Single Triangle)
2D Sweeping Planar Slices
3D Volumetric Orbital Render
Figure 2.1: Visualization of the winding number field of a single triangle. Green tones are positive values, red tones negative. Left: Progressive cross-sectional scans. Right: Path-traced 3D volumetric field. Note how the field smoothly transitions and diminishes with distance.

Discretization: Point Clouds

For Point Clouds lacking explicit connectivity, the surface is treated as a distribution of oriented points. The continuous integration reduces to a discrete sum:

$$\mathcal{W}_{\text{points}}(x) = \frac{1}{4\pi} \sum_{i} \frac{(p_i - x) \cdot n_i}{\|p_i - x\|^3} a_i$$

Where $p_i$ represents the 3D spatial position of point sample $i$, $n_i$ represents its orientation unit normal vector, and $a_i$ denotes the Voronoi surface area associated with that point. This is sometimes called point dipole method as the winding field of the point samples is mathematically equivalent to the electrostatic potential of electric dipoles. The area $a_i$ maps to the magnitude of the electric dipole moment.

With this formula we now can compute the winding field for any oriented point cloud with known surface areas per point. As example we show a single dipole in the following video.

Winding Number Field Representations (Single Dipole)
2D Sweeping Planar Slices
3D Volumetric Orbital Render
Figure 2.2: Visualization of the winding number field of a single dipole. Green tones are positive values, red tones negative. Left: Progressive cross-sectional scans. Right: Path-traced 3D volumetric field.

02. Hierarchical Acceleration & Taylor Expansions

1. The Computational Bottleneck

Evaluating the exact discrete winding number formula requires computing the signed solid angle or dipole potential between every single query location and every surface primitive. This means computing the winding numbers for $M$ query points and $N$ geometric source primitives has a computational complexity of $\mathcal{O}(N M)$. With $M$ and $N$ goint into the millions for real world use cases this computational complexity turns into a serious bottleneck in any pipeline.

To mitigate this, Jacobson et al. (2013) already suggested using a bounding volume hierarchy (BVH) to group together manifold sub parts of a mesh. However, this method's algorithmic performance is highly dependent on clean, continuous mesh connectivity. When confronted with unoriented triangle soups, disconnected components, or broken geometry boundaries, their traversal can break down under worst-case inputs, degrading directly back to a brute-force $\mathcal{O}(NM)$ summation loop. This limitation makes large-scale processing of raw, non-manifold assets a gamble.

2. The Fast Winding Numbers Approximation

Barill et al. (2018) introduced the Fast Winding Numbers framework which solves the complexity issue in an elegant way. This method utilizes the fact that the influence of a geometry element on the winding number quickly decays with distance - similar to gravity. Just like large scale gravity simulations their algorithm is a tree based Barnes-Hut algorithm and builds a spatial hierarchy (such as an octree) to group geometry and their influence on the winding number field into hierarchical clusters.

Instead of computing contribution from all geometry elements for all queries, the collective contribution of distant clusters is evaluated from their center of mass $\tilde{p}$. Fast Winding Numbers approximates the far field clusters with a truncated Taylor series expansion of the Laplace Green's function derivatives, evaluated with respect to node's center of mass $\tilde{p}$.

Point Cloud Expansion Formulation

For a node containing a cluster of $K$ area-weighted point dipoles, Barill et al. (2018) approximate the collective winding number $\tilde{w}$ at an evaluation query coordinate $q$ via the following expansion:

$$\begin{align*}\tilde{w}(q) \approx & \left( \sum_{i=1}^{K} a_i \hat{n}_i \right) \cdot \nabla G(q, \tilde{p}) + \\ & \left( \sum_{i=1}^{K} a_i (p_i - \tilde{p}) \otimes \hat{n}_i \right) \cdot \nabla^2 G(q, \tilde{p}) +\\ & \frac{1}{2} \left( \sum_{i=1}^{K} a_i (p_i - \tilde{p}) \otimes (p_i - \tilde{p}) \otimes \hat{n}_i \right) \cdot \nabla^3 G(q, \tilde{p})\end{align*}$$

Where the operations $\otimes$ and $\cdot$ denote the tensor outer product and inner tensor contraction respectively. Here, $q$ is the 3D query point, $\tilde{p}$ is the precomputed center of mass of the node holding the cluster, $p_i$ is the position of point sample $i$, $\hat{n}_i$ is its oriented unit normal, and $a_i$ is its localized area weight. The terms $\nabla G$, $\nabla^2 G$, and $\nabla^3 G$ represent the gradient (a vector), the Hessian (a symmetric two-tensor), and the third-derivative (a three-tensor) of the standard fundamental Green's function for the Laplace equation, evaluated between the query point and the center of mass. Their formulation is given in the appendix of Barill et al. (2018).

The structural coefficients grouped inside the parentheses are entirely independent of the query position $q$, meaning they can be precomputed for each cluster during tree construction and cached inside each node to be used for all queries.

Triangle Soup Expansion Formulation

This Taylor series framework scales naturally to raw Triangle Soups by substitute the discrete point sums inside the expansion coefficients with continuous surface integrals evaluated across the flat area elements of each triangle $t$ inside the node:

$$\begin{align*}\tilde{w}(q) \approx & \left( \sum_{t=1}^{K} \int_{t} \hat{n}_t \, dA \right) \cdot \nabla G(q, \tilde{p}) +\\ & \left( \sum_{t=1}^{K} \int_{t} (x - \tilde{p}) \otimes \hat{n}_t \, dA \right) \cdot \nabla^2 G(q, \tilde{p}) +\\ & \frac{1}{2} \left( \sum_{t=1}^{K} \int_{t} (x - \tilde{p}) \otimes (x - \tilde{p}) \otimes \hat{n}_t \, dA \right) \cdot \nabla^3 G(q, \tilde{p})\end{align*}$$

Where $\hat{n}_t$ represents the surface normal of triangle $t$, $x$ represents the continuous integration coordinate tracking across the triangle's surface area element $dA$, and $\tilde{p}$ remains the fixed center of mass of the node. The first coefficient maps out the total area-weighted normal across the grouped geometry, projecting the intuition that a distant cluster of triangles can be approximated as a single, representative larger primitive. Crucially, all three of these continuous spatial integrals possess exact, numerically stable closed-form expressions for triangles which are also given in the appendix of Barill et al. (2018).

"Far Field" definition

To compute the winding number of a set of queries the constructed tree is traversed by all queries. The entire contribution of branches that are identified as "Far Field" are computed using the tailor expansion shown above with the coefficients stored in the current node. To determine that for a given query $q$ a node is in the far field, each node stores its center of mass $\tilde{p}$ and the maximum distance $\tilde{r}$ of any contained geometry element to that center of mass. The sub-tree represented by the node is considered "Far Field" w.r.t. $q$ if:

$$\|q-\tilde{r}\|_2 < \beta\cdot\tilde{r}$$

$\beta$ is a hyper parameter that can be set by the user to tune how aggressive the algorithm approximates sub-trees. Larger $\beta$ lead to preciser, but slower results.

Complexity

The expansion can theoretically be extended to higher orders indefinitely. However Barill eg al. found that, in practice, expanding beyond the second order ($p=2$) yields diminishing accuracy returns relative to the computational overhead.

Evaluating a $p$-th order Taylor expansion adds a localized tensor processing complexity scaling at $\mathcal{O}(3^p)$ operations per evaluated node. Assuming a well-formed tree structure, the total amortized complexity encompassing both the structural precomputations and far-field evaluations across $M$ queries on a dataset of $N$ geometric primitives drops to:

$$\mathcal{O}(3^p N \log N + 3^p M \log N)$$

This log-linear complexity allows the Fast Winding Numbers formulation to maintain high processing speeds when evaluating massive meshes, unoriented point clouds, and raw geometric soups containing millions of elements. Winder uses this elegant formulation and ports it to GPU with massively parallel tree construction and traversal algorithms.

03. Differentiability: Motivation & Bottlenecks in Field-Based Optimization

1. Global Gradient Propagation

The core motivation for constructing a differentiable winding number field stems from its mathematical formulation as a global summation. Because the winding potential $\mathcal{W}(x)$ at any point in space integrates the collective influence of all primitives, every single geometry element contributes to the continuous field.

Unlike Signed Distance Fields (SDFs) or discrete voxel grids—which suffer from localized clipping, sharp gradient discontinuities, or zero-gradient dead zones—the winding number landscape provides an uninterrupted, infinitely smooth gradient flow. This is a direct consequence of its harmonic nature away from boundaries: Because the field satisfies $\nabla^2 \mathcal{W} = 0$, it is structurally barred from forming localized coordinate traps, local maxima, or local minima in empty space.

Consequently, an error or loss metric evaluated at any arbitrary query position $x$ can be backpropagated directly to every point-normal or triangle vertex in the scene. This unlocks long-range spatial signals, allowing optimization routines to smoothly steer, deform, and correct geometry from anywhere in the field.

2. Analysis of Existing Frameworks

While several contemporary works leverage this global differentiability, they remain constrained by severe architectural bottlenecks, lack of structural scaling, or specialized application scenarios:

  • Winding Clearness: The framework introduced by Xiao et al. uses the field's global boundary behavior to denoise point clouds. However, its implementation relies on a dense linear least-squares solver (torch.lsqr) to optimize the area-scaled normals. They formulate a dense linear system $Ax=b$ with $A$ having shape $3N\times 3N$, which is solved once in every single optimization step to derive their loss "winding clearness". The used LU decomposition has a time complexity of $\mathcal{O}((3N)^3)$ and memory complexity of $\mathcal{O}((3N)^2)$. This complexity limits the usefulness of the approach to point clouds with only a couple of thousand points.
  • Fast Dipole Sums: This method by Chen et al. successfully applies differentiable fast winding fields to an inverse rendering pipeline, introducing a spatial regularization kernel to smooth out proximity singularities near points. However, the system is limited to its specific neural radiance field use case. It completely lacks gradients for point positions (relying on fixed Structure-from-Motion initializations), provides no support for triangle mesh primitives, and suffers from an unparallelized Octree generation routine with very slow initialization times. While using the Fast Winding Numbers algorithm they only use a truncated zero order Taylor approximation for the far field, which reduces precision by an order of magnitude (RMS).

3. A General-Purpose Winding Number Field

Winder decouples the differentiable winding field from any specific downstream task or rigid solver constraint. It is engineered as an application-agnostic, high-performance GPU engine designed to handle any arbitrary optimization routine involving winding number fields at real-world scale.

By supporting full, simultaneous analytical gradients for point positions, area-scaled normals, and triangle vertices, Winder fully exposes the continuous landscape of the winding field to modern gradient-based solvers.

Instead of relying on slow CPU trees or crushing dense matrix solvers, Winder computes these complex derivatives via a highly parallelized, fully differentiable GPU acceleration structure. By lifting classical complexity constraints, the continuous winding field is transformed from a computationally prohibitive luxury into an efficient, first-class mathematical primitive ready for large-scale computer vision, physical simulation, and deep learning pipelines.

04. Winder: Core Vision, Architectural Philosophy & Hardware Co-Design

1. The Construction-First Paradigm

Existing compute pipelines treating fast winding numbers rely on offline-computed spatial structures. For instance the Dipole Sums implementation uses a tight Octree layout. In their method the geometry is assumed to be static, so tree construction only happens once. This construction is allowed to be slow.

Winder was designed with an inverted priority in mind: spatial hierarchy generation must be virtually free. Rather than chasing a perfect static spatial partition offline, winder prioritizes super fast GPU-optimized spatial partitioning trees for dynamic geometry. Everything else - the precision handling, layout design, and mathematical formulation - is a direct architectural consequence built to exploit this lightning-fast, construction-first engine baseline.

Fast Dipole Sums Build (Octree) 566.2 sec Static offline build executed over 2,559,132 points
Winder Build (GPU Parallel BVH8) 0.0132 sec Online build execution over 2,559,132 points
Note: This extreme build time for Fast Dipole Sums is largely an artifact of the single-threaded implementation rather than an inherent topological limitation of Octrees. Indeed, the original Fast Winding Numbers framework demonstrates a significantly more performant tree-building routine. However, it underscores the build-centric focus of Winder.

2. BVH8 and Cooperative Traversal

To achieve this millisecond build speed, Winder utilizes a parallel Bounding Volume Hierarchy (LBVH). The trick is to sort the underlying geometry by the Morton order (Z-order), which can be done with a radix sort (linear complexity). The tree structure can be built on top of the sorted geometry in parallel, where each thread evaluates adjacent primitives to establish node boundaries independently. This avoids expensive top-down spatial partitioning and guarantees a predictable, linear construction time. The resulting binary tree is then converted into a more traversal friendly, wide 8-ary structure (BVH8), where each node has up to 8 children. This conversion also runs in O(N) linear time using a single, persistent cooperative kernel pass. By utilizing fast on-chip shared memory scans to manage thread work-allocation and executing vectorized 32-byte writes, the implementation minimizes global memory overhead while greedily collapsing the binary levels.

The tree is getting even more shallow because every leaf node contains 32 geometry primitives. This leaf size choice is a deliberate hardware play. Because NVIDIA GPU warps consists of exactly 32 execution lanes, a 32-element leaf structure brings some nice advantages:

  • Zero Intra-Warp Divergence: Whole warps cooperatively traverse down the tree and unroll leaf evaluations simultaneously, completely bypassing execution lane stalls.
  • Hardware Registries & Shuffles: Internal evaluation components utilize lightning-fast register-level data shuffling (__shfl_sync), mitigating costly localized VRAM and shared memory round trips.
  • Shared Data: As all warps cooperatively traverse the tree they can also share- and cooperatively load data from VRAM.

The Structural Penalty: Locking leaf sizes to 32 items to maximize build speed forces a trade-off: In regions with sparse geometry, leaves can become oversized and loose, failing to bound the underlying points as tightly as an offline-built Octree would. This can also be seen in the video with the video below with the big green leaf boxes.

The Work Load: Each thread in a warp is responsible for a different query. Coopertively traversing the tree forces threads whose queries are not even interested in certain subtrees to still go along and help other threads, instead of continuing their own work. We empirically found that with morton sorted queries the cooperation still pays off massively.

3. Hardware Co-design

Because the spatial bounding structure is more loose and has large leafs to favor ultra-fast generation, the evaluation pass may face a slightly higher computational workload. Winder accepts this structural penalty, but does everything to make the traversal blisteringly fast through extreme hardware co-design:

  • float16 Far-Field Math: All calculations in the far-field are done in the float16 loating-point format or bfloat16 for the second order Taylor contributions, leading to high computational throughput thanks to Half2 Arithmetic Functions and Bfloat162 Arighmetic Functions respectively on modern CUDA GPUs.
  • Cache Alignment: Super slim 32-byte BVH8 nodes and the 64 byte tailor coefficients, are loaded cooperatively by entire warps. They are perfectly aligned to L2 cache lines to guarantee maximum memory bandwidth utilization.
  • Second-Order Taylor Expansions: Winder uses second-order Taylor approximations in the far field. Note: The precision of winder is still slightly lower compared to the float32 second-order baseline of the original Fast Winding Numbers due to the lower precision of float16 and bfloat16 computations.
Real BVH8 visualization for benchy.obj
Figure 4.1: Visualization of the the wide 8-ary BVH8 (purple), and the leafs (green) we get from the morton order. Note the leafs can be quite large where the geometry is sparse.

4. Partial Derivatives

To compute the partial derivatives of generalized winding numbers for reverse-mode automatic differentiation, we map an incoming scalar adjoint gradient $g_j = \frac{\partial \mathcal{L}}{\partial w(x_j)} \in \mathbb{R}$ at each query position $x_j$ back to the underlying spatial geometry. Depending on the input representation, the analytical vector gradients ($\mathbb{R}^3$) for the three primitive configurations are defined as:

1. Gradient w.r.t. Scaled Normals (Dipole Moments $m_i = a_i \cdot n_i \in \mathbb{R}^3$): $$\frac{\partial \mathcal{L}}{\partial m_i} = \sum_{j=1}^{M} g_j \nabla G(x_j - p_i)$$

2. Gradient w.r.t. Source Point Positions ($p_i \in \mathbb{R}^3$): $$\frac{\partial \mathcal{L}}{\partial p_i} = -\sum_{j=1}^{M} g_j \mathbf{H}_G(x_j - p_i) m_i$$

3. Gradient w.r.t. Triangle Vertices ($v_k \in \mathbb{R}^3$): $$\frac{\partial \mathcal{L}}{\partial v_k} = \sum_{j=1}^{M} g_j \frac{\partial \Omega_T(x_j)}{\partial v_k}$$

Kernel Definitions

Where the spatial displacement vector is defined as $r = x_j - p_i$ with magnitude $r = \|r\|$, the analytical forms of the 3D potential field kernels are given by:

Free-Space Green's Function Gradient ($\nabla G \in \mathbb{R}^3$): $$\nabla G(r) = -\frac{r}{4\pi r^3}$$

Spatial Hessian Matrix ($\mathbf{H}_G \in \mathbb{R}^{3 \times 3}$): $$\mathbf{H}_G(r) = \frac{1}{4\pi r^5} \left( 3r r^T - r^2 \mathbf{I} \right)$$ where $\mathbf{I}$ is the $3 \times 3$ identity matrix.

Triangle Solid Angle Vertex Gradient ($\frac{\partial \Omega_t}{\partial v_k} \in \mathbb{R}^3$):

For a fixed query point $x_j$, let the relative vertex vectors be $a = v_1 - x_j$, $b = v_2 - x_j$, and $c = v_3 - x_j$, with corresponding unit vectors $\hat{a} = \frac{a}{\|a\|}$, $\hat{b} = \frac{b}{\|b\|}$, and $\hat{c} = \frac{c}{\|c\|}$. The solid angle $\Omega_t$ is defined via the Oosterom–Strackee formula mentioned previously:

$$N = a \cdot (b \times c)$$ $$D = \|a\|\|b\|\|c\| + (a \cdot b)\|c\| + (b \cdot c)\|a\| + (c \cdot a)\|b\|$$ $$\Omega_t(x_j) = 2 \arctan \left( \frac{N}{D} \right)$$

Differentiating $\Omega_t$ with respect to each vertex yields the per-query gradient configuration:

$$\frac{\partial \Omega_t}{\partial v_k} = 2 \frac{D \frac{\partial N}{\partial v_k} - N \frac{\partial D}{\partial v_k}}{N^2 + D^2} \quad \text{for } k \in \{1, 2, 3\}$$

The explicit intermediate vector derivatives for each vertex $v_k$ are evaluated as follows:

Vertex Numerator Derivative ($\frac{\partial N}{\partial v_k}$) Denominator Derivative ($\frac{\partial D}{\partial v_k}$)
$v_1$ $b \times c$ $(\|b\|\|c\| + b \cdot c)\hat{a} + \|c\|b + \|b\|c$
$v_2$ $c \times a$ $(\|c\|\|a\| + c \cdot a)\hat{b} + \|a\|c + \|c\|a$
$v_3$ $a \times b$ $(\|a\|\|b\| + a \cdot b)\hat{c} + \|b\|a + \|a\|b$
Note: These formulations are well-defined everywhere except at singular boundaries (such as directly on the triangle's edges or vertices) where the denominator vanishes ($|D| < 10^{-12}$). In the forward pass, this condition corresponds to a boundary discontinuity where a constant value of $0.5$ is returned; consequently, the gradient is explicitly set to $0$ in this case.

Evaluating these formulations directly requires a dense pairwise computation between all $N$ geometry elements and all $M$ query points, yielding a computational complexity of $\mathcal{O}(NM)$.

However, the underlying field behavior of the backward pass closely mirrors the forward pass. While the forward pass clusters geometry sources because their collective potential field varies smoothly over distance, the backward pass allows us to cluster query points. The adjoint fields, specifically the free-space Green's function gradient $\nabla G$, its spatial Hessian $\mathbf{H}_G$, and the boundary derivative vector field, are smooth fields that decay rapidly with distance to the query position. Consequently, the influence of a distant cluster of query points on any given source primitive can be replaced with a Taylor approximation at the cluster center. This leads to a hierarchical approximation scheme that is very similar to the Fast Winding Numbers scheme in the forward pass.

Truncated Taylor Expansions for Partial Derivatives

Case 1: Scaled Normals (Dipole Moments $m_i$)

Expanding the vector dipole gradient field $\nabla G(x_j - p_i)$ for a single source point $p_i$ relative to the distant query cluster centroid $\tilde{q}$ yields the multi-order tensor series:

$$\begin{align*}\nabla G(x_j - p_i) \approx & \nabla G(\tilde{q} - p_i) +\\ & \mathbf{H}_G(\tilde{q} - p_i)(x_j - \tilde{q}) + \\ & \frac{1}{2} \nabla^3 G(\tilde{q} - p_i) : \left[ (x_j - \tilde{q}) \otimes (x_j - \tilde{q}) \right]\end{align*}$$

Distributing the summation over the query indices $j$ isolates the geometry-independent parameters inside the parentheses:

$$\begin{align*}\frac{\partial \mathcal{L}}{\partial m_i} \approx & \nabla G(\tilde{q} - p_i) \left( \sum_j g_j \right) + \\ & \mathbf{H}_G(\tilde{q} - p_i) \left( \sum_j g_j (x_j - \tilde{q}) \right) +\\ & \frac{1}{2} \nabla^3 G(\tilde{q} - p_i) : \left( \sum_j g_j (x_j - \tilde{q}) \otimes (x_j - \tilde{q}) \right)\end{align*}$$

Note: The target gradient field output is a vector in $\mathbb{R}^3$. Inside the expansion, the first term scales a vector ($\mathbb{R}^3$) by a scalar; the second multiplies a matrix ($\mathbb{R}^{3 \times 3}$) by a vector; the third double-contracts ($:$) a 3rd-order tensor ($\mathbb{R}^{3 \times 3 \times 3}$) against a second-order outer-product matrix ($\mathbb{R}^{3 \times 3}$).

Case 2: Point Positions ($p_i$)

Because the position gradient natively evaluates through the matrix-valued potential Hessian $\mathbf{H}_G(x_j - p_i)$, we construct the Taylor series of this second-order tensor field around the cluster centroid $\tilde{q}$:

$$\begin{align*}\mathbf{H}_G(x_j - p_i) \approx & \mathbf{H}_G(\tilde{q} - p_i) +\\ & \nabla^3 G(\tilde{q} - p_i) \cdot (x_j - \tilde{q}) + \\ & \frac{1}{2} \nabla^4 G(\tilde{q} - p_i) : \left[ (x_j - \tilde{q}) \otimes (x_j - \tilde{q}) \right]\end{align*}$$

Factoring out the invariant source parameters and the dipole vector $m_i$, distributing the scalar loss adjoint $g_j$ isolates our unified query moments inside the contractions:

$$\begin{align*}\frac{\partial \mathcal{L}}{\partial p_i} \approx - \Biggl[ & \mathbf{H}_G(\tilde{q} - p_i) \left( \sum_j g_j \right) + \\ & \nabla^3 G(\tilde{q} - p_i) \cdot \left( \sum_j g_j (x_j - \tilde{q}) \right) + \\ & \frac{1}{2} \nabla^4 G(\tilde{q} - p_i) : \left( \sum_j g_j (x_j - \tilde{q}) \otimes (x_j - \tilde{q}) \right) \Biggr] m_i\end{align*}$$

Note: The bracketed terms resolve to a second-order matrix ($\mathbb{R}^{3 \times 3}$). The first term scales the baseline matrix by a scalar ($G_0$); the second contracts a 3rd-order tensor against a vector ($G_1$) via a single dot product; the third reduces a 4th-order tensor against a matrix ($G_2$) via a double contraction. Post-multiplying the resulting matrix by the vector $m_i \in \mathbb{R}^3$ yields the final derivative vector in $\mathbb{R}^3$.

Case 3: Triangle Vertices ($v_k$)

Denoting the local vertex boundary derivative vector field as $F_k(x_j) = \frac{\partial \Omega_T(x_j)}{\partial v_k} \in \mathbb{R}^3$, we express its Taylor serialization about $\tilde{q}$ directly as:

$$\begin{align*}F_k(x_j) \approx & F_k(\tilde{q}) +\\ & \nabla F_k(\tilde{q})(x_j - \tilde{q}) + \\ & \frac{1}{2} \mathbf{H}_{F_k}(\tilde{q}) : \left[ (x_j - \tilde{q}) \otimes (x_j - \tilde{q}) \right]\end{align*}$$

Distributing the scalar loss adjoint $g_j$ yields the localized structural summation form:

$$\begin{align*}\frac{\partial \mathcal{L}}{\partial v_k} \approx & F_k(\tilde{q}) \left( \sum_j g_j \right) + \\ & \nabla F_k(\tilde{q}) \left( \sum_j g_j (x_j - \tilde{q}) \right) + \\ & \frac{1}{2} \mathbf{H}_{F_k}(\tilde{q}) : \left( \sum_j g_j (x_j - \tilde{q}) \otimes (x_j - \tilde{q}) \right)\end{align*}$$

Note: Because the baseline field $F_k$ maps directly to a vector, its spatial gradient $\nabla F_k$ is a matrix $(\mathbb{R}^{3 \times 3})$ and its Hessian $\mathbf{H}_{F_k}$ is a 3rd-order tensor $(\mathbb{R}^{3 \times 3 \times 3})$. The expression smoothly resolves to $\mathbb{R}^3$ via vector-scalar scaling, matrix-vector multiplication, and a double tensor contraction ($:$) against the matrix-valued displacement outer product.

Unified Query Moments

The critical revelation of this parallel derivation is revealed by inspecting the encapsulated terms in parentheses across all three primitive spaces. These are the components that are entirely independent of the source geometry parameters and can be completely precomputed. Remarkably, the summation configurations match perfectly in every single equation.

All derivatives can be simultaneously computed using the exact same three multi-order Query Moments cached within the hierarchy layout:

0th-Order Moment (Scalar $G_0 \in \mathbb{R}$) $$G_0 = \sum_{j} g_j$$
1st-Order Moment (Vector $G_1 \in \mathbb{R}^3$) $$G_1 = \sum_{j} g_j (x_j - \tilde{q})$$
2nd-Order Moment (Matrix $G_2 \in \mathbb{R}^{3 \times 3}$) $$G_2 = \sum_{j} g_j (x_j - \tilde{q}) \otimes (x_j - \tilde{q})$$

By constructing the tree over the query coordinates on the fly, a single bottom-up parallel pass calculates these moments. During execution, geometry elements act as tree traversers, evaluating field derivatives through local $\mathcal{O}(1)$ contractions against the cached $G_0$, $G_1$, and $G_2$ node blocks. This limits the backward computational pass strictly to a highly efficient execution complexity of $\mathcal{O}((M + N) \log M)$, which natively accounts for the initial $\mathcal{O}(M)$ tree construction overhead.


Engine Architecture Status & Implementation Matrix

Forward Pass Evaluation Pipeline Complete

Reference Implementation (Brute-Force)

Fully functional, templated codebase supporting both point clouds and triangles. Serves as the exact numerical verification baseline for all accelerated tracking paths.

Tree-Accelerated Production Kernels (BVH8)

Highly unified, template-driven pipelines for point clouds and triangles running on an 8-way bounding volume hierarchy. Structural divergences are isolated strictly to the offline tree construction routines. Mesh inputs undergo an automated conversion pass to triangle arrays prior to evaluation.

Backward Pass Differentiation Pipeline In Active Development

Theoretical Adjoint Framework

The mathematical formulation for reverse-mode automatic differentiation via multi-order query moments is established. Adjoint structures map uniformly across all primitive geometries.

Kernel Architecture Design

Currently planning the underlying implementation details for parallel backpropagation. Active tasks focus on determining optimal memory layouts and contraction logic for the dynamic query-tree caching system.

Architecture Layout

Engine Pipeline & Features

High-performance hardware design choices engineered for massive geometric datasets.

05. Interactive Performance Benchmarks - Forward Pass

Compare performance metrics between different data modalities and execution backends using the interactive charts and options below.

Engine Backends & Baselines

  • Brute Force (GPU): Winder’s custom GPU-optimized brute force code path. Serves as the exact ground-truth numerical reference across all modality workflows.
  • Winder: Our production pipeline, evaluating highly optimized accelerated tracking steps across an 8-way bounding volume hierarchy.
  • PyTorch: A raw reference implementation utilizing optimized calculations via compiled execution loops (torch.compile).
  • LibIGL (CPU): The python bindings of the original standard implementation from Fast Winding Numbers for Soups and Meshes (Barill et al., ACM TOG 2018). Note: The current python bindings only expose the mesh variant. Build time and compute time could not be untangled.
  • Fast Dipole Sums: The original reference code implementation for Fast Dipole Sums for Point Clouds (Sharp et al., ACM TOG 2024). ⚠️ Note: Fast Dipole Sums encountered a memory access error during execution on code jaws_attack_benchy.obj and is therefore missing for this method.

Evaluation Workload

Query Dimensionality Performance profiles capture continuous operations mapped across 1,200 temporal frames running over a fixed $256^2$ spatial grid.
Execution Mode All evaluation schemes except for "torch" pass the entire allocation layout of $256^2 \times 1,200$ evaluation targets simultaneously within a single contiguous array kernel invocation to maximize hardware concurrency.

Hardware Environment

GPU NVIDIA RTX PRO 6000 Blackwell Max-Q Workstation Edition
CPU AMD Ryzen Threadripper PRO 5955WX 16-Cores

🎬 Synchronized Grid Comparison

Mesh Track