How to represent an image is an old yet still open topic. Currently, two categories, namely raster image and vector graphic, serve in our daily life. A raster image basically stores color values for all pixels (Compressed formats such as JPG and PNG are not considered for simplicity), and vector graphic is a type of geometric representation. Generally, the predominance of vector graphics is in part due to facility for editing, independence on display resolution, and usually volume reduction of storage, etc. Thus, vector graphic is a potential tool for representing images.

Image meshing, as a class of methods for converting raster images into vector graphics, is a popular branch of 'image vectorization'. With image meshing, many further applications, e.g., image wrapping, image composition, and object editing, become facilitated. The main process of image meshing is to find a (usually triangular) mesh covering the whole input image (This means that no pixel in the image would be left), with each vertex given a color value defined as its located pixel in the input image. Then, color values of other positions on the mesh can be expressed by the vertex colors mentioned above (mostly through linear interpolation). In this way, the input image can be rendered totally by the mesh. Naturally, the goal of image meshing is to find such an optimal mesh, to minimize the difference between the input image and its rendered counterpart as much as possible.

To obtain this optimal mesh, various methods have been proposed, most of which employ the Euclidean distance in the color space to measure the difference. In this paper, a hierarchical optimization method is presented. An example of inputs and outputs is illustrated in Fig. 1. The most difficult point here is not the modeling of the problem, but the solution. As will be shown later, these objectives are usually highly non-linear and non-smooth, and direct use of conventional routines in numerical optimization methods may not achieve the desired results, or can hardly make the variables move even a tiny step in the solution space. In numerical optimization, a global optimal solution to a non-convex problem can rarely be obtained. Generally, initialization plays a key role in influencing the performance of non-convex optimization. Apart from the property of the problem itself, which can hardly be improved, a proper initialization of variables is much easier to achieve. To achieve a better initialization, a hierarchical manner is proposed in this paper. As a preprocessing step, the input image is filtered into several levels, ranging from coarse to fine, by, e.g., the bilateral filter, which preserves meaningful features. In this way, the coarsest level becomes much smoother, and thus may be easier to solve than the original one. Then, the result of a coarser level is recycled and employed as the initial value of its finer descendant. As a result, the performance of the optimizer for the finest level, i.e., the original image, is improved.

In addition to a hierarchical process, the objective function at the stage of geometry optimization is still non-convex. Generally, without using 'greedy' schemes, the remaining choices of global optimization are mostly quasi-Newton methods, e.g., L-BFGS (Liu and Nocedal, 1989), and the result is still highly dependent on the initial values. To further reduce the complexity of the problem, a conversion from a non-convex problem to a convex one is proposed. Since the color space is also a metric space, it can be considered together with the image space, which thus forms a new 5D metric space. In this way, a pixel in the image, as well as its attached color, is mapped into a 5D point in the newly constructed metric space. Therefore, the original non-convex problem is converted into a simple linear least squares one, and can be efficiently solved via some common routines, independent of the initialization. The running time can also be reduced.

The contributions of this paper are summarized as follows:

1. A hierarchical process is proposed, which provides an initial input for a finer level with the result of a coarser level to make the original problem smoother.

2. By combining the color space and image space, a conversion from a non-convex problem to a convex one is proposed to further reduce the complexity of the problem at the stage of geometry optimization.

Recently, many algorithms have been proposed on image meshing. The most representative one, proposed by Adams (2011), is a type of greedy algorithm. The key idea is to simply eliminate redundant vertices iteratively from an initial mesh with a denser distribution of vertices. In each loop, the vertex with the least influence on the result of the whole system was removed, and the edges were rearranged; i.e., vertices were re-connected, based on Delaunay triangulation. To accelerate the algorithm, Yang et al. (2003) employed an error diffusion scheme to generate an initial mesh, instead of the initial mesh used in an adaptive thinning method (Demaret and Iske, 2004), consisting of vertices at all pixels. Although this algorithm is fast, it does not involve the quality of the recovered image when rearranging edges (simply because of using Delaunay triangulation), and thus the quality of the generated mesh is not sufficient.

As an improvement, Xie et al. (2014) took the movement of existing vertices into account, and got rid of the constraints of Delaunay triangulation. In Xie et al. (2014), an initial mesh was optimized by performing geometry and topology optimization alternately. Instead of processing a vertex only once as in Adams (2011), Xie et al. (2014) re-adjusted the position of a vertex by a reasonable number of times, so that the local optimality of that vertex would not be corrupted by later processing of other vertices.

Another type of algorithm involves Xia et al. (2009) and Liao et al. (2012): the authors considered the generated mesh as a surface mesh, with a third coordinate of a vertex, i.e., the color of its located pixel. From this viewpoint, an image can be processed by use of some geometry tools. Note that the three-channel colors are processed separately and consistently. As a result, three meshes with the same connectivity are generated for one given image. Although this conversion gives a good result, the reasons are not clearly stated in either study.

Xia et al. (2009) first constructed an initial dense surface mesh, with some crack holes corresponding to the edge features of the input image, then performed a mesh simplification while maintaining the main edge features of the image, and finally converted it into a Bézier triangular mesh. This algorithm employs a thin-plate spline (TPS) to interpolate colors on the mesh, generating a better result. The shortcomings of this algorithm involve the complex representation of the mesh, and the non-adaptive sampling when calculating the coefficients of the TPS.

Rather than a Bézier triangular mesh, Liao et al. (2012) employed piecewise smooth subdivision surfaces. The initial mesh was the same as the one in Xia et al. (2009), and in the stage of mesh simplification, edge features were maintained at different levels to form a multi-resolution mesh representation. Using these features, one can edit the generated mesh easily, with manipulations such as object deformation and feature enhancement. Although this representation is flexible, it is still somewhat complicated (edge features need to be stored). In addition, there is no standard criterion for the quality of the recovered image for either algorithm proposed by Xia et al. (2009) or Liao et al. (2012). The difference from our algorithm is that these two methods are based mainly on non-global mesh operations, while in ours, a global scheme is employed.

Still, there are many other image meshing methods, such as Demaret et al. (2006), Lecot and Levy (2006), and Swaminarayan and Prasad (2006). These methods assume that the color variation in each region is relatively plain, which can also be considered as being generated by linear interpolation. In another type of algorithm, Sun et al. (2007) and Lai et al. (2009) used gradient mesh to represent images. Their results are good, yet the algorithms depend heavily on image segmentation, as they can deal only with isolated objects with smooth color variance.

Given an input raster image $\boldsymbol{I}$, our goal is to generate a triangular mesh $\boldsymbol{M}=\{\boldsymbol{V}, \boldsymbol{K}\}$ covering the image and representing it as faithfully as possible, where geometry $\boldsymbol{V}$ and topology $\boldsymbol{K}$ are described as follows.

For geometry information, let $\boldsymbol{v}_{i}\in \mathbb{R}^{2}$ ($i = 1, 2, \dots, n_{\mathrm{v}}$) denote a vertex of the triangular mesh, where $n_{\mathrm{v}}$ is the number of vertices, and $\boldsymbol{V}\in \mathbb{R}^{2\times n_{\mathrm{v}}}$, a matrix whose ith column represents vertex $\boldsymbol{v}_{i}$, denotes the set of all vertices on mesh $\boldsymbol{M}$. To recover image $\boldsymbol{I}$, each vertex is attached with a color value, and the colors of other positions on the mesh are linearly interpolated by neighboring vertices, namely all vertices of its belonging triangle. Inspired by Xiong et al. (2014), a dictionary learning framework is introduced for modeling both geometry and topology information.

A sparse matrix $\boldsymbol{K}\in \mathbb{R}^{n_{\mathrm{v}}\times n_{\mathrm{p}}}$ is constructed to denote the connectivity of vertices, whose column $\boldsymbol{\kappa_j}$ denotes the barycentric coordinates of pixel $\boldsymbol{p}_{j} \in \mathbb{R}^{2}$ ($j=1, 2, \dots, n_{\mathrm{p}}$), where $n_{\mathrm{p}}$ is the number of sampling points, i.e., the number of pixels. Note that the sparse matrix $\boldsymbol{K}$ should be constructed such that the mesh can form a manifold. That is, each column of $\boldsymbol{K}$ has at most three non-zero elements, corresponding to up to three vertices of a triangular face, and the summation of the magnitudes of these elements should be one. Thus, our goal is modeled as follows:

(1) $\begin{gathered} ({{\boldsymbol{V}}^*},{{\boldsymbol{K}}^*}) = \mathop {\arg {\text{min}}}\limits_{{\boldsymbol{V}},{\boldsymbol{K}}} E({\boldsymbol{V}},{\boldsymbol{K}}) \hfill \\ {\text{s}}.{\text{t}}.{\mkern 1mu} {{\boldsymbol{\kappa }}_{\boldsymbol{i}}}{_0} \leqslant 3,{\mkern 1mu} {{\boldsymbol{\kappa }}_{\boldsymbol{i}}}{_1} = 1,i = 1,2, \ldots ,{n_{\text{p}}} \hfill \\ \end{gathered} $,

where $E(\boldsymbol{V}, \boldsymbol{K})$ measures the fitting error of the recovered image.

Let $\mathcal{I}\in\{f:\mathbb{R}^{2}\mapsto\mathbb{R}^{3}\}$. Then $\mathcal{I}(\boldsymbol{p})$ denotes the color value attached to pixel $\boldsymbol{p}$, and $E(\boldsymbol{V}, \boldsymbol{K})$ is thus expressed as

(2) $\begin{gathered} E({\boldsymbol{V}},{\boldsymbol{K}}) = {\mkern 1mu} \frac{1}{{{n_{\text{p}}}}}\mathcal{I}({\boldsymbol{P}}) - \mathcal{I}({\boldsymbol{V}}){\boldsymbol{K}}_{\text{F}}^2 \\ = {\mkern 1mu} \frac{1}{{{n_{\text{p}}}}}\sum\limits_i {\mathcal{I}({{\boldsymbol{p}}_i}) - \mathcal{I}({\boldsymbol{V}}){{\boldsymbol{\kappa }}_{\boldsymbol{i}}}_2^2} \end{gathered} $,

where $\|\cdot\|_\mathrm{F}$ is the Frobenius norm, $\boldsymbol{P}\in \mathbb{R}^{2\times n_{\mathrm{p}}}$ is the set of all pixels, whose jth column denotes the jth pixel $\boldsymbol{p}_{j}$, and $\mathcal{I}(\boldsymbol{P}) = \{\mathcal{I}(\boldsymbol{p}_{i})\}_{i=1}^{n_{\mathrm{p}}} \in \mathbb{R}^{3\times n_{\mathrm{p}}}$, $\mathcal{I}(\boldsymbol{V}) = \{\mathcal{I}(\boldsymbol{v}_{i})\}_{i=1}^{n_{\mathrm{v}}} \in \mathbb{R}^{3\times n_{\mathrm{v}}}$.

Note that the sparse matrix $\boldsymbol{K}$ mentioned above is constructed to elegantly combine the geometry and topology information in Eq.(2). In this way, the connectivity of vertices is encoded into a matrix, which makes the objective easy to express. Also, note that the model in this study is actually linear, while a more general higher order model might generate a better visual effect. However, a higher order model can make the whole process more complicated, and entails more running time. Thus, a simple linear model is employed.

Once the problem has been modeled, the next task is to solve it. As shown in Algorithm 1, the pipeline simply involves several steps. As a preprocessing step, a hierarchy with several levels is established for the input image. An initial mesh is then generated, followed by the whole iterative process for all levels. At each level, a geometry and topology alternating optimization is performed. To further simplify the problem, a convexification of the objective function is performed beforehand. Some details are discussed in the following paragraphs.

As shown in Eq. (1), the objective function is highly non-linear and non-smooth, and thus direct use of common routines of optimization algorithms, e.g., L-BFGS, is not quite suitable in this context. Since the body of the algorithm itself can hardly be improved for such a non-convex problem, another path should be sought. For non-convex optimization, a proper initialization is significant. For this, a hierarchical process is proposed. Specifically, the input image is filtered into N levels beforehand. Intuitively, the coarser level is smoother than the finer one, and is supposed to have better performance through the same algorithm. Therefore, as illustrated in Fig. 2, the whole process starts from the coarsest level (level N-1), travels through all other levels of images, ranging from coarser to finer, and finally reaches the original input image (the finest level, level 0). When dealing with level i (i<N-1), the result from level i+1 is used as an initial value.

As for the choice of the filter, our goal is to make the problem smoother, while preserving some necessary features in the image. Thus, general smoothing filters, e.g., the Gaussian filter, are not quite suitable. Instead, filters preserving features are preferable. In our algorithm, the simple bilateral filter is chosen. In fact, other filters, such as the filter proposed in Xu et al. (2011), might have a better result, yet are somewhat time-consuming compared with the bilateral filter.

Before the main process of various levels, an initial mesh for the coarsest level, i.e., level N-1, needs to be generated properly. Although our algorithm has convexified the problem at the stage of geometry optimization, the situation at the stage of topology optimization is not improved, and the whole problem is still somewhat complex. Thus, a totally random initial mesh might generate a undesirable result, and involve a longer running time. Therefore, to achieve better performance and effect, the results generated by Adams (2011) are employed as the very beginning initialization in our algorithm.

To further simplify the problem, a convexification of the objective is performed. Specifically, the position of each vertex $\boldsymbol{v}_{i}\in\mathbb{R}^{2}$ on mesh $\boldsymbol{M}$, together with its attached color value $\mathcal{I}(\boldsymbol{v}_{i})\in\mathbb{R}^{3}$, is mapped into a point $\hat{\boldsymbol{v}}_{i}=(\boldsymbol{v}^\mathrm{T}_{i}, \mathcal{I}(\boldsymbol{v}_{i})^\mathrm{T})^\mathrm{T}\in\mathbb{R}^{5}$. Thus, the planar mesh becomes a 2-manifold in a 5-dimensional (5D) space. Similarly, all pixels of the input image $\boldsymbol{I}$ are converted into such 5D points. In this way, the original non-convex problem is converted into a convex one. Thus, Eq. (2) is converted as follows:

(3)$E\left({\boldsymbol{\hat V}} \right) = \frac{1}{{{\eta _{\text{p}}}}}\left\| {\boldsymbol{\hat P} - \boldsymbol{\hat VK}} \right\|_{\text{F}}^2 = \frac{1}{{{\eta _{\text{p}}}}}\sum\limits_i {\left\| {{{\boldsymbol{\hat p}}_i} - \boldsymbol{\hat V{\kappa _i}}} \right\|_2^2} $,

where $\hat{\boldsymbol{p}}_{i}$ is the ith column of $\hat{\boldsymbol{P}}\in\mathbb{R}^{5\times n_{\mathrm{p}}}$, and $\hat{\boldsymbol{v}}_{i}$ is the ith column of $\hat{\boldsymbol{V}}\in\mathbb{R}^{5\times n_{\mathrm{v}}}$.

The rationale for this conversion is explained as follows. In the original problem, the metrics in the image space and color space are both strongly related to the Euclidian distances. Therefore, the combination of these two spaces is possible in form. On the other hand, it makes sense in geometry. In fact, this is a type of approximation to the original problem. Take a 1D image with one channel as an example (Fig.3). Fig.3a shows the color distance in the original problem, where the segments (red line) connecting the recovered image (slope straight line) and original one (dashed curve) are orthogonal to the image plane (the horizontal line below); in contrast, in Fig.3b, this line is orthogonal to the triangular face of the mesh. When the distances become small enough, these two forms are supposed to generate similar effects. Based on the above explanation, this convexification is proposed and employed in our model.

After the preprocessing, the solution for each level can be solved via common routines. In general, to solve this problem, geometry and topology aspects are usually treated separately and alternately in each iteration, as in Xie et al.(2014). Both sub-optimization stages are discussed below.

In this stage, the topology $\boldsymbol{K}$ of the mesh, i.e., connection of vertices, is fixed, and the only variable concerned is the position matrix of vertices, $\hat{\boldsymbol{V}}$. Since the original problem has been converted into a simple linear least squares one (Section 4.4), the preconditioned linear conjugate gradient method is employed to solve $\hat{\boldsymbol{V}}$. Note that the rows of $\hat{\boldsymbol{V}}$ are solved separately using the same factorization of $\boldsymbol{KK}^\mathrm{T}$, which can be done beforehand to speed up the process.

Given a fixed geometry, i.e., positions of vertices $\hat{\boldsymbol{V}}$, the task of this stage is then to optimize the topology, i.e., the connectivity of vertices $\boldsymbol{K}$. As is known to all, topology operations involve edge flip, edge or face collapse, edge or face split, etc. Since the positions of vertices are constant in this context, which means the number of vertices is not changed, the only operator left is edge flip. To further reduce the total energy, the strategy in Xie et al.(2014) is employed, the process of which is described briefly below. For further details, please refer to Xie et al.(2014).

Before flipping the edges, a suspect set of edges is defined, the element in which is an edge indicating that it might be flipped. At the beginning, all edges excluding boundary edges, i.e., interior edges, are inserted into the set. To sort these suspect edges, a fitting error defined as the summation of fitting errors of two adjacent faces is attached to each suspect as a criterion. Then the suspect set becomes a priority queue. Each time, a suspect is popped up and checked whether to flip or not, depending on whether the total energy declines. If a suspect is flipped, its neighbor edges are also labeled as suspects and pushed into the queue. The process ends when the queue is empty. To avoid possible endless loops, the times visited for each edge are recorded, and edges visited enough times are no longer concerned in this pass.

In our experiments, all images are tested on a PC with AMD Phenom^TM II X4 945 CPU and 8~GB RAM. All the algorithms are implemented in an Arch Linux operating system, using C++ language, with the Surface_mesh library (Sieger and Botsch, 2012) employed to represent the mesh data structure.

The peak signal-to-noise ratio (PSNR) (HuynhThu and Ghanbari, 2008) and running time are used to evaluate our algorithm. Here,

${\text{PSNR}} = 20\log \frac{{{2^\rho } - 1}}{{\sqrt {{\text{MSE}}} }} $,

where $\rho$ is the number of bits or samples in source or target images, and $\mathrm{MSE}$ is the mean squared error between the original image and reconstructed image, which is the total fitting error in this context (Section 4.4). Generally, the larger the PSNR, the better the performance of fitting.

The first parameter to be tuned is the number of hierarchy levels. In general, deeper levels might not produce better performance. To verify this, several images are tested via various levels. In our experiments, we find that using two or three levels usually generates better results. The reason for this is that after filtering enough times, the image is smooth enough, and the influence of initialization on the optimization fades.

Other parameters are described here. In bilateral filtering, the diameter of each pixel neighborhood that is used during filtering is set to 9, and the two parameters used to measure the size of convolution kernels in the color and coordinate space, $\sigma_{\mathrm{Color}}$ and $\sigma_{\mathrm{Space}}$, are both set to 50. As for the number of iterations in each level, we find two iterations give a trade-off between performance and speed.

In this section, two algorithms given in Adams (2011) and Xie et al. (2014) are compared with ours. The images for evaluation are selected diversely, ranging from human to natural plants and animals. The results of PSNR are listed in Table 1. As shown, our algorithm performs better than the other two for most of the test images. Note that the PSNR value obtained here using the algorithm proposed by Xie et al. (2014) has some tiny difference from that given in Xie et al. (2014), due to the implementation. In Xie et al. (2014), the algorithm is for single-channel images, while in this study, it is implemented for color images, for convenience of comparison.

To compare the mesh generated, the wireframes of both initial and final results are shown in Fig. 4. %The initial mesh is shown left, while the optimized counterpart is shown on the right. The mesh from Adams (2011) is a Delaunay triangular mesh, which has less ability than ours to express detailed features. Note that both meshes shown are generated by a 1% sampling rate.

To further illustrate the effect of our algorithm, some images rendered by the three algorithms, and the original ones, are shown in Fig. 5 (see p.39), with some details highlighted and compared.

As demonstrated in Fig. 5, our algorithm has a better visual effect than the other two. For instance, the mouth of the frog in the first image is recovered much better, which can be seen in the enlarged regions in the second row. The competitiveness of our algorithm can also be seen in the edge of the leaf.

Among many factors, the main reason why our algorithm performs better is the choice of the initial value. A good initial value makes the problem much smoother, and thus makes it easier to reach a better local optimality. In addition, our algorithm converts a non-linear problem into a linear least squares one at the stage of geometry optimization, and such convexification of the original problem speeds up the whole process.

The running time of our algorithm is listed in Table 2. Although not in real time, it is already acceptable for daily use. Generally, a good initialization leads to a shorter running time, since it locates near the local optimum in the solution space.

Our algorithm has some flaws. For instance, it is not quite suitable for images with complicated textures. Fig. 6 shows such an example of a failure case. The textures in the water and feather are not well recovered. Since the initial meshes are obtained by bilateral filtering, the results from our algorithm lack some details of complicated textures, while the main features are preserved.

We present a novel algorithm to convert a raster image into a vector graphic represented by a triangular mesh. Compared with other algorithms, the results are improved by establishing a hierarchy of several levels for the input image. In addition, a convexification of the problem in the stage of geometry optimization is performed to further simplify the non-linear problem into a linear least squares one, which can be solved easily via common optimizers.

In the future, we will seek a way to obtain a better initial value. For this, the PatchNet technique (Hu et al., 2013) might be a favorable choice, for our algorithm to achieve a better local optimum. Meanwhile, we will consider using a more general non-linear model to represent images. We believe such a model would render results with better quality.