Inverse Diffusion Curves using Shape Optimization

The inverse diffusion curve problem focuses on automatic creation of diffusion curve images that resemble user provided color fields. This problem is challenging since the 1D curves have a nonlinear and global impact on resulting color fields via a partial differential equation (PDE). We introduce a new approach complementary to previous methods by optimizing curve geometry. In particular, we propose a novel iterative algorithm based on the theory of shape derivatives. The resulting diffusion curves are clean and well-shaped, and the final image closely approximates the input. Our method provides a user-controlled parameter to regularize curve complexity, and generalizes to handle input color fields represented in a variety of formats.


Introduction
Vector graphic images remain invaluable for a broad range of 2D applications because of their resolution independence, compactness of representation, and powerful editability. Recently, diffusion curve images  further improve the expressiveness of vector graphics by providing flexible and easyto-manipulate smooth gradients, and since then inspired a variety of novel applications [Jeschke et al. 2009;Takayama et al. 2010;Sun et al. 2012;Ilbery et al. 2013;Sun et al. 2014]. Defined along the curves, colors are diffused across the image by a Poisson or Laplace reconstruction, and their smoothness can be further controlled by the curve's blurriness through post-processing. Although efficient rendering of diffusion curve images has been well explored, its inverse problem of creating diffusion curves automatically given desired target images remains challenging.
The inverse diffusion curve problem is difficult because even though the curves themselves are 1D, their impact on the final image is nonlinear and global over a 2D domain through a partial differential equation (PDE), Laplace's equation. The geometry of curves largely determines the reconstruction quality. Previous methods have used local heuristics to obtain curve geometry. They place curves at locations indicated by edge detectors applied to the target image ] and its Laplacian or bi-Laplacian [Xie et al. 2014]. While these heuristics work well around sharp edges such as object boundaries, they have difficulty handling color variations in regions that are smooth yet visually rich.
We introduce a new approach to the inverse diffusion curve problem. Complementary to existing methods, our approach solves for curve geometry through a global optimization that takes into account the curves' full impact on the PDE-based color field. To achieve this, we characterize how modifications to a diffusion curve can reduce a global cost function determined by the solution of a PDE with the curves acting as boundary conditions.
Our method is grounded on the theory of shape optimization [Sokolowski and Zolésio 1992]. Given a color field, it computes curve geometry by minimizing a measure of the color reconstruction residual. Starting from an initial set of curves, it iteratively evolves their shapes toward an optimal configuration. Mathematically, the curves are treated as continuous functionals. This allows them to deform arbitrarily, enabling full exploration of possible curve configurations. This iterative process is similar to a surface normal flow: our curve evolution at every iteration is guided by the Fréchet derivative of the residual function with respect to the curve's boundary velocity, leading to an efficient "gradient descent" of the residual. This method is mathematically clean and easy to implement: all computations at each iteration boil down to solving a Laplace equation and a Poisson equation.
Building on our curve placement algorithm, we introduce a complete pipeline for solving the inverse diffusion curve problem. Our method generates curves in a clean and concise way, and the resulting images can accurately capture complex color variations of input color fields (see Figures 1, 14, 15, and 16 as well as the supplementary images).
We demonstrate that our method promises practical applications beyond pixel image vectorization. For instance, it enables automatic rendering of vector graphic images from 3D geometries, analogous to the traditional pixel image rendering. Further, using our algorithm, one can directly transform other formats of vector graphics, such as gradient meshes, into diffusion curve images without rasterizing the input.

Figure 1:
We present a new approach to automatically build diffusion curve images approximating provided color fields (a). Starting from a set of boundary curves (black strokes) indicating color jump discontinuities (b), our method iteratively adds curves (dark yellow strokes) and refines their shapes in an optimized manner (c, d). The resulting image (e) accurately matches the input.

Introduction
Vector graphic images remain invaluable for a broad range of 2D applications because of their resolution independence, compactness of representation, and powerful editability. Recently, diffusion curve images  further improve the expressiveness of vector graphics by providing flexible and easyto-manipulate smooth gradients, and since then inspired a variety of novel applications [Jeschke et al. 2009;Takayama et al. 2010;Sun et al. 2012;Ilbery et al. 2013;Sun et al. 2014]. Defined along the curves, colors are diffused across the image by a Poisson or Laplace reconstruction, and their smoothness can be further controlled by the curve's blurriness through post-processing. Although efficient rendering of diffusion curve images has been well explored, its inverse problem of creating diffusion curves automatically given desired target images remains challenging.
The inverse diffusion curve problem is difficult because even though the curves themselves are 1D, their impact on the final image is nonlinear and global over a 2D domain through a partial differential equation (PDE), Laplace's equation. The geometry of curves largely determines the reconstruction quality. Previous methods have used local heuristics to obtain curve geometry. They place curves at locations indicated by edge detectors applied to the target image ] and its Laplacian or bi-Laplacian [Xie et al. 2014]. While these heuristics work well around sharp edges such as object boundaries, they have difficulty handling color variations in regions that are smooth yet visually rich.
We introduce a new approach to the inverse diffusion curve problem. Complementary to existing methods, our approach solves for curve geometry through a global optimization that takes into account the curves' full impact on the PDE-based color field. To achieve this, we characterize how modifications to a diffusion curve can reduce a global cost function determined by the solution of a PDE with the curves acting as boundary conditions.
Our method is grounded on the theory of shape optimization [Sokolowski and Zolésio 1992]. Given a color field, it computes curve geometry by minimizing a measure of the color reconstruction residual. Starting from an initial set of curves, it iteratively evolves their shapes toward an optimal configuration. Mathematically, the curves are treated as continuous functionals. This allows them to deform arbitrarily, enabling full exploration of possible curve configurations. This iterative process is similar to a surface normal flow: our curve evolution at every iteration is guided by the Fréchet derivative of the residual function with respect to the curve's boundary velocity, leading to an efficient "gradient descent" of the residual. This method is mathematically clean and easy to implement: all computations at each iteration boil down to solving a Laplace equation and a Poisson equation.
Building on our curve placement algorithm, we introduce a complete pipeline for solving the inverse diffusion curve problem. Our method generates curves in a clean and concise way, and the resulting images can accurately capture complex color variations of input color fields (see Figures 1, 14, 15, and 16 as well as the supplementary images).
We demonstrate that our method promises practical applications beyond pixel image vectorization. For instance, it enables automatic rendering of vector graphic images from 3D geometries, analogous to the traditional pixel image rendering. Further, using our algorithm, one can directly transform other formats of vector graphics, such as gradient meshes, into diffusion curve images without rasterizing the input.

Related Work
Diffusion curves ] represent a color field by diffusing the colors defined along control curves over the entire image plane. The diffusion process is described by Laplace's equa-tion solved using a finite volume method. Later, solving Laplace's equation was improved using a multigrid method [Jeschke et al. 2009], triangle mesh interpolation [Pang et al. 2012], Boundary Element method [Sun et al. 2012], 2D ray tracing [Prévost et al. 2014], and Fast Multipole method ].
Inverse diffusion curve problem. Our work focuses on the inverse problem of diffusion curves. Previously, Orzan et al. [2008] proposed to place diffusion curves along edges extracted from input images using the Canny detector [Canny 1986]. Jeschke et al. [2011] introduced a technique to improve curve colorings. Xie et al. [2014] further improved this method by detecting edges in a Laplacian (and/or bi-Laplacian) domain and constructing curves hierarchically. They solve the Laplacian and bi-Laplacian weights using least-squares fitting. In all methods, diffusion curves are placed along the detected edges, and never moved or added in continuous color regions. These methods then rely on optimizing curve coloring for better accuracy.
We introduce a fundamentally complementary solution to the inverse diffusion curve problem. Instead of predetermining curve geometry and optimizing their coloring, we propose doing the opposite by first optimizing the geometry and then determining the coloring accordingly. We demonstrate that with a very simple coloring scheme, our method outperforms prior methods under many situations ( §6.2). Furthermore, our approach accepts input color fields beyond pixel images.
Extensions of diffusion curves. Several methods have been proposed to extend the expressiveness of diffusion curves. Sun et al. [2014] enabled fast diffusion curve cloning and multi-layer composition. Finch et al. [2011] introduced a higher-order notion of smoothness: the colors are defined using a 4th-order linear elliptic PDE rather than a Laplace equation. To accelerate the color evaluation, Boyé et al. [2012] developed a vectorial solver using the Finite Element Method, and Sun et al. [2012] proposed a boundary element based formula, which was later improved in [Ilbery et al. 2013] to handle both Laplacian and bi-Laplacian curves in a unified framework. Higher-order curves offer greater flexibility than the standard diffusion curves, but their inverse problems are more difficult and remain unsolved. In this paper, we focus on the inverse problem for original diffusion curves and discuss potential extension to higher-order domains in §6.3.
Theory and applications of shape optimization. We build our curve optimization on the theoretical foundation of shape optimization [Sokolowski and Zolésio 1992;Haslinger et al. 2003], a subfield of optimal control theory. Mathematically, it solves the problem of finding a bounded set Ω to minimize a continuous functional on Ω. The core idea of shape optimization has been used for image segmentation since the seminal work of [Kass et al. 1988;Mumford and Shah 1989]. It is also related to surface gradient flow widely studied in geometry processing [Schneider and Kobbelt 2001;Crane et al. 2013]. In areas outside of computer graphics, shape optimization has been used to enhance mechanical structures such as airfoils [Mohammadi et al. 2001] and photonic crystals [Burger et al. 2004]. It has also been used in computer vision for image segmentation (e.g., [Herbulot et al. 2006;Jung et al. 2012]). To our knowledge, shape optimization has not yet been applied in vector graphics. In this paper, we solve a shape optimization problem with a PDE constraint ( §5.1), which is significantly more challenging than a conventional shape optimization problem.

Background and Overview
We start by briefly revisiting the mathematical formulation of diffusion curve images. We then present the main focus of this work, the inverse diffusion curve problem, and overview our proposed solution.
Diffusion curve images. In a diffusion curve image, as originally formulated in Jeschke et al. 2009], the color field u is a harmonic function, satisfying a Laplace equation with a Dirichlet boundary condition: where the boundary consists of the entire set of diffusion curves; C and C r specify the colors on the left and right side of each curve, respectively. Typically, both the shapes of the curves and their left-and right-side colors are specified by the user, and the entire color field is uniquely determined by solving the Laplace equation (1).
Since its invention, diffusion curves have been augmented. Orzan et al. [2008] proposed to apply per-pixel blurring to the rasterized image of u, the solution of (1). Finch et al. [2011] further extended to diffuse colors using higher-order elliptic PDEs such as the biharmonic equations.
Inverse diffusion curve problem. While plenty of extensions of the forward diffusion curve problem have been proposed, largely underexplored is the inverse problem, one that computes a set of diffusion curves such that the resulting vector image closely resembles a user-provided 2D color field. In this paper, we address this inverse problem. In particular, we note that the inverse problem involves two subproblems: • Curve geometry. To build a diffusion curve image, one needs to decide where to place the curves (namely, to determine ).
• Curve coloring. Given the curve geometry, the colors on both sides of each curve (namely C and C r ) need to be specified.
As discussed in §2, recent work [Jeschke et al. 2011;Xie et al. 2014] has largely focused on optimizing curve coloring with their geometry predetermined (using edge detection). In contrast, we focus on a complementary problem, the problem of directly optimizing curve geometry. We demonstrate (in §6) that curves with optimized geometries generally yield higher-quality of reconstructions, regardless of the curve coloring schemes.

Overview
Pipeline. We develop a complete pipeline for automatic creation of diffusion curve images ( §4). Figure 2 shows an overview of our pipeline. We take as input a color field I allowing to query for color values at for all x ∈ Ω (where Ω denotes the image domain). Starting with extracting a set of boundary curves ( §4.1) indicating jump discontinuities in I, our method generates a set of curves as "initial guesses" ( §4.2) which are then deformed by our core curve optimization algorithm ( §5) to minimize reconstruction error ( §4.3). Lastly, we post-process the deformed curves ( §4.4) to generate final diffusion curve images. (Section 4.1) Figure 2: Our pipeline. The input to our method is a 2D color field. After obtaining a set of boundary curves indicating color jump discontinuities ( §4.1), our method constructs a set of initial curves ( §4.2) and optimizes their shapes and trajectories ( §4.3 and §5). Finally, we post-process the optimized curves and obtain the resulting diffusion curve image ( §4.4).
Curve optimization. As one of our main contributions, the key component of our pipeline is a curve optimization algorithm that deforms diffusion curves to minimize the reconstruction error. In our algorithm, curves are 1D continuous geometries discretized as polylines, while the cost functional includes an integral over a 2D image region. The optimization problem is constrained by a Laplace equation for color diffusion. Our general idea of solving this optimization problem is using iteration: given a set of curves 0 , we construct a velocity field v so that after deforming 0 according to v for a small time step t, the updated set of curves t leads to improved reconstruction accuracy. This process is repeated until convergence. Conceptually, this routine is similar to gradient descent-in each iteration, a local adjustment is made to obtain a better solution. Our main challenge is to derive a proper form for the derivative of the cost function with respect to the shapes of the curves and use it to construct the velocity field v.
In what follows, we will first provide a complete exposition of our end-to-end pipeline (in §4 and outlined in Figure 2), and then describe in details the algorithms and mathematical derivations for the curve optimization (in §5).

Boundary Curves
Provided an input color field I, we start the pipeline by obtaining a set of boundary curves ∂ Ω indicating the outer boundary and jump discontinuities of I. 1 An example color field and corresponding boundary curves are shown in Figure 3-ab. In practice, we obtain the boundary curves ∂ Ω depending on specific representation of the input color field I: • Pixel images. A common way to represent color fields is using standard pixel images. The boundary curves, however, are not uniquely defined in this case. To obtain these curves in practice, we use Canny edge detection similar to ].
• 3D renderings. If the color field is defined by the rendering of a 3D scene, the boundary curves can be obtained by extracting object contours.
• Other vector formats. For input color fields represented in other vector formats (e.g., gradient mesh), ∂ Ω can be determined directly based on the underlying vector representation (e.g., triangle edges).
Please refer to §6.2 for more details and experimental results on boundary curve computation. 1 The necessity of boundary curves is explained in §5.

Curve Initialization
Desired properties. Similar to gradient descent methods, our curve optimization algorithm takes an initial guess to start with. For ensuring high-quality optimization results, there are a few properties required of the initial curves: 1. Easy to compute. The curve initialization step should not require intense computation: we rely on the optimization step to refine the shapes of these curves.
2. Good coverage. The initial curves should provide a good coverage to the full image domain Ω, so that the optimization is less prone to local optima.
3. Being well-shaped. The initial curves need to be well-shaped. For example, they should have low complexities and not selfintersect or collide with the boundary curves.
To achieve these properties, we use iso-contours of the residual field as the initial curves: In (2), u 0 is given by the diffusion curve image using only the input boundary curves ∂ Ω. These curves can be computed easily from a set of iso-values (Property 1). In addition, as long as the iso-values are distributed properly, the resulting iso-contours will provide a good coverage to the image domain while being well shaped (Properties 2 and 3). For example, to make the initial curves never intersecting with ∂ Ω, we can simply pick strictly positive iso-values as R 0 (∂ Ω; x) = 0 for all x ∈ ∂ Ω.
Our approach. To choose a set of properly distributed isovalues, we start with sampling a set of points in Ω and stacking the residual values (2) at these points into a vector r 0 in ascending order (lines 2 and 3 in Algorithm 1). The resulting vector r 0 provides a picture on the distribution of residuals. We adopt two complementary schemes, global and local, to set iso-values using r 0 , and thereby obtain the iso-contours:

Algorithm 1 Diffusion curve initialization
Require: Color field I (defined on Ω) and boundary curves ∂ Ω 1: procedure CURVEINIT(scheme, ∂ Ω, I, Ω) return iso-contours with iso-values specified in A 11: end procedure Figure 4: Fitting r 0 with 9 components (indicated with purple dots) using a piecewise linear functionr with 2 pieces (m = 1). The value ofr at its internal boundary is selected as iso-value.
• Global. The global scheme constructs a relatively large set of initial curves over the entire domain Ω. Assume that the number of iso-values m is given. Ideally, we would like to find m values such that the consequent iso-contours optimally capture the structure of 2D residual field R 0 . In practice, we solve this problem approximately and rely on our curve optimization algorithm to refine the curves. Particularly, we solve a well-studied 1D problem [Ramer 1972]: to fit a piecewise linear functionr with m+1 pieces that closely describes r 0 (interpreted as a polyline). Then, the values ofr at its m internal piece boundaries are used as iso-values (see Figure 4).
• Local. The local scheme, in contrast to the global one, adds curves locally in regions with high approximation error. In this case, we use only one iso-value determined based on the maximal sampled residual (line 8 of Algorithm 1).
In our curve placement algorithm (detailed in §4.3), we use the global scheme at the beginning to ensure that the initial curves provide a good coverage to the domain Ω (Property 2). Then, the local scheme is applied iteratively to add small sets of curves in high-residual areas. The combination of both schemes offers sufficient approximation accuracy without introducing unnecessarily complex curves (Property 3). We find that m = 2 works well in our experiments.

Curve Placement
Given the initial curves generated by Algorithm 1, our curve optimization algorithm iteratively refines their trajectories to reduce reconstruction errors and finalize curve geometry. We postpone the details of this algorithm (Algorithm 3) and its derivations until §5, but present here the complete curve placement steps (Algorithm 2).
The curve placement process is built on the algorithms of curve initialization and optimization schemes. It takes as input the target color field I defined on domain Ω, the previously obtained boundary curves ∂ Ω, and a tolerance ε 0 on reconstruction error. Based on ∂ Ω, we partition the domain Ω into a number of connected components and process them individually in parallel (line 2 of Algorithm 2). Figure 5 illustrates example boundaries and resulting partitionings. Notice that our approach allows boundary curves to exist inside individual components (e.g., Fig-(  ure 5-b)-these curves will remain fixed throughout the entire pipeline.
For each connected component C, our approach generates diffusion curves via several passes, each of which involves initializing a set of curves (Algorithm 1) and optimizing their shapes (Algorithm 3 in §5). In the first pass, we start with initial curves constructed using the global scheme (lines 5 and 6). After this pass, if the approximation error remains beyond a tolerance ε 0 , additional passes are used in which new curves are initialized using the local scheme (lines 8 to 10). After the error drops below the threshold, we perform a final pass (line 12) in which all curves created in previous passes are optimized together. Finally, we post-process the resulting curves to remove redundant curve segments (line 13 and §4.4). An example of this curve placement process is illustrated in Figure 6.

Curve Post-Processing
Lastly, we post-process the curves returned by Algorithm 2 and generate the final diffusion curve image.
Curve Coloring. Notice that Algorithm 2 returns optimized curve geometry instead of actual diffusion curves. Thus, to turn into a set of diffusion curves, their coloring, namely colors on both sides of each curve, needs to be provided. This corresponds to specifying the values of C and C r in (1).

Algorithm 2 Diffusion curve placement
Require: Color field I (defined on Ω) and boundary curves ∂ Ω 1: procedure CURVEPLACEMENT(∂ Ω, I, Ω, ε 0 ) 2: partition Ω into connected components 3: for each component C do As aforementioned, this curve coloring step is completely orthogonal and complementary to our core technique (Algorithm 2). Thus, in the rest of this paper, we use a simple scheme which directly sample color values on both sides of each curve from the input color field I. That is, for any x ∈ , we set where n and n r respectively denote normal directions pointing left and right side of a point x on a curve (thus, n l = −n r ) and δ is a small positive number that can be set to the size of one pixel when I is represented as a pixel image. Our experiments demonstrate that this simple scheme can yield high-quality results thanks to our optimized curve geometry ( §6.2). In §6.3, we show that more advanced coloring techniques can further improve reconstruction accuracy.
Removing redundant curve segments. As mentioned in §5.4, we represent diffusion curves as polylines consisting of a number of line segments. Some of these segments, however, may be unnecessary. Note that the colors across a line segment are continuous because of the boundary condition (5) on . If the color gradient normal to a segment is also continuous across, then the segment as a boundary has no influence on the solution color field u. A mathematical explanation is in §3 of the supplementary document. Precisely, a normal gradient is continuous when is zero. In practice, we solve u using the Finite Element Method ( §5.4) and check if |d n (x)| at the center point x of each segment is below a threshold. If so, we mark the segment as unnecessary.
Lastly, for each curve output by Algorithm 2, we remove a largest set of connected redundant segments to avoid breaking the curve into many small disconnected components.
To transform the final polyline into a standard diffusion curve made from end-to-end connected Bézier curves, we adopt the Potrace algorithm [Selinger 2003], which was also used in ].
Per-pixel blurring (optional). The curves placed in a smooth color region have continuous color values across the curves. However, since these curves serve as boundaries in the Laplace solve, color gradients may not necessarily remain continuous across curves generated by Algorithm 2. Such gradient discontinuities can sometimes lead to noticeable artifacts [Finch et al. 2011]. Thus, our pipeline includes an optional step following the original framework of diffusion curves ] to perform per-pixel blurring on the rasterized image. The size of blur kernel at each pixel is determined by another Laplace equation: where K 0 (x) gives the desired kernel size along the curves. In particular, we set K 0 (x) = 0 for all x ∈ ∂ Ω since the boundaries and discontinuities in the input color field should never be blurred. For x ∈ , the value d n (x) indicates the magnitude of the gradient domain discontinuity. Thus, we set K 0 (x) = b |d n (x)| a for all x ∈ , where a and b are two global parameters. In our implementation, we set a = 0.2 and b to 5% of the longest axis of Ω's bounding box.
Notice that more advanced curve coloring techniques, such as [Xie et al. 2014], may optimize color gradients across the curves, largely removing gradient discontinuity artifacts. In this case, per-pixel blurring is unnecessary (see §6.3).

Color value (input)
Domain outer boundary (input) Diffusion curves (input/output) Figure 7: Input and output of our curve optimization algorithm. Input: color field I, domain Ω and its outer boundary ∂ Ω, initial diffusion curves ; output: refined curves .

Diffusion Curve Optimization
We now detail the core of our pipeline, the optimization of diffusion curve geometries to approximate a given 2D color field. We first describe an algorithm minimizing the approximation error of diffusion curve images ( §5.1-5.2), and then extend it to balance accuracy against curve length ( §5.3). Lastly, we provide implementation details ( §5.4), followed by the discussions of further extensions ( §5.5).
We introduce Shape Optimization [Sokolowski and Zolésio 1992] to formulate the inverse diffusion curve problem. While building our approach on existing shape optimization concepts and theories ( §5.2), we also develop a new formula for regularizing curve length ( §5.3). Please refer to the supplementary document for complete derivations and a review of related background.
Curve Optimization in a Nutshell. The major steps of our approach are outlined in Algorithm 3. Its input includes the color field I, a 2D closed domain Ω over which I is defined, and a set of initial curves in Ω (Figure 7). In this section, the color field I is treated as a black box, allowing I(x) and ∇I(x) to be evaluated for any x ∈ Ω. Our curve optimization algorithm then iteratively refines the curves by changing their shapes (i.e., the trajectories) and topologies to obtain better approximation. The resulting diffusion curve image consists of the optimized curves and the domain boundary ∂ Ω. During the optimization process, the colors along both sides of these curves (i.e., C and C r of the Laplace equation (1)) are sampled from the given color filed I, and the approximated color value u(x) for all x ∈ Ω is determined according to the equation (1) with the Dirichlet boundary condition u(x) = I(x), ∀x ∈ ∪ ∂ Ω.
We note that rather than sampling color values along the curves, prior methods [Jeschke et al. 2011;Xie et al. 2014] post-optimize color values after the curves are determined. We will discuss the extension of our method to incorporate their post-optimization later ( §5.5) and examine it in our experiments ( §6.3).

PDE-Constrained Optimization Problem
Formally, our iterative curve optimization process minimizes a cost functional defined as the L 2 residual of the color approximation, Here u is the color field determined by diffusion curves. We write as a parameter of R to emphasize the dependence of the residual on through the Dirichlet boundary condition (5). Since u is the solution of the Laplace equation (1), we are concerned with an optimization problem with a PDE constraint, min R(Ω; ) s.t. u satisfies the Laplace eqn. .
PDE-constrained optimization problems are known to be challenging in general [Pinnau and Ulbrich 2008]. In our problem (7),

After timestep
Higher cost Lower cost Figure 8: Curve optimization. Given a set of curves 0 , we construct a velocity field v so that if one deforms 0 according to v, the resulting curves t provide a lower cost.
the optimization variables are the shapes of diffusion curves, that is, the spatial trajectories and topologies of the curves. Ideally, a curve can have an arbitrarily continuous trajectory, and therefore needs to be represented using a continuous functional rather than using individual and discrete parameters. More importantly, the error residual R(Ω; ) depends on the optimization variables (the curves) through the Laplace equation (1) in a complex manner: any local change to the curves has a global impact, one that changes u over the entire domain Ω, which further affects the residual via (6).

Gradient-Descent Solver
We propose a new approach for solving the curve optimization problem (7), following the general spirit of gradient descent.
Starting from a set of initial curves, our approach iteratively decreases the residual (6) by adjusting their shapes. Throughout, a fundamental difficulty we need to address is the computation of the residual's "gradient" with respect to the shapes of the curves, as the conventional gradient in terms of continuous curves is undefined.
We develop our method from the perspective of functional analysis: in each gradient-descent step, we first construct a velocity field v on the curves, specifying v(x) for all x ∈ (Figure 8). We then use v(x) to deform the curves, analogous to a (2D) surface flow in geometry processing [Sethian et al. 2003;Brakke 1992]. In other words, we evolve the curves via a single step of the forward Euler method of integratingẋ = v(x), ∀x ∈ .
In this subsection, we present the details of computing such a v that after deforming the curves accordingly, the residual is guaranteed to decrease (Lines 5-7 of Algorithm 3). Briefly speaking, we will first assume that v(x) is known and analytically express how much the residual would change if the curve is deformed according to v. This analytical expression allows us to formulate the condition of v resulting in a decrease of the residual, and thereby provides us a recipe for computing v.

Algorithm 3 Gradient-descent diffusion curve optimization
Require: initial curves 0 , color field I on Ω with boundary ∂ Ω 1: procedure CURVEOPT( 0 , ∂ Ω, I, Ω) compute v n (x) = −B R (x) using (12) 8: forward-Euler curve advancement, ← + v n t 9: evaluate R using (6), and update its change ∆R 10: end while 11: return current curves 12: end procedure Fréchet derivative as a linear form. Given a domain Ω and a set of initial curves 0 , we consider a general cost functional, where y is continuous on Ω and may depend on the choice of 0 . Our residual (6) takes the form y(x; 0 ) = 1 2 (u(x) − I(x)) 2 and depends on 0 via the Laplace solution u. Assuming a known v, we introduce the Fréchet derivative [Coleman 2012] of C with respect to v. Let t denote the curves evolved according to v after an infinitesimal time period of t, that is, x → x + v(x) t for all x ∈ 0 (Figure 8). The Fréchet derivative of C is a linear form of v satisfying that dC(Ω; 0 ) = lim Conceptually, this derivative measures how quickly the cost functional C changes as we deform the curves using v infinitesimally. According to Hadamard-Zeloésio Structure Theorem [Delfour and Zolésio 2011], such a linear form always exists when Ω, 0 and v are sufficiently regular, which is usually the case in practice. For our cost functional (6), we further reduce the Fréchet derivative into a linear form expressed as a boundary integral where v n (x) := v(x) · n(x) denotes the normal velocity on the curves (Figure 8), Γ 0 = ∂ Ω∪ 0 includes both the domain boundary ∂ Ω and all the inner curves 0 (see Figure 7), and B is another function independent from v but related to the specific integrand y. In the rest of this subsection, we aim to derive a formula to evaluate B(x) for any x ∈ 0 .
Once B is known, setting guarantees a negative derivative value in (9) (assuming that B does not vanish everywhere on 0 ). This provides a formula of constructing v n , which we then apply to deform the curves 0 .
With a sufficiently small timestep size t, the deformed curves t , computed by x + v n (x) t for all x ∈ 0 , is guaranteed by construction to yield a smaller residual value and thus a better approximation of I. Computational Recipe. Shape Optimization Theory has provided a simple recipe of computing B for our particular cost functional (6). Here we simply present the formulas. Please see Appendix A for an outline of the derivation and §3 of the supplementary document for more details.
We first solve the Laplace equation (1) to compute u(x), which in turn allows us to construct a Poisson equation with a Dirichlet boundary condition, Next, the solution p of this equation, together with u, allows the computation of B(x) in a simple form Combining (12) and (10) computes the normal velocity v n , the velocity that can deform the curves 0 and lead to a decrease of the approximation residual (6). This computation is performed at each gradient-decent step, and the optimization process stops when the residual change drops below a threshold ε. Figure 10 illustrates the optimization process with synthetic examples.

Regularizing Curve Complexity
So far, our optimization problem (7) focuses solely on minimizing the L 2 residual (6). However, because the L 2 error along a curve is always zero due to the boundary condition (5), one simple way to yield a very low residual is to use space-filling curves. Indeed, if we start with one curve in a complex color region, it becomes zigzag after running the optimization for many iterations (Figure 9-a). While the numerical residual is low for such curves, their largely increased geometric complexity may be undesirable for certain applications (such as vector graphics editing). Thus, we propose an extension to the cost functional (6), providing users the flexibility to trade approximation accuracy for simpler curves. To this end, we add a regularization term to (6) to penalize the total length of the curves: where α is a user-specified scalar controlling the strength of regularization. It can be shown that similar to (9), the Fréchet derivative of the second term is also a linear form of v n . Let where κ(x) measures the curvature of a point x on the curves. This formula has been used to derive the mean curvature flow [Mantegazza 2011] in geometry processing. It is also a special case of the Fréchet derivative of a general boundary integral (see §1.3 of the supplementary document). Following the derivation of B in §5.2, we obtain the normal velocity for decreas-ingR(Ω; ), that is, v n (x) = −B R (x) − ακ(x). With this slightly different velocity formula, the entire optimization algorithm remains the same as before. In addition, the user is able to control the complexity of resulting curves by adjusting the strength of regularization (Figure 9).

Implementation Details
We now present implementation details of Algorithm 3, wherein two major steps are solving the Laplace equation (1) and the Poisson's equation (11). Both PDEs have Dirichlet boundary conditions defined on the boundary of Ω and the optimized curves (recall (5)). Since we also need to evaluate the domain integral over Ω during the iterations (Line 9 of Algorithm 3), we triangulate the entire domain of Ω and use the Finite Element Method [Zienkiewicz and Morice 1971] for both solves, while other numerical solvers (e.g., the Boundary Element Method) could also be applied.
Finite element discretization. We discretize the boundary and optimized curves into piecewise linear segments, and represent them using polylines. The velocity v n is discretized and stored at every vertex along the polylines. We use the package Triangle [Shewchuk 1996] to triangulate the domain Ω (Line 4 of Algorithm 3). The resulting triangle mesh is then used in the finite element solves. The computation of curves' normal velocity in (12) involves boundary normal derivatives of the finite element solutions (i.e., ∂ p/∂ n and ∂ u/∂ n). We choose the secondorder finite element basis, as it offers higher accuracy especially near the boundary (see §4 of the supplementary document).
Curve tracking. Advancing the curves using the computed normal velocity (i.e., computing t given 0 and v n in Figure 8) is a typical yet nontrivial surface tracking problem. We use a recently developed explicit tracking approach [Brochu and Bridson 2009], which advances the vertices on curve polylines using explicit forward Euler method, and then carefully remeshes the polylines to ensure correct topology changes and a collision-free state.
Timestep size t. To ensure robust curve tracking, we dynamically set the timestep size t for the forward-Euler curve advancement (Line 8 in Algoirthm 3). We start with choosing a t value such that the vertex displacement v n (x) t, ∀x ∈ would not collapse any polyline segment on . This ensures that possible topology changes can be robustly processed. From this starting value, we iteratively halve t until the residual value (after a step of curve deformation) decreases.

Discussions
Measuring geometric complexities. In Equation (13) and the rest of this paper, we use the total length of all diffusion curves to measure their geometric complexity. Depending on specific applications, there may exist other metrics more suitable to user needs. As an example, in §2 of the supplementary document, we discuss another possible measure which can also be incorporated in our curve optimization framework. Coloring schemes. As described at the beginning of this section, given the curve geometry , we specify colors on both sides of each curve by directly sampling color values from the input color field I. Alternatively, prior work [Jeschke et al. 2011;Xie et al. 2014] propose to post-optimize curve colors for better reconstruction accuracy. Our method can easily adopt this approach, post-optimizing the colors after the curves are optimized. We implemented this approach and present the results in §6.3.
Higher-order domains and curves. While our approach focuses on solving the inverse problem of the standard (first-order) diffusion curves, it can be also applied to higher-order domains. For instance, as demonstrated in §6.3, we can feed ∇I instead of I to Algorithm 3 to compute curves offering a higher order of smoothness.
In principle, it is possible to generalize (6) with u(x) directly given by higher-order (e.g., biharmonic diffusion) curves. However, this would dramatically complicate the form of the Fréchet

Figure 10: Synthetic validation of our diffusion curve optimization algorithm. The left-most column contains input color fields given by one closed curve (top) and one open curve (bottom). a varying number of iterations. Our resulting curves precisely match the original.
derivate (9) and thus that of the velocity field (10), causing them significantly more difficult to evaluate numerically. Therefore, we leave the extension of (6) to higher-order curves as future work.

Results
We first ( §6.1) show experimental results demonstrating the validity of our curve optimization algorithm as well as how the regularization behaves in practice. Then, in §6.2, we show reconstructed diffusion curve images using input color fields represented in three forms: pixel image, 3D renderings, and gradient meshes. In addition, we show preliminary results motivating possible future applications.

Experimental Results
Synthetic validations. We design two synthetic tests (Figure 10) to validate our diffusion curve optimization algorithm (Algorithm 3). In these tests, the input color fields I are themselves diffusion curve images with continuous colors. In this case, the optimal is simply the set of diffusion curves used to generate I. Although the shape optimization problem is in general non-convex, our method successfully finds the optimal solutions for both closed and open curves. Please see the accompanying video for curve deformation animations.
Regularizing curve complexity. As discussed in §5.3, our method is able to balance resulting accuracy and curve complexity by varying the strength α of regularization in (13). Figure 11 shows how α influences resulting curves generated with our full pipeline (Algorithm 2). Figure 11-a has simpler curves (due to greater α), but some highlights at the bottom-left of the image are absent. Figure 11-b, on the other hand, provides lower approximation error but at the cost of greater curve complexity vious edge detection based methods Xie et al. 2014]. 2 The parameters for each method are selected such that the resulting curves have approximately identical complexities. Our method outperforms previous ones when handling smoothly varying color fields. Notice that, in the bottom two rows (i.e., Flower and Jade), [Xie et al. 2014] has slightly higher approximation errors (measured in RMSE), because it requires higher curve complexities to work properly in these cases. Please see the supplemental material for more comparisons.
3D renderings. Another kind of color field common to computer graphics applications is renderings of 3D scenes. Our approach can be used to approximate these color fields with diffusion curve images. In this case, we represent I as high-resolution pixel images and obtain the boundary curves ∂ Ω directly using object contours. 3 Figure 15 demonstrates diffusion curve images generated using our method from 3D renderings, where slightly higher curve complexities (compared to Figure 14) are permitted to ensure low reconstruction errors. Our method successfully captures detailed appearances: from glossy surfaces to smooth shadow boundaries.

Gradient meshes.
Our approach can also generate diffusion curve images directly from input color fields I represented using other vector formats. We demonstrate this using input color fields represented as gradient meshes in SVG format [Bah 2011] where each mesh grid is a Coons Patch [Coons 1967]. In this case, the color I(x) and gradient ∇I(x) can be evaluated analytically for any x ∈ Ω, and the boundary curves ∂ Ω are simply the  (resulting from a lower α). In all our results, the complexity is numerically defined as total curve length normalized, so that the longest axis of each image's bounding box has unit length.

Main Results
We now show diffusion curve images generated using our method (Algorithm 2). All our results utilize the per-pixel blurring described in §4.4. Please refer to the supplemental material for unblurred versions.
Theoretically, our approach does not require the input color field I to have any particular representation or discretization. In practice, we demonstrate such flexibility using three types of input: pixel images, 3D renderings, and gradient meshes. The execution time for generating each of these results is summarized in Table 1. The supplemental video contains animations demonstrating the creation process for these results.
Pixel images. One common way to represent a color field I is to use standard pixel images. In this case, I(x) is evaluated using bilinear interpolation, and ∇I(x) using finite difference. As stated in §4.1, we perform edge detection to obtain the boundary curves Conclusion. This paper introduces a novel solution to the inverse diffusion curve problem. The key component of our approach is a curve optimization algorithm that iteratively deforms a set of diffusion curves in a way that the reduction of approximation error is guaranteed. Based upon the core algorithm, we develop a full pipeline that takes an input color field plus a set of boundary curves and produces an image with well-shaped and clean curves that closely matches the input. Our approach offers the generality to take input presented in different formats, which we demonstrate using three types: pixel images, 3D renderings, and gradient meshes.

A Brief Derivation of Equation 12
We now briefly outline the derivation of (12). This derivation has been developed in Shape Optimization Theory. We therefore refer to the supplementary document and Chapter 10.6 of the book [Sokolowski and Zolésio 1992] for a detailed exposition. First, the Fréchet derivative of C with a general integrand function y can be expressed as where y is the so-called shape derivative of y under a given velocity field v, defined as The concept of shape derivative is very much analogous to those in continuum mechanics [Bonet and Wood 1997], which has been widely adopted for creating computer animations. In particular, y is equivalent to the material derivative, measuring the change of y in the undeformed (material) space, while y indicates the derivative value in the deformed space, that is, the change rate of y due to the boundary changes only. See §1.2 of the supplementary document for a rigorous mathematical definition ofẏ.
We notice that because of the Dirichlet boundary condition (5), the approximated color u agrees with the input I for all x on the boundary (i.e. x ∈ Γ 0 ). Therefore, when y(x; 0 ) = (u(x) − I(x)) 2 , the second integral term in (15) vanishes, leaving only the first domain integral term, Here the shape derivative u follows the same definition as in (16); we use the fact that the shape derivatives, just like the conventional ones, satisfy the chain rule.
One can prove that the shape derivative u satisfies another Laplace equation (see §1.5 of the supplementary document), where the Dirichlet boundary condition is determined by the normal derivative of both the provided (i.e., I) and approximated (i.e., u) color fields on the same boundaries. Since the Laplacian operator is self-adjoint and u is used in an integral (17), instead of solving the Laplace equation (18), we solve its adjoint problem (11), whose solution enables us to transform the domain integral (17) into a desired boundary integral (9), because ∂ Ω required by Algorithm 2. Although color discontinuities are not well defined for standard pixel images, our method in practice is robust on the choice of boundary curves. Figure 13 shows three examples with boundary curves detected using Canny detector with three thresholds. Notice how missing boundaries (when increasing the threshold) are handled by additional curves generated by Algorithm 2. All our results for standard pixel input used thresholds between 0.1 and 0.2. Figure 14 contains diffusion curve images reconstructed from pixel input (Figures 3 and 12) using Algorithm 2 as well as previous edge detection based methods Xie et al. 2014]. 2 The parameters for each method are selected such that the resulting curves have approximately identical complexities. Our method outperforms previous ones when handling smoothly varying color fields. Notice that, in the bottom two rows (i.e., Flower and Jade), [Xie et al. 2014] has slightly higher approximation errors (measured in RMSE), because it requires higher curve complexities to work properly in these cases. Please see the supplemental material for more comparisons.
3D renderings. Another kind of color field common to computer graphics applications is renderings of 3D scenes. Our approach can be used to approximate these color fields with diffusion curve images. In this case, we represent I as high-resolution pixel images and obtain the boundary curves ∂ Ω directly using object contours. 3 Figure 15 demonstrates diffusion curve images generated using our method from 3D renderings, where slightly higher curve complexities (compared to Figure 14) are permitted to ensure low reconstruction errors. Our method successfully captures detailed appearances: from glossy surfaces to smooth 2 We thank the authors of [Xie et al. 2014] for confirming the correctness of results in Figure 14 generated with their approach. 3 In this example, we assume all objects to be homogeneous. To handle heterogeneity, ∂ Ω needs to include color jumps across object surfaces.

Input Format Scene Time Scene Time
Pixel images (Figure 14) shadow boundaries.
Gradient meshes. Our approach can also generate diffusion curve images directly from input color fields I represented using other vector formats. We demonstrate this using input color fields represented as gradient meshes in SVG format [Bah 2011] where each mesh grid is a Coons Patch [Coons 1967]. In this case, the color I(x) and gradient ∇I(x) can be evaluated analytically for any x ∈ Ω, and the boundary curves ∂ Ω are simply the mesh boundaries. As shown in Figure 16, our method directly creates diffusion curve images closely approximating the gradient meshes, without having to rasterize the input into pixel formats.

Additional Results
Coloring optimization. As discussed in §5.5, our technique is orthogonal and complementary to coloring optimization techniques which take diffusion curve geometry as input and optimizes colors (and color gradients) at both sides of the curves. These techniques can be used to replace our simple coloring scheme which samples color values directly ( §5.1) and performs per-pixel blurring ( §4.4). Figure 18 demonstrates that our optimized curve geometry can be coupled with the coloring optimization scheme introduced by Xie et al.
[2014] to further reduce reconstruction errors. Notice that, since this technique explicitly optimizes color gradients across each curve, we did not perform per-pixel blurring (as stated in §4.4).
Higher-order domain. As discussed in §5.5, our approach can be applied to higher-order domains for generating curves with higher-order smoothness. For instance, as shown in Figure 19, our pipeline can be applied to color gradients rather than original color values. In other words, given a RGB image I, we can use ∇I, a six-channel image, as the input to Algorithm 2. Given our reconstructed gradient image, we solve an additional least square problem to recover the final image.
However, as observed by Xie et al. [2014], we found that for natural images, solving the optimization at higher-order domains normally does not lead to better approximation accuracy under similar curve complexities. This is because higher-order domains are generally filled with significantly more high-frequency contents that require complex (almost space-filling) curve geometry to accurately reconstruct.
Animated result. Lastly, we show preliminary results to motivate future applications of our approach. Since our method optimizes the shape of diffusion curves iteratively, it is suitable for generating animated results from a sequence of gradually changing input color fields. The basic idea is curve reusing: taking optimized curve geometry from one frame as the initial configuration to "warm start" the next one. Figure 20 and the accompanying video show a proof-of-concept example. The input is the relighting (i.e., the object stays static while the light source moves) of a shiny torus knot. In this case, the boundary curves keep unchanged throughout all frames, and optimized curve geometry from one frame remains valid for all other frames. Previous methods Xie et al. 2014] cannot easily enforce curve coherence across different frames, leading to temporally noisy animations. By modifying the curve initialization step in Algorithm 2 to reuse optimized curve geometry, we are able to accelerate the optimization process by 2.3×, and the resulting animation has lower approximation error and little noise. Please see the supplementary video for full animations.

Limitation and Conclusion
Limitation. Our approach has a few limitations that can inspire future work. First, it requires the color field to be C 0 continuous everywhere except at the given boundaries. Robustly finding clean boundary curves, however, can be challenging. Second, if the color field contains spatially high-frequency features, very fine triangulation may be needed to fully resolve them, slowing down our optimization process.
Conclusion. This paper introduces a novel solution to the inverse diffusion curve problem. The key component of our approach is a curve optimization algorithm that iteratively deforms a set of diffusion curves in a way that the reduction of approximation error is guaranteed. Based upon the core algorithm, we develop a full pipeline that takes an input color field plus a set of boundary curves and produces an image with well-shaped and clean curves that closely matches the input. Our approach offers the generality to take input presented in different formats, which we demonstrate using three types: pixel images, 3D renderings, and gradient meshes.

A Brief Derivation of Equation 12
We now briefly outline the derivation of (12 where y is the so-called shape derivative of y under a given velocity field v, defined as The concept of shape derivative is very much analogous to those in continuum mechanics [Bonet and Wood 1997], which has been widely adopted for creating computer animations. In particular,ẏ is equivalent to the material derivative, measuring the change of y in the undeformed (material) space, while y indicates the derivative value in the deformed space, that is, the change rate of y due to the boundary changes only. See §1.2 of the supplementary document for a rigorous mathematical definition ofẏ.
We notice that because of the Dirichlet boundary condition (5), the approximated color u agrees with the input I for all x on the boundary (i.e. x ∈ Γ 0 ). Therefore, when y(x; 0 ) = (u(x) − I(x)) 2 , the second integral term in (15) Here the shape derivative u follows the same definition as in (16); we use the fact that the shape derivatives, just like the conventional ones, satisfy the chain rule.
One can prove that the shape derivative u satisfies another Laplace equation (see §1.5 of the supplementary document), where the Dirichlet boundary condition is determined by the normal derivative of both the provided (i.e., I) and approximated (i.e., u) color fields on the same boundaries. Since the Laplacian operator is self-adjoint and u is used in an integral (17), instead of solving the Laplace equation (18), we solve its adjoint problem (11), whose solution enables us to transform the domain integral (17) into a desired boundary integral (9), because where the last equality follows the integration by parts and Green's formula. Further, the last two integral terms vanish, following the fact that both p (according to the boundary condition of (11)) and ∇u (according to (18)) are zero on the boundary Γ 0 . Eventually, only the first integral term in the last expression remains, yielding the Fréchet derivative of R as a boundary linear form of v n , that is, dR(Ω; 0 ) = Γ 0 B(x) v n (x) dΓ with B(x) expressed as in (12).

Figure 14:
Comparisons between diffusion curve images generated by our approach and previous methods using pixel images (top row) as input. The parameters are adjusted so that the resulting curves generated by each method have roughly identical complexities. Our approach not only yields lower approximation error (measured in RMSE), but also generates better shaped and relatively simple curves. Figure 14: Comparisons between diffusion curve images generated by our approach and previous methods using pixel images (top row) as input. The parameters are adjusted so that the resulting curves generated by each method have roughly identical complexities. Our approach not only yields lower approximation error (measured in RMSE), but also generates better shaped and relatively simple curves.

Reference
Ours Curves Reference Ours Curves

Twill Knots
Reference Ours Curves Reference Ours Curves

Wobble Chess
Reference Ours Curves Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Apple Tomato
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves

Mango Candle
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries. Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Cornell Box Carved
Reference Ours Curves Reference Ours Curves

Twill Knots
Reference Ours Curves Reference Ours Curves

Wobble Chess
Reference Ours Curves Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Apple Tomato
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves

Mango Candle
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries. Figure 15: Diffusion curve images generated from renderings of 3D scenes using our approach. The boundary curves are obtained using mesh contours extracted from the scene geometries.

Apple Tomato
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves

Mango Candle
Reference Gradient Mesh Ours Curves Reference Gradient Mesh Ours Curves Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries. Figure 16: Diffusion curve images generated from gradient meshes directly (i.e., without rasterizing into pixel images) using our method. The boundary curves are given by the mesh boundaries.