A Survey of Smooth Vector Graphics: Recent Advances in Repr esentation, Creation, Rasterization, and Image Vectorization

The field of smooth vector graphics explores the representation, creation, rasterization, and automatic generation of light-weight image representations, frequently used for scalable image content. Over the past decades, several conceptual approaches on the representation of images with smooth gradients have emerged that each led to separate research threads, including the popular gradient meshes and diffusion curves. As the computational models matured, the mathematical descriptions diverged and article started to focus more narrowly on subproblems, such as on the representation and creation of vector graphics, or the automatic vectorization from raster images. Most of the work concentrated on a specific mathematical model only. With this survey, we describe the established computational models in a consistent notation to spur further knowledge transfer, leveraging the recent advances in each field. We therefore categorize vector graphics article from the last decades based on their underlying mathematical representations as well as on their contribution to the vector graphics content creation pipeline, comprising representation, creation, rasterization, and automatic image vectorization. This survey is meant as an entry point for both artists and researchers. We conclude this survey with an outlook on promising research directions and challenges to overcome in the future.


INTRODUCTION
V ECTOR graphics is a branch of computer graphics describ- ing the representation and synthesis of images based on geometric primitives.Compared to raster graphics, which require to keep per-pixel information, vector graphics provide a more compact representation.In particular, the specification of the image content in vector-based images is resolutionindependent, meaning that they can be scaled and sampled to any desired output resolution.The simplicity of vector graphic representations can be a double-edged sword.While it is easy to create simple color gradients, complex color variations like smooth shadows, caustic effects, and defocus blur are hard to be represented directly through basic gradient functions [101].When the amount of required detail increases, the editing process of vector graphics demands both solid artistic skills and a large amount of working time to obtain complex results.Therefore, vector graphics are more often used to represent stylized images, such as line drawings [33], [37], [100], fonts [92], [107], clipart [30], [38], [59], [114], [139], image triangulations [76], and sketches [12], [117].With the goal of increasing the expressiveness of smooth vector graphics, current research on gradient-based vector graphics can be divided into two directions.One direction is to enhance the artistic control of vector graphics, such that more complicated effects can be constructed and manipulated with less effort.Thereby, the rendering process that converts vector representations into raster images is a crucial component, as high-quality results are expected to be generated at interactive or even real-time speed during editing.The other direction seeks to find the most suitable vector representation for a given image in a (semi-) automatic way with as few primitives as possible, which is called vectorization.The vector representations generated from the two directions are often on opposite ends of the expressiveness spectrum.Usually, complex color variations are hard to represent directly by a single primitive or through simple artistic interaction, hence, the vectorization approach often leads to a richer set of primitives, and the results are therefore more difficult to manipulate directly.The two directions are not completely orthogonal.The advances in the formulation of artist-friendly primitives may inspire vectorization methods, while the vectorized results can act as a basis for further editing and manipulation.Until now, both directions have been studied intensively within their respective vector representations.
Currently, there are two major extensions to simple shapes and gradients: mesh-based and curve-based representations, as explained by Barla and Bousseau[6] in their artist-oriented introduction to gradient art.Mesh-based representations divide the image domain into 2D patches, where color and other attributes are stored at the vertices or in the patch interiors.The research on mesh-based methods can be divided roughly into two directions, aligning with the two aforementioned directions in vector graphics.One path focuses on finding better topological representations of the color variations, such as grouping smaller patches [7], [8], [9], using subdivision-based techniques [84], [144], or adopting irregular shapes like B ezigons [139].These methods mainly concentrate on automatic extraction of a vector-based description from a given image.The other path adds expressiveness by defining additional attributes on the patches, and hence provides more artistic control.For example, the gradient mesh formulation [74], [118] allows users to specify tangents at the patch vertices.In this way, color gradients can be controlled directly.Despite this, mesh-based representations can still be less intuitive as users may need to decompose the domain by themselves.In contrast, the creation process of curve-based methods is closer to free-hand drawing.Curve-based methods use various curves with different attributes like color and normal derivatives as the fundamental primitives.Users can draw the outline of their design and add colors afterwards.The rasterization process of curve-based representations like diffusion curves [101] often requires solving a partial differential equation (PDE) on the entire domain, which inhibits local control.In addition, the solving step may suffer from discretization problems, and an additional blurring step may be required to generate smoother results.Hence, most research in curve-based methods either focused on establishing sound formulations to enhance local control [41], [61], [64], [68], or put effort in developing better rendering techniques to avoid the intrinsic discretization issues [69], [120], [121].As curvebased representations often describe color discontinuities in an image, it can still be difficult to represent either extremely smooth or very detailed regions.
With ongoing research, advances have been made on both fronts, adding new representations, creation tools, rasterization methods, and vectorization algorithms to both meshbased and curve-based methods.Synergies are possible, since the edges of patches are often aligned with the curvilinear features of an image.However, in general, the research threads within both approaches evolved independently with their own mathematical formulations, making it challenging to transfer ideas from one branch to the other.Our goal is not only to present a comprehensive overview of the advances in smooth vector graphics, but also to convey a solid understanding of the underlying mathematical formulations by phrasing them consistently in the same notation.Further, we provide an in-depth discussion of existing algorithms, such that a generalized and unified representation may be developed in the future to unite the advantages of previous methods.

Methodology
Research over the past decades led to a wide range of vector graphics representations, which we categorize into two types: 1) Mesh-based representations decompose an image into 2D patches, where color is interpolated inside each patch.In this survey, we consider this formulation as generalization of simple shapes and gradients.2) Curve-based representations mimic the traditional artistic habits of drawing on paper.The curves usually represent extrema in the color gradient, and the resulting image is often computed by solving a partial differential equation (PDE).For each type of primitive, we classify the research efforts into the following parts of the content creation pipeline: Representation.Research on this topic introduces novel vector graphics representations, which consist of two ingredients: (1) the definition what kind of information is stored on the geometric primitives, such as colors or other features, and ( 2) the underlying mathematical model to fill in void areas between primitives for which appropriate solvers can be chosen or developed.
Creation.Tools for the creation of vector art are often tailored to a particular vector graphics representation.Methods that ease the creation and editing process or that provide more flexibility and richness in the final results are classified into this category.
Rasterization.Given a vector graphics representation, the rasterization computes a rasterized image.Methods in this category seek to improve either the rendering quality or the rendering speed on a desired output device.Currently, solvers are tightly linked to a particular representation, i.e., there is no unified solver for all existing formulations.
Vectorization.For any input image or video, vectorization methods output a vector graphics representation, which resembles the input image or video as closely as possible.When additional information is provided, e.g., when the curve geometry is known [68], [71], the vectorization process may only solve for a part of the model.
As illustrated in Fig. 1, the research topics are internally connected to each other.The creation and editing outputs a vector graphics representation, which is then rendered to a raster image or a video (i.e., an image sequence).On the other hand, vectorization acts as an inverse process that converts rasterized image(s) into a vector graphics representation.

Structure of the Survey
The survey is structured with respect to the categorization proposed in Section 1.1.For each vector graphics primitive, we describe state-of-art methods in representation, creation, rasterization, and vectorization.In Section 2, we introduce the mathematical notation used in this paper, and we discuss related mathematical foundations.We then concentrate on Fig. 1.The connections between the main research topics in vector graphics: a vector graphics representation is defined by its underlying mathematical model, which can be created or edited accordingly.The rasterization process generates raster image(s) from a vector graphics representation, while vectorization describes the inverse problem.
mesh-based representations in Section 3, including simple shapes and gradients.Subsequently, we thoroughly analyze curve-based vector graphics representations in Section 4. Lastly, we conclude the survey in Section 5, where we discuss the connections between the formulations and cover possible avenues for future work.

Notation
In order to unify the mathematical formulations in this paper, we adopt the same notations for all equations.Scalars and scalar-valued functions are expressed with italic letters, such as s or sðtÞ.Scalar-valued basis functions are denoted with capital italic letters, for example BðtÞ.Vectors and vector-valued functions are expressed with bold lowercase letters, like v or vðtÞ.All vectors are laid out as columnvectors.Matrices are denoted by capitalized bold letters, such as M. To make it clear whether a symbol denotes a constant or a function, we always list functions with their arguments, for example xðtÞ.When needed, discrete values are indexed with lower scripts, such as x i .Note that user parameters are treated as constants in the expressions and are therefore not explicitly listed in the function arguments.If a derivative is a user-defined constant, we use @ u f i;j to abbreviate the derivative at a discrete location @fðu;vÞ @u ju¼u i ;v¼v j .

B ezier Curves
B ezier curves are widely used as the fundamental representation of curve-based methods, as well as patch boundaries in mesh-based methods.A B ezier curve xðtÞ : ½0; 1 !R 2 is: where B n i ðtÞ are Bernstein basis functions and b i are B ezier control points.The first and last control point are always interpolated, i.e., xð0Þ ¼ b 0 and xð1Þ ¼ b n .A cubic B ezier curve has a polynomial degree of n ¼ 3 and is specified by 4 control points b 0 ; . ..; b 3 2 R 2 .From the tangent _ xðtÞ ¼ dxðtÞ dt , the curve normal nðtÞ : ½0; 1 !R 2 is defined as: When interpolating multiple points, it is numerically more suitable to join multiple curves, as described next.

Spline Curves
Two (or more) curves x 1 ðtÞ : ½t 0 ; t 1 !R 2 and x 2 ðtÞ : ½t 1 ; t 2 !R 2 can be joined at the end points with suitable continuity constraints to form a B ezier spline curve xðtÞ : ½t 0 ; t 2 !R 2 .For example, C 0 -continuity requires matching end points For a cubic curve, the hardest constraint is C 2 -continuity, i.e., C 1 and matching acceleration vectors given parameterization (t 0 , t 1 , t 2 ) requires the solution of a linear problem with boundary conditions to find a unique solution.When selecting natural end conditions € x 1 ðt 0 Þ ¼ € x 2 ðt 2 Þ ¼ 0, the resulting curve minimizes the approximate bending energy R t 2 t 0 k€ xðtÞk 2 dt.In analogy to the extension of this minimizer to surfaces, this curve is sometimes also referred to as thin plate spline.The solution of a global linear problem can be avoided if a continuity constraint is given up.For example, Catmull-Rom and Hermite splines are only C 1 continuous cubic curves that can be computed locally.Switching from Bernstein basis functions B n i ðtÞ to B-spline basis functions leads to B-splines, which are equivalent to B ezier splines in terms of the function space.The control points of B-splines are called de Boor points.Rational B-splines are an extension that includes a weighting term for each de Boor point, making it possible to also represent conic sections.We refer to Farin [36] for a comprehensive introduction to curve design.

MESH-BASED METHODS
We begin our overview of vector graphics representations with mesh-based methods.An overview of existing methods is presented in Table 1, which follows the categorization from Section 1.1.We start with the different mesh representations.

Representation
Mesh-based methods divide the image domain into nonoverlapping 2D patches across which colors are interpolated.The representation is mainly concerned with the placement and connectivity of patches and determines how color is interpolated.The patch shape can be triangular, rectangular, or even irregular, while color and other attributes are stored at vertices or in the patch interior.Fig. 2 provides an overview.Next, we categorize the methods based on the domain decomposition.

Triangular Meshes
Triangle-based methods interpolate a function fðu; v; wÞ across each triangle with u; v; w 2 R being barycentric coordinates with u þ v þ w ¼ 1.Given discrete control values f i;j;k and associated basis functions F n i;j;k ðu; v; wÞ of degree n with i; j; k 2 f0; . ..; ng and i þ j þ k ¼ n, the continuous function fðu; v; wÞ is: Using Bernstein polynomials B n i;j;k ðu; v; wÞ ¼ n! i! j! k! u i v j w k as basis functions F n i;j;k ðu; v; wÞ gives rise to B ezier triangles and the special case of n ¼ 1 corresponds to linear barycentric interpolation.
The usage of triangulations or triangular patches to represent images has a long history [27], [29], [116], especially in image vectorization.An overview of basic triangle primitives is shown in Fig. 3. Basic triangulations with piecewise constant colors as in Fig. 3a require too many triangles to depict complex color variations [132].Further, C 0 continuity across patch boundaries may become noticeable at higher resolution for smooth regions [84].One way to better align the features is to use curved edges around the triangular patches.Xia et al. [132] proposed a representation using non-overlapping triangular patches with curved edges, see Fig. 3b.Each patch is modeled as a 2D B ezier triangle as specified in Eq. (3), such that every boundary is a 2D cubic B ezier curve as in Eq. ( 1).The color of interior B ezier control points is fitted by thin plate splines.The boundary B ezier curves are aligned with color discontinuities for a given image, which enables larger color variation at the boundaries and smooth colors inside the patches.More recently, Xiao et al. [135] used quadratic B ezier curves as the patch edges.The control points, which are placed on the edges, are computed through an optimization process.However, neither thin plate spline fitting nor energy minimization guarantees C 1 continuity at non-edge patch boundaries [84], [135].Alternatively, subdivision-based methods can be used to tackle the continuity problem, as illustrated in Fig. 3c.Liao et al. [84] improved the continuity property within and across the triangular patches, such that the boundary curves are C 2 continuous.Other than regions with discontinuity features (for example across edges that are intended to be sharp), the final color is no less than C 1 continuous [84].The patches are constructed using a piecewise smooth subdivision scheme.The color attributes are stored at the vertices and are optimized such that the color function over the entire image is C 1 continuous.Zhou et al. [144] also used subdivision-based triangular meshes as their fundamental primitives.They used cubic B-splines to describe the curvilinear features in an image.The triangular meshes are constructed and refined similarly as in Liao et al. [84].The color attributes are stored at the vertices.To model sharp edges, they used two nearby curves with separate colors.In addition, a sharpness factor is stored at the edges to control edge smoothness.Generally, the subdivision-based approaches may generate an excess amount of primitives, and can be  strongly dependent on the feature extraction output [146].Recently, Zhu et al. [146] proposed to use quadratic triangle configuration B-splines (TCBs) [19], [112] to faithfully depict complex color variations while remaining C 1 continuous.An image can be treated as a surface in 5D (geometry and color) and can therefore be approximated and refined as a TCB-spline surface.

Rectangular Meshes
Rectangular meshes interpolate functions fðu; vÞ over a u; v-parameter domain from given control points f i;j with i 2 f0; . ..; ng and j 2 f0; . ..; mg by computing tensor products with the basis functions F n i ðtÞ and F m j ðtÞ: Price and Barret [105] used bi-cubic tensor product surfaces, also known as B ezier patches, where the basis functions F n i ðtÞ are the Bernstein basis functions B n i ðtÞ from Eq. ( 1) and n ¼ m ¼ 3. Thus, every patch has four rows and columns of B ezier control points x i;j , each representing a cubic B ezier curve.Colors c i;j are stored at the control points and are likewise interpolated using Eq. ( 4).
Although the bi-cubic B ezier patch formulation is simple and straightforward, it still lacks control over the color variations from the editing perspective.One alternative is to use gradient meshes.This representation was first introduced by Adobe Illustrator for manually creating complex color variations, although the actual implementations were proprietary.In 2007, Sun et al. [118] published a formal definition of gradient meshes that are composed of Ferguson patches [40].A Ferguson patch is a bi-cubic tensor-product surface following Eq.( 4), which is specified through four corners f i 0 j 0 with i 0 ; j 0 2 f0; 3g and u; v-tangent handles at each of the four corners @ u f i 0 j 0 and @ v f i 0 j 0 see Fig. 4a.By definition, the second-order mixed partials @ u @ v f i 0 ;j 0 are set to zero at the corners, which is referred to as zero twist, cf.later Section 3.3.2.A gradient mesh is composed by aligning Ferguson patches along mesh-lines, see Fig. 4b.Each vertex q of the gradient mesh's mesh-lines stores: Position.The 2D position x q ¼ ðx q ; y q Þ T of the vertex.
Tangents.Every corner has four u; v tangents, shared by the four neighboring patches, see Fig. 4b.The four tangents are constructed by the two directions @ u x q and @ v x q while the opposite directions are scaled by factors a u and a v .
Color.Each vertex has an RGB value c q ¼ ðr q ; g q ; b q Þ T .Color Derivatives.For every vertex q, its color derivatives @ u c q , @ v c q along the u; v tangents are estimated per color channel with monotonic cubic spline interpolation [131].For example, in Fig. 4b, the color derivative in the u direction can be analytically calculated from the cubic splines between vertices q À 1; q and q; q þ 1.The color derivative for the v direction can be computed analogously.
One limitation in Sun et al. [118] is that regions with holes are hard to represent.Thus, Lai et al. [74] proposed a topologypreserving gradient mesh representation, which models a cut in a surface xðu; vÞ by duplicating consecutive horizontal grid lines, see Fig. 4c.Let xðuÞ :¼ xðu; v c Þ be the horizontal curve with u 2 ½u 0 ; u 1 and v c ¼ const, along which the surface is cut.The two formerly adjacent patches receive distinct boundary curves xðuÞ and e xðuÞ, requiring ) at the end points to form a closed loop.If there is more than one hole in the mesh, several cuts can be performed accordingly.By splitting a topology-preserving mesh [74] with n holes along the horizontal grid lines, this representation can be transformed into n þ 1 gradient meshes [118].The two representations were adopted by many following works [127], [130], [136], [137].Later, Barendrecht et al. [5] used cubic rectangular patches and extended the primitives to support local refinement.

Irregular Meshes
Swaminarayan and Prasa [123] used polygonal patches with constant color, see Fig. 5a, where the edges of the polygons align with curvilinear features of an input image.Apart from polygonal patches, constant colors and linear gradients have been used in other shapes, such as areas of homogeneous features [43], see Fig. 5b.More recently, Yang et al. [139] used B ezigons (closed regions bounded by B ezier curves) to depict clipart images.The color in each region is represented by a spatially-varying color function.
Later, Hettinga et al. [53] investigated several schemes to interpolate the values and gradients of color from boundary conditions, including the approach for topologically unrestricted gradient meshes by Lieng et al. [86], the cubic mean value coordinates by Li et al. [81], and two newly-proposed interpolation methods based on the Gregory generalized B ezier patch and the Charrot-Gregory corner interpolator.The extension to gradient meshes of arbitrary topology by Lieng et al. [86] is illustrated in Fig. 6, where each patch can have an arbitrary number of vertices and edges.Each control point has an assigned color value and color gradients, and the edges can be set to different discontinuity levels (C 1 smooth, C 0 crease, or C À1 discontinuous).

Creation
While simple shapes and gradients are easy to control, vector-based representations are more difficult to create and edit when the amount of required detail increases.In the following, we cover manual and data-driven approaches, as well as alternative recoloring and texturing methods, leading to representations as illustrated in Fig. 7.

Manual Methods
Simple shapes are directly specified by shape parameters, e.g., position and radius of a circle.In Adobe Illustrator, gradient meshes are created manually by selecting an object and specifying the number of vertices per row and column.Users then need to choose a color for each grid point.Sun et al. [118] enabled users to create meshes in a similar way.After specifying the objects and dividing each object into 4 components, each sub-object is uniformly split parametrically.Users can click to add or remove mesh lines.The position of a vertex can be dragged or typed in, the four derivatives can be dragged through tangent handles, and the color can be picked from a color palette.Xia et al. [132] allowed users to select patches of interest and to edit colors using thin plate splines.Similarly, Liao et al. [84] let users select and edit either an individual vertex or a feature, or even a local region by defining a radius around a selected vertex or feature.After selecting the center of deformation, the geometry can be refined through distortion minimization [63].Users can change the color of the chosen vertex or feature by assigning a new color or by blending the new color with the old color.Sverja et al. [122] edited irregular meshes by adding faces to the mesh interactively, and by splitting or collapsing edges of the mesh.Verstraaten and Kosinka [126] further let users add sharp color transitions in this formulation.Barendrecht et al. [5] added branching such that the mesh is not restricted to a rectangular grid.They also extended sharp color transitions such that a control point can be assigned with up to 4 different colors.Subdivision-based methods provide users with multi-resolution editing [147].Liao et al. [84] and Zhou et al. [144] let users deform the mesh at a desired resolution level.
Lieng et al. [86] enabled users to change colors at a higher resolution, while the geometry is manipulated at a coarser level.

Color Transfer, Recoloring, and Textures
Rather than relying on manual color editing, color transfer methods have been proposed to ease color assignment.Xiao et al. [136] proposed an example-based color transfer method for gradient meshes.Given a reference image, the algorithm transforms the colors of a gradient mesh using a linear operator constructed by PCA-based color transfer algorithms [1], [134].The method assumes Gaussian color distributions in both the reference image and the gradient mesh, but may fail to preserve image details like textures after color transfer [136], [137].Later, they proposed a refined method using an optimization-based linear operator [137].The color transfer process is modeled by minimizing the difference in color distribution, which also enables adding image gradients as additional constraints such that fine details can be preserved after color transfer.Instead of using a reference image, Wan et al. [127] developed a scribble-based method to recolor a gradient mesh.Given an input mesh and scribbles by the user, an auxiliary mesh called control net is constructed.The control net structure can be computed for both classic gradient meshes [118] or topology-preserving gradient meshes [74], which are then recolored with extended chrominance blending [140] and are in the end converted back to a gradient mesh representation.
High-frequency image details often lead to a large number of patches when converted to a mesh-based representation, and can therefore be a burden for future editing.Methods have been presented to add texture details in different ways.Hettinga et al. [55] added Perlin, Worley and Gabor noise.Baksteen et al. [3] used mesh colors to add textures to gradient meshes, which is another option for 3D texture mapping in 2D [142].In this representation, colors are defined at both the vertices, and at uniformly sampled points along the edges and inside the patches to represent complex color variations with a limited number of patches.

Rasterization
The rasterization process determines how colors are computed within each patch, which may be constant, a linear gradient, a combination of simple gradients, or the result of a fitting or interpolation from boundary colors.The rasterization process of mesh-based representations can be roughly divided into subdivision-based methods [84], [86], [144] that may use adaptive GPU-accelerated tessellation [52], and interpolation-based methods [81].In this section, we discuss approaches for the various mesh formulations.

Triangular Meshes
Xia et al. [132] rendered triangular B ezier patches in three steps.First, they recursively subdivided each patch into 4 smaller B ezier patches until triangles were smaller than a pixel.Second, they determined the associated patch and its barycentric coordinates for each pixel.Third, the color of the pixel was computed by thin plate spline fitting using the barycentric coordinates.Liao et al. [84] accelerated the rendering by a GPU-based patch subdivision.Zhu et al. [146] used TCB-spline interpolation for image reconstruction.

Rectangular Meshes
Across rectangular patches, colors are interpolated using tensor products.Using Bernstein basis functions B n i ðtÞ from Eq. ( 1) with n ¼ m ¼ 3 in Eq. ( 4) gives rise to bi-cubic B ezier patches.Instead of specifying 16 control points f 0;0 , f 1;0 , . .., f 3;3 for the color function, the same tensor product surface fðu; vÞ can be expressed equivalently using only the corner control points and their first-order and second-order mixed partial derivatives, cf.Sun et al. [118]: with the monomial basis functions tðtÞ ¼ ð1; t; t 2 ; t 3 Þ T and Since the Ferguson patch formulation specifies only the corner control points and its first-order partial derivatives, the mixed second-order partial derivatives in Eq. ( 6) are usually set to 0, such that @ u @ v f 0;0 ¼ @ u @ v f 0;3 ¼ @ u @ v f 3;0 ¼ @ u @ v f 3;3 ¼ 0. By substituting positions x q and tangents @ u x q , @ v x q or colors c q and tangents @ u c q , @ v c q , respectively, into f q and its tangents @ u f q , @ v f q in Eq. ( 6) allows for the interpolation of positions and colors.Equivalently, Barendrecht et al. [5] used cubic Hermite boundary curves.

Irregular Meshes
SVGenie [9] and SVGWave [7] both use polygonal mesh representations, for which a constant color or a simple color gradient can be assigned directly.Lecot and L evy [77] approximate the color in each polygonal region with a linear combination of constant, linear and quadratic basis functions for each color channel.Later, more advanced techniques have been developed to interpolate polygons.Li et al. [81] used cubic mean values to interpolate the colors and gradients at the vertices, which was developed from the mean value property for biharmonic functions.Beatson et al. [10] then extended the domain to piecewise quadratic boundaries, where they proposed the Hermite mean value interpolation technique for polygons.Hettinga et al. [53] further adjusted the interpolation schemes from generalized B ezier patches [125] and Charrot-Gregory corner interpolation patches [22] (both are multi-sided parametric patches) such that the interpolation scheme can be used for 2D gradient meshes.By using the algorithm from Chiyokura and Kimura [25], they achieved at least G 1 continuity for the color.Lieng et al. [86] proposed an interpolation method based on Catmull-Clark subdivisions for rendering gradient meshes with arbitrary topology [21], [98].

Vectorization
In contrast to the manual creation of meshes, a vectorization process decomposes a given image into 2D patches.Common starting points are image segmentations and triangulations based on edge detections.When complex color models are used, such as gradient meshes with color tangents [118] or B ezigons with parametric edges [139], optimization techniques that minimize a reconstruction error are commonly used.Apart from the generic mesh types discussed in the following, several learning-based techniques have been proposed to vectorize line drawings [33], [49], sketches or manga [35], [117], and floorplans [91].

Triangular Meshes
Triangular patches are often extracted in two steps.The first step determines the mesh geometry, then the second step extracts or optimizes the color attributes at the vertices or within the patches.Xia et al. [132] performed a triangulation at pixel resolution as their first step.In order to keep image features, subpixels are inserted at edges.By treating each color channel as a height field, every pixel is lifted to a vertex in a 3D domain, and every subpixel is transformed into two vertices with the same 2D projection but with distinct color values to avoid blurry edges.Afterwards, the 3D meshes are simplified collectively regarding the quadric error metric [47] and are projected to the 2D image space.Non-overlapping triangular B ezier patches are fit to the triangulated domain through an optimization process.Colors are fit to the B ezier patches using thin plate spline fitting [103].Liao et al. [84] performed the vectorization in a hierarchical manner.First, a Canny edge detector [18] is used to extract image features.In order to obtain region boundaries, GrabCut [109] is applied for image segmentation.Second, an initial mesh is constructed, where each pixel is initialized with a vertex, and additional vertices are inserted at the extracted features if needed.A retriangulation step is performed locally.Third, the dense mesh is simplified using quadric error metrics [47].The color at the vertices is globally optimized.Zhou et al. [144] developed the vectorization step similarly for their subdivision-based patches.First, curvilinear features are extracted using contour detection [2] and corner detection [96].Cubic B-splines are then fit to the detected edges.Second, a constrained Delaunay triangulation is carried out to obtain the initial mesh.Colors are assigned to the subdivision patches through an optimization process similar to Hoppe et al. [58].The color at each vertex is then optimized through a least squares minimization.For places where the color fitting error exceeds a threshold, new vertices are added, and a constrained Delaunay triangulation is further applied for refinement.More recently, Hettinga et al. [51] fit cubic B ezier triangles to an input image, and mapped textures as mesh colors through an optimization process.Later, they extended the vectorization framework to enhance user control and feature extraction [54].

Rectangular Meshes
Early methods like Price and Barret [105] vectorize an image in two steps.First, the mesh geometry is constructed.For this, a graph cut [16] is used to extract the object contour, which is guided by the user.Once an object is detected, curvature is used to locate four corner points on the object.Then, axes can be added to the object by conducting a least cost search to find a minimum cost path.Mesh lines are added similarly by subdividing the object and recursively performing the search.Each mesh cell is then fit by a B ezier patch.Second, color attributes are assigned to the B ezier control points by sampling colors from the input image.As gradient meshes have a richer set of attributes, the geometry and the attributes are often optimized together in a single minimization step.Sun et al. [118] required a manual initialization in which users first decompose the input image into multiple sub-objects, and then divided the boundary of each sub-object into 4 segments to which cubic B ezier splines were fitted.Users can interactively subdivide the patches further.The manual result is then used as initialization to a non-linear least-squares minimization to fit Ferguson patches to the input image, which includes the color values and color derivatives.The non-linear optimization is solved with the Levenberg-Marquardt algorithm [79], [99].Users can guide the optimization further by drawing preferred directions of the mesh lines.Lai et al. [74] extended the representation to topology-preserving gradient meshes, and proposed a fully automatic vectorization method for this representation.They first segmented the image into objects, using a graph-based method for coarse segmentation [39], and refined the boundary with Grab-Cut [109].
The objects are converted into triangular meshes by performing a constrained Delaunay triangulation on the samples from a saliency map with error diffusion [67].This is computed by applying a compass filter on the image.The triangular mesh is then mapped into a rectangular domain, and is reparameterized for local refinement.Color are sampled from the underlying image and for each control point, color is computed through interpolation.The derivatives of color are calculated using monotonic cubic interpolation [131].Wei et al. [130] developed a faster method by adopting quadrangulation techniques for 3D meshes that are often used in geometry processing tasks.

Irregular Meshes
The Data Dependent Triangulation method (DDT [32]) in SVGenie [9] and the Wavelet Based Triangulation method (WBT [78]) in SVGWave [7] both decompose the image domain into triangles and group the triangles into polygonal patches based on similarity.SVGenie [9] seeks for a locally optimal triangulation iteratively.Based on Yu et al. [141], each pixel is first divided into two triangles by connecting the lowest cost diagonal.For each convex polygon composed of adjacent triangles, a look-ahead step is performed to swap polygon edges if necessary.The triangulation is further optimized by minimizing a cost function iteratively.SVGWave [7] performs a wavelet transformation hierarchically.In each level, the triangulation can be computed with the wavelet coefficient and is refined iteratively.Swaminarayan and Prasad [123] first used an edge detector to extract edge pixel chains and performed a constrained Delaunay triangulation such that the edges of the triangles align with the extracted edges.Each triangle is assigned a constant color by sampling several points within the triangle, which constitutes a so-called trixel.Neighboring trixels are grouped together to form polygons with respect to their proximity and continuity.Triangulation may also be used as an intermediate step to accelerate vectorization.Lecot and L evy [77] developed a vectorization technique to compute closed regions bounded by cubic splines.The algorithm partitions the domain into a Voronoi diagram.In every region, the color is estimated with a quadratic polynomial.In a basic setup, the algorithm can be performed at pixel-level.Each pixel grows into a region in a greedy way and the polynomial coefficients are updated iteratively.In order to accelerate the vectorization process, instead of applying the algorithm on each pixel, an intermediate triangulated structure of trixels is constructed.Yang et al. [139] proposed a method to vectorize clipart images into B ezigons.A B ezigon is a closed B ezier curve with a color function defined in the interior.The initial B ezigons can be extracted automatically using image segmentation [39] and a curve fitting method [113], or they can be created by hand.In addition to the reconstruction error, the optimization seeks to minimize failure cases, such as self-intersection and angle-variation.He et al. [50] introduced a perception-aware vectorization for quantized raster images.

CURVE-BASED METHODS
In their seminal work, Elder et al. [34] used edges with attributes (intensity, blur scale, and gradient direction) to encode images, which has spurred further research as discussed throughout this section.A comprehensive overview based on the taxonomy in Section 1.1 is presented in Table 2.

Representation
First, we concentrate on the representations that have been developed to describe images by means of curves.An overview of the evolution is provided in Fig. 8.

Basic Formulation
Inspired from Elder et al. [34], Orzan et al. [101] proposed a new vector graphics primitive called diffusion curves, which uses B ezier curves as geometric primitives with colors defined on the left and right sides of the curves.The final image is constructed by diffusing the colors from both sides of the curves, which allows for the modeling of color discontinuities.The diffusion process may be followed by a blurring step to obtain smooth edges.Fig. 9 illustrates the three components of a diffusion curve, as defined below: A Cubic B ezier Spline.The geometric curves are specified as splines, formed from cubic B ezier curves xðtÞ : R !R 2 , see Eq. (1).
Color Control Points.On each side of the curve xðtÞ, a color is specified by the user, formally denoted using functions c l ðtÞ and c r ðtÞ.The colors are piecewise linearly interpolated from a set of color control points, placed along the curve.Each set contains at least two colors to define the color at the start (t ¼ 0) and end (t ¼ 1) of the curve.
Blur Control Points.As an additional degree of freedom, Orzan et al. [101] introduced a smoothness parameter sðtÞ along the curve, which controls the filter size of a Gaussian smoothing carried out in a post-process.As with the colors, the function sðtÞ is piecewise linearly interpolated from a set of blur control points.
The final image uðxÞ is then defined as the smoothest image that meets both the color gradient wðxÞ along the diffusion curves as well as specific colors along the image boundary, formally solved for via: Eq. ( 7) requires the x-and y-partial derivatives of the color field @wðxÞ @x : R 2 !R 3 and @wðxÞ @y : R 2 !R 3 that are zero everywhere, except along the curves xðtÞ: @wðxÞ @x jx¼xðtÞ ¼ ðc l ðtÞ À c r ðtÞÞ Á n x ðtÞ @wðxÞ @y jx¼xðtÞ ¼ ðc l ðtÞ À c r ðtÞÞ Á n y ðtÞ (10) which requires the curve normal components n x ðtÞ and n y ðtÞ, cf.Eq. ( 2).Note that in regions away from the diffusion curves, we have @wðxÞ @x ¼ @wðxÞ @y ¼ 0 to achieve smoothness.Eqs. ( 7)-( 8) were solved on a discrete image.To place the colors on the left and right of the diffusion curves, Orzan et al. [101] moved the color sources with an offset of a few pixels away from the curve geometry.It is worth noting that smooth color transitions across diffusion curves are difficult to represent directly.Hence a blurring post-process is essential for this model.Later, several variations have been proposed to enrich the modeling abilities of diffusion curves.Based on how the diffusion process is formulated, we categorize these models into the following types.

Harmonic Equation
Jeschke et al. [69] generalized the formulation of Orzan et al. [101] by treating diffusion curves as domain boundary curves with colors c l ðtÞ and c r ðtÞ, such that there is no need to define the gradient fields @wðxÞ @x and @wðxÞ @y as in Eq. ( 7) anymore.The diffusion process can then be formulated as a Laplacian equation, which is also called a harmonic equation: In this formulation, the final image is smooth everywhere except at the boundaries, which now includes not only the rectangular image boundaries but also the diffusion curves.As in Eq. ( 8), Dirichlet boundary conditions are used in Eq. ( 12).This formulation has become the most popular one for rendering methods [69], [102], [120], vectorization techniques [94], [143] and editing tools [71], [93].
One problem of the harmonic formulation is the lack of control over the color gradients.The formulation currently seeks smoothness inside the domain but not across diffusion curves, as illustrated in Fig. 10.This may generate tent-like artifacts around the curves.Jeschke [68] addressed the problem by using a linear combination of harmonic functions: where u 1 ðxÞ and u 2 ðxÞ are the solutions to two individual Laplacian equations defined according to Eq. ( 11) on the same B ezier curve boundary with different color functions g 1 and g 2 for the Dirichlet boundary conditions.Thereby, w l ðxÞ is an interpolation function that spatially blends between u 1 ðxÞ and u 2 ðxÞ.The function w l ðxÞ is computed by a Poisson equation and a parameter l controls the blending speed between u 1 ðxÞ and u 2 ðxÞ.Likewise, the subsequent blurring is adapted and optimized.When u 1 ¼ u 2 , this reduces to Eq. ( 7).Both the harmonic model and the composition model require a subsequent blurring step, though in cases when the shape geometry and color control points are optimized, this process can sometimes be omitted [143].

Biharmonic Equation
The basic harmonic formulation lacks control over smoothness and direction of color gradients.In addition, with each added curve the derivative continuity is broken [41].One way to address this problem is to specify constraints on higher-order derivatives.Finch et al. [41] proposed a bi-harmonic formulation based on thin plate splines, which seeks solutions that are harmonic in their Laplacian domain: Here, gðxÞ is a Cauchy boundary condition function, which can be considered as a combination of Dirichlet and Neumann boundary conditions.The boundary curves are constrained in both value and gradient of color.Later, Illbery Fig. 10.The harmonic equation only enforces smoothness away from the curves' constraints, therefore artifacts may occur near the curves.Images adapted from [6].et al. [64] generalized the formulation of biharmonic equations, such that harmonic diffusion curves can be considered as a special case of the biharmonic formulation.Since directional color derivatives can be controlled directly for biharmonic diffusion curves, there is no need for a subsequent blurring process.However, the biharmonic formulation may bring other artifacts.One problem is that extrapolation can create new extrema at unconstrained positions, especially at colinear control points [15], [41].This can be solved by manually adjusting the control point locations [41] or by applying a non-linear optimization [66].Another problem is that when the curve geometry is scaled, the diffused color range will scale with the geometry, which may result in color oversaturation [68].Later methods use other formulations such as a composition of Laplacian equations (Eq.( 13)) [68] or a combination of Laplacian and Poisson equations (Eq.( 17)) [61] to avoid this artifact.

Poisson Equation
Orzan's representation [101] and the harmonic equation in Eq. ( 11) can be further generalized to a Poisson equation: Therefore, when fðxÞ ¼ 0, Eq. ( 17) is reduced to the harmonic model specified in Eq. ( 11).The initial representation in Eq. ( 7) can also be considered a special case of Eq. ( 17) when setting fðxðtÞÞ ¼ ðc l ðtÞ À c r ðtÞÞ Á ðn x ðtÞ þ n y ðtÞÞ on the curves and fðxÞ ¼ 0 elsewhere.By manipulating the expression of fðxÞ and gðxÞ, more control over the diffusion process can be achieved.Bezerra et al. [11] solved Eq. ( 17) linearly, adding additional diffusion constraints to the system, such that the strength and direction of the diffusion process could be controlled.Later, Hou et al. [61] set fðxÞ as a piecewise constant function in order to control color gradients locally and globally.They used harmonic diffusion curves to depict color variants in an image, and added Poisson curves or Poisson regions controlled by fðxÞ for controlling shading details.Same as with the harmonic equations, gðxÞ is used to specify Dirichlet boundary conditions for the curve colors.
If fðxÞ ¼ a @uðx;tÞ @t , Eq. ( 17) is a standard heat equation, with time t and thermal diffusivity a, i.e., the rate of heat transfer.By introducing a time t, the result image is extended to uðx; tÞ : R 2 Â R !R 3 .Lin et al. [88] modified the heat equation to control diffusion strength and direction via a scalarvalued diffusion coefficient cðxÞ : R 2 !½0; 1: Eq. ( 20) is similar to previous methods in that gðxÞ is zero everywhere except along the boundary curves where it denotes the Dirichlet boundary condition.At time t ¼ 0, the initial state of the image uðx; 0Þ contains the undiffused colored curves.We can reorganize Eq. ( 19) to a Poisson equation when cðxÞ 6 ¼ 0: r 2 uðx; tÞ ¼ @uðx;tÞ @t À ruðxÞ Á rcðxÞ À gðxÞ cðxÞ In this case, the function fðxÞ in Eq. ( 17) is simply the right hand side of Eq. ( 21).The multi-layer method contains at least one color layer defined by the boundary curves gðxÞ.
The function cðxÞ is used to model the strength and direction layers that control the diffusion process.Instead of diffusing to an equilibrium, Li et al. [83] used the heat equation to only diffuse for a fixed time interval.Here, the scalar-valued opacity f : R 2 Â R ! ½0; 1 is diffused: fðx; 0Þ ¼ gðxÞ As before, cðxÞ is the diffusion coefficient and controls the diffusion width.Eq. ( 22) can be reorganized similarly to Eq. ( 21) to form a Poisson equation fðx; tÞ ¼ @fðx;tÞ @t Á 1 cðxÞ when cðxÞ 6 ¼ 0. The scalar-valued opacity f can be further modified to model different stroke effects like watercoloring and oil painting.The opacity distribution is combined with the curve colors to determine the final image uðxÞ.
By extending the formulation from Laplacian equations to Poisson equations, a significantly larger solution space is introduced, providing more control over the curves [60].An additional post-processing step is not needed for the Poisson models as it can be added directly in the diffusion process by manipulating the diffusion coefficients [88].

Ray Tracing
The diffusion can be approximated by 2D ray tracing, similar to the final gathering concept in global illumination.The colored curves can be considered as virtual light sources [14]: where the pixel color uðxÞ is determined by the weighted radiance integral over all directions.For a ray starting from a 2D point x in direction u, the closest intersection point of the ray with a diffusion curve is x i ðx; uÞ, and Lðx i ðx; uÞÞ is its radiance.The normalized weight w ¼ w c Á w d is determined by a curve importance w c that models diffusion barriers [11], [14], and a distance weighting w d ðx i ; xÞ ¼ kx i À xk Àp .When no blur is taking place and the rays are not occluded, the ray tracing formulation is equivalent to mean value interpolation, which is a solution to the Laplacian equation [14], [42], [72].This formulation utilizes shaders for flexible color control and avoids the rasterization problems [14] of the original representation [101].It was later adopted by Pr evost et al. [104] and further extended by Lieng [87] to diffusion points.

Shading Curves
Shading curves [87] are an alternative for the modeling of shading and lighting effects.A shading curve is composed of a B-spline curve geometry and shading profiles at the left and right sides of the curve.During creation, users manually draw areas with the curves in constant tone, and each area is filled with constant color, resulting in a piecewise flat image.The shading effects are then added using shading profiles attached to the left and right sides of the curves, which are rendered as meshes rather than curves, making this a combination of curve-based and mesh-based methods.

Creation
Research on methods to enhance artistic control can be roughly divided into two directions.The first stream studies how users can better interact with existing features, such that they can create or edit primitives more easily.The second stream seeks to enrich the feature set of a certain model to create more expressive results.The fundamental functionality in diffusion curve frameworks is the manipulation of curve control points.Users can manipulate the curve geometry, as well as add, delete, duplicate, or change the values of attributes.In this section, we discuss methods and techniques that further improve artistic control on top of the basic editing functionality for different curve-based models, and we additionally group them into contributions to user interaction or feature enhancement wherever applicable.Fig. 11 displays results of different frameworks.

Harmonic Equation
User Interaction.Sun et al. [119] proposed a multi-touch sketching interface for creating and editing diffusion curves.Jeschke et al. [68], [71] adopted a "click and drag" interaction.Sun et al. [120] enabled seamless cloning, where users can copy-paste control curves within a region.Lu et al. [94] used depth information to group curves semantically into objects, which are edited collectively.Later, they enabled users to add arbitrary handles to diffusion curves, which achieves global shape manipulation through linear blending deformation [93].Bao and Fu [4] developed a scribble-based interface for editing diffusion curves, where users draw scribbles and colors are automatically assigned to curve points such that the rendered result aligns with the input.Feature Enhancement.As diffusion always strives for smooth regions, Jeschke et al. [71] proposed to add texture parameters in a procedural way for creating finer details using Gabor noise.The noise parameters were automatically estimated to fit an input image with associated diffusion curves.This led to the application of stippling and hatching processes.Sun et al. [120] used Gaussian radial basis functions as diffusion points in their fast multipole representation of the harmonic equation.This enabled users to add non-harmonic color fields to the image.

Biharmonic Equation
User Interaction.The biharmonic formulation enables users to change not only the color values, but also the color derivatives in the curve normal direction.Finch et al. [41] allowed users to sketch curves freely and let them add color constraints at different points.Curves from vectorization are often too complicated to be edited directly.Xie et al. [138] improved this with their hierarchically extracted curves.Users can view the multi-scale curves in order to manipulate curves in different resolutions.
Feature Enhancement.Finch et al.'s [41] representation enables user to add point values, critical points and five curve types, including value curves, tear and crease curves, and contour and slope curves, which specify different boundary conditions.Boy e et al. [15] generalized the formulation from Finch et al. [41] and support diffusion, barrier, tear, crease, crease-value, value, value-slope, and slope curves, as well as nil curves.Ilbery et al. [64] used sharpprofiles and smooth-profiles to control colors at boundaries and in the interior.

Poisson Equation
User Interaction.Orzan et al. [101] let users create diffusion curves manually by sketching.Zooming and panning is interactive when only diffusing at low resolution in the visible viewport.Lin et al. [89] provided an implementation for mobile devices.In order to make the selection of curves easier, Lin et al. [88] labeled the diffusion curves during the rasterization process, such that when any pixel that belongs to a curve is selected, the entire curve can be immediately chosen.To achieve seamless cloning, Hou et al. [61] added Laplacian constraints to the places where the curves intersect with each other, such that users do not have to manually break diffusion curves into short disjoint segments.
Feature Enhancement.Be zerra et al. [11] added further constraints, including barrier curves that do not emit colors but block colors from other curves, they supported anisotropic diffusion, and smoothly blended between regions via soft constraints.Lin et al. [88] achieved this similarly by solving the diffusion process with diffusion coefficients.Further, they blended multiple layers and added diffusion point sources.Hou et al. [61] introduced Poisson curves and Poisson regions, which enable effects like core shadows, specular highlights, halos, and translucency.They separated hues and tones explicitly so that users could edit hue or tone locally, giving the possibility of editing specular highlights.Li et al. [83] focused on generating stylized images like watercolor or oil paint by introducing different strokes styles.

Ray Tracing
User Interaction.Lieng et al. [85] control how colors are propagated locally around a diffusion point by drawing line constraints for the triangulation near the point.
Feature Enhancement.Using shaders, Bowers et al. [14] added gradient fill and texture to diffusion curve images, and provided parameters for the influence w c and the distance-weighting w d in Eq. (24).Pr evost et al. [104] have both curve-domain shaders (attached along the curves, similar to Bowers et al. [14]) and image-domain shaders (defined in the whole image domain) to avoid undersampling due to sparse interpolation.A multi-layer scheme is implemented in the framework, such that users can add local curves (for example to create specular highlights) and clone objects in the scene.Similar to Pr evost et al. [104], Lieng et al. [85] added local diffusion points in a particular image layer.

Rasterization
Curve-based vector graphics are frequently formulated as solutions to partial differential equations.Hence, two components are essential in the rasterization process.The discretization of the curves and the domain determines the rendering quality.Especially for diffusion curves, a small discretization error may result in large shifts after the diffusion process [69].The solving scheme determines the speed of the rendering result.In this section, we discuss the contributions made to both components.

Harmonic Equation: Discretization
Based on the type of grid that the solvers are applied on, we divide the discretization step into the following types.
Regular Grid.Jeschke et al. [69] proposed a robust rasterization scheme.The curves are first discretized into smaller line segments, and a Voronoi diagram is constructed [57].The color of each pixel is determined by its closest point on the curves, which sets the initial condition.The initial guess and a distance map are used for the diffusion.By diffusing from an entire image (the initial condition) instead of from discrete curves, the discretization error is reduced.Although at pixels near curves, strobing artifacts may still occur [69].
Triangular Grid.Pang et al. [102] rendered diffusion curve images on a triangulated domain.They first discretized diffusion curves and the image boundaries into line segments, then they triangulated the domain using a constrained Delaunay triangulation.Instead of solving at each pixel, they solved for every vertex.The triangulated space was also used as intermediate representation in FEM solvers [15].Boy e et al. [15] applied a constrained Delaunay triangulation, and used quadratic Lagrange polynomials as basis functions.Zhao et al. [143] used a similar triangulation method from the Triangle package [115] with secondorder basis functions.
Irregular Grid.Dai et al. [26] used image segmentation to avoid the discretization of individual curves.Their curves are extracted from a set of superpixels of an input image.The left and right pixels adjacent to a superpixel boundary are considered as the left and right color sources of a curve.Therefore, overlapping of color sources would not occur during the diffusion process.

Harmonic Equation: Solving Scheme
Closed-Form.Sun et al. [121] derived a closed-form solution for closed diffusion curves using Green's functions: with x 0 2 @V and where Gðx; and where the curve normal n points outwards the domain, Gðx; x 0 Þ is Green's function, and G n ðx; x 0 Þ is its normal derivative.In this formulation, every point in the image can be evaluated directly as boundary integral in Eq. ( 26).In practice, Eq. ( 25) is evaluated with boundary element methods (BEM) by sampling points on each diffusion curve.Adaptive sampling methods select a random position on the curve, and iteratively add more samples such that the approximation error is below a threshold.Sun et al. [120] reformulated Eq. ( 25) to a fast multipole representation.Using a hierarchical lattice on the image domain, they calculated boundary integrals hierarchically only on curves in adjacent lattice cells to speed up the computation.Iterative Scheme.The harmonic equation can be solved iteratively using Jacobi relaxations.In each iteration k, the discrete pixel colors u i;j of pixel ði; jÞ are update d: Solving this for every pixel is slow when the required resolution is large.Hence, Orzan et al. [101] used a multi-grid solver, see later in Section 4.3.5.Jeschke et al. [69] designed a GPU-based variable stencil size solver such that a bigger step size can be used.After rasterization, a closest point map and an initial guess of the final image are generated and fed into the solver.In this setup, the Laplacian is not computed as an average of its 4 neighboring pixels, but can be treated as an average of its surrounding circle according to the mean value theorem [42], [69].The circle radius can be increased as long as boundary curves are not crossed, which is determined by the closest point map.Alternatively, Dai et al. [26] formulated the diffusion process as a random walk process.Direct Solvers.Pang et al. [102] solved the diffusion process by directly evaluating the colors on a triangulated domain.For curve nodes, colors are defined on the curve.For other triangle vertices, colors are interpolated using mean value coordinates (MVC) [42] from the colored curve nodes.MVC algorithms assume that nodes are surrounded by a closed polygon.Therefore, an additional visibility test is performed to select the nearest curve nodes around each non-curve node before evaluating the MVC.Colors inside each triangle can then be computed simply using barycentric interpolation.Alternatively, Boy e et al. [15] and Zhao et al. [143] used the finite-element method to solve the harmonic equation linearly in the weak formulation.

Biharmonic Equation: Discretization
In addition to curve geometry and colors, other attributes used in biharmonic formulations also need to be discretized.The primitives can be rasterized onto a uniform grid or they can be discretized onto different elements, using for example finite elements (FEM) or boundary elements (BEM).
Regular Grid.The formulation of Finch et al. [41] contains a set of critical points and feature curves.They rasterize each critical point with respect to the bilinear weights of its four neighbor pixels.The curves are rasterized by computing the intersections of the curves and the edges of the pixels.Since the feature curves are used to control local gradients, these constraints are added at the intersections.
Triangular Grid.Boy e et al. [15] discretized the domain into non-overlapping quadratic triangular patches for both the harmonic and biharmonic formulation.However, two additional conditions need to be fulfilled for biharmonic equations.First, the weak formulation of biharmonic equations requires C 1 continuity across conforming elements.Therefore, unlike the harmonic formulation, quadratic Lagrange patches with C 0 continuity cannot be directly used as FEM elements.Boy e et al. [15] use a third-order non-conforming element (Fraeijs de Veubeke (FV)) [28] such that C 1 continuity is not needed for convergence.Second, biharmonic equations enable gradient constraints, which can be attached or interpolated at the vertices and edges when using FV as FEM elements.
Irregular Grid.Ilbery et al. [64] discretized curves into line segments through recursive subdivision with a threshold on curve straightness.The constraints on color values and color derivatives in normal direction are computed by interpolation for each line segment and are fed into a BEM solver for line by line evaluation.

Biharmonic Equation: Solving Scheme
Closed-Form.Weber et al. [129] and Ilbery et al. [64] derived a closed-form solution using Green's functions: where x 0 2 @V is a boundary point.The Green's function G H and its normal derivative G H n are the fundamental solutions to a harmonic equation and are the same as in Eq. ( 26).The Green's function G B and G B n for biharmonic equations are: In order to improve computation efficiency, a curve-aware upsampling scheme may be used.As the line segment field becomes smoother when the distance between a point and a line segment increases, its value can be approximated through interpolation.In a multi-level formulation, most computation is done on a coarse grid.Iterative Scheme.Finch et al. [41] solved on a smaller grid and upsampled the result using a hierarchical strategy [13].When solving on a coarse domain, they adopt a direct solver (see below).This solution is then upsampled to a finer level, and Gauss-Seidel iterations are performed to refine the solution at that level.
Direct Solvers.Finch et al. [41] used a direct solver for small images.They structured their biharmonic formulation as a linear system Ax ¼ b.The final image is solved by calculating the banded-diagonal Cholesky decomposition (CD) of the matrix A. In the cases when not sufficiently many value constraints are added, such that A is positive definite rather than positive semi-definite, a small regularization weight is added to obtain a unique solution.Boy e et al. [15] used the weak formulation for biharmonic equations from Lascaux and Lessant [75].The final image can be approximated as a linear system, which is solved directly.

Poisson Equation: Discretization
Regular Grid.The original diffusion curve formulation [101] requires to rasterize the color sources, the color derivatives across the curves, and an additional blur source, see Fig. 12.The color curves are rasterized with an offset in the normal direction to avoid overlapping in the rasterized results.However, overlaps can still occur, especially at high curvatures, thin structures, and intersections [101].
Triangular Grid.Zhao et al. [143] used a finite element method to solve the harmonic equation for their diffusion curves.In addition, their method requires solving a Poisson equation during the vectorization process for the cost function.They discretized the domain similarly as described for the harmonic formulation, and the velocity attribute (for vectorization) is discretized along polylines.
Irregular Grid.Hou et al. [60] first decomposed the image domain into disjoint sub-regions, such that each subdomain is bounded by a closed diffusion curve.For each sub-domain, they discretized the primitives and constraints on a quad-tree, which was used for their closed-form solver.

Poisson Equation: Solving Scheme
The solving process of the general Poisson formulation in Eq. ( 17) varies for different definitions of fðxÞ and gðxÞ.
Closed-Form.Like for harmonic and biharmonic functions, Green's identities can be used to derive a closed-form solution to Poisson equations.Hou et al. [60] used Green's third identities to structure the solution to their formulation: When r 2 uðx 0 Þ ¼ 0, integral 1 vanishes, and the equation is reduced to Eq. ( 25): the closed-form solution to a harmonic equation.The definition of G and G n are the same as in Eq. ( 26).Intuitively, Eq. ( 30) can be directly evaluated as was done for the harmonic formulation.However, integral 1 is defined on the entire domain, which is expensive to compute.Hou et al. [60] discretized the domain onto a quad-tree and derived an approximation for each cell.Li et al. [83] used a Fourier transform to solve their heat equation formulation in Eqs. ( 22)-( 23).Therefore, for any 2D point x ¼ ðx; yÞ, the scalar-valued opacity f at time t can be evaluated as: The closed-form solution provides a direct result for any time t in the diffusion process.The calculated result is then rasterized as a whole with respect to the desired resolution.Iterative Scheme.Similar to Eq. ( 27), the original formulation from Orzan et al. [101] can be solved with Jacobi relaxations where for iteration k, the image u kþ1 is: on curve u k iÀ1;j þu k iþ1;j þu k i;jÀ1 þu k i;jþ1 ÀdivðrwÞ i;j 4 otherwise 8 <

:
The divergence of the gradient field divðrwÞ is To achieve interactivity, they used a multi-grid approach [17], [48] and recursively decreased the domain resolution to solve on coarser domains, and upsampled until the original resolution is reached.Similarly, Lin et al. [88] solved their Poisson formulation in Eq. ( 21) for each color layer iteratively.Direct Solvers.Bezerra et al. [11] reformulated the Poisson equation in Eq. ( 7) as a constrained optimization to add more control over the diffusion.Due to the existence of nonlocal constraints (e.g., two closed areas can be linked to enable diffusion between them), the multi-grid solver generates unsatisfying results and can be inefficient.They used a global linear solver instead.Alternatively, finite element methods can be used to solve harmonic equations by solving a linear system iteratively or directly, cf.Section 4.3.1.Zhao et al. [143] used a second-order basis to solve the weak form of a Poisson equation linearly.

Ray Tracing: Discretization
Currently, two domains have been studied for ray tracing: a regular grid or a triangulated domain.For regular grids, the ray starts from each point in the domain, while the ray originates from each triangle vertex in a triangulated domain.
Regular Grid.Bowers et al. [14] subdivided the curves into a set of line segments for their ray tracing algorithm.They built a uniform acceleration grid with a grid size 2 ffiffiffi n p Â 2 ffiffiffi n p for n line segments, such that an increase in the number of curves would not cause a drastic decrease in performance.Triangular Grid.Pr evost et al. [104] solved the ray tracing process on a triangulated domain.Similarly, they discretized curves into line segments, and triangulated the image domain with a constrained Delaunay triangulation.Lieng [85] used the same concept for the diffusion points extension.For every triangle, evaluation points are sampled to compute the final result.For pixels close to a diffusion point, evaluation points receive the diffusion point color.

Ray Tracing: Solving Scheme
The ray tracing formulation is solved stochastically.Bowers et al. [14] performed parallel ray tracing on the GPU for each pixel.The angles of the rays are determined through stratified sampling.Simple adaptive sampling can also be applied for acceleration.This can be done by first dividing the image into a set of blocks and performing ray tracing with fewer rays.During this process, the shortest intersection distance is stored for each block.If this distance is above a threshold, ray tracing is further applied to avoid undersampling.Pr evost et al. [104] performed ray tracing on each triangle vertex.Unlike Pang et al. [102] who used barycentric interpolation, cubic interpolation is performed to avoid mach banding effects.Lieng et al. [85] performed ray tracing similar to Pr evost et al. [104].To add local diffusion points, they render a Catmull-Clark subdivision surface [98].The final result is composed by blending the local and global results linearly.Recently, novel mesh-less stochastic PDE solvers [110] have been introduced for vector graphics [106] to avoid discretization.

Anti-Aliasing
Anti-aliasing techniques have been added in some of the methods to reduce aliasing artifacts.Bowers et al. [14] used a classic technique in their ray tracing formulation, such that the origin of each ray is set as a 2 Â 2 jittered pattern.Pang et al. [102] performed aliasing reduction on curve boundaries by using an accumulation buffer.Sun et al. [121] proposed a random access solver that treats the image plane as a continuous domain.They evaluate the integral over a small rectangular region instead of at each individual point.This formulation, however, may cause a performance drop for the fast multipole formulation [120].Instead, Sun et al. [120] used supersampling with 9 samples per pixel.

Vectorization
Current research on the extraction of curve-based representations from images is mainly focused on diffusion curves, especially the original, harmonic and biharmonic formulations.Since the original formulation and the harmonic formulation extract the same attributes, we cover them together.

Basic Formulation and Harmonic Equation
The simplest harmonic formulation comprises curve geometry, color control points, and usually blur control points.
Geometry.Diffusion curves are based on the idea that edges can represent the majority of color variations in an image [73], [97], [101].Therefore, edge detection often serves in the initialization.Orzan et al. [101] used a Canny detector in the Gaussian scale space (the collection of Gaussian blurred images) on the input image, which captures local maxima in the gradient domain for different scales.A scale-space analysis [34], [90] is performed to select the best scale for each edge.This approach was used by several follow-up works [82], [94], [143] as starting point, too.However, as discussed by Orzan et al. [101], it is hard to detect edges in smooth regions, and these areas are often extracted as several smaller edges, which necessitates a smoothing post-process to represent smooth regions.
Instead of using a Canny detector, Dai et al. [26] segmented the input image into several superpixels, which were then merged into a single superpixel map with k-labeled regions through bipartite graph partitioning.The curves were later extracted as superpixel boundaries.Smooth regions were segmented into several superpixels, and hence the reblurring process was unnecessary.Zhao et al. [143] used iso-contours of the difference between the reconstructed image and the original image as the initial curves in an optimization step.The curves were iteratively refined through gradient descent until the reconstruction error is small enough.In addition to pixel images, the approach can also be used to generate diffusion curve representations for other color fields such as 3D renderings and gradient meshes, where the initial reconstruction is created using object contours or triangle edges.Lu et al. [94] used an additional depth image to improve the object contour detection.A Canny detector is used on the multi-scale color image, while an enhanced cartoon edge detector [24] is applied on the depth map.The edges of both images are combined to generate diffusion curves, where coarse edges in depth maps represent object contours and edges in the color map represent color details.Li et al. [82] extended the vectorization process from images to videos.In each frame, a Canny detector is used and curves are classified as salient and non-salient curves based on the average gradient along the curves.Non-salient curves have a smaller average gradient and are the ones that cause the flickering artifacts in videos.By assuming that each curve is a rigid body, non-salient curves are propagated from one frame to another frame using optical flow.The curves are eventually transformed to B ezier curves.
Color.Given a curve geometry, the color control points are computed next, e.g., by sampling color points on the left and right sides of the curves in the original image [143].After uniformly sampling color points, Orzan et al. [101] simplified the piecewise linear color curves using Douglas-Peucker [31] with color differences measured in the perceptually uniform CIE L*a*b color space.Jeschke et al. [71] first computed a dense set of color control points using a geometric heuristic and removed outliers similar to Orzan et al. [101].A fitting process is then expressed as a linear system Ax ¼ b, where x are the color points, b are the image colors, and A is the influence of control points onto the input image.The influence matrix A is computed during a single diffusion process and the optimal color points are calculated using linear least squares.Later, Jeschke et al. [68] further optimized the computation of the influence matrix, leading to an easier implementation and a more accurate reconstruction.The new scheme keeps a fixed set of the most influential color control points for each pixel.To achieve temporal coherence in videos, Li et al. [82] sampled colors on the left and right side of salient curves as an initial guess, and iteratively refined color points by minimizing a color energy.
Blur.Harmonic equations often require a smoothing step to achieve smooth edges, cf.Section 4.1.2.Thus, blur control points have to be extracted during the inverse process.Orzan et al. [101] computed Gaussian blurred images of the input image and selected the best scale to represent the edges through a scale-space analysis.This ideal scale is used as blur value for the edge pixels.The blur control points are then sampled in a way similar to the color control points.Jeschke et al. [68] first estimated blur values locally at each piecewise linear curve segment.The blur values are then fit to blur control points using a least squares optimization.
Noise.In addition to blur control points, Jeschke et al.
[71] estimated a Gaussian distribution of Gabor noise parameters to mimic textures in an input image.

Biharmonic Equation
The inverse problem for biharmonic equations is harder than for harmonic equations, as more constraints are required.So far, only Xie et al. [138] proposed a vectorization method.
Geometry.Xie et al. [138] have shown that color variations in the Laplacian domain cannot be properly described by edges extracted in the gradient domain.To solve this, they applied a Canny detector on the Laplacian and bi-Laplacian domain of the input image hierarchically.This multi-scale representation avoids the smoothing process and extends the application to biharmonic diffusion curves.As derived by Ilbery et al. [64], Laplacian curves can be considered a special case of bi-Laplacian curves.Therefore, a naı ¨ve approach is to directly solve for all curves in the bi-Laplacian domain.However, it is faster to solve for Laplacian curves than for bi-Laplacian curves.Therefore, Xie et al. [138] first extracted curves in the Laplacian domain, and used a voting method to determine if a curve is a Laplacian curve.Then, they extracted curves in the bi-Laplacian domain and subtracted the unclassified curves extracted in the Laplacian domain.The subtracted bi-Laplacian curves and the classified Laplacian curves are the output diffusion curves.
Attributes.Since Green's formulation from Sun et al. [121] and Ilbery et al. [64] is adopted, the curve color, the normal derivative of color, the Laplacian along curves, and the normal derivative of the Laplacian are solved with a direct linear system.Xie et al. [138] reduced the size of the system matrix by linearly interpolating most of the weights, using a variable stencil size similar to Jeschke et al. [69], solving for both Laplacian and bi-Laplacian together, and by using B ezier splines to share control points.

CONCLUSION AND OUTLOOK
We provided a comprehensive overview of state-of-art methods for the representation, creation, rasterization, and vectorization of mesh-based and curve-based vector graphics.Although mesh-and curve-based methods have conceptual differences in their underlying mathematical representations, there are potential synergies in their rasterization and vectorization processes, as well as in the interaction techniques that enable artistic control.Next, we point out a number of opportunities for future work.
3D Vector Graphics.This survey concentrated on meshbased and curve-based methods for smooth 2D vector graphics.It is mentioned in several works [11], [14], [26], [69], [70] that their representations may be extended to 3D (2D animation or 3D image) or even 4D (3D + time).Compared to 3D images, 2D animations or vector videos have been studied more intensively.One popular commercial tool for creating and editing vector graphics animations is Adobe Flash, which enables users to manipulate with simple shapes and gradients.Video vectorization algorithms have also been proposed for both curve-based [82] and mesh-based [128] representations.On the other hand, the inverse processes, such as the creation and rasterization of vector videos, are less researched.Takayama et al. [124] extended diffusion curves in 2D to diffusion surfaces in 3D.Similarly, Lu et al. [65] proposed a tool to create 2D paths that can be joined to form meshes.However, current results are still much less expressive compared to the achievements in 2D images, and so far the concept of 4D vector graphics remains untouched to our best knowledge.
Vectorization.The vectorization process is a popular research topic in mesh-based and curve-based methods.Both representations use similar steps.For many meshbased and all the curve-based representations, current algorithms follow a two-step scheme.The first step extracts the geometry from an input raster image, where edge and corner detection methods are usually adopted to obtain the image features.Mesh-based methods will then perform a patch fitting or triangulation process to locate the patch vertices, while curve-based methods usually fit a parametric curve or spline to the detected features.It is also common that the geometry may be further optimized through a remeshing or shape optimization process to reduce the construction error [84], [144].The second step is then to extract the attributes of the geometry.Both representations sample colors from the input image.At present, there is only one method designed for biharmonic curves, while all the other methods focused on the original or harmonic representations.Further, there is no vectorization method for the Poisson representation yet.Contrarily, to solve for more attributes, the gradient mesh representation computes both the geometry and the attributes in one optimization step.
Artistic Degrees of Freedom.At present, users may freely edit and control color and shape through control points, which is a rather fine-grained user control.In raster graphics, plenty of research went into the artistic editing of appearance, lighting and materials [111].It is well imaginable that similarly more such high-level editing tools will be developed for vector art.Especially in stylized images, where the shading might not be physically plausible, relighting and appearance editing are challenging tasks.
Sketch-Based Systems.Currently, mainly curve-based methods support free-hand drawing of the curve geometry.The input and adjustment of attributes are more restricted to a select-and-change format for both curve-and mesh-based representations.It would be interesting to extend sketchbased input for all formulations, such that attributes can be added in a more intuitive way.In addition, the ideas from sketch-based systems and applications may also be transferred to vector graphics [62], [133].
Data-Driven Methods.So far, data-driven methods have been intensively studied for generating new vector-based images or vectorizing a given image, mainly for simple mesh-based representations, especially in the SVG format.For example, DeepSVG [20] presented a hierarchical transformer-based network to generate a set of drawing commands that can be rendered into SVG images.SVG-VAE [92] generated similar drawing instructions through a variational autoencoder.Later, Reddy et al. [107] proposed a generative model that does not require direct supervision through a differentiable rasterizer [80].So far, the training data are often restricted to stylized images [107], [108].Although this has been extended to photorealistic images [95] recently, the more complex vector representations like diffusion curves and gradient meshes are still not covered.One reason is the lack of data, as these representations are less available to artists and automatic vectorization results often lead to too many primitives.It is then not easy to generate sufficient ground truth data that has a good balance between image expressiveness and low number of primitives.Another reason might be the richness of the attributes in these representations, as more parameters need to be properly learned to generate convincing results, which also adds complexity in the training model.
Standardization.The SVG format has become industry standard for the representation of basic vector graphics primitives.For the many curve-based and mesh-based formulations, no such default format is available yet.Standardization will be an important step on the road to a wider adoption, because content creation tools and rendering APIs will need a common basis for communication.
Synergies.Curve-based and mesh-based methods have both their individual strengths and weaknesses.While editing tools are more advanced for curve-based methods, vectorization is more matured for mesh-based techniques.First hybrid techniques have already been proposed.Lieung et al. [87] introduced the idea of shading curves with which users can draw lines on flat images, while the shading effects are rendered as Catmull-Clark surfaces along the curves.Chen et al. [23] embedded image features into parametric patches and solved for the color through thin plate spline interpolation to enable real-time vectorization.In the future, it will be rewarding to develop new mathematical models that unify curve-based and mesh-based methods to leverage the benefits of both formulations.
Applications.Vector graphics have been adapted in other areas in visual computing, such as for textures [70], [121] or height fields [56].The compactness of diffusion curve representations can also be used as an intermediate representation for image compression in data-driven methods.For example Poisson vector graphics has been used in the color transfer for faces and portraits [44], [45].With more advances in the representation, creation, rasterization, and image vectorization in the future, we look forward to more exciting applications of vector graphics in other domains.

Fig. 11 .
Fig. 11.Images created with different curve-based frameworks, along with visualizations of input curves and control points.

TABLE 1 An
Overview of Existing Mesh-Based Methods

TABLE 2 An
Overview of the Existing Curve-Based Methods 1660 IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. 30, NO. 3, MARCH 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
TIAN AND G € UNTHER: SURVEY OF SMOOTH VECTOR GRAPHICS: RECENT ADVANCES IN REPRESENTATION, CREATION,... 1665 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.