Euler’s Elastica Based Cartoon-Smooth-Texture Image Decomposition

Roy Y. He¹¹1Department of Mathematics, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong. Email: [email protected]. Hao Liu²²2Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Kowloon, Hong Kong. Email: [email protected].

Abstract

We propose a novel model for decomposing grayscale images into three distinct components: the structural part, representing sharp boundaries and regions with strong light-to-dark transitions; the smooth part, capturing soft shadows and shades; and the oscillatory part, characterizing textures and noise. To capture the homogeneous structures, we introduce a combination of $L^{0}$ -gradient and curvature regularization on level lines. This new regularization term enforces strong sparsity on the image gradient while reducing the undesirable staircase effects as well as preserving the geometry of contours. For the smoothly varying component, we utilize the $L^{2}$ -norm of the Laplacian that favors isotropic smoothness. To capture the oscillation, we use the inverse Sobolev seminorm. To solve the associated minimization problem, we design an efficient operator-splitting algorithm. Our algorithm effectively addresses the challenging non-convex non-smooth problem by separating it into sub-problems. Each sub-problem can be solved either directly using closed-form solutions or efficiently using the Fast Fourier Transform (FFT). We provide systematic experiments, including ablation and comparison studies, to analyze our model’s behaviors and demonstrate its effectiveness as well as efficiency.

1 Introduction

In many tasks of image processing, such as denoising, segmentation, compression, and pattern analysis, a given image $f:\Omega\to\mathbb{R}$ defined on a compact domain $\Omega\subset\mathbb{R}^{2}$ can be written as a sum of several components, each of which represents a distinctive feature of $f$ . For example, the well known Rudin-Osher-Fatemi (ROF) model [47] separates a noisy image $f=u+n$ into its noise-free part $u$ and the noisy remainder $n$ by considering

\displaystyle\min_{(u,v)\in\text{BV}(\Omega)\times L^{2}(\Omega)/f=u+n}\left[\frac{1}{2\lambda}\|n\|^{2}_{L^{2}(\Omega)}+\int_{\Omega}|\nabla u|\right].

(1)

Here $\int_{\Omega}|\nabla u|$ is the total variation (TV)-norm of a bounded variation function $u\in\text{BV}(\Omega)=\{u\in L^{1}(\Omega):\int_{\Omega}|\nabla u|<+\infty\}$ , and $\lambda>0$ is a weight parameter. Various extensions of the ROF model are proposed in the literature mainly from two perspectives. The first branch of works is motivated by the fact that $L^{2}$ -norm in (1) can fail to capture noise [41, 3], thus leaving visible oscillations in $u$ . Another class of works is based on the observation that TV-norm leads to various artifacts including staircase effects, and loss of geometry as well as contrast [12, 1, 38, 62]. See Figure 1 (e) for an illustration.

Refer to caption — Figure 1: Motivation for curvature regularuzation on level lines. (a) Clean image. (b) Noisy image ( $\sigma=120/255$ ). (c) Denoised image by $L^{0}$ -gradient. (d) Denoised image by the MC penalty function [59]. (e) Denoised image by TV. (f) Denoised image by Euler’s elastica [42]. The noise is captured by $H^{-1}$ -norm approximating the Meyer’s $G$ -norm. Observe that $L^{0}$ -gradient and MC penalty function produce irregular boundaries, and TV-norm creates staircase effects. These artifacts are reduced in (f) when incorporating the curvature regularization on the level lines.

In his seminal work [41], Meyer suggested to characterize the oscillating patterns of an image including textures and noise by the dual space of the closure of the Schwartz class in $\text{BV}(\mathbb{R}^{2})$ . Later, Aubert and Aujol [3] refined the construction by restricting to the case of open and bounded domain $\Omega$ . They proposed the $G(\Omega)$ -space as a subspace of the dual of the Sobolev space $W^{1,1}_{0}(\Omega)$ , and formulated the Meyer’s model as

\displaystyle\min_{(u,v)\in\text{BV}(\Omega)\times G(\Omega)/f=u+n}\left[\frac{1}{\lambda}\|v\|_{G(\Omega)}+\int_{\Omega}|\nabla u|\right],

(2)

where $\|\cdot\|_{G(\Omega)}$ stands for the $G$ -norm. A key property of $G$ -norm is that, oscillating functions converging to $0$ in $G(\Omega)$ does not have to converge strongly in $L^{2}(\Omega)$ [3], which justifies the benefit of using $G$ -norm in place of $L^{2}$ -norm in the ROF model for noise. As $G$ -norm is related to $L^{\infty}$ -norm, it is computationally challenging to solve (2). Vese and Osher [52] proposed the first practical solution by considering

\displaystyle\min_{(u,\mathbf{g})\in\text{BV}(\Omega)\times\left(L^{p}(\Omega)\right)^{2}}\left[\int_{\Omega}|\nabla u|+\frac{1}{2\lambda}\int_{\Omega}|f-u-\partial_{x}g_{1}-\partial_{y}g_{2}|^{2}\,d\mathbf{x}+\mu\|\mathbf{g}\|_{L^{p}(\Omega)}\right],

(3)

where $d\mathbf{x}$ stands for the Lebesgue measure in $\mathbb{R}^{2}$ , $\lambda>0$ and $\mu>0$ are weight parameters, $p\geq 1$ is an integer index, and $\|\mathbf{g}\|_{L^{p}(\Omega)}=\|\sqrt{g_{1}^{2}+g_{2}^{2}}\|_{L^{p}(\Omega)}$ for $\mathbf{g}=(g_{1},g_{2})$ . They showed that as $\lambda\to 0$ and $p\to\infty$ , $\partial_{x}g_{1}+\partial_{y}g_{2}$ belongs to $G(\Omega)$ . They derived the Euler-Lagrange equation related to (3) and solved the problem after discretization in the case of $p=1$ . In [45], a simplified version of (3) when $p=2$ was studied, which leads to using $H^{-1}(\Omega)$ seminorm in place of Meyer’s $G$ -norm, where $H^{-1}(\Omega)$ is the dual of $H_{0}^{1}(\Omega)$ . Approaching from a different perspective, Aujol et al. [4] proposed to tackle Meyer’s problem by restricting the function space $G(\Omega)$ to a convex subset $G_{\mu}(\Omega)=\{v\in G(\Omega)\leavevmode\nobreak\ |\leavevmode\nobreak\ \|v\|_{G(\Omega)}\leq\mu\}$ for some $\mu>0$ . Analogous to Chambolle’s method [11] for (1), a projection algorithm was proposed in [4] to solve their model. In [53], Wen et al. considered the Legendre-Fenchel duality of the Meyer’s $G$ -norm and proposed a primal-dual algorithm for the cartoon-texture decomposition. In the literature, there are many other types of regularizations for capturing texture and noise. For instance, the Bounded Mean Oscillation (BMO) space was considered in [31, 22], the negative Hilbert-Sobolev space $H^{-s}$ was applied in [33], and a general family of Hilbert spaces were studied in [6]. With a more refined decomposition, Aujol and Chambolle [5] suggested to separate a digital image $f$ as a sum of three components: the geometric part $u\in\text{BV}(\Omega)$ , the texture part $v\in G_{\mu}(\Omega)$ , and the noise component $w$ in the Besov space [41] consisting of functions whose wavelet coefficients are in $\ell^{\infty}$ . A fast filter method as well as a nice summary of variational models for cartoon $+$ texture decomposition is available in [10]. These variational decomposition frameworks have also been extended and applied to color images [7, 32, 20, 54], RGB-d images [15], images on manifold [55], segmented image decomposition [48], and image registration [14].

To address the unsatisfying artifacts of TV-norm [12, 62, 50] in the ROF model, much efforts are devoted to constructing new regularizers that incorporate higher-order features or non-convex terms. The Euler’s elastica energy [42, 39, 49, 1]

\displaystyle\mathcal{R}_{\text{Euler}}(v;a,b)=\int_{\Omega}\left(a+b\left(\nabla\cdot\frac{\nabla v}{|\nabla v|}\right)^{2}\right)\left|\nabla v\right|\,d\mathbf{x}

(4)

has been used for disocclusion [42], inpainting [49], and denoising [1, 58], where $a,b>0$ are weight parameters. See Figure 1 (f) for the reduction of staircase effects by the Euler’s elastica. In [38], Lysaker et al. proposed a regularizer based on the Frobenius norm of Hessian

\displaystyle\mathcal{R}_{H}(v)=\|\partial^{2}v\|_{L^{1}(\Omega)}=\int_{\Omega}\sqrt{v_{xx}^{2}+v_{xy}^{2}+v_{yx}^{2}+v_{yy}^{2}}\,d\mathbf{x}.

Regarding the image as a manifold in the Euclidean space with higher dimension, Zhu et al. [62] proposed to use the $L^{1}$ -norm of the mean curvature as the regularizer

\displaystyle\mathcal{R}_{\text{mean}}(v)=\int_{\Omega}\left|\nabla\cdot\left(\frac{\nabla v}{\sqrt{1+|\nabla v|^{2}}}\right)\right|\,d\mathbf{x}.

From the same geometric point of view, Liu et al. [35] considered the Gaussian curvature

\displaystyle\mathcal{R}_{\text{Gauss}}(v)=\int_{\Omega}\frac{\text{det}\begin{bmatrix}v_{xx}&v_{xy}\\ v_{yx}&v_{yy}\end{bmatrix}}{(1+|\nabla v|^{2})^{3/2}}\,d\mathbf{x}

to reduce image noise while preserving geometric features. Other regularizers substituting the TV-norm include the Total Generalized Variation (TGV) [60, 56], total variation with adaptive exponent [9], non-local TV [23], locally modified TV [61], and Higher Degree TV (HDTV) [28].

The aforementioned two streams of research focus on improving different parts of the ROF model (1), and there are only a few works combining them. In [13], Chan et al. extended [5] by considering the inf-convolution of TV-norm and the squared $L^{2}$ -norm of the Laplacian and arrived at the following model

	$\displaystyle\min_{(v,w,n_{1},n_{2})\in\text{BV}(\Omega)\times H^{2}(\Omega)\times G_{\mu}(\Omega)\times B_{1,1}^{1}}$	$\displaystyle\bigg{[}\int_{\Omega}\|\nabla v\|+\frac{\alpha}{2}\\|\Delta w\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+J^{}\left(\frac{n_{1}}{\mu}\right)+B^{}\left(\frac{n_{2}}{\delta}\right)+\frac{1}{2\lambda}\int_{\Omega}\|f-v-w-n_{1}-n_{2}\|^{2}\,d\mathbf{x}\bigg{]}\;.$		(5)

Here $\alpha>0$ , $J^{*}\left(\frac{n_{1}}{\mu}\right)$ and $B^{*}\left(\frac{n_{2}}{\delta}\right)$ are the Legendre-Fenchel conjugate of the $G$ -norm and the $B_{-1,\infty}^{\infty}$ -norm, which enforce boundedness $\|n_{1}\|_{G(\Omega)}\leq\mu$ and $\|n_{2}\|_{B_{-1,\infty}^{\infty}}\leq\delta$ , respectively, and $B_{-1,\infty}^{\infty}$ is the dual to the homogeneous Besov space $B_{1,1}^{1}$ [41]. In (5), $v$ corresponds to the image’s piecewise constant part, and $w$ can be associated with the smoothly varying part such as soft lighting and blurry shadows. The $n_{1}$ component corresponds to the textures, and the $n_{2}$ component captures the noise. Recently, Huska et al. [29] proposed another three-component decomposition model, where the image $f=v+w+n$ consists of a cartoon part $v$ that captures the piecewise constant structures of the image, a smooth part $w$ that models the uneven lighting conditions in the scene, and an oscillatory component $n$ that includes textures and noise. Specifically, their model is

\displaystyle\min_{(v,w,n)}\left[\beta_{1}\|\nabla v\|_{0}+\beta_{2}\|\partial^{2}w\|_{L^{2}(\Omega)}^{2}+\beta_{3}\|n\|^{4}_{H^{-1}(\Omega)}+\frac{1}{2\lambda}\int_{\Omega}|f-v-w-n|^{2}\,d\mathbf{x}\right],

(6)

where $\|\nabla v\|_{0}$ is the $L^{0}$ -gradient energy [57] that counts the pixels with non-zero gradient, $\partial^{2}w(\mathbf{x})\in\mathbb{R}^{2\times 2}$ for each $\mathbf{x}\in\Omega$ is the Hessian matrix, $\|\partial^{2}w\|_{L^{2}(\Omega)}^{2}=\left(\int_{\Omega}(w_{xx}^{2}+2w^{2}_{xy}+w^{2}_{yy})\,d\mathbf{x}\right)^{2}$ , $\|\cdot\|_{H^{-1}(\Omega)}$ is the seminorm of the negative Sobolev space $H^{-1}(\Omega)$ as an approximation of the $G$ -norm, and $\beta_{i}$ , $i=1,2,3$ are positive weight parameters. For numerical computation, they replaced $\|\nabla v\|_{0}$ by a modified minimax concave (MC) penalty function [59]. Note that the regularization on $n$ is the fourth power of the $H^{-1}(\Omega)$ seminorm, which was heuristically found to be useful in distinguishing $w$ and $n$ [29]. One of the major drawbacks of solely using the $L^{0}$ -gradient regularization is visualized in Figure 1. Since $L^{0}$ -gradient regularization is global and heights of boundary jumps do not affect its value, it can lead to strong staircase effects and loss of geometry [12].

In this paper, we propose a new three-component decomposition model to address the issue discussed above. Our model decomposes any grayscale image into a structural part, a smooth part, and an oscillatory part. For the structural part, our model enforces strong sparsity on the gradient while reducing the presence of staircase effects and encouraging shape-preserving boundary smoothing. We achieve this by combining the $L^{0}$ -gradient energy with the Euler’s elastica (4) which controls the magnitude of level-line curvature. For the smooth part, instead of using the Frobenius norm of the Hessian as in (6), we substitute it with squared $L^{2}$ -norm of the Laplacian, $\|\Delta w\|_{L^{2}(\Omega)}^{2}$ , to encourage isotropic smoothness. It can be shown that $\|\Delta w\|^{2}_{L^{2}(\Omega)}=\|\partial^{2}w\|^{2}_{L^{2}(\Omega)}$ for compactly supported smooth functions [13]. Moreover, we prove that they have the same critical points under certain boundary conditions (see Appendix B). For the oscillatory part, we use the squared $H^{-1}(\Omega)$ semi-norm as a computationally practical substitute for the Meyer’s $G$ -norm. It yields not only effective extraction of the oscillations similarly to the 4-th power in (6) but also significant improvement on the computational efficiency.

Figure 2 illustrates the general performance of our model. The structural part, $v^{*}$ , contains definite contours and shades with sharp transition; the smooth part, $w^{*}$ , captures soft intensity variation on the cheeks and chin; and the oscillatory part, $n^{*}$ , includes noise and textures such as those of the hair. In case the input image does not contain fine-scale textures, the composition $u^{*}=v^{*}+w^{*}$ gives a denoised version of $f$ where the noisy perturbation is removed and sharp boundaries are preserved. The level lines of $u^{*}$ by our model are more regular compared to those without the curvature regularization. This shows that the proposed model effectively reduces the staircase effects caused by the $L^{0}$ -gradient regularization.

Our model involves both non-convex and nonlinear high-order terms, thus designing an effective numerical scheme poses a particular challenge. For the non-convex $L^{0}$ -gradient minimization, various algorithms and approximations are employed in the literature, such as hard-thresholding [57, 18], region-fusion [43], and penalty approximation [29]. In [44], the $L^{0}$ -gradient projection was considered and solved by using a mixed $\ell^{0,1}$ pseudo-norm. Another nontrivial part of our proposed model is the presence of the nonlinear high-order term for the curvatures of level lines. A classical approach is to derive the first variation and then discretize the associated Euler-Lagrange equation [49]. With higher efficiency, various Operator-Splitting Methods (OSMs) and Augmented Lagrangian Methods (ALMs) [24, 26, 50, 35, 27, 29, 37, 17, 58] have been proposed to handle variational models which contain curvature regularization. OSMs are commonly derived from associating an initial value problem with the formal optimality condition, and ALMs are usually developed via the augmented Lagrangian. These approaches are deeply related and sometimes equivalent [24, 26, 21]. The main advantage of these paradigms come with appropriately splitting the original problem into several sub-problems, each of which can be efficiently handled. Recently, connections between OSMs and neural networks are investigated in [30, 51, 34].

To efficiently complete the decomposition, we propose a fast operator-splitting scheme that leverages the Fast Fourier Transform (FFT) and explicit formulas for the involved sub-problems. For an image of size $N\times M$ , each iteration of the proposed algorithm has a computational complexity of $\mathcal{O}\left(MN(\log M+\log N)\right)$ . In comparison, the computational complexity in [29] is $\mathcal{O}((MN)^{3})$ . We present various examples to validate our method. Specifically, we conduct ablation studies to verify the combination of regularization terms in our method. We design experiments to investigate the impact of regularization parameters. We also demonstrate the effectiveness and efficiency of our approach by comparing it to state-of-the-art models for three-component-decomposition.

Our main contributions in this work are

•

We propose a new model for decomposing a grayscale image into its structural part, smooth part, and oscillatory part. To achieve satisfying results, we incorporate curvature regularization on level lines along with the sparsity-inducing $L^{0}$ -gradient. Our model effectively reduces the undesirable staircase effects and irregular boundary curves caused by the geometry-agnostic $L^{0}$ -gradient, yielding more visually appealing components.
•

We develop a new operator-splitting scheme to tackle the non-convex non-smooth optimization problem associated with our model. The proposed algorithm is computationally efficient. Each sub-problem is directly solved or numerically addressed using FFT techniques. Remarkably, when compared to several state-of-the-art methods, our algorithm achieves a speedup of over 10 times.
•

We present systematic ablation and comparison studies to justify the effectiveness and efficiency of our proposed method. We also numerically investigate the effects of model parameters and discuss the problem of parameter selection.

The paper is organized as follows. In Section 2, we propose our model. In particular, we introduce our model accompanied by mathematical remarks in Section 2.1; and reformulate it in Section 2.2 for the subsequent algorithmic developments. In Section 3, we propose our operator-splitting scheme to address the non-convex non-smooth functional minimization problem associated with the reformulated model. In Section 4, we describe the discretization and implementation details. In Section 5, we present various experiments to justify the proposed model and compare it with other methods from different aspects. We conclude the paper in Section 6.

2 Proposed Cartoon-Smooth-Texture Decomposition Model

In this section, we propose a new model for cartoon-smooth-texture image decomposition integrating Euler’s elastica. We present a formal version of our model in Section 2.1 where each component is characterized by a functional. Considering the difficulty of the problem, we introduce a reformulation of our model in Section 2.1 and derive a formal optimality condition, which paves our way towards developing an efficient operator splitting scheme in Section 3.

2.1 Proposed model

Suppose $\Omega\subset\mathbb{R}^{2}$ is a bounded open domain with a Lipschitz boundary. We assume that a grayscale image $f:\Omega\to\mathbb{R}$ is a sum of three components with distinctive features. In particular, we consider $f$ to be composed by a structural part $v^{*}:\Omega\to\mathbb{R}$ that defines objects with sharp boundaries or regions outlined by strong transitions between the light and the dark; a smooth part $w^{*}:\Omega\to\mathbb{R}$ that characterizes the soft intensity variation, such as shadows and uneven lighting; and an oscillatory part $n^{*}:\Omega\to\mathbb{R}$ that captures fine-scale textures and noise. Our objective is to extract these three components from $f$ , and we propose to achieve this by considering the following problem

\displaystyle\{v^{*},w^{*},n^{*}\}

\displaystyle\in\operatorname*{arg\,min}_{\begin{subarray}{c}v,w,n\\ v+w+n=f\end{subarray}}\left[\alpha_{0}\|\nabla v\|_{0}+\alpha_{\text{curv}}\int_{\Omega}\left(\nabla\cdot\frac{v}{|\nabla v|}\right)^{2}\left|\nabla v\right|\,d\mathbf{x}+\alpha_{w}\|\Delta w\|_{L^{2}(\Omega)}^{2}+\alpha_{n}\|n\|^{2}_{H^{-1}(\Omega)}\right].

(7)

In our proposed model (7), each component of the given image $f$ is characterized by different functionals. Particularly, for the structural part, we associate it with the $L^{0}$ -gradient energy and the total squared curvature of the level-lines; for the smooth part, we consider the squared $L^{2}$ -norm of the Laplacian; and for the oscillatory part, we use the squared inverse Sobolev norm. Each of these terms is weighted by a non-negative parameter, which controls various properties of the reconstructed components.

We design functionals in (7) for better visual effects and computational efficiency. The sparsity-inducing $L^{0}$ -gradient energy enforces the gradient of $v$ to be piecewise constant. Compared to the well-known TV regularization, it punishes discontinuities independent of their jump sizes, yielding results with higher contrast; this also implies stronger stair-casing effects. To reduce such artifacts in $v$ , our second term regularizes the TV-weighted squared curvatures of its level-lines. It combines high-order derivative information of the image, which has been shown to be an effective strategy for improving the visual effects of the reconstruction [13]. The squared $L^{2}$ -norm of the Laplacian encourages isotropic smoothness of $w$ , modeling the gradual changes of illumination over the scene. Its effect is similar to that of the squared $L^{2}$ -norm of the second-order derivatives, i.e., $\|\partial^{2}w\|_{L^{2}(\Omega)}^{2}$ used in (6). Under certain boundary conditions, we show that they are equivalent in the sense of having the same critical points (see Appendix B), which generalizes the observation in [13]. Using this term, we gain significant computational efficiency thanks to the simpler first variation (see Section 3). As for the oscillatory part, we use the squared inverse Sobolev norm

\|n\|_{H^{-1}(\Omega)}^{2}=\|\nabla\Delta^{-1}n\|_{L^{2}(\Omega)}^{2}

as a practical approximation for the Meyer’s $G$ -norm.

Remark 1.

We note that some regularizations in our model (7) are rather formal. For concreteness, it is natural to assume that $w\in H^{2}(\Omega)$ and $n\in H^{-1}(\Omega)$ ; thus there is no difficulty in understanding the regularizations on the smooth and oscillatory parts. Since a classical setting for piecewise constant images is the space of bounded variation functions $\text{BV}(\Omega)$ or de Giorgi and Ambrosio’s $\text{SBV}(\Omega)$ [16], both terms associated with the structural part call for interpretations. For the $L^{0}$ -gradient energy $\|\nabla v\|_{0}$ , it is often understood in a discrete sense [57], then the energy is deduced from a counting measure supported on grid points. An alternative is to consider a sufficiently smooth approximation such as the minimax concave penalty function [59, 29]. For the Euler’s elastica part, the major technical problem is that level-lines of BV functions generally lack the sufficient regularity to define curvature in the classical sense. One strategy is to restrict $v$ to $C^{2}(\Omega)$ [40]; another approach is to define the weak curvature for BV functions after mollification [49]; a different paradigm involves considering the solution of the $p$ -Laplacian equation as $p\to 1^{+}$ ; yet another method relies on the pairing theory [2] for distributions [46].

Remark 2.

No matter which mathematical definitions are used for interpretation, our model (7) is challenging to analyze for the existence of minimizers. The main difficulty is caused by the $L^{0}$ -gradient energy and the Euler’s elastica energy. If we use the TV-norm for regularizing $v$ instead, the resulting model becomes convex, and the existence property can be derived similarly to [13]. With an additional restriction on the mean value of $w$ , then the uniqueness can be established. If we replace the $L^{0}$ -gradient energy by the TV-norm, we recover the ordinary Euler’s elastica [39, 40, 49] for the $v$ component, and the related analysis on semi-continuity in the context of $\Gamma$ -convergence can be found in [8].

Remark 3.

In our proposed model (7), both of the regularization on the curvature of level-lines of $v^{*}$ and the regularization on the Laplacian of $w^{*}$ encourage smoothness. However, they emphasize different types of smoothing. The curvature regularization on $v^{*}$ admits a geometric prior that prefers straight level lines, whereas the regularization on $w^{*}$ has no preferences in diffusion directions, i.e., it is isotropic.

There exist functions $w$ defined on $\Omega\subset\mathbb{R}^{2}$ whose curvatures of level lines are constantly zero but not harmonic. For example, consider $w(x_{1},x_{2})=x_{1}^{2}$ defined over $(x_{1},x_{2})\in\Omega=[-1,1]\times[-1,1]$ . There also exist harmonic functions whose level lines have non-zero curvatures, such as $v(x_{1},x_{2})=x_{1}^{2}-x_{2}^{2}$ . For general images, the amount of smoothness in a given image can be transferred between $v^{*}$ and $w^{*}$ depending on the ratio $\alpha_{\text{curv}}/\alpha_{w}$ . We investigate this special property of our model in Section 5 using real images.

We emphasize that the $v^{*}$ component is also regularized by the $L^{0}$ -gradient, which is crucial in distinguishing the structural features in $v^{*}$ from the soft shading and blurry shapes represented by $w^{*}$ . Consider a Gaussian blurred smooth function $f_{\varepsilon}(x_{1},x_{2})=G_{\varepsilon}*h(x_{1},x_{2})$ defined on $[0,1]\times[0,1]$ where $h(x_{1},x_{2})=1$ if $0\leq x_{1}\leq 0.5,0\leq x_{2}\leq 1$ , and $h(x_{1},x_{2})=0$ otherwise, $G_{\varepsilon}$ is the Gaussian kernel with standard deviation $\varepsilon>0$ , and $*$ is the convolution. The level lines of $f_{\varepsilon}$ remain straight and $\|\Delta f\|_{L^{2}(\Omega)}^{2}=\mathcal{O}(\varepsilon^{-4})>0$ which decreases in $\varepsilon$ . If $v^{*}$ is only penalized by the curvature of its level line and we ignore the oscillatory component for simplicity, then the optimal decomposition is $v^{*}=f_{\varepsilon}$ and $w^{*}=0$ for any choice of $\varepsilon$ . This is contradictory to the visual perception of homogeneous structures with clear boundaries. To consider the $L^{0}$ -gradient, we first note that the Gaussian kernel has a compact support in the discrete sense. In particular, for a fixed threshold $\gamma>0$ , we note that

\frac{1}{\varepsilon\sqrt{2\pi}}\exp\left(-\frac{x_{1}^{2}}{2\varepsilon^{2}}\right)>\gamma\;

only if

-\sqrt{-2\varepsilon^{2}\ln(\gamma\varepsilon\sqrt{2\pi})}\leq x_{1}\leq\sqrt{-2\varepsilon^{2}\ln(\gamma\varepsilon\sqrt{2\pi})}\;.

Hence, $\|\nabla f_{\varepsilon}\|_{0}\approx 2\sqrt{-2\varepsilon^{2}\ln(\gamma\varepsilon\sqrt{2\pi})}$ , which increases when $\varepsilon$ is bigger, i.e., $f_{\varepsilon}$ is blurrier. With appropriate choices of the parameter $\alpha_{0}$ and $\alpha_{n}$ , the optimal decomposition will be $v^{*}=0$ and $w^{*}=f_{\varepsilon}$ if $\varepsilon$ is sufficiently large, and the optimal decomposition becomes $v^{*}=f_{\varepsilon}$ and $w^{*}=0$ if $\varepsilon$ is sufficiently small. By considering the $L^{0}$ -gradient of $v^{*}$ component, the decomposition is compatible with visual perception.

2.2 Model reformulation and formal optimality condition

As discussed above, Model (7) is non-trivial to analyze and solve. To arrive at a practical form that can be numerically addressed, we first write the constraint as a quadratic penalization term to convert (7) as follows

	$\displaystyle\min_{\begin{subarray}{c}v\in H^{1}(\Omega),w\in H^{2}(\Omega),n\in H^{-1}(\Omega)\end{subarray}}\bigg{[}$	$\displaystyle\alpha_{0}\\|\nabla v\\|_{0}+\alpha_{\text{curv}}\int_{\Omega}\left(\nabla\cdot\frac{v}{\|\nabla v\|}\right)^{2}\left\|\nabla v\right\|\,d\mathbf{x}+\alpha_{w}\\|\Delta w\\|_{L^{2}(\Omega)}^{2}+\alpha_{n}\\|n\\|^{2}_{H^{-1}(\Omega)}$
		$\displaystyle+\frac{1}{2}\\|f-(v+w+n)\\|_{L^{2}(\Omega)}^{2}\bigg{]}.$		(8)

We then resolve the high-order and non-linear term associated with the curvature by introducing auxiliary variables $\mathbf{q}=\nabla v$ and $\boldsymbol{\mu}=\frac{\nabla v}{|\nabla v|}$ , which satisfy the obvious constraints $\mathbf{q}\cdot\boldsymbol{\mu}=|\mathbf{q}|$ and $|\boldsymbol{\mu}|\leq 1$ as in [17]. Moreover, if we let $\mathbf{g}=(g_{1},g_{2})=\nabla\Delta^{-1}v\in L^{2}(\Omega)\times L^{2}(\Omega)$ as in [45], then (8) is equivalent to

\displaystyle\begin{cases}\displaystyle\min_{\begin{subarray}{c}v\in H^{1}(\Omega)\\ \mathbf{q}\in(L^{2}(\Omega))^{2},\boldsymbol{\mu}\in(H^{1}(\Omega))^{2},\\ w\in H^{2}(\Omega),\mathbf{g}\in(H^{1}(\Omega))^{2}\end{subarray}}\bigg{[}\alpha_{0}\|\mathbf{q}\|_{0}+\alpha_{\text{curv}}\int_{\Omega}\left(\nabla\cdot\boldsymbol{\mu}\right)^{2}\left|\mathbf{q}\right|\,d\mathbf{x}+\alpha_{w}\|\Delta w\|_{L^{2}(\Omega)}^{2}+\alpha_{n}\|\mathbf{g}\|^{2}_{L^{2}(\Omega)}\\ \hskip 56.9055pt+\frac{1}{2}\|f-(v+w+\nabla\cdot\mathbf{g})\|_{L^{2}(\Omega)}^{2}\bigg{]},\\ \mathbf{q}=\nabla v,\\ \boldsymbol{\mu}=\frac{\nabla v}{|\nabla v|},\\ |\boldsymbol{\mu}|\leq 1.\end{cases}

(9)

We note that the problem can be further simplified by observing the following characterization of minimizers of Model (8), whose proof is in Appendix A.

Proposition 2.1.

Let $(v^{*},w^{*},n^{*})$ be a minimizer of (8), then we have

\displaystyle\int_{\Omega}v^{*}+w^{*}+n^{*}\,d\mathbf{x}=\int_{\Omega}f\,d\mathbf{x}.

(10)

This implies that if $(v^{*},\mathbf{q}^{*},w^{*},\mathbf{g}^{*})$ is a minimizer of Model (9), $v^{*}$ is in fact uniquely determined by the solution of the following problem

\displaystyle\begin{cases}\nabla^{2}v=\nabla\cdot\mathbf{q}^{*}&\mbox{ in }\Omega,\\ (\nabla v-\mathbf{q}^{*})\cdot\mathbf{n}=0&\mbox{ on }\partial\Omega,\\ \int_{\Omega}v\,d\mathbf{x}=\int_{\Omega}f-w^{*}-\nabla\cdot\mathbf{g}^{*}\,d\mathbf{x}.\end{cases}

(11)

where $\mathbf{n}$ denotes the outward normal direction along $\partial\Omega$ . Therefore, we consider an approximation of (9) as

	$\displaystyle\min_{\begin{subarray}{c}\mathbf{q}\in(L^{2}(\Omega))^{2},\boldsymbol{\mu}\in(H^{1}(\Omega))^{2},\\ w\in H^{2}(\Omega),\mathbf{g}\in(H^{1}(\Omega))^{2}\end{subarray}}\bigg{[}$	$\displaystyle\alpha_{0}\\|\mathbf{q}\\|_{0}+\alpha_{\rm curv}\int_{\Omega}\left\|\nabla\cdot\boldsymbol{\mu}\right\|^{2}\left\|\mathbf{q}\right\|\,d\mathbf{x}+\alpha_{w}\int_{\Omega}\|\nabla^{2}w\|^{2}d\mathbf{x}$
		$\displaystyle+\alpha_{n}\int_{\Omega}\|\mathbf{g}\|^{2}d\mathbf{x}+\frac{1}{2}\int_{\Omega}\|f-(v_{\mathbf{q},w,\mathbf{g}}+w+\nabla\cdot\mathbf{g})\|^{2}d\mathbf{x}+I_{\Sigma_{f}}(\mathbf{q},w,\mathbf{g})+I_{S}(\mathbf{q},\boldsymbol{\mu})\bigg{]},$		(12)

where $v_{\mathbf{q},w,\mathbf{g}}$ is uniquely determined by $\mathbf{q},w,\mathbf{g}$ via (11); and the indicator functions

\displaystyle I_{\Sigma_{f}}(\mathbf{q},w,\mathbf{g})=\begin{cases}0&\mbox{ if }(\mathbf{q},w,\mathbf{g})\in\Sigma_{f},\\ +\infty&\mbox{ otherwise},\end{cases}\quad I_{S}(\mathbf{q},\boldsymbol{\mu})=\begin{cases}0&\mbox{ if }(\mathbf{q},\boldsymbol{\mu})\in S,\\ +\infty&\mbox{ otherwise}.\end{cases}

with

	$\displaystyle\Sigma_{f}=\bigg{\{}(\mathbf{q},w,\mathbf{g})\in(L^{2}(\Omega))^{2}\times L^{2}(\Omega)\times(H^{1}(\Omega))^{2}\leavevmode\nobreak\ \|\leavevmode\nobreak\$
	$\displaystyle\hskip 85.35826pt\exists v\in H^{1}(\Omega)\leavevmode\nobreak\ \text{such that}\leavevmode\nobreak\ \mathbf{q}=\nabla v\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \int_{\Omega}(v+w+\nabla\cdot\mathbf{g}-f)\,d\mathbf{x}=0\bigg{\}},$
	$\displaystyle S=\left\{(\mathbf{q},\boldsymbol{\mu})\in(L^{2}(\Omega))^{2}\times(L^{2}(\Omega))^{2}\leavevmode\nobreak\ \|\leavevmode\nobreak\ \mathbf{q}\cdot\boldsymbol{\mu}=\|\mathbf{q}\|,\|\boldsymbol{\mu}\|\leq 1\right\}.$

Suppose we find $\{\mathbf{q}^{*},\boldsymbol{\mu}^{*},w^{*},\mathbf{g}^{*}\}$ as a minimizer of (12), then we have $v^{*}=v_{\mathbf{q}^{*},w^{*},\mathbf{g}^{*}}$ for the structural part, $w^{*}$ for the smooth part, and $n^{*}=\nabla\cdot\mathbf{g}^{*}$ for the oscillatory part. Moreover, if we denote

	$\displaystyle J_{0}(\mathbf{q},w,\mathbf{g})=\frac{1}{2}\int_{\Omega}\|f-(v_{\mathbf{q},w,\mathbf{g}}+w+\nabla\cdot\mathbf{g})\|^{2}d\mathbf{x},$
	$\displaystyle J_{1}(\mathbf{q})=\alpha_{0}\\|\mathbf{q}\\|_{0},$
	$\displaystyle J_{2}(\mathbf{q},\boldsymbol{\mu})=\alpha_{\rm curv}\int_{\Omega}\left\|\nabla\cdot\boldsymbol{\mu}\right\|^{2}\left\|\mathbf{q}\right\|\,d\mathbf{x},$
	$\displaystyle J_{3}(w)=\alpha_{w}\int_{\Omega}\|\nabla^{2}w\|^{2}d\mathbf{x},$
	$\displaystyle J_{4}(\mathbf{g})=\alpha_{n}\int_{\Omega}\|\mathbf{g}\|^{2}d\mathbf{x},$

the minimizer $\{\mathbf{q}^{*},\boldsymbol{\mu}^{*},w^{*},\mathbf{g}^{*}\}$ satisfies the formal optimality condition

\displaystyle\begin{cases}\partial_{\mathbf{q}}J_{1}(\mathbf{q}^{*})+\partial_{\mathbf{q}}J_{2}(\mathbf{q}^{*},\boldsymbol{\mu}^{*})+\partial_{\mathbf{q}}J_{0}(\mathbf{q}^{*},w^{*},\mathbf{g}^{*})+\partial_{\mathbf{q}}I_{\Sigma_{f}}(\mathbf{q}^{*},w^{*},\mathbf{g}^{*})+\partial_{\mathbf{q}}I_{S}(\mathbf{q}^{*},\boldsymbol{\mu}^{*})\ni 0,\\ \partial_{\boldsymbol{\mu}}J_{2}(\mathbf{q}^{*},\boldsymbol{\mu}^{*})+\partial_{\boldsymbol{\mu}}I_{S}(\mathbf{q}^{*},\boldsymbol{\mu}^{*})\ni 0,\\ \partial_{w}J_{3}(w^{*})+\partial_{w}J_{0}(\mathbf{q}^{*},w^{*},\mathbf{g}^{*})+\partial_{w}I_{\Sigma_{f}}(\mathbf{q}^{*},w^{*},\mathbf{g}^{*})\ni 0,\\ \partial_{\mathbf{g}}J_{4}(\mathbf{g}^{*})+\partial_{\mathbf{g}}J_{0}(\mathbf{q}^{*},w^{*},\mathbf{g}^{*})+\partial_{\mathbf{g}}I_{\Sigma_{f}}(\mathbf{q}^{*},w^{*},\mathbf{g}^{*})\ni 0.\end{cases}

(13)

where $\partial$ is the subdifferential.

3 Operator-Splitting Algorithm for the Proposed Model

In this section, we propose an efficient algorithm for Model (12) based on an operator-splitting framework. In Section 2.2, we associate the model with a gradient flow based on its formal optimality condition. In Section 3.1, we present a novel operator-splitting scheme to solve the reformulated model. This non-convex non-smooth minimization problem is decomposed into several subproblems, whose solutions are described in Section 3.2-3.4. In Section 3.5, we adapt the proposed algorithm to periodic boundary condition. We summarize the proposed scheme in Algorithm 1.

3.1 Proposed scheme

We associate (13) with the following initial value problem in analogous to the gradient flow

\displaystyle\begin{cases}\frac{\partial\mathbf{p}}{\partial t}+\partial_{\mathbf{q}}J_{1}(\mathbf{p})+\partial_{\mathbf{q}}J_{2}(\mathbf{p},\boldsymbol{\lambda})+\partial_{\mathbf{q}}J_{0}(\mathbf{p},r,\mathbf{s})+\partial_{\mathbf{q}}I_{\Sigma_{f}}(\mathbf{p})+\partial_{\mathbf{q}}I_{S}(\mathbf{p},\boldsymbol{\lambda})\ni 0,\\ \gamma_{1}\frac{\partial\boldsymbol{\lambda}}{\partial t}+\partial_{\boldsymbol{\mu}}J_{2}(\mathbf{p},\boldsymbol{\lambda})+\partial_{\boldsymbol{\mu}}I_{S}(\mathbf{p},\boldsymbol{\lambda})\ni 0,\\ \gamma_{2}\frac{\partial r}{\partial t}+\partial_{w}J_{3}(r)+\partial_{w}J_{0}(\mathbf{p},r,\mathbf{s})=0,\\ \gamma_{3}\frac{\partial\mathbf{s}}{\partial t}+\partial_{\mathbf{g}}J_{4}(\mathbf{s})+\partial_{\mathbf{g}}J_{0}(\mathbf{p},w,\mathbf{s})=0,\\ \{\mathbf{p}(0),\boldsymbol{\lambda}(0),r(0),\mathbf{s}(0)\}=\{\mathbf{p}_{0},\boldsymbol{\lambda}_{0},r_{0},\mathbf{s}_{0}\},\end{cases}

(14)

where $\{\mathbf{p}_{0},\boldsymbol{\lambda}_{0},r_{0},\mathbf{s}_{0}\}$ is a given initial condition, and $\gamma_{1},\gamma_{2},\gamma_{3}>0$ are parameters that control the evolution speed of $\mathbf{p},\boldsymbol{\lambda},r,\mathbf{s}$ , respectively. We note that the steady state solution of (14) satisfies the optimality condition (13) for Model (12). To avoid confusion, we distinguish the time-dependent variables in (14) from the optimization variables in (12), and Table 1 shows the correspondences.

	Correspondence
General variables	$\mathbf{q}$	$\boldsymbol{\mu}$	$w$	$\mathbf{g}$
PDE solutions	$\mathbf{p}$	$\boldsymbol{\lambda}$	$r$	$\mathbf{s}$

Table 1: Notation of general variables and PDE solutions and their correspondences.

We find the steady state solution of (14) by employing the operator-splitting framework [26, 25]. Our proposed scheme proceeds as follows:

Initialization

\displaystyle\text{Initialize}\leavevmode\nobreak\ \{\mathbf{p}^{0},\boldsymbol{\lambda}^{0},r^{0},\mathbf{s}^{0}\}=\{\mathbf{p}_{0},\boldsymbol{\lambda}_{0},r_{0},\mathbf{s}_{0}\}.

Fix some $\tau>0$ as the time step. For $k=0,1,2,\dots$ , denote $t^{k}=k\tau$ for $k=0,1,2\dots$ ; and we update $\{\mathbf{p}^{k},\boldsymbol{\lambda}^{k},r^{k},\mathbf{s}^{k}\}\rightarrow\{\mathbf{p}^{k+1},\boldsymbol{\lambda}^{k+1},r^{k+1},\mathbf{s}^{k+1}\}$ via four fractional steps:

Fractional step 1: Solve

\displaystyle\begin{cases}\begin{cases}\frac{\partial\mathbf{p}}{\partial t}+\partial_{\mathbf{q}}J_{1}(\mathbf{p})\ni 0,\\ \gamma_{1}\frac{\partial\boldsymbol{\lambda}}{\partial t}=0,\\ \gamma_{2}\frac{\partial r}{\partial t}=0,\\ \gamma_{3}\frac{\partial\mathbf{s}}{\partial t}=0,\end{cases}\mbox{ in }\Omega\times(t^{k},t^{k+1}),\\ \{\mathbf{p}(t^{k}),\boldsymbol{\lambda}(t^{k}),r(t^{k}),\mathbf{s}(t^{k})\}=\{\mathbf{p}^{k},\boldsymbol{\lambda}^{k},r^{k},\mathbf{s}^{k}\},\end{cases}

and set

\displaystyle\{\mathbf{p}^{k+1/4},\boldsymbol{\lambda}^{k+1/4},r^{k+1/4},\mathbf{s}^{k+1/4}\}=\{\mathbf{p}(t^{k+1}),\boldsymbol{\lambda}(t^{k+1}),r(t^{k+1}),\mathbf{s}(t^{k+1})\}.

Fractional step 2: Solve

\displaystyle\begin{cases}\begin{cases}\frac{\partial\mathbf{p}}{\partial t}+\partial_{\mathbf{q}}J_{2}(\mathbf{p},\boldsymbol{\lambda})=0,\\ \gamma_{1}\frac{\partial\boldsymbol{\lambda}}{\partial t}+\partial_{\boldsymbol{\mu}}J_{2}(\mathbf{p},\boldsymbol{\lambda})\ni 0,\\ \gamma_{2}\frac{\partial r}{\partial t}=0,\\ \gamma_{3}\frac{\partial\mathbf{s}}{\partial t}=0,\end{cases}\mbox{ in }\Omega\times(t^{k},t^{k+1}),\\ \{\mathbf{p}(t^{k}),\boldsymbol{\lambda}(t^{k}),r(t^{k}),\mathbf{s}(t^{k})\}=\{\mathbf{p}^{k+1/4},\boldsymbol{\lambda}^{k+1/4},r^{k+1/4},\mathbf{s}^{k+1/4}\},\end{cases}

and set

\displaystyle\{\mathbf{p}^{k+2/4},\boldsymbol{\lambda}^{k+2/4},r^{k+2/4},\mathbf{s}^{k+2/4}\}=\{\mathbf{p}(t^{k+1}),\boldsymbol{\lambda}(t^{k+1}),r(t^{k+1}),\mathbf{s}(t^{k+1})\}.

Fractional step 3: Solve

\displaystyle\begin{cases}\begin{cases}\frac{\partial\mathbf{p}}{\partial t}+\partial_{\mathbf{q}}I_{S}(\mathbf{p},\boldsymbol{\lambda})\ni 0,\\ \gamma_{1}\frac{\partial\boldsymbol{\lambda}}{\partial t}+\partial_{\boldsymbol{\mu}}I_{S}(\mathbf{p},\boldsymbol{\lambda})\ni 0,\\ \gamma_{2}\frac{\partial r}{\partial t}=0,\\ \gamma_{3}\frac{\partial\mathbf{s}}{\partial t}=0,\end{cases}\mbox{ in }\Omega\times(t^{k},t^{k+1}),\\ \{\mathbf{p}(t^{k}),\boldsymbol{\lambda}(t^{k}),r(t^{k}),\mathbf{s}(t^{k})\}=\{\mathbf{p}^{k+2/4},\boldsymbol{\lambda}^{k+2/4},r^{k+2/4},\mathbf{s}^{k+2/4}\},\end{cases}

and set

\displaystyle\{\mathbf{p}^{k+3/4},\boldsymbol{\lambda}^{k+3/4},r^{k+3/4},\mathbf{s}^{k+3/4}\}=\{\mathbf{p}(t^{k+1}),\boldsymbol{\lambda}(t^{k+1}),r(t^{k+1}),\mathbf{s}(t^{k+1})\}.

Fractional step 4: Solve

\displaystyle\begin{cases}\begin{cases}\frac{\partial\mathbf{p}}{\partial t}+\partial_{\mathbf{q}}J_{0}(\mathbf{p},r,\mathbf{s})+\partial_{\mathbf{q}}I_{\Sigma_{f}}(\mathbf{p},r,\mathbf{s})\ni 0,\\ \gamma_{1}\frac{\partial\boldsymbol{\lambda}}{\partial t}=0,\\ \gamma_{2}\frac{\partial r}{\partial t}+\partial_{w}J_{0}(\mathbf{p},r,\mathbf{s})+\partial_{w}J_{3}(r)+\partial_{w}I_{\Sigma_{f}}(\mathbf{p},r,\mathbf{s})\ni 0,\\ \gamma_{3}\frac{\partial\mathbf{s}}{\partial t}+\partial_{\mathbf{g}}J_{0}(\mathbf{p},r,\mathbf{s})+\partial_{\mathbf{g}}J_{4}(\mathbf{s})+\partial_{\mathbf{g}}I_{\Sigma_{f}}(\mathbf{p},r,\mathbf{s})\ni 0,\end{cases}\mbox{ in }\Omega\times(t^{k},t^{k+1}),\\ \{\mathbf{p}(t^{k}),\boldsymbol{\lambda}(t^{k}),r(t^{k}),\mathbf{s}(t^{k})\}=\{\mathbf{p}^{k+3/4},\boldsymbol{\lambda}^{k+3/4},r^{k+3/4},\mathbf{s}^{k+3/4}\},\end{cases}

and set

\displaystyle\{\mathbf{p}^{k+1},\boldsymbol{\lambda}^{k+1},r^{k+1},\mathbf{s}^{k+1}\}=\{\mathbf{p}(t^{k+1}),\boldsymbol{\lambda}(t^{k+1}),r(t^{k+1}),\mathbf{s}(t^{k+1})\}.

We time discretize the systems above by the Marchuk-Yanenko scheme and get

	$\displaystyle\begin{cases}\frac{\mathbf{p}^{k+1/4}-\mathbf{p}^{k}}{\tau}+\partial_{\mathbf{q}}J_{1}(\mathbf{p}^{k+1/4})\ni 0,\\ \gamma_{1}\frac{\boldsymbol{\lambda}^{k+1/4}-\boldsymbol{\lambda}^{k}}{\tau}=0,\\ \gamma_{2}\frac{r^{k+1/4}-r^{k}}{\tau}=0,\\ \gamma_{3}\frac{\mathbf{s}^{k+1/4}-\mathbf{s}^{k}}{\tau}=0,\end{cases}$		(15)
	$\displaystyle\begin{cases}\frac{\mathbf{p}^{k+2/4}-\mathbf{p}^{k+1/4}}{\tau}+\partial_{\mathbf{q}}J_{2}(\mathbf{p}^{k+2/4},\boldsymbol{\lambda}^{k+1/4})\ni 0,\\ \gamma_{1}\frac{\boldsymbol{\lambda}^{k+2/4}-\boldsymbol{\lambda}^{k+1/4}}{\tau}+\partial_{\boldsymbol{\mu}}J_{2}(\mathbf{p}^{k+2/4},\boldsymbol{\lambda}^{k+2/4})\ni 0,\\ \gamma_{2}\frac{r^{k+2/4}-r^{k+1/4}}{\tau}=0,\\ \gamma_{3}\frac{\mathbf{s}^{k+2/4}-\mathbf{s}^{k+1/4}}{\tau}=0,\end{cases}$		(16)
	$\displaystyle\begin{cases}\frac{\mathbf{p}^{k+3/4}-\mathbf{p}^{k+2/4}}{\tau}+\partial_{\mathbf{q}}I_{S}(\mathbf{p}^{k+3/4},\boldsymbol{\lambda}^{k+3/4})\ni 0,\\ \gamma_{1}\frac{\boldsymbol{\lambda}^{k+3/4}-\boldsymbol{\lambda}^{k+2/4}}{\tau}+\partial_{\boldsymbol{\mu}}I_{S}(\mathbf{p}^{k+3/4},\boldsymbol{\lambda}^{k+3/4})\ni 0,\\ \gamma_{2}\frac{r^{k+3/4}-r^{k+2/4}}{\tau}=0,\\ \gamma_{3}\frac{\mathbf{s}^{k+3/4}-\mathbf{s}^{k+2/4}}{\tau}=0,\end{cases}$		(17)
	$\displaystyle\begin{cases}\frac{\mathbf{p}^{k+1}-\mathbf{p}^{k+3/4}}{\tau}+\partial_{\mathbf{q}}J_{0}(\mathbf{p}^{k+1},r^{k+1},\mathbf{s}^{k+1})+\partial_{\mathbf{q}}I_{\Sigma_{f}}(\mathbf{p}^{k+1},r^{k+1},\mathbf{s}^{k+1})\ni 0,\\ \gamma_{1}\frac{\boldsymbol{\lambda}^{k+1}-\boldsymbol{\lambda}^{k+3/4}}{\tau}=0,\\ \gamma_{2}\frac{r^{k+1}-r^{k+3/4}}{\tau}+\partial_{w}J_{0}(\mathbf{p}^{k+1},r^{k+1},\mathbf{s}^{k+1})+\partial_{w}J_{3}(r^{k+1})+\partial_{w}I_{\Sigma_{f}}(\mathbf{p}^{k+1},r^{k+1},\mathbf{s}^{k+1})\ni 0,\\ \gamma_{3}\frac{\mathbf{s}^{k+1}-\mathbf{s}^{k+3/4}}{\tau}+\partial_{\mathbf{g}}J_{0}(\mathbf{p}^{k+1},r^{k+1},\mathbf{s}^{k+1})+\partial_{\mathbf{g}}J_{4}(\mathbf{s}^{k+1})+\partial_{\mathbf{g}}I_{\Sigma_{f}}(\mathbf{p}^{k+1},r^{k+1},\mathbf{s}^{k+1})\ni 0.\end{cases}$		(18)

for $k=0,1,2,\dots$ . We note that each system above corresponds to a subproblem which can be solved explicitly or approximated efficiently. In the following sections, we describe the details.

3.2 On the solution to sub-problem (15) and (16)

Given $\mathbf{p}^{k}$ for $k=0,1,2,\dots$ , the vector field $\mathbf{p}^{k+1/4}$ solves

\displaystyle\mathbf{p}^{k+1/4}=\operatorname*{arg\,min}_{\mathbf{q}\in(L^{2}(\Omega))^{2}}\left[\frac{1}{2}\int_{\Omega}|\mathbf{q}-\mathbf{p}^{k}|^{2}d\mathbf{x}+\tau\alpha_{0}\|\mathbf{q}\|_{0}\right]\;.

This problem has the closed-form solution [57]

\displaystyle\mathbf{p}^{k+1/4}=\begin{cases}\mathbf{0}&\mbox{ if }|\mathbf{p}^{k}|^{2}\leq\tau\alpha_{0},\\ \mathbf{p}^{k}&\mbox{ otherwise}.\end{cases}

(19)

The first equation in (16) is the Euler-Lagrange equation of the following minimization problem

\displaystyle\mathbf{p}^{k+2/4}=\operatorname*{arg\,min}_{\mathbf{q}\in(L^{2}(\Omega))^{2}}\left[\frac{1}{2}\int_{\Omega}|\mathbf{q}-\mathbf{p}^{k+1/4}|^{2}d\mathbf{x}+\tau\alpha_{\rm curv}\int_{\Omega}\left|\nabla\cdot\boldsymbol{\lambda}^{k+1/4}\right|^{2}|\mathbf{q}|d\mathbf{x}\right]\;,

whose closed-form solution is given as

\displaystyle\mathbf{p}^{k+2/4}=\max\left\{0,1-\frac{\tau\alpha_{\rm curv}\left|\nabla\cdot\boldsymbol{\lambda}^{k+1/4}\right|^{2}}{|\mathbf{p}^{k+1/4}|}\right\}\mathbf{p}^{k+1/4}.

(20)

For the second equation of (16), we note that $\boldsymbol{\lambda}^{k+2/4}$ solves

\displaystyle\boldsymbol{\lambda}^{k+2/4}=\operatorname*{arg\,min}_{\boldsymbol{\mu}\in(H^{1}(\Omega))^{2}}\left[\frac{\gamma_{1}}{2}\int_{\Omega}|\boldsymbol{\mu}-\boldsymbol{\lambda}^{k+1/4}|^{2}d\mathbf{x}+\tau\alpha_{\rm curv}\int_{\Omega}\left|\nabla\cdot\boldsymbol{\mu}\right|^{2}|\mathbf{p}^{k+2/4}|d\mathbf{x}\right],

(21)

and the corresponding Euler-Lagrange equation is

\displaystyle\begin{cases}\gamma_{1}\frac{\boldsymbol{\lambda}^{k+2/4}-\boldsymbol{\lambda}^{k+1/4}}{\tau}-2\alpha_{\rm curv}\nabla(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4})=0\mbox{ in }\Omega,\\ |\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4}=0\mbox{ on }\partial\Omega.\end{cases}

(22)

3.3 On the solution to sub-problem (17)

We update $\mathbf{p}^{k+3/4}$ and $\boldsymbol{\lambda}^{k+3/4}$ using the method proposed in [17, Section 3.5]. Note that $\mathbf{p}^{k+3/4}$ and $\boldsymbol{\lambda}^{k+3/4}$ solve

\displaystyle(\mathbf{p}^{k+3/4},\boldsymbol{\lambda}^{k+3/4})=\operatorname*{arg\,min}_{(\mathbf{q},\boldsymbol{\mu})\in S}\left[\frac{1}{2}\int_{\Omega}|\mathbf{q}-\mathbf{p}^{k+2/4}|^{2}d\mathbf{x}+\frac{\gamma_{1}}{2}\int_{\Omega}|\boldsymbol{\mu}-\boldsymbol{\lambda}^{k+2/4}|^{2}d\mathbf{x}\right].

(23)

Problem (23) can be solved in a pointwise manner. For each given $\mathbf{x}$ , $(\mathbf{p}^{k+3/4}(\mathbf{x}),\boldsymbol{\lambda}^{k+3/4}(\mathbf{x}))$ satisfies

\displaystyle(\mathbf{p}^{k+3/4}(\mathbf{x}),\boldsymbol{\lambda}^{k+3/4}(\mathbf{x}))=\operatorname*{arg\,min}_{(\mathbf{q}(\mathbf{x}),\boldsymbol{\mu}(\mathbf{x}))\in\chi}\left[|\mathbf{q}(\mathbf{x})-\mathbf{q}^{k+2/4}(\mathbf{x})|^{2}+\gamma_{1}|\boldsymbol{\mu}(\mathbf{x})-\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})|^{2}\right]

(24)

with $\chi=\{(\mathbf{a},\mathbf{b})\in\mathbb{R}^{2}:\mathbf{a}\cdot\mathbf{b}=|\mathbf{a}|,|\mathbf{b}|\leq 1\}$ . Define

\displaystyle\chi_{1}=\{(\mathbf{a},\mathbf{b})\in\mathbb{R}^{2}:|\mathbf{a}|=0,|\mathbf{b}|\leq 1\},\;\text{and}\;\chi_{2}=\{(\mathbf{a},\mathbf{b})\in\mathbb{R}^{2}:|\mathbf{a}|\neq 0,\mathbf{a}\cdot\mathbf{b}=|\mathbf{q}|,|\mathbf{b}|=1\}.

We have $\chi=\chi_{1}\cup\chi_{2}$ . We next discuss the solution of (24) on $\chi_{1}$ and $\chi_{2}$ separately.

Over $\chi_{1}$ , denote the solution to (24) by $(\widehat{\mathbf{p}}^{k+3/4}(\mathbf{x}),\widehat{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x}))$ . Problem (24) reduces to

\displaystyle\widehat{\boldsymbol{\lambda}}(\mathbf{x})^{k+3/4}=\operatorname*{arg\,min}_{\boldsymbol{\mu}(\mathbf{x})\in\mathbb{R}^{2}:|\boldsymbol{\mu}(\mathbf{x})|\leq 1}|\boldsymbol{\mu}(\mathbf{x})-\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})|^{2}\;,

and $\widehat{\mathbf{p}}^{k+3/4}(\mathbf{x})=\mathbf{0}$ . The closed-form solution for $\widehat{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x})$ is given as

\displaystyle\widehat{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x})=\frac{\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})}{\max\{1,|\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})|\}}.

(25)

Over $\chi_{2}$ , denote the solution to (24) by $(\widetilde{\mathbf{p}}^{k+3/4}(\mathbf{x}),\widetilde{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x}))$ , and it can be written as

\displaystyle(\widetilde{\mathbf{p}}^{k+3/4}(\mathbf{x}),\widetilde{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x}))=\operatorname*{arg\,min}_{\begin{subarray}{c}(\mathbf{q}(\mathbf{x}),\boldsymbol{\mu}(\mathbf{x}))\in\mathbb{R}^{2}\times\mathbb{R}^{2},\\ \mathbf{q}(\mathbf{x})\neq\mathbf{0},\mathbf{q}(\mathbf{x})\cdot\boldsymbol{\mu}(\mathbf{x})=|\mathbf{q}(\mathbf{x})|,\\ |\boldsymbol{\mu}(\mathbf{x})|=1\end{subarray}}\left[|\mathbf{q}(\mathbf{x})-\mathbf{p}^{k+2/4}(\mathbf{x})|^{2}+\gamma_{1}|\boldsymbol{\mu}(\mathbf{x})-\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})|^{2}\right].

(26)

Denote $\theta=|\mathbf{q}(\mathbf{x})|$ , $\mathbf{y}(\mathbf{x})=\mathbf{p}^{k+2/4}(\mathbf{x})$ , and $\mathbf{w}(\mathbf{x})=\boldsymbol{\mu}^{k+2/4}(\mathbf{x})$ for each $\mathbf{x}\in\Omega$ . Problem (26) is equivalent to solving

\displaystyle\min_{(\theta,\mathbf{q}(\mathbf{x}))\in\mathbb{R}\times\mathbb{R}^{2},\theta>0,|\mathbf{q}(\mathbf{x})|=1}\left[\theta^{2}-2\theta\mathbf{q}(\mathbf{x})\cdot\mathbf{y}(\mathbf{x})-2\gamma_{1}\mathbf{q}(\mathbf{x})\cdot\mathbf{w}(\mathbf{x})\right]\;.

(27)

For a fixed $\theta$ , the solution $\mathbf{q}^{*}(\mathbf{x},\theta)$ of (27) is given by

\displaystyle\mathbf{q}^{*}(\mathbf{x},\theta)=\frac{\theta\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x})}{|\theta\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x})|}\;.

(28)

Substituting (28) to (27) leads to

\displaystyle\min_{\theta>0}\left[\theta^{2}-2|\theta\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x})|\right].

(29)

We solve (29) by a fixed point method. Denote

\displaystyle E(\mathbf{x},\theta)=\theta^{2}-2|\theta\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x})|\;,

and thus

\displaystyle\frac{dE}{d\theta}(\mathbf{x},\theta)=2\theta-\frac{2\mathbf{y}(\mathbf{x})\cdot(\theta\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x}))}{|\theta\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x})|}.

We initialize $\theta^{0}=|\mathbf{y}(\mathbf{x})|$ . For $m=0,1,2,\dots$ , we update $\theta^{m+1}$ from $\theta^{m}$ via

\displaystyle\theta^{m+1}=\max\left\{0,\frac{\mathbf{y}(\mathbf{x})\cdot(\theta^{m}\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x}))}{|\theta^{m}\mathbf{y}(\mathbf{x})+\gamma_{1}\mathbf{w}(\mathbf{x})|}\right\}

(30)

until $|\theta^{m+1}-\theta^{m}|<\varepsilon$ for some small $\varepsilon>0$ .

Denote the converged $\theta$ by $\theta^{*}(\mathbf{x})$ , then we compute

\displaystyle\widetilde{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x})=\frac{\theta^{*}(\mathbf{x})\mathbf{p}^{k+2/4}(\mathbf{x})+\gamma_{1}\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})}{|\theta^{*}(\mathbf{x})\mathbf{p}^{k+2/4}(\mathbf{x})+\gamma_{1}\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})|},\quad\widetilde{\mathbf{p}}^{k+3/4}(\mathbf{x})=\theta^{*}\widetilde{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x})\;.

(31)

Finally, we define

\displaystyle G(\mathbf{q},\boldsymbol{\mu};\mathbf{x})=|\mathbf{q}(\mathbf{x})-\mathbf{q}^{k+2/4}(\mathbf{x})|^{2}+\gamma_{1}|\boldsymbol{\mu}(\mathbf{x})-\boldsymbol{\lambda}^{k+2/4}(\mathbf{x})|^{2},

and compare the minimizers in $\chi_{1}$ and $\chi_{2}$ to obtain $(\mathbf{p}^{k+3/4}(\mathbf{x}),\boldsymbol{\lambda}^{k+3/4}(\mathbf{x}))$ as

\displaystyle(\mathbf{p}^{k+3/4}(\mathbf{x}),\boldsymbol{\lambda}^{k+3/4}(\mathbf{x}))=\begin{cases}(\widehat{\mathbf{p}}^{k+3/4}(\mathbf{x}),\widehat{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x}))&\mbox{ if }G(\widehat{\mathbf{p}}^{k+3/4},\widehat{\boldsymbol{\lambda}}^{k+3/4};\mathbf{x})\\ &\hskip 28.45274pt\leq G(\widetilde{\mathbf{p}}^{k+3/4},\widetilde{\boldsymbol{\lambda}}^{k+3/4};\mathbf{x}),\\ (\widetilde{\mathbf{p}}^{k+3/4}(\mathbf{x}),\widetilde{\boldsymbol{\lambda}}^{k+3/4}(\mathbf{x}))&\mbox{ otherwise}.\end{cases}

(32)

3.4 On the solution to sub-problem (18)

In (18), $(\mathbf{p}^{k+1},r^{k+1},\mathbf{s}^{k+1})$ solves

	$\displaystyle\min_{\begin{subarray}{c}(\mathbf{q},w,\mathbf{g})\in\Sigma_{f},w\in H^{2}(\Omega),\\ \mathbf{g}\in(H^{1}(\Omega))^{2}\end{subarray}}$	$\displaystyle\bigg{[}\frac{1}{2}\int_{\Omega}\|\mathbf{q}-\mathbf{p}^{k+3/4}\|^{2}d\mathbf{x}+\frac{\gamma_{2}}{2}\int_{\Omega}\|w-r^{k+3/4}\|^{2}d\mathbf{x}+\frac{\gamma_{3}}{2}\int_{\Omega}\|\mathbf{g}-\mathbf{s}^{k+3/4}\|^{2}d\mathbf{x}$
		$\displaystyle+\tau\alpha_{w}\int_{\Omega}\|\nabla^{2}w\|^{2}d\mathbf{x}+\tau\alpha_{n}\int_{\Omega}\|\mathbf{g}\|^{2}d\mathbf{x}+\frac{\tau}{2}\int_{\Omega}\|f-(v_{\mathbf{q},w,\mathbf{g}}+w+\nabla\cdot\mathbf{g})\|^{2}d\mathbf{x}\bigg{]},$		(33)

or equivalently

\displaystyle\begin{cases}(v^{k+1},r^{k+1},\mathbf{s}^{k+1})=\displaystyle\operatorname*{arg\,min}_{\begin{subarray}{c}v\in H^{1}(\Omega),w\in H^{2}(\Omega),\\ \mathbf{g}\in(H^{1}(\Omega))^{2}\end{subarray}}\bigg{[}\frac{1}{2}\int_{\Omega}|\nabla v-\mathbf{p}^{k+3/4}|^{2}d\mathbf{x}+\frac{\gamma_{2}}{2}\int_{\Omega}|w-r^{k+3/4}|^{2}d\mathbf{x}\\ \hskip 113.81102pt+\frac{\gamma_{3}}{2}\displaystyle\int_{\Omega}|\mathbf{g}-\mathbf{s}^{k+3/4}|^{2}d\mathbf{x}+\tau\alpha_{w}\int_{\Omega}|\nabla^{2}w|^{2}d\mathbf{x}+\tau\alpha_{n}\int_{\Omega}|\mathbf{g}|^{2}d\mathbf{x}\\ \hskip 113.81102pt+\frac{\tau}{2}\displaystyle\int_{\Omega}|f-(v+w+\nabla\cdot\mathbf{g})|^{2}d\mathbf{x}\bigg{]},\\ \mathbf{p}^{k+1}=\nabla v^{k+1}.\end{cases}

(34)

The optimality condition for $(v^{k+1},r^{k+1},\mathbf{s}^{k+1})$ is

\displaystyle\begin{cases}\begin{cases}\frac{-\nabla^{2}v^{k+1}+\nabla\cdot\mathbf{p}^{k+3/4}}{\tau}+v^{k+1}-(f-r^{k+1}-\nabla\cdot\mathbf{s}^{k+1})=0&\mbox{ in }\Omega,\\ (\nabla v^{k+1}-\mathbf{p}^{k+3/4})\cdot\mathbf{n}=0&\mbox{ on }\partial\Omega,\end{cases}\\ \begin{cases}\gamma_{2}\frac{r^{k+1}-r^{k+3/4}}{\tau}+r^{k+1}-(f-v^{k+1}-\nabla\cdot\mathbf{s}^{k+1})+2\alpha_{w}(\nabla^{4}r^{k+1})=0&\mbox{ in }\Omega,\\ \nabla^{2}r^{k+1}=0,\ \nabla^{3}r^{k+1}\cdot\mathbf{n}=0&\mbox{ on }\partial\Omega,\end{cases}\\ \begin{cases}\gamma_{3}\frac{\mathbf{s}^{k+1}-\mathbf{s}^{k+3/4}}{\tau}-\nabla(\nabla\cdot\mathbf{s}^{k+1})+2\alpha_{n}\mathbf{s}^{k+1}+\nabla(f-v^{k+1}-r^{k+1})=0&\mbox{ in }\Omega,\\ f-(v^{k+1}+r^{k+1}+\nabla\cdot\mathbf{s}^{k+1})=0&\mbox{ on }\partial\Omega.\end{cases}\end{cases}

(35)

The proposed operator splitting scheme is summarized in Algorithm 1.

0: Gray image

f

; model parameters:

\alpha_{0},\alpha_{\text{curv}},\alpha_{w},\alpha_{n}

; algorithmic parameters: time step

\tau

, evolution speed factor

\gamma_{i}

i=1,2,3

, threshold

\rho

for termination, maximal number of iteration

\text{Iter}_{\max}

1: for

k=1,\dots,\text{Iter}_{\max}

2: Update

\mathbf{p}^{k}\to\mathbf{p}^{k+1/4}

by (19), and set

\boldsymbol{\lambda}^{k+1/4}=\boldsymbol{\lambda}^{k}

r^{k+1/4}=r^{k}

\mathbf{s}^{k+1/4}=\mathbf{s}^{k}

3: Update

\mathbf{p}^{n+1/4}\to\mathbf{p}^{k+2/4}

by (20),

\boldsymbol{\lambda}^{k+1/4}\to\boldsymbol{\lambda}^{k+2/4}

by solving (22), and set

r^{k+2/4}=r^{k+1/4}

\mathbf{s}^{k+2/4}=\mathbf{s}^{k+1/4}

4: Update

\mathbf{p}^{k+2/4}\to\mathbf{p}^{k+3/4}

and

\boldsymbol{\lambda}^{k+2/4}\to\boldsymbol{\lambda}^{k+3/4}

by (32), and set

r^{k+3/4}=r^{k+2/4}

\mathbf{s}^{k+3/4}=\mathbf{s}^{k+2/4}

5: Compute

v^{k+1}

r^{k+1}

, and

\mathbf{s}^{k+1}

by solving (35), and set

\mathbf{p}^{k+1}=\nabla v^{k+1}

n^{k+1}=\nabla\cdot\mathbf{s}^{k+1}

\boldsymbol{\lambda}^{k+1}=\boldsymbol{\lambda}^{k+3/4}

6: if

\max\{\frac{\|r^{k+1}-r^{k}\|_{2}}{\|r^{k}\|_{2}},\frac{\|v^{k+1}-v^{n}\|_{2}}{\|v^{k}\|_{2}}\}<\rho

then

7: break.

8: end if

9: end for

10: return the converged homogeneous structure

v^{*}

, smooth part

w^{*}

, and oscillatory part

n^{*}

, such that

f

is approximated by

v^{*}+w^{*}+n^{*}

Algorithm 1 Proposed operator-splitting algorithm

3.5 Adapt to periodic boundary conditions for efficiency

We note that while (22) and (35) contain high-order derivatives, they are linear. Therefore, we can leverage the Fast Fourier Transform (FFT) to convert these PDE systems to pointwise algebraic systems, for which we can deduce explicit solutions. To garnish such computational efficiency, we extend our model (7) to periodic data, and this requires slight modifications on some boundary conditions.

Without loss of generality, consider a rectangular domain $\Omega=[0,L_{1}]\times[0,L_{2}]$ and define the spaces

\displaystyle H^{k}_{P}(\Omega)=\{v\in H^{k}(\Omega):v(0,y)=v(L_{1},y),\ v(:,0)=v(:,L_{2})\}

(36)

for $k=1,2$ . We revise the set $\Sigma_{f}$ to

	$\displaystyle\Sigma_{f}=\big{\{}(\mathbf{q},w,\mathbf{g})\in($	$\displaystyle L^{2}(\Omega))^{2}\times L^{2}(\Omega)\times(H^{1}(\Omega))^{2}\|$
		$\displaystyle\exists v\in H_{P}^{1}(\Omega)\leavevmode\nobreak\ \text{such that}\leavevmode\nobreak\ \mathbf{q}=\nabla v\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \int_{\Omega}(v+w+\nabla\cdot\mathbf{g}-f)\,d\mathbf{x}=0\big{\}}.$		(37)

Problem (12) and (11) are replaced by

	$\displaystyle\min_{\begin{subarray}{c}\mathbf{q}\in(H_{P}^{1}(\Omega))^{2},\boldsymbol{\mu}\in(H^{1}_{P}(\Omega))^{2},\\ w\in H^{2}_{P}(\Omega),\mathbf{g}\in(H^{1}_{P}(\Omega))^{2}\end{subarray}}\bigg{[}$	$\displaystyle\alpha_{0}\int_{\Omega}\\|\mathbf{q}\\|_{0}\,d\mathbf{x}+\alpha_{\rm curv}\int_{\Omega}\left\|\nabla\cdot\boldsymbol{\mu}\right\|^{2}\left\|\mathbf{q}\right\|\,d\mathbf{x}+\alpha_{w}\int_{\Omega}\|\nabla^{2}w\|^{2}d\mathbf{x}$
		$\displaystyle+\alpha_{n}\int_{\Omega}\|\mathbf{g}\|^{2}d\mathbf{x}+\frac{1}{2}\int_{\Omega}\|f-(v_{\mathbf{q}}+w+\nabla\cdot\mathbf{g})\|^{2}d\mathbf{x}+I_{\Sigma_{f}}(\mathbf{q},w,\mathbf{g})+I_{S}(\mathbf{q},\boldsymbol{\mu})\bigg{]}$		(38)

and

\displaystyle\begin{cases}\nabla^{2}v_{\mathbf{q},w,\mathbf{g}}=\nabla\cdot\mathbf{q}&\mbox{ in }\Omega,\\ \int_{\Omega}v_{\mathbf{q},w,\mathbf{g}}\,d\mathbf{x}=\int_{\Omega}f-w-\nabla\cdot\mathbf{g}\,d\mathbf{x},\\ v_{\mathbf{q},w,\mathbf{g}}(0,x_{2})=v_{\mathbf{q},w,\mathbf{g}}(L_{1},x_{2})&\mbox{ for }0\leq x_{2}\leq L_{2},\\ v_{\mathbf{q},w,\mathbf{g}}(x_{1},0)=v_{\mathbf{q},w,\mathbf{g}}(x_{1},L_{2})&\mbox{ for }0\leq x_{1}\leq L_{1},\\ \left(\frac{\partial v_{\mathbf{q},w,\mathbf{g}}}{\partial x_{1}}-q_{1}\right)(0,x_{2})=\left(\frac{\partial v_{\mathbf{q},w,\mathbf{g}}}{\partial x_{1}}-q_{1}\right)(L_{1},x_{2})&\mbox{ for }0\leq x_{2}\leq L_{2},\\ \left(\frac{\partial v_{\mathbf{q},w,\mathbf{g}}}{\partial x_{2}}-q_{2}\right)(x_{1},0)=\left(\frac{\partial v_{\mathbf{q},w,\mathbf{g}}}{\partial x_{2}}-q_{2}\right)(x_{1},L_{2})&\mbox{ for }0\leq x_{1}\leq L_{1}.\end{cases}

(39)

For subproblem solvers, we replace (21) and (22) by

\displaystyle\boldsymbol{\lambda}^{k+2/4}=\operatorname*{arg\,min}_{\boldsymbol{\mu}\in(H_{P}^{1}(\Omega))^{2}}\left[\frac{\gamma_{1}}{2}\int_{\Omega}|\boldsymbol{\mu}-\boldsymbol{\lambda}^{k+1/4}|_{2}^{2}d\mathbf{x}+\tau\alpha_{w}\int_{\Omega}\left|\nabla\cdot\boldsymbol{\mu}\right|^{2}|\mathbf{p}^{k+2/4}|d\mathbf{x}\right]

(40)

and

\displaystyle\begin{cases}\gamma_{1}\frac{\boldsymbol{\lambda}^{k+2/4}-\boldsymbol{\lambda}^{k+1/4}}{\tau}-2\alpha_{w}\nabla(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4})=0&\mbox{ in }\Omega,\\ \left(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4}\right)(0,x_{2})=\left(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4}\right)(L_{1},x_{2})&\mbox{ for }0\leq x_{2}\leq L_{2},\\ \left(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4}\right)(x_{1},0)=\left(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4}\right)(x_{1},L_{1})&\mbox{ for }0\leq x_{1}\leq L_{1},\end{cases}

(41)

respectively.

Replace (33), (34) and (35) by

	$\displaystyle\min_{\begin{subarray}{c}(\mathbf{q},w,\mathbf{g})\in\Sigma_{f},w\in H_{P}^{2}(\Omega),\\ \mathbf{g}\in(H_{P}^{1}(\Omega))^{2}\end{subarray}}$	$\displaystyle\bigg{[}\frac{1}{2}\int_{\Omega}\|\mathbf{q}-\mathbf{p}^{k+3/4}\|^{2}d\mathbf{x}+\frac{\gamma_{2}}{2}\int_{\Omega}\|w-r^{k+3/4}\|^{2}d\mathbf{x}+\frac{\gamma_{3}}{2}\int_{\Omega}\|\mathbf{g}-\mathbf{s}^{k+3/4}\|^{2}d\mathbf{x}$
		$\displaystyle+\tau\alpha_{w}\int_{\Omega}\|\nabla^{2}w\|^{2}d\mathbf{x}+\tau\alpha_{n}\int_{\Omega}\|\mathbf{g}\|^{2}d\mathbf{x}+\frac{\tau}{2}\int_{\Omega}\|f-(v_{\mathbf{q},w,\mathbf{g}}+w+\nabla\cdot\mathbf{g})\|^{2}d\mathbf{x}\bigg{]},$		(42)

\displaystyle\begin{cases}(v^{k+1},r^{k+1},\mathbf{s}^{k+1})=\displaystyle\operatorname*{arg\,min}_{\begin{subarray}{c}v\in H_{P}^{1}(\Omega),w\in H_{P}^{2}(\Omega),\\ \mathbf{g}\in(H_{P}^{1}(\Omega))^{2}\end{subarray}}\bigg{[}\frac{1}{2}\int_{\Omega}|\nabla v-\mathbf{p}^{k+3/4}|^{2}d\mathbf{x}+\frac{\gamma_{2}}{2}\int_{\Omega}|w-r^{k+3/4}|^{2}d\mathbf{x}\\ \hskip 113.81102pt+\frac{\gamma_{3}}{2}\displaystyle\int_{\Omega}|\mathbf{g}-\mathbf{s}^{k+3/4}|^{2}d\mathbf{x}+\tau\alpha_{w}\int_{\Omega}|\nabla^{2}w|^{2}d\mathbf{x}+\tau\alpha_{n}\int_{\Omega}|\mathbf{g}|^{2}d\mathbf{x}\\ \hskip 113.81102pt+\frac{\tau}{2}\displaystyle\int_{\Omega}|f-(v+w+\nabla\cdot\mathbf{g})|^{2}d\mathbf{x}\bigg{]},\\ \mathbf{p}^{k+1}=\nabla v^{k+1},\end{cases}

(43)

and

\displaystyle\begin{cases}\begin{cases}\frac{-\nabla^{2}v^{k+1}+\nabla\cdot\mathbf{p}^{k+3/4}}{\tau}+v^{k+1}-(f-r^{k+1}-\nabla\cdot\mathbf{s}^{k+1})=0,&\mbox{ in }\Omega,\\ \left(\frac{\partial v^{k+1}}{\partial x_{1}}-p_{1}^{k+3/4}\right)(0,x_{2})=\left(\frac{\partial v^{k+1}}{\partial x_{1}}-p_{1}^{k+3/4}\right)(L_{1},x_{2})&\mbox{ for }0\leq x_{2}\leq L_{2},\\ \left(\frac{\partial v^{k+1}}{\partial x_{2}}-p_{2}^{k+3/4}\right)(x_{1},0)=\left(\frac{\partial v^{k+1}}{\partial x_{2}}-p_{2}^{k+3/4}\right)(x_{1},L_{2})&\mbox{ for }0\leq x_{1}\leq L_{1},\end{cases}\\ \begin{cases}\gamma_{2}\frac{r^{k+1}-r^{k+3/4}}{\tau}+r^{k+1}-(f-v^{k+1}-\nabla\cdot\mathbf{s}^{k+1})+2\alpha_{w}(\nabla^{4}r^{k+1})=0&\mbox{ in }\Omega,\\ (\nabla^{2}r^{k+1})(0,x_{2})=(\nabla^{2}r^{k+1})(L_{1},x_{2}),\ \left(\frac{\partial(\nabla^{2}r^{k+1})}{\partial x_{1}}\right)(0,x_{2})=\left(\frac{\partial(\nabla^{2}r^{k+1})}{\partial x_{1}}\right)(L_{1},x_{2})&\mbox{ for }0\leq x_{2}\leq L_{2},\\ (\nabla^{2}r^{k+1})(x_{1},0)=(\nabla^{2}r^{k+1})(x_{1},L_{2}),\ \left(\frac{\partial(\nabla^{2}r^{k+1})}{\partial x_{2}}\right)(x_{1},0)=\left(\frac{\partial(\nabla^{2}r^{k+1})}{\partial x_{2}}\right)(x_{1},L_{2})&\mbox{ for }0\leq x_{1}\leq L_{1},\end{cases}\\ \begin{cases}\gamma_{3}\frac{\mathbf{s}^{k+1}-\mathbf{s}^{k+3/4}}{\tau}-\nabla(\nabla\cdot\mathbf{s}^{k+1})+2\alpha_{n}\mathbf{s}^{k+1}+\nabla(f-v^{k+1}-r^{k+1})=0&\mbox{ in }\Omega,\\ (f-(v^{k+1}+r^{k+1}+\nabla\cdot\mathbf{s}^{k+1}))(0,x_{2})=(f-(v^{k+1}+r^{k+1}+\nabla\cdot\mathbf{s}^{k+1}))(L_{1},x_{2})&\mbox{ for }0\leq x_{2}\leq L_{2},\\ (f-(v^{k+1}+r^{k+1}+\nabla\cdot\mathbf{s}^{k+1}))(x_{1},0)=(f-(v^{k+1}+r^{k+1}+\nabla\cdot\mathbf{s}^{k+1}))(x_{1},L_{2})&\mbox{ for }0\leq x_{1}\leq L_{1},\end{cases}\end{cases}

(44)

respectively, with $\mathbf{p}=(p_{1},p_{2})$ . In the rest of this paper, we always assume these adaptations and consider the periodic boundary condition.

4 Numerical Discretization

In this section, we provide implementation details for solving the proposed minimization problem.

4.1 Synopsis

Let $\Omega$ be a rectangular domain with $M\times N$ pixels. Denote the two spatial directions by $x_{1},x_{2}$ with step size $\Delta x_{1}=\Delta x_{2}=h$ for some $h>0$ . For an image $f$ defined on $\Omega$ , we denote $f(i,j)=f(ih,jh)$ . Assume that $f$ satisfies the periodic boundary condition. We define the backward ( $-$ ) and forward ( $+$ ) approximation for $\partial f/\partial x_{1}$ and $\partial f/\partial x_{2}$ as

	$\displaystyle\partial_{1}^{-}f(i,j)=\begin{cases}\frac{f(i,j)-f(i-1,j)}{h}&\ 1<i\leq M,\\ \frac{f(1,j)-f(M,j)}{h}&\ i=1,\end{cases}$	$\displaystyle\partial_{1}^{+}f(i,j)=\begin{cases}\frac{f(i+1,j)-f(i,j)}{h}&\ 1\leq i<M,\\ \frac{f(1,j)-f(M,j)}{h}&\ i=M,\end{cases}$
	$\displaystyle\partial_{2}^{-}f(i,j)=\begin{cases}\frac{f(i,j)-f(i,j-1)}{h}&\ 1<i\leq N,\\ \frac{f(i,1)-f(i,N)}{h}&\ i=1,\end{cases}$	$\displaystyle\partial_{2}^{+}f(i,j)=\begin{cases}\frac{f(i,j+1)-f(i,j)}{h}&\ 1\leq i<N,\\ \frac{f(i,1)-f(i,N)}{h}&\ i=N.\end{cases}$

For any scalar valued function $f$ and vector valued function $\mathbf{p}=(p_{1},p_{2})$ , the backward ( $-$ ) and forward ( $+$ ) approximation for gradient and divergence are defined as

\displaystyle\nabla^{\pm}f(i,j)=(\partial^{\pm}_{1}f(i,j),\partial^{\pm}_{2}f(i,j)),\quad\text{div}^{\pm}\mathbf{p}(i,j)=\partial_{1}^{\pm}p_{1}(i,j)+\partial^{\pm}_{2}p_{2}(i,j).

The discretized Laplacian is defined as

\displaystyle\nabla^{2}f(i,j)=\text{div}^{-}(\nabla^{+}f(i,j))=\partial_{1}^{-}\partial_{1}^{+}f(i,j)+\partial_{2}^{-}\partial_{2}^{+}f(i,j)\;,

which recovers the central difference scheme.

Define the shifting and identity operator as

\displaystyle(\mathcal{S}_{1}^{\pm}f)(i,j)=f(i\pm 1,j),\quad(\mathcal{S}_{2}^{\pm}f)(i,j)=f(i,j\pm 1),\quad(\mathcal{I}f)(i,j)=f(i,j)\;,

where the periodic boundary condition is assumed. Denote the discrete Fourier transform and its inverse by $\mathcal{F}$ and $\mathcal{F}^{-1}$ , respectively. We have

\displaystyle\mathcal{F}(\mathcal{S}_{1}^{\pm}f)(i,j)=e^{\pm 2\pi\sqrt{-1}(i-1)/M}\mathcal{F}(f)(i,j)\;,\text{and}\leavevmode\nobreak\ \mathcal{F}(\mathcal{S}_{2}^{\pm}f)(i,j)=e^{\pm 2\pi\sqrt{-1}(j-1)/N}\mathcal{F}(f)(i,j)\;.

We use $\text{Real}(\cdot)$ to denote the real part of a complex quantity.

4.2 Discrete analogue of $\mathbf{p}^{k+1/4}$ , $\mathbf{p}^{k+2/4}$ and $\boldsymbol{\lambda}^{k+2/4}$

For $\mathbf{p}^{k+1/4}$ , we discretize (19) as

\displaystyle\mathbf{p}^{k+1/4}(i,j)=\begin{cases}\mathbf{0}&\mbox{ if }(p_{1}^{k}(i,j))^{2}+(p_{2}^{k}(i,j))^{2}\leq\alpha_{0}\tau,\\ \mathbf{p}^{k}(i,j)&\mbox{ otherwise},\end{cases}

(45)

and the computational cost is $\mathcal{O}(MN)$ .

According to (20), $\mathbf{p}^{k+2/4}$ is updated pointwisely. We thus have

\displaystyle\mathbf{p}^{k+2/4}(i,j)=\max\left\{0,1-\frac{\tau\alpha_{w}\left|\text{div}^{-}\boldsymbol{\lambda}^{k+1/4}(i,j)\right|^{2}}{\sqrt{(p_{1}^{k+1/4}(i,j))^{2}+(p_{2}^{k+2/4}(i,j))^{2}}}\right\}\mathbf{p}^{k+1/4}(i,j)\;,

(46)

whose computational cost is $\mathcal{O}(MN)$ .

To compute $\boldsymbol{\lambda}^{k+2/4}=(\lambda_{1}^{k+2/4},\lambda_{2}^{k+2/4})$ , the updating formula (41) can be rewritten as

\displaystyle\gamma_{1}\boldsymbol{\lambda}^{k+2/4}-2\tau\alpha_{w}\nabla(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+2/4})=\gamma_{1}\boldsymbol{\lambda}^{k+1/4}\;.

(47)

Instead of directly solving (47), we use the frozen coefficient method [50, 27, 36] and consider the following approximation

\displaystyle\gamma_{1}\boldsymbol{\lambda}^{k+2/4}-c\nabla^{2}\boldsymbol{\lambda}^{k+2/4}=\gamma_{1}\boldsymbol{\lambda}^{k+1/4}-c\nabla^{2}\boldsymbol{\lambda}^{k+1/4}+2\tau\alpha_{w}\nabla(|\mathbf{p}^{k+2/4}|\nabla\cdot\boldsymbol{\lambda}^{k+1/4})

(48)

for some small coefficient $c>0$ . In our experiment, we set $c=1\times 10^{-9}$ . We discretize (48) as

\displaystyle\gamma_{1}\boldsymbol{\lambda}^{k+2/4}-c\nabla^{+}(\text{div}^{-}\boldsymbol{\lambda}^{k+2/4})=\gamma_{1}\boldsymbol{\lambda}^{k+1/4}-c\nabla^{+}(\text{div}^{-}\boldsymbol{\lambda}^{k+1/4})+2\tau\alpha_{w}\nabla^{+}(|\mathbf{p}^{k+2/4}|(\text{div}^{-}\boldsymbol{\lambda}^{k+1/4})).

(49)

Equation (49) can be written in a matrix form as

\displaystyle\begin{pmatrix}\gamma_{1}-c\partial_{1}^{+}\partial_{1}^{-}&-c\partial_{1}^{+}\partial_{2}^{-}\\ -c\partial_{2}^{+}\partial_{1}^{-}&\gamma_{1}-c\partial_{2}^{+}\partial_{2}^{-}\end{pmatrix}\begin{pmatrix}\lambda^{k+2/4}_{1}\\ \lambda^{k+2/4}_{2}\end{pmatrix}=\begin{pmatrix}b_{1}\\ b_{2}\end{pmatrix}

(50)

with

\displaystyle\begin{cases}&b_{1}=\gamma_{1}\lambda_{1}^{k+1/4}-c\partial_{1}^{+}(\text{div}^{-}\boldsymbol{\lambda}^{k+1/4})+2\tau\alpha_{w}\partial_{1}^{+}(|\mathbf{p}^{k+2/4}|(\text{div}^{-}\boldsymbol{\lambda}^{k+1/4})),\\ &b_{2}=\gamma_{1}\lambda_{2}^{k+1/4}-c\partial_{2}^{+}(\text{div}^{-}\boldsymbol{\lambda}^{k+1/4})+2\tau\alpha_{w}\partial_{2}^{+}(|\mathbf{p}^{k+2/4}|(\text{div}^{-}\boldsymbol{\lambda}^{k+1/4})).\end{cases}

System (50) is equivalent to

\displaystyle\begin{pmatrix}\gamma_{1}-c(\mathcal{S}_{1}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{1}^{-})/h^{2}&-c(\mathcal{S}_{1}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{2}^{-})/h^{2}\\ -c(\mathcal{S}_{2}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{1}^{-})/h^{2}&\gamma_{1}-c(\mathcal{S}_{2}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{2}^{-})/h^{2}\end{pmatrix}\begin{pmatrix}\lambda^{k+2/4}_{1}\\ \lambda^{k+2/4}_{2}\end{pmatrix}=\begin{pmatrix}b_{1}\\ b_{2}\end{pmatrix}.

(51)

Applying discrete Fourier transform on both sides of (51) gives

\displaystyle\begin{pmatrix}a_{11}&a_{12}\\ a_{21}&a_{22}\end{pmatrix}\mathcal{F}\begin{pmatrix}\mu^{k+2/4}_{1}\\ \mu^{k+2/4}_{2}\end{pmatrix}=\mathcal{F}\begin{pmatrix}b_{1}\\ b_{2}\end{pmatrix},

(52)

where

	$\displaystyle a_{11}(i,j)=\gamma_{1}-2c(\cos\zeta_{i}-1)/h^{2},\ a_{22}(i,j)=\gamma_{1}-2c(\cos\eta_{j}-1)/h^{2},$
	$\displaystyle a_{12}(i,j)=c(\cos\zeta_{i}+\sqrt{-1}\sin\zeta_{i}-1)(\cos\eta_{j}-\sqrt{-1}\sin\eta_{j}-1)/h^{2},$
	$\displaystyle a_{21}(i,j)=c(\cos\eta_{j}+\sqrt{-1}\sin\eta_{j}-1)(\cos\zeta_{i}-\sqrt{-1}\sin\zeta_{i}-1)/h^{2},$

with $\zeta_{i}=2\pi(i-1)/M,\ \eta_{j}=2\pi(j-1)/N$ for $i=1,...,M$ and $j=1,...,N$ . We thus update $\boldsymbol{\lambda}^{k+2/4}$ via

\displaystyle\begin{pmatrix}\lambda_{1}^{k+2/4}\\ \lambda_{2}^{k+2/4}\end{pmatrix}=\text{Real}\left(\mathcal{F}^{-1}\left[\frac{1}{a_{11}a_{22}-a_{12}a_{21}}\begin{pmatrix}a_{11}\mathcal{F}(b_{1})-a_{12}\mathcal{F}(b_{2})\\ -a_{21}\mathcal{F}(b_{1})+a_{22}\mathcal{F}(b_{2})\end{pmatrix}\right]\right)\;.

(53)

Since the linear system is solved explicitly under the forward and inverse FFT, the computational cost is $\mathcal{O}(MN(\log M+\log N))$ .

4.3 Discrete analogue of $\mathbf{p}^{k+3/4}$ , $\boldsymbol{\lambda}^{k+3/4}$

According to Section 3.3, $\mathbf{p}^{k+3/4}$ and $\boldsymbol{\lambda}^{k+3/4}$ can be computed pointwisely. On $\chi_{1}$ , for each $i,j$ , we set $\widehat{\mathbf{p}}^{k+3/4}(i,j)=\mathbf{0}$ and discretize (25) as

\displaystyle\widehat{\boldsymbol{\lambda}}^{k+3/4}(i,j)=\frac{\boldsymbol{\lambda}^{k+2/4}(i,j)}{\max\left\{1,\sqrt{|\boldsymbol{\lambda}^{k+2/4}(i,j)|^{2}}\right\}}.

(54)

On $\chi_{2}$ , we discretize (31) as

\displaystyle\widetilde{\boldsymbol{\lambda}}^{k+3/4}(i,j)=\frac{\theta^{*}\mathbf{p}^{k+2/4}(i,j)+\gamma_{1}\boldsymbol{\lambda}^{k+2/4}(i,j)}{|\theta^{*}\mathbf{p}^{k+2/4}(i,j)+\gamma_{1}\boldsymbol{\lambda}^{k+2/4}(i,j)|},\quad\widetilde{\mathbf{p}}^{k+3/4}(i,j)=\theta^{*}\widetilde{\boldsymbol{\lambda}}^{k+3/4}(i,j)\;,

where $\theta^{*}$ is computed by (30) with $\mathbf{y}(\mathbf{x})=\mathbf{p}^{k+2/4}(i,j),\mathbf{w}(\mathbf{x})=\boldsymbol{\lambda}^{k+2/4}(i,j)$ . Functions $\mathbf{p}^{k+3/4}$ and $\boldsymbol{\lambda}^{k+3/4}$ are then updated using (32).

4.4 Discrete analogue of $\mathbf{p}^{k+1},r^{k+1}$ and $\mathbf{s}^{k+1}$

Reorganizing the system (44), we get a linear system in $v^{k+1},r^{k+1}$ and $\mathbf{s}^{k+1}=(s_{1}^{k+1},s_{2}^{k+1})$

\displaystyle\begin{cases}(-\nabla^{2}+\tau I)v^{k+1}+\tau r^{k+1}+\tau\nabla\cdot\mathbf{s}^{k+1}=-\nabla\cdot\mathbf{p}^{k+3/4}+\tau f,\\ \tau v^{k+1}+[(\gamma_{2}+\tau)I+2\tau\alpha_{w}\nabla^{4}]r^{k+1}+\tau\nabla\cdot\mathbf{s}^{k+1}=\gamma_{2}r^{k+3/4}+\tau f,\\ -\tau\nabla v^{k+1}-\tau\nabla r^{k+1}+[(\gamma_{3}+2\tau\alpha_{n})I-\tau\nabla\text{div}]\mathbf{s}^{k+1}=\gamma_{3}\mathbf{s}^{k+3/4}-\tau\nabla f.\end{cases}

(55)

Upon discretization, we obtain a linear system

\displaystyle\mathbf{A}\begin{bmatrix}v^{k+1}\\ r^{k+1}\\ s_{1}^{k+1}\\ s_{2}^{k+1}\end{bmatrix}=\begin{bmatrix}b^{k}_{1}\\ b^{k}_{2}\\ b^{k}_{3}\\ b^{k}_{4}\end{bmatrix}\leavevmode\nobreak\ \text{with}\leavevmode\nobreak\ \begin{cases}&b^{k}_{1}=-\text{div}^{-}\mathbf{p}^{k+3/4}+\tau f,\\ &b^{k}_{2}=\gamma_{2}r^{k+3/4}+\tau f,\\ &b^{k}_{3}=\gamma_{3}s_{1}^{k+3/4}-\tau\partial_{1}^{+}f,\\ &b^{k}_{4}=\gamma_{3}s_{2}^{k+3/4}-\tau\partial_{2}^{+}f,\end{cases}

where

	$\displaystyle\mathbf{A}=$
	$\displaystyle\begin{bmatrix}\begin{aligned} &-[(\mathcal{I}-\mathcal{S}_{1}^{-})(\mathcal{S}_{1}^{+}-\mathcal{I})/h^{2}\\ &+(\mathcal{I}-\mathcal{S}_{2}^{-})(\mathcal{S}_{2}^{+}-\mathcal{I})/h^{2}]\\ &+\tau\mathcal{I}\end{aligned}&\tau&\tau(\mathcal{I}-\mathcal{S}_{1}^{-})/h&\tau(\mathcal{I}-\mathcal{S}_{2}^{-})/h\\ \tau&(\gamma_{2}+\tau)\mathcal{I}+2\tau\alpha_{w}\mathcal{P}&\tau(\mathcal{I}-\mathcal{S}_{1}^{-})/h&\tau(\mathcal{I}-\mathcal{S}_{2}^{-})/h\\ -\tau(\mathcal{S}_{1}^{+}-\mathcal{I})/h&-\tau(\mathcal{S}_{1}^{+}-\mathcal{I})/h&\begin{aligned} &\gamma_{3}+2\tau\alpha_{n}\\ &-\tau(\mathcal{S}_{1}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{1}^{-})/h^{2}\end{aligned}&-\tau(\mathcal{S}_{1}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{2}^{-})/h^{2}\\ -\tau(\mathcal{S}_{2}^{+}-\mathcal{I})/h&-\tau(\mathcal{S}_{2}^{+}-\mathcal{I})/h&-\tau(\mathcal{S}_{2}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{1}^{-})/h^{2}&\begin{aligned} &\gamma_{3}+2\tau\alpha_{n}\\ &-\tau(\mathcal{S}_{2}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{2}^{-})/h^{2}\end{aligned}\end{bmatrix}$

with

	$\displaystyle\mathcal{P}=$	$\displaystyle\frac{1}{h^{4}}\Big{[}((\mathcal{I}-\mathcal{S}_{1}^{-})(\mathcal{S}_{1}^{+}-\mathcal{I}))^{2}+(\mathcal{I}-\mathcal{S}_{2}^{-})(\mathcal{S}_{2}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{1}^{-})(\mathcal{S}_{1}^{+}-\mathcal{I})$
		$\displaystyle\quad+(\mathcal{I}-\mathcal{S}_{1}^{-})(\mathcal{S}_{1}^{+}-\mathcal{I})(\mathcal{I}-\mathcal{S}_{2}^{-})(\mathcal{S}_{2}^{+}-\mathcal{I})+((\mathcal{I}-\mathcal{S}_{2}^{-})(\mathcal{S}_{2}^{+}-\mathcal{I}))^{2}\Big{]}.$

Note that $\mathcal{A}$ does not depend on the iteration, thus it can be constructed beforehand and applied repeatedly in each iteration. Applying the discrete Fourier transform $\mathcal{F}$ , we get

\displaystyle\mathbf{D}\ \mathcal{F}\left(\begin{bmatrix}v^{k+1}\\ r^{k+1}\\ s_{1}^{k+1}\\ s_{2}^{k+1}\end{bmatrix}\right)=\mathcal{F}\left(\begin{bmatrix}b^{k}_{1}\\ b^{k}_{2}\\ b^{k}_{3}\\ b^{k}_{4}\end{bmatrix}\right)

(56)

with $\mathbf{D}\in\mathbb{R}^{4\times 4}$ whose $(l,m)$ -th entries $d_{lm}$ , $1\leq l,m\leq 4$ are listed below

	$\displaystyle d_{11}(i,j)=(\tau-2(\cos\zeta_{i}-1)/h^{2}-2(\cos\eta_{j}-1)/h^{2}),\quad d_{12}=\tau,$
	$\displaystyle d_{13}(i,j)=\tau[1-(\cos\zeta_{i}-\sqrt{-1}\sin\zeta_{i})]/h,\quad d_{14}(i,j)=\tau[1-(\cos\eta_{j}-\sqrt{-1}\sin\eta_{j})]/h,$
	$\displaystyle d_{21}=\tau,\quad d_{22}(i,j)=\gamma_{2}+\tau+2\tau\alpha_{w}\frac{4}{h^{4}}\left[(\cos\zeta_{i}-1)^{2}+2(\cos\zeta_{i}-1)(\cos\eta_{j}-1)+(\cos\eta_{j}-1)^{2}\right],$
	$\displaystyle d_{23}=d_{13},\quad d_{24}=d_{14},$
	$\displaystyle d_{31}(i,j)=d_{32}(i,j)=-\tau(\cos\zeta_{i}+\sqrt{-1}\sin\zeta_{i}-1)/h,$
	$\displaystyle d_{33}(i,j)=\gamma_{3}+2\tau\alpha_{n}-2\tau(\cos\zeta_{i}-1)/h^{2},$
	$\displaystyle d_{34}(i,j)=\tau(\cos\zeta_{i}+\sqrt{-1}\sin\zeta_{i}-1)(\cos\eta_{j}-\sqrt{-1}\sin\eta_{j}-1)/h^{2},$
	$\displaystyle d_{41}(i,j)=d_{42}(i,j)=-\tau(\cos\eta_{j}+\sqrt{-1}\sin\eta_{j}-1)/h,$
	$\displaystyle d_{43}(i,j)=\tau(\cos\eta_{j}+\sqrt{-1}\sin\eta_{j}-1)(\cos\zeta_{i}-\sqrt{-1}\sin\zeta_{i}-1)/h^{2},$
	$\displaystyle d_{44}(i,j)=\gamma_{3}+2\tau\alpha_{n}-2\tau(\cos\eta_{j}-1)/h^{2}.$

To stabilize the computation, we add a small constant $\kappa=1\times 10^{-9}$ on the diagonal of $\mathbf{D}$ . The inverse $\mathbf{D}^{-1}$ can be explicitly computed using the Cramer’s rule. We can efficiently get $v^{k+1},r^{k+1}$ and $\mathbf{s}^{k+1}$ by inverse FFT

\displaystyle\begin{bmatrix}v^{k+1}\\ r^{k+1}\\ s_{1}^{k+1}\\ s_{2}^{k+1}\end{bmatrix}=\text{Real}\left(\mathcal{F}^{-1}\left[\mathbf{D}^{-1}\mathcal{F}\left(\begin{bmatrix}b^{k}_{1}\\ b^{k}_{2}\\ b^{k}_{3}\\ b^{k}_{4}\end{bmatrix}\right)\right]\right)\;.

(57)

Similarly to (53), as we obtain the updated variables by explicit formula, the computational cost is still $\mathcal{O}(MN(\log M+\log N))$ . Overall, the computational cost of a single iteration of our proposed algorithm that updates $(\mathbf{p}^{k},\boldsymbol{\lambda}^{k},r^{k},\mathbf{s}^{k})$ to $(\mathbf{p}^{k+1},\boldsymbol{\lambda}^{k+1},r^{k+1},\mathbf{s}^{k+1})$ is in the order of $\mathcal{O}(MN(\log M+\log N))$ .

5 Numerical Experiments

In this section, we present various experiments to show the behaviors of the proposed model, analyze the effects of model parameters, and justify its effectiveness in the task of decomposition and denoising by ablation and comparison studies. To evaluate the image quality, we employ the PSNR metric

\displaystyle\text{PSNR}(u_{\text{ref}},\widehat{u})=10\log_{10}\left(\frac{1}{\|\widehat{u}-u_{\text{ref}}\|_{2}}\right)\;,

where $u_{\text{ref}}$ is the reference image, $\widehat{u}$ is the corresponding approximation, and their intensity values are scaled to $[0,1]$ . We test our method on various synthetic images as well as grayscale photos. For model parameters $\alpha_{0},\alpha_{\text{curv}},\alpha_{w}$ and $\alpha_{n}$ , we discuss their selections in Section 5.5 and fine-tune them according to the image contents. For algorithmic parameters, we fix $\tau=0.1,\gamma_{1}=0.01,\gamma_{2}=20$ , and $\gamma_{3}=1$ . We set the global constants for numerical stability $c=1\times 10^{-9}$ and $\kappa=1\times 10^{-9}$ ; here $c$ is the frozen coefficient in (48), and $\kappa$ is added to the diagonal of $\mathbf{D}$ in (56) for stability. We take the maximal number of iterations $\text{Iter}_{\max}=1000$ , and fix the threshold for the terminating condition $\rho=1\times 10^{-6}$ , i.e., the algorithm is terminated if $\max\left\{\frac{\|r^{k+1}-r^{k}\|_{2}}{\|r^{k}\|_{2}},\frac{\|v^{k+1}-v^{n}\|_{2}}{\|v^{k}\|_{2}}\right\}<\rho$ . For the default initial conditions, we fix $v_{0}=0.001\times\widetilde{f}+0.999\times 0.5$ where $\widetilde{f}$ is the input $f$ convolved with a standard Gaussian kernel, $\mathbf{p}_{0}=\nabla^{+}v_{0}$ , $r_{0}=f-v_{0}$ , and $\boldsymbol{\lambda}_{0}=\mathbf{s}_{0}=\mathbf{0}$ . To reduce the boundary artifacts, we pad the input with 30 pixels on each side in a symmetric manner. The algorithm is implemented with MATLAB 2022a in the environment of Windows 11, and all the experiments are conducted in a laptop equipped with 12th Gen Intel(R) Core(TM) i7-12800H, 2400Mhz, and 14 cores.

5.1 General performances

We illustrate the general performances of the proposed model using two sets of experiments. In Figure 3, we show the decomposition results for synthetic images generated by convex combinations of piecewise constant components and smoothly varying shades. In Figure 4, we show the results for grayscale photos with richer details. For both sets of examples, Gaussian noise with intensity $\sigma=20/255$ is added.

In the first row of Figure 3, a white cross in the dark is shed by a soft light at the center. The proposed model successfully extracts the cross $v^{*}$ where the intensity variation due to the lighting is removed, and the important geometric features such as sharp corners as well as straight boundaries are accurately reconstructed. The lighting is captured by the $w^{*}$ component, and it is well separated from $v^{*}$ . The $n^{*}$ component corresponds to the noise whose standard deviation is close to the true noise level. In the second row, we test the proposed model on a combination of silhouettes of ellipses and a circular wave of lighting. Under the non-uniform lighting, the ellipses are correctly identified in $v^{*}$ and filled with nearly constant intensity. The $w^{*}$ component shows the correct wave pattern of the lighting. In the third row, we have a shaded globe, where the source of the light is fixed at the top left corner. The proposed method correctly assigns in the $v^{*}$ component the silhouettes of the globe where boundaries of continents are sharp. Noticeably, the proposed model keeps the soft shadow characterizing the volumetric feature of a ball in the $w^{*}$ component with the correct direction of lighting. In the last row, we have a challenging case where a binarized portrait of Einstein is blended with multiple blurry shades of the mass-energy equation. The $w^{*}$ component is more complicated than those in the previous examples. Our method successfully extracts the soft contents with complex shapes and leaves a relatively clean configuration for the $v^{*}$ component.

In Figure 4, we apply the proposed model to photo-realistic grayscale images. In all examples, we observe that contents with sharp boundaries are captured by $v^{*}$ ; soft shadows and illumination variation are modeled by $w^{*}$ ; and the $n^{*}$ keeps fine textures plus noise. Noticeably, in the first row, the rectangular reflectances on the dome are separated from the soft light and shadow defining the volume of the dome. In the second row, the foreground with crisp boundaries is distinguished from the background with blurry shapes.

The examples in Figure 3 and 4 justify the effectiveness of the proposed method. Thanks to the curvature regularization, the extracted $v^{*}$ component preserves the defining geometric features of the homogeneous part in the given input image, such as sharp boundaries and corners. The extracted $w^{*}$ component consistently matches the soft intensity variation caused by lighting and shadows. The extracted $n^{*}$ component holds the oscillatory residuals including noise and fine-scale textures. We note that when the input images do not contain rich textures, the $n^{*}$ components rendered by the proposed model do not contain structural details such as boundaries. It allows us to use the summation of $v^{*}$ and $w^{*}$ , which is denoted by $u^{*}$ , to reconstruct the clean image. This is justified by the PSNR values reported in Figure 3 and 4.

5.2 Ablation study

To justify the combination of regularizations in our model (7), we compare it with the following settings:

Model I. An image $f$ is decomposed into two parts: cartoon part $u^{*}$ and oscillatory part $n^{*}$

\displaystyle\{u^{*},n^{*}\}=\operatorname*{arg\,min}_{u,n}\left[\alpha_{0}\|\nabla u\|_{0}+\alpha_{n}\|n\|^{2}_{H^{-1}(\Omega)}+\frac{1}{2}\int_{\Omega}\left|f-(u+n)\right|^{2}\,d\mathbf{x}\right]\;.

(58)

Model II. An image $f$ is decomposed into three parts: cartoon part $v^{*}$ , smooth part $w^{*}$ , and oscillatory part $n^{*}$

\displaystyle\{v^{*},w^{*},n^{*}\}=\operatorname*{arg\,min}_{v,w,n}\left[\alpha_{0}\|\nabla v\|_{0}+\alpha_{w}\|\Delta w\|_{L^{2}(\Omega)}^{2}+\alpha_{n}\|n\|^{2}_{H^{-1}(\Omega)}+\frac{1}{2}\int_{\Omega}\left|f-(v+w+n)\right|^{2}\,d\mathbf{x}\right]\;.

(59)

Model III. An image $f$ is decomposed into two parts: homogeneous structure $u^{*}$ and oscillatory part $n^{*}$

\displaystyle\{u^{*},n^{*}\}=\operatorname*{arg\,min}_{u,n}\left[\alpha_{0}\|\nabla u\|_{0}+\alpha_{\text{curv}}\int_{\Omega}\left(\nabla\cdot\frac{\nabla u}{|\nabla u|}\right)^{2}\left|\nabla u\right|\,d\mathbf{x}+\alpha_{n}\|n\|^{2}_{H^{-1}(\Omega)}+\frac{1}{2}\int_{\Omega}\left|f-(u+n)\right|^{2}\,d\mathbf{x}\right]\;.

(60)

Compared to the proposed model, there is no smoothly varying component in Model I, and the regularity of $u^{*}$ is mainly affected by the $L^{0}$ -gradient. In Model II, the noise-less component $u^{*}$ is decomposed into a cartoon part $v^{*}$ and a smooth part $w^{*}$ , the same as in the proposed model; but there is no curvature regularization on $v^{*}$ . In Model III, we impose both $L^{0}$ -gradient and curvature regularization on the noise-free approximation $u^{*}$ , but there is no component that models the smooth shading.

Figure 5 shows the $u^{*}$ components learned from the various settings above. In column (a), we present a noisy grayscale image of a squared ring, and the regions inside red boxes are zoomed-in and shown in the subsequent rows. We added Gaussian noise with intensity $\sigma=50/255$ on the input. In column (b), we show the $u^{*}$ component reconstructed by Model I (58). As was expected, using the regularization $\|\nabla u\|_{0}$ enforces strong sparsity on the gradient of reconstructed $u^{*}$ component and yields sharp boundaries. However, analogous to the TV-norm regularization, using only $L^{0}$ -gradient as regularization creates strong staircase effects. Moreover, $L^{0}$ -gradient is local, yet larger-scale geometric information is necessary to render better results. The boundary in the third row shows undesirable irregularities. The smooth region is better reconstructed in Model II, where an extra component $w^{*}$ is introduced. However, due to the lack of geometric information in $L^{0}$ -gradient, the boundary lines still present zig-zags. In Model III, the curvature regularization on $u^{*}$ is used, which imposes geometric priors on the level lines of $v^{*}$ . The straight boundary lines are better preserved compared to Model I and Model II. Without the smooth part $w^{*}$ , the intensity transition in the last row is not as smooth as in Model II. Combining the benefits of geometric regularization on the level lines and adding an extra smoothly varying component, our proposed model outperforms the others. Observe that the transition is smooth in the second row indicated by the red contour curves, and the boundary line is straight in the third row. These reconstructed features match the underlying image shown in the last column.

From this set of experiments, we conclude that $L^{0}$ -gradient regularization induces strong sparsity on the gradient which causes significant staircase artifacts and irregular boundaries. These undesirable phenomenons are effectively mitigated by the curvature regularization on the level lines. In order to recover smoothly varying regions while preserving sharp boundaries, our experiments show that adding an extra smooth part in the decomposition is important.

5.3 Effect of parameter $\alpha_{w}$ and $\alpha_{\text{curv}}$

In our proposed model, both $v^{*}$ and $w^{*}$ favor regularity. Particularly, $v^{*}$ prefers shapes with definite and regular boundaries; whereas $w^{*}$ tends to be blurry. We demonstrate that image content can be distributed and converted between these two components, yielding different decomposition results.

First, we study the effect of the model parameter $\alpha_{w}$ associated with the term $\int_{\Omega}\|\nabla^{2}w\|_{2}^{2}\,d\mathbf{x}$ that encourages isotropic smoothness of the $w^{*}$ component. Intuitively, larger values of $\alpha_{w}$ induce smoother $w^{*}$ where image details are smeared. If $\alpha_{w}$ is too large, regions with smoothly varying intensity and sufficiently regular level lines can be captured by the $v^{*}$ component.

We illustrate this conversion between $w^{*}$ and $v^{*}$ affected by $\alpha_{w}$ in Figure 6. Given a noisy image ( $\sigma=10/255$ ), we apply the proposed model with $\alpha_{w}$ increasing from $10$ to $1000$ while fixing the other parameters. For each value of $\alpha_{w}$ , we show the scaled $v^{*}$ and $w^{*}$ components in the first and second row, respectively, and in the third row, we show the respective sum $u^{*}=v^{*}+w^{*}$ . Observe that with a wide range of values of $\alpha_{w}$ , the noise-free reconstruction $u^{*}$ remains visually similar, and the PSNR values are stable. When $\alpha_{w}$ is small, most of the intensity variation is captured by $w^{*}$ , while $v^{*}$ only captures objects’ large-scale contours. Although hair, face, eyes, and many recognizable shapes are visible in $w^{*}$ , these distinct parts are blended with each other smoothly without sharp boundaries. As $\alpha_{w}$ increases, we see that $w^{*}$ becomes more blurry, and $v^{*}$ shows richer intensity variations such as those in the background. Note that these smoothly varying components in $v^{*}$ have definite boundaries. For example, the zoom-ins of the red boxed region in the fourth row show that intensity transitions become sharper when $\alpha_{w}$ increases. Other softer smooth parts without distinctive boundaries, such as the shades on the cheeks and chin, are captured by the $w^{*}$ component.

Second, we show the effect of the parameter $\alpha_{\text{curv}}$ which is closely related to the regularity of the level lines in the $v^{*}$ component. Similar to the previous study on the effect of $\alpha_{w}$ , by changing $\alpha_{\text{curv}}$ , smoothly varying components can transfer between $w^{*}$ and $v^{*}$ . In particular, if $\alpha_{\text{curv}}$ is large, certain regions with smooth intensity variation can be pushed into $w^{*}$ for lower total loss.

Figure 7 presents results of the proposed model on a noisy image $(\sigma=10/255)$ with varying values of $\alpha_{\text{curv}}$ when the other regularization parameters are fixed. When $\alpha_{\text{curv}}$ increases from $0.01$ to $10$ , we observe that soft shades on peppers in $v^{*}$ are suppressed while the boundaries are preserved. Meanwhile, the $w^{*}$ component gains more variation in the intensity. If we keep increasing $\alpha_{\text{curv}}$ above $0.1$ , the intensity variation near the boundaries and the highlights due to reflection slowly fade away in $v^{*}$ . In contrast, shapes of peppers are more apparent in $w^{*}$ . Since the sharp gradients in $v^{*}$ are reduced and the intensity variation in $w^{*}$ is intrinsically smooth, the reconstructed $u^{*}$ component becomes slightly blurry. This transition is also reflected by the decreasing PSNR. In Figure 7, we also show the level lines of $u^{*}$ in a zoomed-in region. The effect of $\alpha_{\text{curv}}$ on $u^{*}$ is obvious. Increasing $\alpha_{\text{curv}}$ yields more regular level lines, and the staircase effect due to $L^{0}$ -gradient regularization is thus ameliorated.

5.4 Comparison on cartoon+smooth+oscillation decomposition

We compare the proposed model with models proposed by Chan et al. [13] (CEP) and Huska et al. [29] (HKLM), which are stated in (5) and (6), respectively. In both the proposed model and HKLM, an input image is decomposed into three components corresponding to cartoon-like part $v^{*}$ , smooth part $w^{*}$ , and oscillatory part $n^{*}$ . For the four-component model (5), in addition to $v^{*}$ and $w^{*}$ , we take $n^{*}=n_{1}^{*}+n_{2}^{*}$ , i.e., the oscillatory part is the summation of the texture and noise.

In Figure 8, we show the decomposition results from CEP, HKLM, and DOSH when applied to a noisy synthetic image ( $\sigma=50/255$ ), respectively. For CEP [13] (see also (5)), we set $\lambda=0.5,\alpha=10$ , $\mu=2$ and the wavelet threshold $\delta$ is determined using Stein’s Unbiased Risk Estimate (SURE) [19]. For HKLM [29] (see also (6)), the authors proposed an automatic scheme for tuning the regularization parameters. However, we find that they fail to yield good decomposition results in general for our test images. To reach good PSNR values in $u^{*}$ , after experimenting with different ranges of parameters, we manually set $\beta_{1}=1\times 10^{-4}$ , $\beta_{2}=0.5$ , $\beta_{3}=5\times 10^{-6}$ , and the MC penalty (an approximation of the $L^{0}$ -norm, see [29] for details) parameter $a=30$ . For the proposed model, we use $\alpha_{0}=2\times 10^{-2}$ , $\alpha_{\text{curv}}=0.15$ , $\alpha_{w}=50$ and $\alpha_{n}=1\times 10^{-3}$ . We also present the summations of $v^{*}$ and $w^{*}$ as the noise-less approximations for $f_{\text{noise}}$ and evaluate their respective PSNR values. In this example with heavy noise, both HKLM and our proposed model successfully separate the piecewise constant part from $f_{\text{noise}}$ with strong sparsity of gradient thanks to the $L^{0}$ -gradient regularization, and CEP with TV-norm fails to produce clean $v^{*}$ part. Moreover, we observe that, with a heavy noise, regularization using only $L^{0}$ -gradient as in HKLM yields irregular boundaries. In contrast, with the additional regularization on the curvature of level-set curves of $v^{*}$ , our proposed model renders more satisfying results.

In Figure 9, we compare these methods using a grayscale photo of butterfly with noise ( $\sigma=10/255$ ) and show the zoom-in comparisons in Figure 10. For CEP (see also (5)), we set $\lambda=5,\alpha=10,\mu=2$ and $\delta$ determined from SURE [19]. For HKLM (see also (6)), we use $\beta_{1}=1\times 10^{-4}$ , $\beta_{2}=3$ , $\beta_{3}=5\times 10^{-6}$ , and $a=30$ for the MC penalty. For the proposed method, we use $\alpha_{0}=1\times 10^{-3}$ , $\alpha_{\text{curv}}=0.05$ , $\alpha_{w}=30$ and $\alpha_{n}=10$ . All methods separate the structural part including patterns on the wing in the front, the smooth part including the blurry leaves in the back, as well as the oscillatory component such as noise. Compared to CEP, both HKML and our proposed model produce better separations between the foreground details and background variations. Since HKLM uses only $L^{0}$ -gradient for regularizing the $v^{*}$ component, we observe that $u^{*}$ of HKML has the sharpest boundaries compared to CEP and the proposed model in Figure 10. This also yields irregular curves and strong staircase effects in (d), which are effectively suppressed in CEP in (c) and the proposed model in (e). Since CEP uses TV-norm and the proposed model uses the stronger sparsity-inducing $L^{0}$ -gradient, the $v^{*}$ component and consequently the $u^{*}$ component resulted from the proposed model are clearer than those from CEP. Moreover, for this example, the noise-less approximation $u^{*}$ produced by the proposed model has the highest PSNR value compared to the other two. We conclude that our proposed model produces more visually appealing decomposition results with better quality.

Our proposed model is also more efficient than CEP and HKLM. For the image of size $512\times 512$ in Figure 9, CEP takes around 1655 seconds, HKLM takes over 4000 seconds, whereas our model only requires around 180 seconds. We note that HKLM requires to solve a dense linear system of size $(4NM)\times(4NM)$ , which costs around $\mathcal{O}((NM)^{3})$ . As discussed in Section 4, by leveraging FFT, the heaviest computation in our model requires solving a linear system of size $4\times 4$ for each pixel in parallel, and the solution can be explicitly expressed. Overall, the computational cost of each iteration of our method is in the order or $\mathcal{O}(MN(\log M+\log N))$ . This is similar for CEP, which involves both FFT and the discrete wavelet transform.

5.5 Selection of model parameters

There are four model parameters in the proposed model: $\alpha_{0},\alpha_{\text{curv}},\alpha_{w}$ and $\alpha_{n}$ , and it is generally difficult to determine their optimal values a priori. We demonstrate proper ranges for these model parameters that render good quality for $u^{*}$ in terms of PSNR metric.

We take house, pepper, and blonde in Figure 11 with increasing levels of details. In the first set of experiments, we add Gaussian noise ( $\sigma=10/255$ ) and apply the proposed model to them with different combinations of $\alpha_{\text{curv}}$ and $\alpha_{w}$ while fixing $\alpha_{0}=1\times 10^{-3}$ and $\alpha_{n}=10$ from heuristics. Figure 12 shows the PSNR contours of the $u^{*}$ components learned from the test images. We observe that, for $\alpha_{w}<1$ , the PSNR values are relatively low due to lack of smoothness, and they increase as $\alpha_{w}$ grows. When $1\leq\alpha_{w}\leq 10^{3}$ , PSNR values become relatively insensitive to the choice of $\alpha_{w}$ . As for $\alpha_{\text{curv}}$ , we observe that the optimal values depend on image contents. Specifically, if the images are relatively simple, such as house, $\alpha_{\text{curv}}$ around $0.1\sim 1$ yields slightly better quality. As the images contain richer details such as those in pepper and blonde, the maximal PSNR values are obtained when $\alpha_{\text{curv}}$ takes smaller values around $10^{-3}\sim 10^{-2}$ . We suggest fixing $\alpha_{\text{curv}}$ in the order of $0.1$ and $\alpha_{w}$ in the order of $10$ . In case fine-tuning is needed, a proper range for $\alpha_{\text{curv}}$ is $10^{-3}\sim 1$ , where smaller values are suitable for complex images, and larger values for simple ones. As for $\alpha_{w}$ , keeping it above $1$ is a good practice. Investigations of their effects on the decomposition are collected in Section 5.3.

We conduct our second set of experiments to search for appropriate $\alpha_{n}$ while fixing $\alpha_{\text{curv}}=0.1$ and $\alpha_{w}=50$ . Since $\alpha_{n}$ is directly associated with the noise level, we show in Figure 13 the PSNR values of $u^{*}$ components learned from test images by the proposed model, where both $\alpha_{n}$ and the noise level $\sigma$ vary. We observe that for noise levels below $20$ , the optimal choices of $\alpha_{n}$ stably remain around $0.1\sim 10$ . When the noise level grows above $20$ , we observe that the optimal choice of $\alpha_{n}$ has a clear dependence on $\sigma$ . In particular, smaller $\alpha_{n}$ is needed for larger $\sigma$ to capture oscillations due to noises. We could fit $\alpha_{n}$ in log-scale with $\sigma>20$ using a quadratic polynomial; instead, we suggest a strategy for images with noise $\sigma<100$ with better generalizability. We roughly divide the levels of noise into three regimes: low noise ( $\sigma<20$ ), medium noise ( $20\leq\sigma<60$ ), and high noise ( $60\leq\sigma<100$ ). For images with low noise, use $\alpha_{n}=10$ ; for images with medium noise, use $\alpha_{n}=1\times 10^{-2}$ ; and for those with heavy noise, use $\alpha_{n}=1\times 10^{-4}$ .

With $\alpha_{\text{curv}}$ , $\alpha_{w}$ , and $\alpha_{n}$ determined above, we proceed to tune $\alpha_{0}$ . Using various values of $\alpha_{0}$ , we record in Figure 14 the corresponding PSNR values of the $u^{*}$ components for the three test images with different levels of noise. Specifically, we add Gaussian noises with $\sigma=10/255$ in (a), $\sigma=40/255$ in (b), and $\sigma=80/255$ in (c). As discussed above, we use $\alpha_{n}=10$ for (a), $\alpha_{n}=1\times 10^{-2}$ for (b), and $\alpha_{n}=1\times 10^{-4}$ for (c); and we fix $\alpha_{\text{curv}}=0.1$ and $\alpha_{w}=50$ . For the tested images, when the noise is mild as in (a) and (b), we observe two plateaus in PSNR values of $u^{*}$ components connected by decreasing phases between $1\times 10^{-3}<\alpha_{0}<1\times 10^{-2}$ . The transitions from high PSNR to low PSNR regime are smooth, and they correspond to removing rich image details such as textures. However, when the noise level is high as in (c), we observe that as $\alpha_{0}$ increases, the PSNR values improve slightly. This shows that $L^{0}$ -gradient can contribute to perturbation reduction in case of heavy noise. For general images, we suggest using $\alpha_{0}=1\times 10^{-3}$ as the default.

6 Conclusion

In this paper, we propose a novel model for decomposing grayscale images into three components: the structural part, the smooth part, and the oscillatory part, each of which has a distinctive visual feature. To take advantage of the strong sparsity induced by the $L^{0}$ -gradient regularization while avoiding the undesirable staircase effect, our model integrates the curvature regularization on the level lines of the structural part. As a result, the reconstructed structural part effectively preserves sharp boundaries enclosing homogeneous regions. To solve the associated energy minimization problem, we develop an efficient numerical algorithm based on operator splitting. The original non-convex non-smooth problem is separated into four sub-problems. Each of them either admits a closed-form solution or enjoys a form suitable for applying FFT. We present various experiments to justify the proposed model. In particular, we show that our model successfully separates underlying patterns from soft shadows even when both of them have complicated geometries with fine details. We also justify the regularization terms included in the proposed model and illustrate the effects of parameters. We compare our model with CEP [13] and HKLM [29] for three-component decomposition and show that our model achieves visually appealing results and outstanding computational efficiency. To amend the loss of fine-scale details, a natural extension is to model noise and textures separately. The proposed method can also be extended to color image decomposition when applied channel-wise. To avoid potential artifacts, we leave it as a future work to investigate appropriate regularizers that leverage correlations among color channels.

Acknowledgement

We would like to thank Dr. Martin Huska at University of Bologna for kindly sharing the code of [29]. We would also like to thank Prof. Sung Ha Kang at Georgia Institute of Technology for inspiring discussions on high-order regularizations.

References

[1] L. Ambrosio and S. Masnou. A direct variational approach to a problem arising in image reconstruction. Interfaces and Free Boundaries, 5(1):63–81, 2003.
[2] G. Anzellotti. Pairings between measures and bounded functions and compensated compactness. Annali di Matematica pura ed applicata, 135:293–318, 1983.
[3] G. Aubert and J.-F. Aujol. Modeling very oscillating signals. Application to image processing. Applied Mathematics and Optimization, 51:163–182, 2005.
[4] J.-F. Aujol, G. Aubert, L. Blanc-Féraud, and A. Chambolle. Image decomposition into a bounded variation component and an oscillating component. Journal of Mathematical Imaging and Vision, 22:71–88, 2005.
[5] J.-F. Aujol and A. Chambolle. Dual norms and image decomposition models. International Journal of Computer Vision, 63:85–104, 2005.
[6] J.-F. Aujol, G. Gilboa, T. Chan, and S. Osher. Structure-texture image decomposition—modeling, algorithms, and parameter selection. International Journal of Computer Vision, 67:111–136, 2006.
[7] J.-F. Aujol and S. H. Kang. Color image decomposition and restoration. Journal of Visual Communication and Image Representation, 17(4):916–928, 2006.
[8] G. Bellettini et al. Variational approximation of functionals with curvatures and related properties. Journal of Convex Analysis, 4:91–108, 1997.
[9] P. Blomgren, T. F. Chan, P. Mulet, and C.-K. Wong. Total variation image restoration: numerical methods and extensions. In Proceedings of International Conference on Image Processing, volume 3, pages 384–387. IEEE, 1997.
[10] A. Buades, T. M. Le, J.-M. Morel, and L. A. Vese. Fast cartoon+texture image filters. IEEE Transactions on Image Processing, 19(8):1978–1986, 2010.
[11] A. Chambolle. An algorithm for total variation minimization and applications. Journal of Mathematical Imaging and Vision, 20:89–97, 2004.
[12] T. Chan, S. Esedoglu, F. Park, A. Yip, et al. Recent developments in total variation image restoration. Mathematical Models of Computer Vision, 17(2):17–31, 2005.
[13] T. F. Chan, S. Esedoglu, and F. E. Park. Image decomposition combining staircase reduction and texture extraction. Journal of Visual Communication and Image Representation, 18(6):464–486, 2007.
[14] K. Chen and H. Han. Three-stage approach for 2d/3d diffeomorphic multi-modality image registration with textural control: diffeomorphic multi-modality image registration. SIAM Journal on Imaging Sciences, 2024.
[15] Q. Chen and V. Koltun. A simple model for intrinsic image decomposition with depth cues. In Proceedings of the IEEE International Conference on Computer Vision, pages 241–248, 2013.
[16] E. De Giorgi and L. Ambrosio. New functionals in calculus of variations. In Nonsmooth Optimization and Related Topics, pages 49–59. Springer, 1989.
[17] L.-J. Deng, R. Glowinski, and X.-C. Tai. A new operator splitting method for the Euler elastica model for image smoothing. SIAM Journal on Imaging Sciences, 12(2):1190–1230, 2019.
[18] Q. Ding, T. Niu, X. Zhang, and Y. Long. Image-domain multimaterial decomposition for dual-energy CT based on prior information of material images. Medical physics, 45(8):3614–3626, 2018.
[19] D. L. Donoho and I. M. Johnstone. Adapting to unknown smoothness via wavelet shrinkage. Journal of the American Statistical Association, 90(432):1200–1224, 1995.
[20] V. Duval, J.-F. Aujol, and L. Vese. Projected gradient based color image decomposition. In Scale Space and Variational Methods in Computer Vision: Second International Conference, SSVM 2009, Voss, Norway, June 1-5, 2009. Proceedings 2, pages 295–306. Springer, 2009.
[21] E. Esser. Applications of lagrangian-based alternating direction methods and connections to split Bregman. CAM report, 9:31, 2009.
[22] J. B. Garnett, P. W. Jones, T. M. Le, and L. A. Vese. Modeling oscillatory components with the homogeneous spaces $\dot{\text{BMO}}^{-\alpha}$ and $\dot{W}^{-\alpha,p^{*}}$ . Pure and Applied Mathematics Quarterly, 7(2):275–318, 2011.
[23] G. Gilboa and S. Osher. Nonlocal operators with applications to image processing. Multiscale Modeling & Simulation, 7(3):1005–1028, 2009.
[24] R. Glowinski and P. Le Tallec. Augmented Lagrangian and operator-splitting methods in nonlinear mechanics. SIAM, 1989.
[25] R. Glowinski, S. J. Osher, and W. Yin. Splitting methods in communication, imaging, science, and engineering. Springer, 2017.
[26] R. Glowinski, T.-W. Pan, and X.-C. Tai. Some facts about operator-splitting and alternating direction methods. Splitting Methods in Communication, Imaging, Science, and Engineering, pages 19–94, 2016.
[27] Y. He, S. H. Kang, and H. Liu. Curvature regularized surface reconstruction from point clouds. SIAM Journal on Imaging Sciences, 13(4):1834–1859, 2020.
[28] Y. Hu and M. Jacob. Higher degree total variation (HDTV) regularization for image recovery. IEEE Transactions on Image Processing, 21(5):2559–2571, 2012.
[29] M. Huska, S. H. Kang, A. Lanza, and S. Morigi. A variational approach to additive image decomposition into structure, harmonic, and oscillatory components. SIAM Journal on Imaging Sciences, 14(4):1749–1789, 2021.
[30] Y. Lan, Z. Li, J. Sun, and Y. Xiang. Dosnet as a non-black-box pde solver: When deep learning meets operator splitting. Journal of Computational Physics, 491:112343, 2023.
[31] T. M. Le and L. A. Vese. Image decomposition using total variation and div(BMO). Multiscale Modeling & Simulation, 4(2):390–423, 2005.
[32] N. Le Bihan and S. Sangwine. Color image decomposition using quaternion singular value decomposition. In 2003 International Conference on Visual Information Engineering VIE 2003, pages 113–116. IET, 2003.
[33] L. H. Lieu and L. A. Vese. Image restoration and decomposition via bounded total variation and negative Hilbert-Sobolev spaces. Applied Mathematics and Optimization, 58:167–193, 2008.
[34] H. Liu, J. Liu, R. Chan, and X.-C. Tai. Double-well net for image segmentation. arXiv preprint arXiv:2401.00456, 2023.
[35] H. Liu, X.-C. Tai, and R. Glowinski. An operator-splitting method for the Gaussian curvature regularization model with applications to surface smoothing and imaging. SIAM Journal on Scientific Computing, 44(2):A935–A963, 2022.
[36] H. Liu, X.-C. Tai, R. Kimmel, and R. Glowinski. A color elastica model for vector-valued image regularization. SIAM Journal on Imaging Sciences, 14(2):717–748, 2021.
[37] H. Liu, X.-C. Tai, R. Kimmel, and R. Glowinski. Elastica models for color image regularization. SIAM Journal on Imaging Sciences, 16(1):461–500, 2023.
[38] M. Lysaker, A. Lundervold, and X.-C. Tai. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Transactions on Image Processing, 12(12):1579–1590, 2003.
[39] S. Masnou and J.-M. Morel. Level lines based disocclusion. In Proceedings 1998 International Conference on Image Processing. ICIP98 (Cat. No. 98CB36269), pages 259–263. IEEE, 1998.
[40] S. Masnou and J.-M. Morel. On a variational theory of image amodal completion. Rendiconti del Seminario Matematico della Universita di Padova, 116:211–252, 2006.
[41] Y. Meyer. Oscillating patterns in image processing and nonlinear evolution equations: the fifteenth Dean Jacqueline B. Lewis memorial lectures, volume 22. American Mathematical Soc., 2001.
[42] D. Mumford. Elastica and computer vision. In Algebraic Geometry and its Applications: Collections of Papers from Shreeram S. Abhyankar’s 60th Birthday Conference, pages 491–506. Springer, 1994.
[43] R. M. Nguyen and M. S. Brown. Fast and effective $L_{0}$ gradient minimization by region fusion. In Proceedings of the IEEE International Conference on Computer Vision, pages 208–216, 2015.
[44] S. Ono. $L_{0}$ gradient projection. IEEE Transactions on Image Processing, 26(4):1554–1564, 2017.
[45] S. Osher, A. Solé, and L. Vese. Image decomposition and restoration using total variation minimization and the $H^{-1}$ norm. Multiscale Modeling & Simulation, 1(3):349–370, 2003.
[46] M. T. Pimenta and M. Montenegro. Existence of a BV solution for a mean curvature equation. Journal of Differential Equations, 299:51–64, 2021.
[47] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena, 60(1-4):259–268, 1992.
[48] J. Shen. Piecewise $h^{-1}+h^{0}+h^{1}$ images and the mumford-shah-sobolevmodel for segmented image decomposition. Applied Mathematics Research Express, 2005(4):143–167, 2005.
[49] J. Shen, S. H. Kang, and T. F. Chan. Euler’s elastica and curvature-based inpainting. SIAM Journal on Applied Mathematics, 63(2):564–592, 2003.
[50] X.-C. Tai, J. Hahn, and G. J. Chung. A fast algorithm for Euler’s elastica model using augmented lagrangian method. SIAM Journal on Imaging Sciences, 4(1):313–344, 2011.
[51] X.-C. Tai, H. Liu, and R. Chan. PottsMGNet: A mathematical explanation of encoder-decoder based neural networks. SIAM Journal on Imaging Sciences, 17(1):540–594, 2024.
[52] L. A. Vese and S. J. Osher. Modeling textures with total variation minimization and oscillating patterns in image processing. Journal of Scientific Computing, 19:553–572, 2003.
[53] Y.-W. Wen, H.-W. Sun, and M. K. Ng. A primal-dual method for the meyer model of cartoon and texture decomposition. Numerical Linear Algebra with Applications, 26(2):e2224, 2019.
[54] Y.-W. Wen, M. Zhao, and M. Ng. Cartoon and texture decomposition for color image in opponent color space. Applied Mathematics and Computation, 414:126654, 2022.
[55] X. Wu, J. Zheng, C. Wu, and Y. Cai. Variational structure–texture image decomposition on manifolds. Signal Processing, 93(7):1773–1784, 2013.
[56] J. Xu, X. Feng, Y. Hao, and Y. Han. Image decomposition and staircase effect reduction based on total generalized variation. Journal of Systems Engineering and Electronics, 25(1):168–174, 2014.
[57] L. Xu, C. Lu, Y. Xu, and J. Jia. Image smoothing via $L_{0}$ gradient minimization. In Proceedings of the 2011 SIGGRAPH Asia Conference, pages 1–12, 2011.
[58] M. Yashtini and S. H. Kang. Alternating direction method of multiplier for Euler’s elastica-based denoising. In Scale Space and Variational Methods in Computer Vision: 5th International Conference, SSVM 2015, Lège-Cap Ferret, France, May 31-June 4, 2015, Proceedings 5, pages 690–701. Springer, 2015.
[59] C.-H. Zhang. Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, 38:894–942, 2010.
[60] H. Zhang, L. Tang, Z. Fang, C. Xiang, and C. Li. Nonconvex and nonsmooth total generalized variation model for image restoration. Signal Processing, 143:69–85, 2018.
[61] W. Zhu. A first-order image denoising model for staircase reduction. Advances in Computational Mathematics, 45(5):3217–3239, 2019.
[62] W. Zhu and T. Chan. Image denoising using mean curvature of image surface. SIAM Journal on Imaging Sciences, 5(1):1–32, 2012.

Appendix A Proof of Proposition 2.1

Proof.

The proof is analogous to [17]. Denote

	$\displaystyle E_{1}(v,w,n)=\alpha_{0}\\|\nabla v\\|_{0}+\alpha_{\text{curv}}\int_{\Omega}\left(\nabla\cdot\frac{\nabla v}{\|\nabla v\|}\right)^{2}\,\left\|Dv\right\|\,d\mathbf{x}+\alpha_{w}\\|\Delta w\\|_{L^{2}(\Omega)}^{2}+\alpha_{n}\\|n\\|_{H^{-1}(\Omega)}^{2},$
	$\displaystyle E_{2}(v,w,n)=\frac{1}{2}\int_{\Omega}\|f-(v+w+n)\|_{2}^{2}d\mathbf{x}.$

Let $(v^{*},w^{*},n^{*})$ be a minimizer of Model (8). Consider the tuple $(v^{*}+c_{1},w^{*}+c_{2},n^{*})$ , where $c_{1},c_{2}\in\mathbb{R}$ are scalar fields. Notice that

\displaystyle E_{1}(v^{*},w^{*},n^{*})=E_{1}(v^{*}+c_{1},w^{*}+c_{2},n^{*})\;,

and

	$\displaystyle E_{2}(v^{}+c_{1},w^{}+c_{2},n^{*})=$	$\displaystyle\frac{1}{2}\int_{\Omega}\|f-(v^{}+c_{1}+w^{}+c_{2}+n^{*})\|^{2}d\mathbf{x}$
	$\displaystyle=$	$\displaystyle E_{2}(v^{},w^{},n^{})+(c_{1}+c_{2})\int_{\Omega}(v^{}+w^{}+n^{}-f)d\mathbf{x}+\frac{\|\Omega\|}{2}(c_{1}+c_{2})^{2},$		(61)

where $|\Omega|$ is the area of $\Omega$ . The function on the right-hand-side of (61) is quadratic in $(c_{1}+c_{2})$ , which takes its minimal value when

\displaystyle c_{1}+c_{2}=\frac{\int_{\Omega}(f-(v^{*}+w^{*}+n^{*}))d\mathbf{x}}{|\Omega|}.

(62)

If $\int_{\Omega}(f-(v^{*}+w^{*}+n^{*}))d\mathbf{x}\neq 0$ , it would imply that

\displaystyle E_{1}(v^{*},w^{*},n^{*})+E_{2}(v^{*},w^{*},n^{*})>E_{1}(u^{*}+c_{1},w^{*}+c_{2},n^{*})+E_{2}(u^{*}+c_{1},w^{*}+c_{2},n^{*}),

which contradicts to the assumption that $(v^{*},w^{*},n^{*})$ is a minimizer of Model (8). We thus conclude that $\int_{\Omega}v^{*}+w^{*}+n^{*}d\mathbf{x}=\int_{\Omega}fd\mathbf{x}$ . ∎

Appendix B Equivalence between Hessian regularizer and Laplacian regularizer

We compare the following two regularizers

\displaystyle E_{H}(w)=\frac{1}{2}\|\partial^{2}w\|_{L^{2}(\Omega)}^{2}\quad\mbox{ and }\quad E_{L}(w)=\frac{1}{2}\|\Delta w\|_{L^{2}(\Omega)}^{2}.

Functional $E_{H}$ is used in [29] and $E_{L}$ is used in this paper as a regularizer for the $w$ component. It was observed in [13] that $E_{H}=E_{L}$ for compactly supported smooth functions. We show that more general boundary conditions, $E_{H}$ and $E_{L}$ have the same set of critical points. To prove that, it is sufficient to show that their variations concerning $w$ have the same form. For simplicity, we assume $\Omega=[0,L_{1}]\times[0,L_{2}]$ for some $L_{1},L_{2}>0$ .

The functional $E_{H}$ can be written as

\displaystyle E_{H}(w)=\frac{1}{2}\int_{\Omega}\left(\frac{\partial^{2}w}{\partial x_{1}^{2}}\right)^{2}+2\left(\frac{\partial^{2}w}{\partial x_{1}\partial x_{2}}\right)^{2}+\left(\frac{\partial^{2}w}{\partial x_{2}^{2}}\right)^{2}d\mathbf{x}.

The Fréchet derivative of $E_{H}$ with respect to $w$ in the direction $h$ is

$\displaystyle\left\langle\frac{\partial E_{H}}{\partial w},h\right\rangle=$	$\displaystyle\int_{\Omega}\left(\frac{\partial^{4}w}{\partial x_{1}^{4}}+2\frac{\partial^{4}w}{\partial x_{1}^{2}\partial x_{2}^{2}}+\frac{\partial^{4}w}{\partial x_{2}^{4}}\right)h\ d\mathbf{x}$
	$\displaystyle-\left.\left\langle\frac{\partial^{2}w}{\partial x_{1}^{2}},\frac{\partial h}{\partial x_{1}}\right\rangle\right\|_{x_{1}=0}^{x_{1}=L_{1}}+\left.\left\langle\frac{\partial^{3}w}{\partial x_{1}^{3}},h\right\rangle\right\|_{x_{1}=0}^{x_{1}=L_{1}}-\left.\left\langle\frac{\partial^{2}w}{\partial x_{2}^{2}},\frac{\partial h}{\partial x_{2}}\right\rangle\right\|_{x_{2}=0}^{x_{2}=L_{2}}+\left.\left\langle\frac{\partial^{3}w}{\partial x_{2}^{3}},h\right\rangle\right\|_{x_{2}=0}^{x_{2}=L_{2}}$
	$\displaystyle-\left.\left\langle\frac{\partial^{2}w}{\partial x_{1}\partial x_{2}},\frac{\partial h}{\partial x_{2}}\right\rangle\right\|_{x_{1}=0}^{x_{1}=L_{1}}+\left.\left\langle\frac{\partial^{3}w}{\partial x_{1}^{2}\partial x_{2}},h\right\rangle\right\|_{x_{2}=0}^{x_{2}=L_{2}}$
	$\displaystyle-\left.\left\langle\frac{\partial^{2}w}{\partial x_{1}\partial x_{2}},\frac{\partial h}{\partial x_{1}}\right\rangle\right\|_{x_{2}=0}^{x_{2}=L_{2}}+\left.\left\langle\frac{\partial^{3}w}{\partial x_{1}\partial x_{2}^{2}},h\right\rangle\right\|_{x_{1}=0}^{x_{1}=L_{1}}.$	(63)

The functional $E_{L}$ can be written as

\displaystyle E_{L}(w)=\frac{1}{2}\int_{\Omega}\left(\frac{\partial^{2}w}{\partial x_{1}^{2}}\right)^{2}+2\left(\frac{\partial^{2}w}{\partial x_{1}^{2}}\frac{\partial^{2}w}{\partial x_{2}^{2}}\right)+\left(\frac{\partial^{2}w}{\partial x_{2}^{2}}\right)^{2}d\mathbf{x}.

The Fréchet derivative of $E_{L}$ with respect to $w$ in the direction $h$ is

$\displaystyle\left\langle\frac{\partial E_{L}}{\partial w},h\right\rangle=$	$\displaystyle\int_{\Omega}\left(\frac{\partial^{4}w}{\partial x_{1}^{4}}+2\frac{\partial^{4}w}{\partial x_{1}^{2}\partial x_{2}^{2}}+\frac{\partial^{4}w}{\partial x_{2}^{4}}\right)h\ d\mathbf{x}$
	$\displaystyle-\left.\left\langle\left(\frac{\partial^{2}w}{\partial x_{1}^{2}}+\frac{\partial^{2}w}{\partial x_{2}^{2}}\right),\frac{\partial h}{\partial x_{1}}\right\rangle\right\|_{x_{1}=0}^{x_{1}=L_{1}}+\left.\left\langle\left(\frac{\partial^{3}w}{\partial x_{1}^{3}}+\frac{\partial^{3}w}{\partial x_{1}^{2}\partial x_{2}}\right),h\right\rangle\right\|_{x_{1}=0}^{x_{1}=L_{1}}$
	$\displaystyle-\left.\left\langle\left(\frac{\partial^{2}w}{\partial x_{2}^{2}}+\frac{\partial^{2}w}{\partial x_{2}^{2}}\right),\frac{\partial h}{\partial x_{2}}\right\rangle\right\|_{x_{2}=0}^{x_{2}=L_{2}}+\left.\left\langle\left(\frac{\partial^{3}w}{\partial x_{2}^{3}}+\frac{\partial^{3}w}{\partial x_{1}\partial x_{2}^{2}}\right),h\right\rangle\right\|_{x_{2}=0}^{x_{2}=L_{2}}.$	(64)

By using the boundary conditon

\displaystyle\frac{\partial^{2}w}{\partial x_{1}^{2}}=\frac{\partial^{2}w}{\partial x_{2}^{2}}=\frac{\partial^{3}w}{\partial x_{1}^{3}}=\frac{\partial^{3}w}{\partial x_{2}^{3}}=\frac{\partial^{3}w}{\partial x_{1}^{2}\partial x_{2}}=\frac{\partial^{3}w}{\partial x_{1}\partial x_{2}^{2}}=0\mbox{ on }\partial\Omega,

all boundary integrals in (63) and (64) vanish. We get the same Fréchet derivative for both $E_{H}$ and $E_{L}$

\displaystyle\left\langle\frac{\partial E_{H}}{\partial w},h\right\rangle=\left\langle\frac{\partial E_{L}}{\partial w},h\right\rangle=\int_{\Omega}\left(\frac{\partial^{4}w}{\partial x_{1}^{4}}+2\frac{\partial^{4}w}{\partial x_{1}^{2}\partial x_{2}^{2}}+\frac{\partial^{4}w}{\partial x_{2}^{4}}\right)h\ d\mathbf{x}.

Therefore, $E_{H}$ and $E_{L}$ have the same optimality condition and thus the same set of critical points.

	$\displaystyle J_{0}(\mathbf{q},w,\mathbf{g})=\frac{1}{2}\int_{\Omega}\|f-(v_{\mathbf{q},w,\mathbf{g}}+w+\nabla\cdot\mathbf{g})\|^{2}d\mathbf{x},$
	$\displaystyle J_{1}(\mathbf{q})=\alpha_{0}\\|\mathbf{q}\\|_{0},$
	$\displaystyle J_{2}(\mathbf{q},\boldsymbol{\mu})=\alpha_{\rm curv}\int_{\Omega}\left\|\nabla\cdot\boldsymbol{\mu}\right\|^{2}\left\|\mathbf{q}\right\|\,d\mathbf{x},$
	$\displaystyle J_{3}(w)=\alpha_{w}\int_{\Omega}\|\nabla^{2}w\|^{2}d\mathbf{x},$
	$\displaystyle J_{4}(\mathbf{g})=\alpha_{n}\int_{\Omega}\|\mathbf{g}\|^{2}d\mathbf{x},$

$f_{\text{noise}}$	$u^{*}$	$n^{*}_{\text{scaled}}$

PSNR $=22.10$	PSNR $=43.60$	STD $=19.69/255$

PSNR $=22.10$	PSNR $=35.87$	STD $=19.97/255$

PSNR $=22.12$	PSNR $=36.44$	STD $=19.82/255$

PSNR $=22.12$	PSNR $=25.69$	STD $=23.19/255$

$f_{\text{noise}}$	$u^{*}$	$n^{*}_{\text{scaled}}$

PSNR $=22.11$	PSNR $=28.50$	STD $=19.87/255$

PSNR $=22.11$	PSNR $=30.57$	STD $=19.24/255$

$\alpha_{w}$	10	50	100	500	1000
$v_{\text{scaled}}^{*}$
$w_{\text{scaled}}^{*}$
$u^{*}$
Zoom-in $u^{*}$
PSNR( $u^{*}$ )	$34.58$	34.59	$34.53$	$34.34$	$34.28$

$\alpha_{\text{curv}}$	0.01	0.1	1	5	10
$v_{\text{scaled}}^{*}$
$w_{\text{scaled}}^{*}$
$u^{*}$
Zoom-in $u^{*}$
PSNR( $u^{*}$ )	$33.06$	$32.77$	$32.05$	$31.19$	$30.67$

(a)	(b)	(c)	(d)	(e)	(f)
Input	Model I	Model II	Model III	Proposed	True



PSNR $(u^{*})$	$32.31$	$34.64$	$35.03$	$\mathbf{35.66}$

Euler’s Elastica Based Cartoon-Smooth-Texture Image Decomposition

Abstract

1 Introduction

2 Proposed Cartoon-Smooth-Texture Decomposition Model

2.1 Proposed model

Remark 1.

Remark 2.

Remark 3.

2.2 Model reformulation and formal optimality condition

Proposition 2.1.

3 Operator-Splitting Algorithm for the Proposed Model

3.1 Proposed scheme

3.2 On the solution to sub-problem (15) and (16)

3.3 On the solution to sub-problem (17)

3.4 On the solution to sub-problem (18)

3.5 Adapt to periodic boundary conditions for efficiency

4 Numerical Discretization

4.1 Synopsis

4.2 Discrete analogue of 𝐩k+1/4\mathbf{p}^{k+1/4}, 𝐩k+2/4\mathbf{p}^{k+2/4} and 𝝀k+2/4\boldsymbol{\lambda}^{k+2/4}

4.3 Discrete analogue of 𝐩k+3/4\mathbf{p}^{k+3/4}, 𝝀k+3/4\boldsymbol{\lambda}^{k+3/4}

4.4 Discrete analogue of 𝐩k+1,rk+1\mathbf{p}^{k+1},r^{k+1} and 𝐬k+1\mathbf{s}^{k+1}

5 Numerical Experiments

5.1 General performances

5.2 Ablation study

5.3 Effect of parameter αw\alpha_{w} and αcurv\alpha_{\text{curv}}

5.4 Comparison on cartoon+smooth+oscillation decomposition

5.5 Selection of model parameters

6 Conclusion

Acknowledgement

References

Appendix A Proof of Proposition 2.1

Proof.

Appendix B Equivalence between Hessian regularizer and Laplacian regularizer

4.2 Discrete analogue of $\mathbf{p}^{k+1/4}$ , $\mathbf{p}^{k+2/4}$ and $\boldsymbol{\lambda}^{k+2/4}$

4.3 Discrete analogue of $\mathbf{p}^{k+3/4}$ , $\boldsymbol{\lambda}^{k+3/4}$

4.4 Discrete analogue of $\mathbf{p}^{k+1},r^{k+1}$ and $\mathbf{s}^{k+1}$

5.3 Effect of parameter $\alpha_{w}$ and $\alpha_{\text{curv}}$