A Color Elastica Model for Vector-Valued Image Regularization

In image processing and computer vision, a fundamental question related to denoising and enhancing the quality of given images is what is the appropriate metric. In the literature, image regularization for gray-scale images has been extensively studied. One class of such models takes advantage of the Euler elastica energy defined by [19, p.63]

where $v:\Omega\rightarrow\mathbb{R}^{+}$ is a gray-scale image given as a function on a bounded domain $\Omega\subset\mathbb{R}^{2}$, $d\mathbf{x}=dx_{1}dx_{2}$, with $x_{1},x_{2}$ the coordinates of the generic point $\mathbf{x}$ of $\Omega$, while $a$ and $b$ are two positive scalar model parameters.

Euler’s elastica has a wide range of uses in image processing, such as image denoising [21, 38, 47], image segmentation [2, 11, 46, 49], image inpainting [35, 38, 45], and image segmentation with depth [13, 27, 48], to name just a few. An important image denoising model, incorporating the Euler’s elastica energy, is

where $f:\Omega\rightarrow\mathbb{R}^{+}$ is the input image we would like to enhance and denoise. To find the minimizer of (2), one class of methods is based on the augmented Lagrangian method. In [38], an augmented Lagrangian method was proposed for image denoising, inpainting and zooming. Based on the method in [38], [12] suggested a fast augmented Lagrangian method and [47] proposed a linearized augmented Lagarangian method. A split-Bregman method was suggested in [45] to solve a linearized elastica model introduced in [1]. Recently, an almost ‘parameter free’ operator-splitting method was presented in [9]; this method is efficient, robust, and has less parameters to adjust compared to augmented Lagrangian based methods.

A color image is a vector-valued signal which can be represented, for example, by three RGB channels or four CMYK channels. While there are many papers and models that deal with gray-scale image regularization, most of them can not be trivially extended to handle color images when taking the geometric properties of the color image into consideration. Inspired by an early paper [10] that focused on regularization functionals on vector-valued images, [34] proposed an anisotropic diffusion framework and [5] introduced a color TV norm for vector-valued images. Another generalization of the scalar TV, known as the vectorial total variation, was proposed in [33] and studied in [20]. A total curvature model was proposed in [39], in which the color TV term in [5] is replaced by the sum of the curvatures of each channel. A generic anisotropic diffusion framework which unifies several PDE based methods was studied in [40]. The Beltrami framework was proposed in [36], it considers the image as a two-dimensional manifold embedded in a five-dimensional space-color $(x,y,R,G,B)$ space. The image is regularized by minimizing the Polyakov action [28], a surface area related functional. Evolving the image according to the Euler–Lagrange equation of the functional gives rise to the Beltrami flow. To accelerate the convergence rate of the Beltrami flow, a fixed-point method was used in [3] and a vector-extrapolation method was explored in [29]. In [37], the authors used a short-time kernel of the Beltrami flow to selectively smooth two-dimensional images on manifolds. Later on, a semi-implicit operator-splitting method [30] and an ALM based method [31] were proposed to minimize the Polyakov action for processing images. In [50], the authors used a primal-dual projection gradient method to solve a simplified Beltrami functional. A fidelity-Beltrami-sparsity model was proposed in [42], in which a sparsity penalty term was added to the Polyakov action. A generalized gradient operator for vector bundles was introduced in [4], which is applied in vector-valued image regularization. In [6], the Beltrami framework has been applied in active contour for image segmentation. Finally, in [32, 43] the Beltrami framework was extended to deal with patch based geometric flows, or anisotropic diffusion in patch space.

Most existing models for vector-valued image regularization contain only first order derivatives, like the Polyakov action [36], color TV [5] and vectorial TV [20]. This maybe insufficient when considering image formation models containing geometric properties, like the curvature as captured in the single channel Euler’s elastica model. The Polyakov action was shown to be a natural measure which minimization leads to selective smoothing filters for color images. To enrich the variety of image formation models, here, we propose a color elastica model which incorporates the Polyakov action and a Laplace–Beltrami operator acting on the image channels, as a way to regularize color images. While the Polyakov action measures the area of the image manifold, the surface elastica can be computed by applying the Laplace–Beltrami operator on the color channels. It describes a second order geometric structure that allows for more flexibility of the color image manifold and could relax the over-restrictive area minimization. The new model is a natural extension of the Euler elastica model (2) to color images, which takes the geometrical relation between channels into account. Compared to the Polyakov action model, the suggested model is more challenging to minimize since the Laplace–Beltrami term involves high order nonlinear terms. Here, we use an operator-splitting method to find the minimizer of the proposed model. To decouple the non-linearity, three vector-valued and matrix-valued variables are introduced. Then, finding the minimizer of the proposed model can be shown to be equivalent to solving a time-dependent PDE system until the steady state is reached. The PDE system is time-discretized by the operator-splitting method such that each sub-problem can be solved efficiently. To the best of our knowledge, this is the first numerical method to solve the color elastica model. Unlike deep learning methods, the suggested model does not require training, though the model itself could be used as an unsupervised loss in general settings. Numerical experiments show that the proposed model is effective in selectively smoothing color images, while keeping sharp edges.

While the model we propose may not outperform some existing methods like the BM3D [8] and some other deep learning methods [25, 44, 24] on vector-valued image regularization, it manifests properties which have simple geometric interpretations of images. The proposed model provides a new perspective on how to regularize vector-valued images and can be used in many other image processing applications, like image inpainting and segmentation, or as a semi-supervised regularization loss in network based denoising and deblurring. In addition, the model and algorithm proposed here could be also important for many other applications with vector-valued functions that need to be regularized. Many inverse problems, surface processing and image construction problems need this.

The rest of this paper is structured as follows: In Section 2, we briefly review the Polyakov action and the Beltrami framework. The color elastica model is introduced in Section 3. We derive the operator-splitting scheme in Section 4 and describe its finite difference discretization in Section 5. In Section 6, we demonstrate the efficiency, and robustness of the proposed method by systematic numerical experiments. We conclude the paper in Section 7.

In this section we briefly review the Riemannian geometry relevant to the Polyakov action as applied to color images. Since the Polyakov action is adopted from high energy physics, we temporarily follow the Einstein’s summation notation convention in this section. We first introduce the metric on the Riemannian manifold. Denote a two-dimensional Riemannian manifold embedded in $\mathds{R}^{d}$ by $\Sigma$. Let $(\sigma^{1},\sigma^{2})$ be a local coordinate system of $\Sigma$ and $F(\sigma^{1},\sigma^{2})=(F^{1}(\sigma^{1},\sigma^{2}),...,F^{d}(\sigma^{1},\sigma^{2}))$ be the embedding of $\Sigma$ into $\mathds{R}^{d}$. Given a metric $\{g_{ij}\}$ for $i,j=1,2$, of $\Sigma$ which is a symmetric positive definite matrix-valued metric tensor, the squared geodesic distance $(ds)^{2}$ on $\Sigma$ between two points $(\sigma^{1},\sigma^{2})$ and $(\sigma^{1}+d\sigma^{1},\sigma^{2}+d\sigma^{2})$ is defined by,

Assume $F$ embeds $\Sigma$ into a $d$-dimensional Riemannian manifold $M$ with metric $\{h_{ij}\},i,j=1,...,d$. If $\{h_{ij}\}$ is known, we can deduce the metric of $\Sigma$ by the pullback,

where $\partial_{i}$ denotes the partial derivative with respect to the $i^{\mbox{th}}$ coordinate, $\partial_{i}\equiv\partial/\partial\sigma^{i}$.

Consider a color image $\mathbf{f}(\sigma^{1},\sigma^{2})=(f_{1},f_{2},f_{3})$ where $(\sigma^{1},\sigma^{2})$ denote the spatial coordinates in the image domain, and $f_{1},f_{2},f_{3}$ denote the RGB color components. $\mathbf{f}$ can be considered as a two-dimensional surface embedded in the five-dimensional space-feature space: $F=(\sqrt{\alpha}\sigma^{1},\sqrt{\alpha}\sigma^{2},f_{1},f_{2},f_{3})$, where $\alpha>0$ is a scalar parameter controlling the ratio between the space and color. Assuming the spatial and color spaces to be Euclidean, we have $h_{ij}=\delta_{ij}$ on the space-feature space, where $\delta_{ij}$ is the Kronecker delta, $\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ otherwise. Set $\sigma^{1}=x_{1},\sigma^{2}=x_{2}$. According to (3), the metric on the image manifold is

Denote $(\Sigma,g)$ the image manifold and its metric, and $(M,h)$ its space-feature space and the metric on it. The weight of the mapping $F:\Sigma\rightarrow M$ is measured by

for $i,j=1,...,\dim(\Sigma)$ and $\mu,\nu=1,...,\dim(M)$, where $g=\det(\mathbf{G})$ and $\|dX\|_{g,h}^{2}=\partial_{i}F^{\mu}\partial_{j}F^{\nu}g^{ij}h_{\mu\nu}$ with $F=(F^{1},...,F^{\dim(M)})$ and $\{g^{ij}\}$ being the inverse of the image metric $\{g_{ij}\}$. With $m=2,h_{ij}=\delta_{ij}$, the functional is known as the Polyakov action [28] in String Theory.

Substituting $h_{ij}=\delta_{ij}$ and (4) into (5) gives rise to

Minimizing (6) by variational gradient descent, yields the Beltrami flow [23],

where $\nabla F^{k}$ is a column vector by convention. The term on the right hand side of (7) is known as the Laplace–Beltrami operator applied to $F^{k}$.

We assume the domain $\Omega=[0,L_{1}]\times[0,L_{2}]$ to be a rectangle. Periodic boundary condition is a common condition in image processing. In this paper, we assume all variables (scalar-valued, vector-valued and matrix-valued) have the periodic boundary condition. Let $\mathbf{f}=(f_{1},...,f_{d})^{T}$ be the given $d$-dimensional signal to be smoothed, in which each component $f_{k}:\Omega\rightarrow\mathds{R}^{+}$ represents a channel. Based on the Polyakov action (6), we define the color elastica model as,

where $\mathbf{v}=(v_{1},...,v_{d})^{T}:\Omega\rightarrow\mathds{R}^{d}$, $\beta,\eta>0$ are weight parameters, $\mathcal{V}$ denotes the Sobolev space with periodic boundary condition $\mathcal{H}^{2}_{P}(\Omega)$ on $\Omega$:

Similar to (4), let $\mathbf{G}=(g_{kr})_{1\leq k,r\leq 2}$ be a symmetric positive-definite metric matrix depending on $\mathbf{v}$, where,

for some $\alpha>0$. We denote $g=\det\mathbf{G}$. $\Delta_{g}$ in (8) is the Laplace–Beltrami operator associated to $\mathbf{G}$,

Here, we consider the case $d=3$ for color images: $\mathbf{f}=(f_{1},f_{2},f_{3})^{T}$ whose components represent the RGB channels. Then $\mathbf{G}$ is the metric of the surface $F(x_{1},x_{2})=(\sqrt{\alpha}x_{1},\sqrt{\alpha}x_{2},v_{1},v_{2},v_{3}$ induced by the denoised image $\mathbf{v}$.

The energy in (8) is nonlinear and is difficult to minimize, since it contains the determinate and inverse of a Hessian matrix. To simplify the nonlinearity, we introduce three vector-valued and matrix-valued variables. For $k=1,2,3$ and $r=1,2$, let us denote by $q_{kr}$ the real valued function $\frac{\partial v_{k}}{\partial x_{r}}$ and by $\mathbf{q}$ the $3\times 2$ Jacobian matrix

Denote $\mathbf{q}_{k}=\begin{pmatrix}q_{k1}&q_{k2}\end{pmatrix},k=1,2,3$. We introduce the $3\times 2$ matrix $\bm{\mu}=\sqrt{g}\mathbf{q}\mathbf{G}^{-1}$ with

for each row of $\bm{\mu}$. We immediately have $\mathbf{q}_{k}=\frac{1}{\sqrt{g}}\bm{\mu}_{k}\mathbf{G}$. We shall use $\mathbf{M}(\mathbf{q})$ to be the matrix-valued function defined by

and denote $\det\mathbf{M}(\mathbf{q})$ by $m(\mathbf{q})$, that is,

Define $\mathbf{v}_{\mathbf{q}}=\{(v_{\mathbf{q}})_{k}\}_{k=1}^{3}$ to be the solution of

Denote

We define the sets $\Sigma_{f}$ and $S_{\mathbf{p}}$ as

and their corresponding indicator functions as

If $(\mathbf{p},\bm{\lambda},\mathbf{G})$ is the minimizer of the following energy

then $\mathbf{v}$ solving (14) minimizes (8).

In (17), we rewrite the Laplace–Beltrami term as $\nabla\cdot\bm{\mu}_{k}$ where $\bm{\mu}$ linearly depends on $\mathbf{q}$. The metric $\sqrt{g}=\det\mathbf{M}(\mathbf{q})$ is also a function of $\mathbf{q}$. In the fidelity term, $\mathbf{v}$ which solves (14) can be uniquely determined by $\mathbf{q}$. Thus the functional only depends on $\bm{\mu}$ and $\mathbf{q}$ under some constraints on their relations. We then remove the constraints by incorporating the functional with the two indicator function $I_{\Sigma_{f}}(\mathbf{q})$ and $I_{S_{\mathbf{G}}}(\mathbf{q},\bm{\mu})$. The resulting problem (17) is an unconstrained optimization problem with respect to $\bm{\mu}$ and $\mathbf{q}$ only.

In this section, we derive the operator-splitting scheme to solve (17). We refer the readers to [18] for a complete introduction of the operator-splitting method and to [17] for its applications in image processing.