1 Introduction↩︎

Non-Uniform Rational B-Splines (NURBS) [1] are the standard curve and surface representation in CAD and manufacturing. Every major CAD kernel—SolidWorks, CATIA, Rhino, Siemens NX—and exchange formats such as STEP (ISO 10303) and IGES use NURBS, and isogeometric analysis [2] builds directly on them for finite-element simulation. Their advantage over polynomial B-splines is well established: rational weights and non-uniform knots let NURBS exactly represent conic sections—circles, ellipses, parabolic arcs—that polynomial splines can only approximate [1], [3], while also providing per-control-point shape modulation for tighter local geometric control.

As pixel-based image generation models improve, differentiable rendering offers a promising way to convert these pixel outputs into physical, manufacturable vector formats. Because splines are compact and resolution-independent, they are highly suited for such tasks. Berio  [4], for example, showed that optimizing long smoothing splines can produce high-quality artistic abstractions that can be reproduced with robotic pen plotters. However, spline representations remain only partially explored in differentiable vector graphics. Current differentiable 2D renderers [5], [6] are restricted to B’ezier curves, while Berio  [4] support B-splines but do not introduce rational weights or non-uniform knots. Consequently, differentiable renderers cannot faithfully represent shapes that CAD systems natively support—such as exact circular arcs and conic profiles—and remain disconnected from the NURBS-based modeling infrastructure that industry has relied on for decades.

We address this with NURBS Splatting. We adaptively sample isotropic Gaussians along NURBS curves—positions, widths, and analytical derivatives computed via the Cox–de Boor recursion—and inside closed regions via an SDF-modulated grid, then composite them through a tile-based Gaussian splatting rasterizer. The pipeline is fully differentiable: image-space losses drive gradients back to control points, rational weights, knot intervals, stroke widths, colors, and opacities. For filled closed curves, grid-step annealing progressively refines spatial resolution during optimization.

The main contributions of this work are as follows:

• The first differentiable renderer for NURBS curves in 2D image space, supporting joint optimization of rational weights and non-uniform knot vectors.

• A Gaussian splatting formulation for both open stroke contours and filled closed regions, combining adaptive arc-length sampling with SDF-modulated interior fill and grid-step annealing.

• Evaluation on calligraphy reconstruction, image vectorization, and neural image abstraction showing that rational splines improve both quality and speed over polynomial baselines.

2 Related Work↩︎

2.1 Differentiable Rendering↩︎

DiffVG [5] pioneered differentiable rasterization, which enabled gradient-based optimization of vector images. A variety of applications have since been explored, including image vectorization [7]–[10], text-guided drawing synthesis fransCLIPDrawExploringTexttoDrawing2022?, vinkerCLIPassoSemanticallyAwareObject2022?, text-to-SVG generation via score distillation [11]–[13], semantic typography [14], collage generation [15], animation generation [16], sketch refinement [17], . However, the DiffVG approach relies on per-pixel intersection tests, making high-resolution optimization slow and memory-intensive.

An alternative line of work uses Gaussian primitives for rendering. In 3D vision, Kerbl  [18] propose 3D Gaussian Splatting for novel-view synthesis. BézierGS [19] and SplineGS [20] further extend this paradigm with spline-based motion modeling. The idea has since been adapted to 2D domains. Huang  [21] introduce 2D Gaussian Splatting for radiance-field reconstruction, collapsing the 3D volume into oriented 2D Gaussian disks to enforce view-consistent geometry. Likewise, Zhang  [22] present GaussianImage, encoding images directly as collections of 2D Gaussians to achieve ultra-fast rendering and compression. In the context of vector graphics, Liu  [6] propose Bézier Splatting, which samples 2D Gaussian kernels along Bézier curves and rasterizes via a Gaussian-splatting pipeline. These works demonstrate that Gaussian splatting is a powerful differentiable primitive for both 3D scene reconstruction and 2D vector graphics rendering.

Notably, all the aforementioned methods employ either polynomial Bézier curves or polynomial B-splines as the underlying geometric primitives. To the best of our knowledge, no existing differentiable rendering framework supports rational spline primitives for 2D image-space optimization.

2.2 Non-Uniform Rational B-Splines↩︎

NURBS [1] generalize polynomial B-splines by introducing rational weights and non-uniform knot vectors, allowing exact representation of conic sections, analytic curves, and flexible free-form shapes. They have become the de facto industry standard for modeling curves and surfaces in CAD and graphics. While classical algorithms like the Cox–de Boor recursion [23], deBoorCalculatingBSplines1972? evaluate and sample points for rendering, many modern professional tools employ GPU curve rasterization based on the implicit form of rational Bézier curves obtained from NURBS via knot insertion [24]. Beyond geometric modeling, NURBS have found broader scientific impact through isogeometric analysis [2], which unifies CAD geometry and finite-element analysis by using NURBS basis functions directly as the computational basis.

Recent work has begun to integrate NURBS into differentiable and learning-based pipelines. Deva Prasad [25] introduce NURBS-Diff, a differentiable NURBS layer that provides forward and backward passes through the NURBS evaluation process, enabling tasks such as curve and surface fitting, surface offsetting, . Worchel and Alexa [26] propose differentiable rendering of parametric geometry, including B-spline surfaces, by directly deriving screen-space gradients for parametric primitives. Tojo  [27] optimize B-spline curves via differentiable rendering for fabricable 3D wire art. On the generative side, SplineGen [28] leverages a generative model to fit B-spline surfaces to unorganized point clouds, NeuroNURBS [29] learns compact NURBS surface representations for 3D solids, and BrepGen [30] employs diffusion models to generate B-rep CAD models whose faces are NURBS surfaces. Despite this growing interest, existing differentiable and generative NURBS methods predominantly target 3D geometry fitting or CAD reconstruction rather than image-based rendering.

3 Method↩︎

Our pipeline is illustrated in 1. We first define NURBS curve primitives with learnable control points, rational weights, and knot vectors (3.1). Isotropic Gaussians are then adaptively sampled along curve contours (3.2) or within closed filled regions via an SDF-modulated grid (3.3). These Gaussians are composited through a tile-based splatting rasterizer to produce a differentiable image (3.4), and image-space losses drive gradients back to all curve parameters (3.5). Although we focus on NURBS, our Gaussian sampling formulation is generic and readily extends to other parametric primitives.

3.1 NURBS Primitives↩︎

Each curve in our system is a NURBS curve of degree \(p\) (order \(k{=}p{+}1\)) with \(n+1\) control points: \[\mathbf{C}(u) = \frac{\sum_{i=0}^{n} N_{i,p}(u)\, w_i \, \mathbf{P}_i}{\sum_{i=0}^{n} N_{i,p}(u)\, w_i}=\frac{\mathbf{S}(u)}{W(u)}, \quad u \in [u_0, u_m],\] where \(\{\mathbf{P}_i \in \mathbb{R}^2\}\) are control points, \(\{w_i \in \mathbb{R}^+\}\) are rational weights, and \(\{N_{i,p}(u)\}\) are B-spline basis functions defined over a non-decreasing knot vector \(\mathbf{U} = \{u_0, u_1, \dots, u_m\}\). The basis functions \(N_{i,p}(u)\) are defined by the Cox–de Boor recursion [23], deBoorCalculatingBSplines1972?, which guarantees compact support and \(C^{p-1}\) continuity.

The curve derivative, needed for arc-length estimation and regularization, follows by the quotient rule: \[\label{eq:curve95derivative} \mathbf{C}'(u) = \frac{ W(u) \sum_{i=0}^{n} N'_{i,p}(u)\, w_i \, \mathbf{P}_i - \mathbf{S}(u) \sum_{i=0}^{n} N'_{i,p}(u)\, w_i }{W^2(u)},\tag{1}\] where the basis-function derivatives \(N'_{i,p}\) are also obtainable via the Cox–de Boor recursion. The arc length used in 3.2 can be approximated by numerical quadrature: \[\label{eq:arc95length} L = \int_{u_p}^{u_{m-p}} \|\mathbf{C}'(u)\|_2 \, du \approx \sum_{j=1}^{N_L} \|\mathbf{C}'(u_j)\|_2 \, \mathrm{\Delta} u.\tag{2}\]

Rather than optimizing control points and knots directly, we leverage key points \(\{\mathbf{K}_i\}_{i=1}^{n_k}\) as the primary learnable parameters and derive control points and knot vectors from them.

For open curves, the control points coincide with the key points, and the knot vector is clamped at both endpoints with multiplicity \(k\). Only the \(n_k - k + 1\) interior knot intervals \(\{\tau_j > 0\}\) are stored as learnable parameters; the full knot vector is reconstructed by prepending and appending \(p\) zero-intervals (to enforce endpoint multiplicity) and taking a cumulative sum, which inherently guarantees monotonicity.

For closed curves, the control polygon is formed by cyclically wrapping the last \(\lceil p/2 \rceil\) and first \(\lfloor p/2 \rfloor\) key points to produce \(n_k + p\) control points. We store \(n_k\) core knot intervals; the periodic knot vector is then built by wrapping the last \(p\) core intervals before and the first \(p\) after, followed by a cumulative sum.

Each curve also carries an RGB color \(\mathbf{c} \in [0,1]^3\) and an opacity \(o \in [0,1]\). We extend the control points to \(\mathbb{R}^3\) by appending a third coordinate encoding local stroke width, so evaluating \(\mathbf{C}(u)\) simultaneously yields position and width at parameter \(u\).

3.2 Gaussian Sampling on Contours↩︎

Given the knot vector \(\mathbf{U}\) and control points derived from the key points (3.1), we adaptively place isotropic Gaussians along each curve with density proportional to arc length.

We first approximate the arc length \(L\) via 2 using \(N_L\) uniformly spaced parameter samples, then set the number of Gaussians as \(M = \lceil \delta_\text{c} \cdot L \rceil\), where \(\delta_\text{c}\) is a user-specified contour density. The \(M\) parameter values \(\{u_j\}_{j=1}^{M}\) are distributed uniformly within the valid knot domain \([u_p, u_{m-p}]\).

To evaluate each sample efficiently, we exploit the compact support of B-spline basis functions: at parameter value \(u_j\), only the \(p{+}1\) basis functions whose support contains \(u_j\) are non-zero. The Gaussian center is therefore computed via a span-local formulation: \[\boldsymbol{\mu}_j = \frac{ \sum_{l=0}^{p} N_{i_j-p+l,p}(u_j)\, w_{i_j-p+l}\, \mathbf{P}_{i_j-p+l} }{ \sum_{l=0}^{p} N_{i_j-p+l,p}(u_j)\, w_{i_j-p+l} },\] where \(i_j\) is the knot span index satisfying \(u_{i_j} \le u_j < u_{i_j+1}\). This reduces evaluation cost from \(O(Mn)\) to \(O(Mp)\) per curve. The local stroke width \(s_j\) is read from the third coordinate of \(\mathbf{C}(u_j)\); higher-order analytical derivatives, used for regularization (3.5), are obtained through the same local recursion.

Each sample is assigned an isotropic Gaussian whose scale depends solely on the local stroke width: \[\sigma_{j,x} = \sigma_{j,y} = s_j / \rho, \qquad \theta_j = 0, \qquad o_j = o,\] where \(\rho\) is a global scale ratio controlling overlap between neighboring Gaussians and \(o\) is the curve-level opacity from 3.1. Unlike Bézier Splatting [6], which aligns anisotropic Gaussians to the curve tangent, this isotropic formulation eliminates rotation estimation and avoids visual artifacts—blurry edges at large stroke widths and spikes at high-curvature regions. Because \(M\) scales with arc length, Gaussian centers remain sufficiently dense that circular kernels provide smooth, uniform coverage without directional alignment.

For open curves, we replicate the first and last Gaussian centers \(\kappa\) times, where \(\kappa\) is a fixed number of times proportional to \(\delta_\text{c}\), ensuring opaque, well-defined endpoints without additional parameters.

3.3 Gaussian Sampling in Regions↩︎

For closed curves marked as filled, we replace contour Gaussians with a dense set of interior Gaussian primitives whose opacity is modulated by a differentiable signed distance field (SDF).

We first sample \(\lceil\delta_\text{b} \cdot L\rceil\) boundary points \(\{\mathbf{b}_j\}\) along the closed curve in pixel coordinates, where \(\delta_\text{b}\) is a user-specified boundary density and \(L\) is the arc length from 3.2. These boundary points define the polyline against which the SDF is evaluated.

Next, we compute the axis-aligned bounding box of the boundary with symmetric padding proportional to its extent, clipped to the canvas. A uniform 2D grid at step size \(h\) is placed inside this padded box, yielding \(Q = n_x \times n_y\) query points \(\{\mathbf{q}_j\}\) in pixel coordinates. A larger \(h\) produces a coarser grid with fewer Gaussians, trading density for speed; we anneal \(h\) during optimization to accelerate early convergence while preserving final detail (4.2).

For each query point \(\mathbf{q}_j\), we generate an isotropic Gaussian primitive: \[\boldsymbol{\mu}_j = \mathbf{q}_j, \qquad \theta_j = 0, \qquad \sigma_{j,x} = \sigma_{j,y} = \frac{\eta \cdot h}{\rho},\] where \(\eta\) is a coverage factor (we use \(\eta{=}1.5\)) that ensures sufficient overlap between neighboring Gaussians, and \(\rho\) is the global scale ratio from 3.2. The base opacity is the product of the curve-level opacity \(o\) and a sigmoid soft mask derived from the SDF: \[o_j = o \cdot \mathrm{sigmoid}\!\left(\frac{\gamma}{h} \cdot \mathrm{sdf}(\mathbf{q}_j)\right),\] where \(\gamma\) controls the sharpness of the boundary transition, and \(\mathrm{sdf}(\mathbf{q}_j)\) returns positive values inside the region and negative values outside, computed with respect to the boundary polyline \(\{\mathbf{b}_j\}\). Detailed algorithmic steps for computing the SDF are provided in the supplementary material.

Each Gaussian inherits its color from the parent curve. This pixel-space SDF formulation scales to large filled shapes while remaining fully differentiable with respect to boundary point positions.

3.4 Differentiable Splatting Rendering↩︎

Once Gaussian primitives are sampled from NURBS curves, we rasterize them using the tile-based Gaussian splatting rasterizer of GaussianImage [22].

The rasterizer sorts the 2D Gaussians by their stacking order (depth) and composites them front-to-back for each pixel. The rendered color \(\mathcal{I}(\mathbf{x}_n)\) at pixel \(\mathbf{x}_n\) is obtained by alpha blending the contributions of overlapping Gaussians: \[\mathcal{I}(\mathbf{x}_n) = \sum_{i \in \mathcal{N}} \mathbf{c}_i \, \alpha_i(\mathbf{x}_n) \prod_{j=1}^{i-1} \bigl(1-\alpha_j(\mathbf{x}_n)\bigr),\] where \(\mathcal{N}\) is the depth-ordered set of Gaussians covering pixel \(\mathbf{x}_n\), \(\mathbf{c}_i\) is the color of the \(i\)-th Gaussian, and \(\alpha_i(\mathbf{x}_n)\) is its opacity at that location. Any residual transmittance is composited against a predefined background color.

Since all sampled Gaussians are isotropic with \(\sigma_{i,x} = \sigma_{i,y} \eqqcolon \sigma_i\) and \(\theta_i = 0\) ([sec:gaussian_contour_sampling] [sec:interior_filling]), the covariance reduces to \(\boldsymbol{\Sigma}_i = \sigma_i^2 \mathbf{I}\) and the per-pixel Gaussian opacity simplifies to \[\alpha_i(\mathbf{x}_n) = o_i \exp\!\biggl(-\frac{\|\mathbf{x}_n - \boldsymbol{\mu}_i\|^2}{2\,\sigma_i^2}\biggr),\] where \(o_i\) is the base opacity and \(\boldsymbol{\mu}_i\) is the Gaussian center. This formulation enables efficient backpropagation of gradients from the image loss through the Gaussian parameters to the NURBS control points and weights.

3.5 Optimization↩︎

The overall objective combines an image-space reconstruction loss with geometric regularization: \[\mathcal{L} = \mathcal{L}_{\text{img}}(\mathcal{I}, \mathcal{I}^{\ast}) + \lambda_{\text{geo}} \, \mathcal{L}_{\text{geo}},\] where \(\mathcal{I}\) and \(\mathcal{I}^{\ast}\) are the rendered and target images, \(\mathcal{L}_{\text{img}}\) is an image-based loss, and \(\mathcal{L}_{\text{geo}}\) enforces regularity on the underlying NURBS curves. Gradients propagate from pixels through the Gaussian parameters to NURBS control points and weights via the differentiable pipeline above; the specific choice of \(\mathcal{L}_{\text{img}}\) is application-dependent. Below we describe the two regularization terms shared across all applications; additional task-specific losses are introduced in 4.

Derivative Loss, the squared magnitude of the \(r\)-th derivative of a parametric curve, is a classical smoothing energy in curve fairing [31]: the second derivative relates to curvature, and the third captures its rate of change. We penalize this energy averaged over the parameter domain: \[\mathcal{L}_{\text{deriv}}^r(\mathbf{C}) = \frac{1}{u_{\max} - u_{\min}}\int_{u_{\min}}^{u_{\max}} \left\| \frac{d^r \mathbf{C}(u)}{d u^r} \right\|_2^2 \, du,\] where \(r\) denotes the derivative order. We evaluate this integral using analytical derivatives obtained from 1 during contour sampling (3.2), avoiding the numerical errors that arise from finite-difference approximations. The total loss averages \(\mathcal{L}_{\text{deriv}}^r\) over all curves; following Berio  [4], we set \(r=3\).

Bounding Box Loss prevents curves from drifting outside the target region. We penalize points that violate the image bounding box: \[\mathcal{L}_{\text{bbox}}=\sum_i \varphi(\mathbf{b}_{\min} - \mathbf{P}_i) + \varphi(\mathbf{P}_i - \mathbf{b}_{\max}),\] where \(\mathbf{b}_{\min}\) and \(\mathbf{b}_{\max}\) are the bounding-box corners, \(\mathbf{P}_i\) ranges over all sampled curve positions, and \(\varphi\) is a soft penalty (ReLU or Softplus) that activates only when a point lies outside the permitted region.

4 Applications↩︎

We demonstrate the versatility of NURBS Splatting on three representative vector-graphics tasks: calligraphy reconstruction (4.1), layer-wise image vectorization (4.2), and neural image abstraction (4.3).

4.1 Calligraphy Reconstruction↩︎

In calligraphy, stroke trajectories and their varying widths encode both style and meaning. Recovering these as editable vector strokes from raster images is useful for font design, digital archiving, and robotic handwriting reproduction.

Most calligraphy generation methods produce raster images via GANs [32], [33] or diffusion models [34], discarding stroke structure entirely. Methods that output vector outlines as Bézier sequences [35]–[37] describe glyphs as closed regions, not as strokes with varying width. Classical stroke extraction [38], [39] and the segmentation method of StrokeStyles [40] do model individual strokes but lack end-to-end differentiable optimization in image space.

Our framework represents each stroke as a NURBS curve with a per-control-point width channel. For initialization, we use the line vectorization method of Magne and Sorkine-Hornung [41]: a medial axis is extracted from the binarized image via the Voronoi diagram method [42], a CNN classifies intersections to produce ordered sample points along the axis, and NURBS curves are fit to these points with constant initial width.

The curves are then optimized end-to-end against the raster target with MSE loss and geometric regularization: \[\label{eq:calligraphy-loss} \mathcal{L} = \lambda_{\text{MSE}} \mathcal{L}_{\text{MSE}} + \lambda_{\text{deriv}} \mathcal{L}_{\text{deriv}}^3 + \lambda_{\text{bbox}} \mathcal{L}_{\text{bbox}}.\tag{3}\] Control point positions, NURBS weights, knot intervals, and per-control-point widths are updated jointly, recovering both stroke trajectories and thickness variations. Results are presented in 5.2.

4.2 Layer-wise Image Vectorization↩︎

Ma  [7] propose LIVE, a layer-wise vectorization approach that greedily adds closed Bézier paths one at a time, each initialized at the region of highest reconstruction error and optimized via DiffVG [5]. A distance-weighted loss modulates per-pixel gradients based on proximity to the newly added shape boundary, concentrating optimization effort where it matters most. We adapt LIVE to our framework, replacing the DiffVG rendering backend with NURBS Splatting while introducing a grid-step annealing strategy for filled shapes.

To accelerate convergence without sacrificing fidelity, we anneal the fill grid step \(h\) from a coarse value (\(h_{\max}{=}4\)) to a fine one (\(h_{\min}{=}1\)) using a cosine schedule between iteration fractions \(t_b{=}0.1\) and \(t_e{=}0.8\): \[h(t) = h_{\max} + (h_{\min} - h_{\max})\,\frac{1 + \cos(\pi(1 - \phi))}{2}, \quad t/T \in (t_b, t_e),\] where \(\phi = (t/T - t_b)/(t_e - t_b)\) and \(h\) is clamped outside the range. The coarse grid reduces the number of primitives in the early stages, thereby stabilizing large-scale deformations (2).

The overall loss for each round combines a pixel-wise reconstruction term with geometric regularization: \[\mathcal{L} = \mathcal{L}_{\text{img}} + \lambda_{\text{deriv}} \mathcal{L}_{\text{deriv}}^3 + \lambda_{\text{xing}} \mathcal{L}_{\text{xing}},\] where \(\mathcal{L}_{\text{img}}\) is the distance-weighted MSE loss, \(\mathcal{L}_{\text{deriv}}^3\) is the third-order derivative smoothing loss, and \(\mathcal{L}_{\text{xing}}\) is a self-crossing penalty adapted from LIVE that penalizes control-point configurations whose consecutive segments reverse winding direction. Results are presented in 5.3.

4.3 Neural Image Abstraction↩︎

Berio  [4] apply long smoothing B-splines to generate abstracted vector images from raster inputs. Because NURBS generalize polynomial B-splines, these applications transfer naturally to our framework while gaining rational-weight flexibility. We demonstrate two representative tasks: single-stroke area filling and diffusion-guided image abstraction.

The Single-Stroke Area Filling task aims to produce a single NURBS curve that fills a colored shape given its bitmap image. We initialize a long NURBS curve from a TSP route [43] and optimize it with an MSE coverage loss between the rendered and target images. Lowering the target opacity reduces curve density in the interior. To enable semantic stylization, we add a term \(\mathcal{L}_{\text{style}}\) based on the patch-wise directional CLIP loss of Kwon and Ye [44], which aligns encoded features of the rendered curves with those of a reference style image. Geometric regularization via \(\mathcal{L}_{\text{deriv}}^3\) and \(\mathcal{L}_{\text{bbox}}\) enforces smoothness and prevents control-point drift. Representative results appear in 3.

For the Single-Stroke Image Abstraction task, we couple our renderer with a diffusion-based Score Distillation Sampling (SDS) pipeline. Following Berio  [4], the pipeline uses ControlNet [45] with Canny edge detection and IP-Adapter [46], guided by the text prompt “A black and white ink drawing.” The SDS gradient with respect to NURBS parameters \(\boldsymbol{\psi}\) is \[\nabla_{\boldsymbol{\psi}} \mathcal{L}_{\text{SDS}}(\boldsymbol{\psi}) = \mathbb{E}_{t} \left[ w(t) \left(\hat{\boldsymbol{\epsilon}}(x_t, t, y) - \boldsymbol{\epsilon}\right) \frac{\partial g(\boldsymbol{\psi})}{\partial \boldsymbol{\psi}}\right],\] where \(\hat{\boldsymbol{\epsilon}}(x_t, t, y)\) is the predicted denoising direction of the frozen diffusion model for a latent \(x_t\) at time step \(t\), \(\boldsymbol{\epsilon}\) is the sampled noise, and \(w(t)\) is a weighting function dependent on \(t\). The final loss combines SDS with geometric regularization \(\mathcal{L}_{\text{deriv}}^3\) and \(\mathcal{L}_{\text{bbox}}\). As shown in 4, our method achieves \(1.27\times\) faster optimization than Berio  [4] at comparable stylization quality—even though differentiable rendering is not the primary bottleneck of the SDS pipeline and NURBS evaluation is more complex than polynomial B-spline evaluation.

5 Experiments↩︎

We implement our NURBS Splatting framework in PyTorch. We show general comparisons on rendering efficiency, quality, and editability in 5.1; experiments for our calligraphy reconstruction application are detailed in 5.2; experiments for LIVE are in 5.3. All experiments are run on a single NVIDIA RTX 4090. Code and dataset are available at https://github.com/AnicoderAndy/nurbs-splatting.

5.1 Efficiency, Quality, and Editability↩︎

Our method shows improved efficiency compared to Berio  [4] when the number of control points scales, as it avoids matrix multiplication to convert the splines. However, the superiority is not seen when rendering with simple geometric primitives, comparing our method with Berio  [4] and Bézier Splatting (BS) [6]. Details are provided in the supplementary material.

As for rendering quality, we qualitatively compare our method with BS which uses anisotropic Gaussian sampling on three diagnostic cases (5). For thick strokes, BS exhibits severely blurred edges. In the high-curvature case, BS produces spike artifacts due to the discontinuity and instability of Gaussian rotation evaluation via atan2. In the filled-region case, BS generates spikes at region junctions caused by incorrect curve interpolation at the shared endpoints of paired Bézier curves. Our method resolves all three artifacts through purely isotropic Gaussian sampling, consistently yielding sharp boundaries and smooth geometry.

A major advantage of vector graphics is their editability. The control points of optimized NURBS curves from our method provide a clean structure allowing manual adjustment. We present a visualization in the supplementary material.

5.2 Calligraphy Reconstruction↩︎

We construct a calligraphy benchmark of 192 characters: all 92 Japanese Kana and the first 100 characters from General Standard Chinese Characters [47]. Each character is rasterized from existing fonts [48], [49] to \(512{\times}512\) images.

Since, to the best of our knowledge, no existing method performs stroke-level reconstruction of calligraphy, we implement a baseline based on Berio  [4], which optimizes uniform B-splines. Both methods use identical initialization via medial-axis extraction [41]. We evaluate reconstruction quality using standard pixel-level measures (MSE, PSNR, SSIM) and geometry-aware metrics (symmetric Hausdorff distance and F1 score on binarized stroke masks).

1 shows the quantitative results averaged over 191 characters. At the contour density \(\delta_\text{c}{=}18\), our full method reduces MSE by \(33\%\), improves PSNR by \(+1.6\) dB, and lowers the Hausdorff distance by \(30\%\), while running \(1.2{\times}\) faster. The speed gain stems from the avoidance of the matrix multiplication required to convert B-splines into Bézier curves. Gains are most pronounced on Japanese Kana (\(+\)​2.10 dB PSNR), where frequent sharp bends and hooks benefit from non-uniform knot spacing that concentrates basis-function support in high-curvature regions. 6 presents qualitative comparisons on representative characters, showing that our method more faithfully preserves stroke trajectories and width variation than the baseline.

The default configuration drops slightly in SSIM because isotropic Gaussians produce softer boundary transitions than the analytic edge rendering of DiffVG. Increasing the sampling density to \(\delta_\text{c}{=}30\) sharpens these boundaries, recovering the best SSIM (\(\boldsymbol{0.9824}\)) and achieving the highest F1 (\(\boldsymbol{0.9743}\)) at only a modest increase in optimization time.

Comparing only uniform polynomial B-splines—the Baseline against our “w/o w & k” variant—isolates the effect of the rendering backend: splatting yields a \(2.1{\times}\) speedup (\(5.12\) s vs.\(10.65\) s) while simultaneously improving PSNR by \(+1.2\) dB and Hausdorff by \(26\%\). This confirms that the Gaussian rasterization pipeline is consistently beneficial even without NURBS-specific parameters.

Ablation studies (bottom three rows of 1) isolate the individual contributions of learnable weights and knots. Removing rational weights (“w/o weights”) incurs the largest single-component degradation (\(-0.69\) dB PSNR, \(+1.01\) px Hausdorff relative to the full model), indicating that per-control-point weighting is especially valuable for modulating local curvature. Removing knots (“w/o knots”) has a smaller effect on reconstruction quality but provides the largest time saving (\(5.63\) s). Jointly dropping both (“w/o w & k”) still comfortably outperforms the baseline on all metrics except SSIM, while reducing optimization time to roughly half that of the baseline (\(5.12\) s)—offering a strong quality–speed trade-off for applications where rational parameterization is unnecessary. Importantly, the benefits of weights and knots are complementary: the full model improves \(+0.25\) dB over the next-best single-ablation variant (“w/o knots”), confirming that rational weights and non-uniform knots address distinct geometric limitations.

5.3 Layer-wise Image Vectorization↩︎

Since the original dataset used in LIVE [7] is not publicly available, we construct a new evaluation benchmark using Noto Emoji v1 [50], following their protocol. We select 100 representative emojis and rasterize them at a resolution of 480\(\times\)​480. We compare our method against the original LIVE [7], which uses DiffVG [5] as its rendering backend.

In each round, both methods initialize a new closed path with 12 key points as a circle and optimize for 300 iterations; 10 rounds are run in total, yielding 10 paths per image. We report standard reconstruction metrics: MSE, SSIM, and LPIPS [51], and measure the total optimization time per image.

2 summarizes quantitative results averaged over 100 emojis. All our variants substantially outperform the LIVE baseline on every metric while also requiring less optimization time.

Our full method with the grid-step annealing strategy reduces MSE by \(32\%\), improves SSIM by \(+0.017\), and lowers LPIPS by \(41\%\), while achieving a \(1.4{\times}\) speedup. Disabling step annealing (“w/o step”) further improves reconstruction quality at the cost of longer optimization (\(527\) s), yet remains \(11\%\) faster than LIVE, indicating that our optimization framework remains computationally efficient even without the annealing strategy. Removing rational weights and knots (“w/o w & k”) still outperforms the baseline on all metrics—reducing MSE by \(28\%\), improving SSIM by \(+0.016\), and lowering LPIPS by \(40\%\)—while achieving the fastest optimization at \(354\) s (\(1.7{\times}\) speedup), albeit with a smaller quality margin than the full model.

Qualitatively, 7 demonstrates that our method is more robust to complex topologies, particularly for thin interleaving structures in facial expressions, where the DiffVG-based LIVE baseline often exhibits artifacts such as missing details and geometric distortions.

6 Conclusion↩︎

We presented NURBS Splatting, the first differentiable rendering framework capable of optimizing rational B-splines in 2D image space. By representing curves and filled regions as continuous Gaussian fields, our method unifies vector graphics rendering with gradient-based optimization, allowing joint refinement of control points, rational weights, and non-uniform knot vectors. Experiments on calligraphy reconstruction and image vectorization demonstrate consistent improvements over polynomial baselines, achieving higher reconstruction fidelity. Crucially, the ability to optimize rational weights introduces geometric flexibility beyond prior differentiable vector graphics methods, narrowing the gap between learning-based reconstruction and CAD-standard representations.

Our current region-filling strategy relies on a heuristic SDF-based grid, which provides limited advantages over scanline methods for regular geometries. In addition, although splatting enables efficient rendering, the number of Gaussian primitives can grow with curve complexity, potentially increasing memory usage for highly detailed illustrations compared to implicit representations. Future work will focus on extending the framework to 3D NURBS surfaces to support native CAD geometry, developing a more efficient analytic region-filling algorithm, and exploring integration with generative models to synthesize editable vector designs directly from text or images.

Acknowledgements↩︎

We thank the anonymous reviewers for their valuable feedback and suggestions. This work was funded by the National Natural Science Foundation of China, No. 62076090; the Huxiang Youth Talent Support Program, Hunan Province, China, No. 2020RC3014; and the Natural Science Foundation of Hunan Province, China, No. 2022JJ30173.

7 Code and Data Availability↩︎

Code will be made available upon publication at https://github.com/AnicoderAndy/nurbs-splatting. All fonts used to construct our dataset are freely available online. However, due to licensing restrictions, some of these fonts cannot be redistributed directly. Instead, we will release the source code used for dataset construction.

8 Signed Distance Field Evaluation↩︎

The unsigned distance is the minimum Euclidean distance from \(\mathbf{q}\) to the nearest boundary segment: \[d(\mathbf{q}) = \min_{i} \left\| \mathbf{q} - \mathrm{proj}_{\overline{\mathbf{b}_i \mathbf{b}_{i+1}}}(\mathbf{q}) \right\|_2,\] where \(\mathrm{proj}_{\overline{\mathbf{b}_i \mathbf{b}_{i+1}}}(\mathbf{q})\) denotes the closest point on segment \(\overline{\mathbf{b}_i \mathbf{b}_{i+1}}\) to \(\mathbf{q}\), computed by clamping the projection parameter to \([0,1]\). The sign is determined by the winding number \(\mathrm{wn}(\mathbf{q})\), which is detached from the computation graph to prevent gradient flow through the discrete sign: \[\begin{gather} \mathrm{wn}(\mathbf{q}) = \frac{1}{2\pi} \sum_{i=1}^{M} \angle \left( \mathbf{b}_i - \mathbf{q},\; \mathbf{b}_{i+1} - \mathbf{q} \right),\\ \mathrm{sdf}(\mathbf{q}) = \mathop{\mathrm{sgn}}\!\big(\mathrm{wn}(\mathbf{q})\big) \cdot d(\mathbf{q}), \end{gather}\] where \(\mathop{\mathrm{sgn}}\left(\mathrm{wn}\right) = +1\) when \(\mathrm{wn} > 0.5\) (inside) and \(-1\) otherwise. This formulation yields positive values inside the region and negative values outside, and the smooth unsigned distance provides well-behaved gradients near the boundary.

9 Symbols and Settings↩︎

We utilize numerous mathematical symbols throughout the main paper and supplementary material. For clarity, we summarize their definitions and default settings in 1. In this table, “Calligraphy Reconstruction” is abbreviated as “Cal.”, and “Layer-wise Image Vectorization” is abbreviated as “LIVE”.

10 Experiments↩︎

We provide additional details on the efficiency test (10.1), the editability demonstration (10.2), the experimental setup and results for our calligraphy reconstruction experiment (10.3) and the LIVE experiment (10.4), as well as further discussion on our comparison with Bézier Splatting (10.5).

10.1 Efficiency Test↩︎

We conduct an efficiency test to evaluate the computational performance of our method compared to the baseline. 4 shows the comparison of forward/backward passes and the total time for each method under different resolution settings. It is shown that our method is slightly slower than the baseline methods at lower resolutions. Based on the profiling results, we hypothesize that the performance gap primarily stems from the evaluation of NURBS, which appears to dominate the overall runtime. Compared with methods based on Bézier curves or uniform B-splines, NURBS evaluation is inherently more computationally demanding because of its higher computational complexity, making it a likely performance bottleneck. Additionally, Bézier Splatting exhibits superior scalability as the number of geometric primitives increases. This advantage may be attributed to its representation of all curves as a single fixed-shape tensor, which enables efficient batched GPU computation but requires all curves to share the same geometric structure, such as an identical number of control points.

Nevertheless, our method demonstrates promising scalability as the number of control points increases, outperforming Berio  [4] in terms of efficiency. As shown in Fig. 8, our method exhibits a more favorable runtime trend as the number of control points increases. We attribute this advantage to the fact that Berio first convert B-spline segments into Bézier curves through Bézier extraction (implemented as a sequence of local matrix operations) before evaluation, introducing additional preprocessing overhead. In contrast, our method operates directly on the original NURBS representation, eliminating the need for such conversions. Consequently, although Bézier extraction has linear complexity with respect to the number of control points, its additional preprocessing cost accumulates as curve complexity increases. By avoiding this conversion entirely, our method scales more gracefully with the number of control points.

10.2 Editability Demonstration↩︎

9 shows the control points of the NURBS curves optimized for the calligraphy reconstruction task. Since mainstream vector graphics editors generally do not support NURBS editing, providing an editing demonstration in existing software is impractical. Instead, we visualize the optimized control-point layout, which remains compact and interpretable for manual adjustment.

10.3 Calligraphy Reconstruction↩︎

We extract initial curves from the input calligraphy images using SLDVec [41]. Because SLDVec generates a dense set of \(n_\text{SLD}\) points per stroke, we downsample each stroke to \(\max(30, 0.01 n_\text{SLD})\) control points to improve optimization efficiency. The dataset consists of 192 characters, one of which failed the extraction process and was therefore excluded from all experiments.

Both our method and the Berio  [4] baseline optimize degree \(p{=}5\) curves for 150 Adam iterations. We apply cosine learning-rate annealing, decaying the rate to \(10\%\) of its initial value. For our method, we set the learning rates to 1.5 for control-point positions, 0.2 for stroke widths, and 0.1 for both rational weights and knot intervals. We clamp the rational weights to \([0.01, 10]\) and knot intervals to \([0, 2]\). The baseline uses the same position learning rate but a slightly higher width learning rate (0.3), and it does not optimize rational weights or knots.

We weight our objective function components with \(\lambda_\text{MSE}{=}1\), \(\lambda_\text{deriv}{=}1\), and \(\lambda_\text{bbox}{=}5\). The baseline uses the DiffVG rasterizer [5] with independently tuned loss weights (\(\lambda_\text{MSE}{=}20\), \(\lambda_\text{deriv}{=}3\), \(\lambda_\text{bbox}{=}10\)).

We evaluate reconstruction quality using pixel-level and stroke-geometry metrics. To measure pixel-level fidelity, we compute MSE, PSNR, and SSIM between the rendered output and the ground-truth raster image. To evaluate stroke geometry, we binarize the rendered and target images at an intensity threshold of \(0.5\) to extract stroke masks. We then calculate the Hausdorff distance between the edge pixels of these masks, as well as the F1 score.

5 reports the detailed quantitative results for this application.

10.4 Layer-wise Image Vectorization↩︎

In our Layer-wise Image Vectorization (LIVE) experiments, we replace the cubic Bézier paths from the original LIVE framework [7] with our closed, filled NURBS curves. For the LIVE baseline, we retain all of its default hyperparameter settings with the sole exception of the initial circle radius, which we set to 16 pixels. Our method also uses an initial circle radius of 16 pixels, with our specific hyperparameter settings and optimization schedule detailed below.

Each of our curves is initialized as a circle using 12 key points and a degree of \(p{=}3\). For the optimization process, we train each path addition for 300 Adam iterations. We apply learning rate decay with a ratio of 0.4 over the course of the iterations. The initial learning rates are set to 2.0 for control-point positions, and 0.01 for colors, rational weights, knot intervals, and opacities. In addition to the standard reconstruction loss, our objective function incorporates a self-crossing penalty formulated in the original LIVE method with weight \(\lambda_\text{xing}{=}0.01\). Note that the third-order derivative smoothing loss is not activated in this experiment because the curve degree is \(p{=}3\).

To accelerate the differentiable filling process while preserving fine boundary details, we employ grid-step annealing. Within each 300-iteration round, the fill grid step \(h\) is cosine-annealed from \(h_{\max}{=}4.0\) to \(h_{\min}{=}1.0\) between \(10\%\) and \(80\%\) of the iterations (, \(t_b{=}0.1\) and \(t_e{=}0.8\)).

10.4.0.1 Discussion on Performance.

Our method achieves both higher efficiency and better reconstruction quality than the LIVE baseline. Readers may notice an apparent inconsistency: the results in 10.1 indicate that the Berio baseline (based on DiffVG [5]) is faster than our method, whereas in this experiment our approach outperforms the baseline (also based on DiffVG) in both speed and reconstruction quality. We believe that two factors contribute to this difference. (a) We employ a customized udf_weight loss computation tailored to our framework. Unlike the LIVE baseline, which evaluates the SDF on the CPU, our implementation performs this computation on the GPU, resulting in substantially higher efficiency. (b) Our step grid annealing strategy partially alleviates the excessive Gaussian generation caused by region filling, thereby improving both optimization efficiency and reconstruction performance.

10.5 Discussion on Bézier Splatting↩︎

Bézier Splatting (BS) [6] demonstrates excellent efficiency for vector graphic reconstruction and, to the best of our knowledge, remains the current state-of-the-art method in terms of computational performance. Nevertheless, its representation is less flexible than ours in several aspects. First, it focuses on vectorization tasks and is temporarily not well-suited for other applications. Second, its paired Bézier representation requires each filled region to consist of pairs of curves. Boundaries with an odd number of segments (, a triangle) therefore require an additional dummy point to complete the pairing, introducing an unnecessary degree of freedom that may adversely affect gradient-based optimization. Third, as discussed in 10.1, BS requires all curves in a scene to share the same geometric configuration, such as the number of control points, whereas our method allows heterogeneous curve representations within a single scene.

We have included a quantitative comparison of computational performance between our method and BS in 10.1, as well as a qualitative comparison of rendering quality in Sec. 5.1 of the main paper. However, we did not include quantitative comparisons on downstream tasks. The reasons are explained below.

Calligraphy reconstruction. The publicly released BS codebase requires all Bézier curves in a scene to share a fixed number of control points. Our calligraphy benchmark, however, contains strokes extracted via SLDVec [41] with varying numbers of control points per stroke, reflecting the inherent complexity differences among strokes within a single character. This constraint makes it infeasible to run BS on our dataset without substantial modification of its optimization pipeline. Moreover, BS operates on cubic Bézier curves only and does not support B-splines, varying stroke widths, or long splines—all of which are central to our calligraphy formulation. These architectural differences preclude an apples-to-apples comparison under identical experimental conditions.

Nonetheless, to provide a partial reference, we note that an early prototype of our framework employed the same contour-sampling strategy as BS (, anisotropic Gaussians aligned to curve tangents). We test this prototype on our calligraphy benchmark and the results are included in 5. As shown in the table, while the anisotropic sampling strategy achieves slightly faster optimization due to fewer splats, it significantly degrades reconstruction quality, with lower PSNR and higher Hausdorff distances compared to our final formulation. Furthermore, as shown in 10, the BS prototype completely fails to preserve the intricate stroke boundaries and varying widths of the original calligraphy.

Layer-wise image vectorization. The BS authors have not released their code for the layer-wise vectorization application described in their paper. Consequently, no compatible implementation is available for us to evaluate on our LIVE benchmark.

Summary. Our qualitative comparison in the main paper, alongside the quantitative evaluation against the BS prototype (BS Proto.) in this supplement, demonstrates clear rendering and reconstruction advantages of isotropic over anisotropic Gaussian sampling. While a direct quantitative comparison with the official BS codebase on downstream tasks is precluded by practical incompatibilities, our proxy experiment confirms that anisotropic sampling significantly degrades reconstruction fidelity. We emphasize that our primary comparisons—against Berio  [4] for calligraphy and against LIVE [7] for vectorization—are conducted under fair, controlled conditions with identical initialization and training iterations.

References↩︎

---

1. Corresponding author.↩︎

2. \(W,H\) are the width and height of the image.↩︎

3. \(\delta_\text{c}=10\) is sufficient to capture curve details and reduce blurry edges in our tests. However, we use \(\delta_\text{c}=18\) for better quality in calligraphy reconstruction.↩︎