# Adaptive Image Restoration via a Relaxed Regularization of Mean Curvature

Authors: Mingxi Ma; Jun Zhang; Chengzhi Deng; Zhaoyang Liu; Yuanyun Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, a new relaxation model based on mean curvature for adaptive image restoration is proposed. To solve the problem efficiently, an alternating direction method of multipliers (ADMMs) is proposed. Furthermore, a rigorous convergence theory of the proposed algorithm is established. We also give the complexity analysis of our proposed method. Experimental results are provided to demonstrate the effectiveness and efficiency of the proposed method over a state-of-the-art method on synthetic and natural images.


---

## Body

## 1. Introduction

Image restoration has always been a research hotspot in the field of image processing. To solve this inverse problem, one of the most popular research directions is applying variational regularization methods. It mainly includes total variation (TV) regularization [1–3], nonlocal regularization [4], sparse regularization [5], higher-order regularization based on higher-order derivatives [6–9], and fractional-order regularization based on fractional-order derivatives [10, 11].As a classic regularization method, TV regularization has been successfully applied to image restoration, image segmentation, image reconstruction, image decomposition, and phase recovery. When using it for image denoising, the edge-preserving effect is better. For the TV-based image denoising problem, researchers have proposed a variety of fast and efficient numerical algorithms, such as dual algorithm [12], primal-dual algorithm [13], split Bregman method [14], augmented Lagrangian method [15, 16], alternating direction method of multipliers (ADMMs) [17], proximity algorithm [18], domain decomposition method [19], and nonlinear multigrid method [3]. However, after denoising with TV regularization, the image often exhibits staircase effects in the flat area. Second, this type of regularization model will cause image contrast loss.To reduce the staircase effects, a number of higher-order models have been developed [20, 21]. Among them, let us specially note the mean curvature model [9, 22] and Euler’s elastica model [23–27], which are nonconvex so that the numerical algorithms converge to the local optimal solution. The mean curvature-based image denoising model (called the MC model) can be formulated as follows:(1)minu∫Ω∇·∇u1+∇u2dxdy+12λ∫Ωu−f2dxdy,where u and f represent the restored image and the noisy image, ∇ is the gradient operator, the weighting parameter λ>0 controls the amount of denoising, and ∇·∇u/1+∇u2 is the mean curvature of the image surface ϕx,y,z=ux,y−z=0. As mentioned in [28], the MC model preserves image contrast better for cleaned images than the TV model. In addition, this model handles corners better for synthetic images than the TV and Euler’s elastica models. In view of a number of valuable properties of the MC model, designing fast iterative algorithms for its solution has important theoretical and practical significance. There are some different augmented Lagrangian methods for solving the MC model. In [28], original unconstrained problem (1) was transformed into a constrained optimization problem and the saddle-points of the corresponding augmented Lagrangian functional are obtained from solving several subproblems alternatively. Some subproblems had closed form solutions while others were solved by fast Fourier transform (FFT). Subsequently, in [29], the solutions of two subproblems were obtained by a Gauss–Seidel method. In addition to the reduction of computational cost, this approach can be applied to computational domains with nonperiodic boundaries. Besides, a fast linearized augmented Lagrangian method was proposed to further improve the convergence rates [30]. In this method, all subproblems have closed form solutions. However, due to the nonconvexity of the MC model, it is still very hard to find suitable model parameter and penalization parameters.By choosing the regularization parameter based on the inverse gradient, in [31], Thanh et al. proposed some adaptive image restoration methods. Besides, in [32, 33], the authors proposed the adaptive TVp regularization and the adaptive weighted TVp regularization for image denoising, respectively. In [34], a novel adaptive image denoising method was proposed, which combines a Tchebichef moment-based sparse regularizer with an adaptive steerable total variation regularizer.By minimizing theL1 norm of the gradient of the unit normal, Duan et al. proposed a novel variational model as follows [8]:(2)minu∫Ω∇∇u1+∇u2dxdy+12λ∫Ωu−f2dxdy,and proved this model with the ability to preserve image contrast and sharp edges during denoising. Due to the high nonlinearity and nonconvexity of its regularization term, the relaxation form of model (2) was studied. Namely, they proposed a spatially adaptive hybrid regularization-based minimization problem, which can be given as(3)minu∫Ωαx,y∇u+βx,y∇2udxdy+12λ∫Ωu−f2dxdy,where operator ∇2 represents the second-order gradient, αx,y=∇1/1+∇u2 and βx,y=1/1+∇u2. For solving this relaxation model, the ADMM-based numerical algorithm was presented and its convergence was proved [8, 35]. ADMM is a fast method for solving many regularization problems. It has been widely applied in many fields, such as hyperspectral image unmixing [36], classification [37], and fusion [38].Inspired by the successful application of the relaxation technique, in this paper, we present a new form for the MC model (1) and propose its relaxation model. Moreover, we propose an efficient ADMM with guaranteed convergence. Compared with the ADMM-based algorithm for solving model (3), our method is faster.The organization of the paper is as follows: In Section2, the new form of the MC model (1) is presented, and then its geometric properties and its relaxation form together with connections to the existing models are introduced. Section 3 details an efficient ADMM for solving the proposed relaxation model. The convergence and complexity of the proposed algorithm are also studied in this section. Experimental results on synthetic and natural images are presented in Section 4 to show the effectiveness and efficiency of our proposed method. Finally, Section 5 concludes the paper.


## 2. Adaptive TV and Laplacian Regularization Model

In this section, we first give a new form for the MC model (1). Then, we introduce its geometric properties for image denoising. Finally, we propose a relaxation version using a spatially adaptive TV and Laplacian regularization and analyze the connections to the existing models.


### 2.1. New Form for the MC Model

Let us first review the MC model.For a given imagef:Ω⟶ℝ where Ω⊂ℝ2 is the image domain, the corresponding image surface is x,y,fx,y or z=fx,y. The aim of denoising is to find a piecewise smooth surface z=ux,y to approximate it. Meanwhile, we would like to preserve the geometric features, such as edges of objects and corners. We also expect to preserve image contrasts. To this end, the L1 norm of mean curvature of the image surface was used in [9]. In the following, to introduce the mean curvature of surface z=ux,y, a level set function ϕx,y,z=ux,y−z is considered. Obviously, its zero level set corresponds to the surface. Then, the mean curvature of this surface can be defined as(4)Hu=∇·∇ϕ∇ϕ=∇·∇u1+∇u2.Using theL1 norm of mean curvature as a regularization term, the MC model was proposed as follows:(5)minu∫Ω∇·∇u1+∇u2dxdy+12λ∫Ωu−f2dxdy.Since(6)∇·∇u1+∇u2=∇u·∇11+∇u2+11+∇u2Δu,we can obtain a new form of the MC model which is(7)minu∫Ω∇u·∇11+∇u2+11+∇u2Δudxdy+12λ∫Ωu−f2dxdy.Because this model is nonconvex, high order, and highly nonlinear, its complete study is nontrivial. However, one can get some insight into its features.


### 2.2. Features of Model (7)

In this subsection, the features of model (7) will be introduced. One can refer to [9] for more details.If∇u≪1 and Δu≪1 for each point in a region, the evolution equation corresponding to the gradient descent flow of model (7) behaves like the biharmonic heat equation in that region. Owing to the property that the biharmonic heat equation can damp high frequency signals quickly, model (7) can remove noise efficiently.In order to illustrate why model (7) can preserve image contrasts, one may consider a simple image f=hχB0,Rx,y defined on Ω=−2R,2R×−2R,2R with χ being the characteristic function, B0,R being an open disk in ℝ2 centered at the origin with radius R, and h>0. Since f is radial symmetric, and one can approximate it by a sequence of smooth radial symmetric function gn. Referring to [9], the properties of preserving image contrasts and sharp edges about model (7) are guaranteed by the following lemma and theorem:Lemma 1.
Assuming thatf=hχB0,Rx,y is an image defined as above, we have(8)∫ΩHfdxdy=4πR,where Hf=∇f·∇1/1+∇f2+1+∇f2Δf.The integral indicates that the regularizer of model (7) does not rely on heights of signals. With the conclusion of this lemma, one can further obtain the following theorem:Theorem 1.
Letf=hχB0,Rx,y be an image defined as above. Then, there exists a constant C>0 such that if λ<C, f is a minimizer of model (7).This theorem demonstrates that model (7) is able to preserve image contrasts, once λ is chosen to be small enough. In contrast, a loss of image contrast may appear while using the ROF model, for whatever small λ. Besides, this theorem illustrates that the model can keep sharp edges as the ROF model.To illustrate why the model can preserve object corners, one needs to consider a different imagef=hχE defined on Ω=−R,R×−R,R with E=0,R×0,R [9]. Similar results are presented in the following lemma and theorem:Lemma 2.
Letf=hχE be an image defined on a rectangle Ω=−R,R×−R,R with E=0,R×0,R. Then,(9)∫ΩHfdxdy=4R.Theorem 2.
Letf=hχE be defined as the above lemma, then there exists a constant C>0 such that if λ<C, f is a minimizer of model (7).Theorem2 indicates another key feature of model (7) on preserving object edges once λ is small enough.


### 2.3. Proposed Relaxation Model and Connections to the Previous Models

Although the new form of the MC model has good properties in preserving image contrasts, sharp edges, and object corners, it is still nonconvex and highly nonlinear. Motivated by the success of relaxation in [8], by relaxing the regularization term, we reformulate (7) into a spatially adaptive TV and Laplacian regularization model, which can be expressed as(10)minu∫Ωαx,y∇u+βx,yΔudxdy+12λ∫Ωu−f2dxdy,where(11)αx,y=∇11+∇u2,βx,y=11+∇u2.It is straightforward to see that model (10) is just the relaxation version of the MC model (1). In [39], an effective model which combines the total variation with the Laplacian regularizer (TVL) was introduced:(12)minuμ1∇u1+μ2Δu1+12KAu−d22,where K and A represent the blurring operator and the Abel transform, respectively. μ1 and μ2 are two constants. The regularization of our model (10) can be seen as an adaptive form of this TVL regularization. Replacing Δ in (10) by ∇2, it becomes the spatially adaptive first- and second-order regularization model (3) proposed in [8]. Compared with the isotropic model, the corresponding anisotropic model is faster when it achieves almost the same denoising effect. Therefore, in order to improve the computational complexity, it is reasonable to use the Laplacian operator Δ to replace the isotropic second-order gradient operator ∇2.


## 2.1. New Form for the MC Model

Let us first review the MC model.For a given imagef:Ω⟶ℝ where Ω⊂ℝ2 is the image domain, the corresponding image surface is x,y,fx,y or z=fx,y. The aim of denoising is to find a piecewise smooth surface z=ux,y to approximate it. Meanwhile, we would like to preserve the geometric features, such as edges of objects and corners. We also expect to preserve image contrasts. To this end, the L1 norm of mean curvature of the image surface was used in [9]. In the following, to introduce the mean curvature of surface z=ux,y, a level set function ϕx,y,z=ux,y−z is considered. Obviously, its zero level set corresponds to the surface. Then, the mean curvature of this surface can be defined as(4)Hu=∇·∇ϕ∇ϕ=∇·∇u1+∇u2.Using theL1 norm of mean curvature as a regularization term, the MC model was proposed as follows:(5)minu∫Ω∇·∇u1+∇u2dxdy+12λ∫Ωu−f2dxdy.Since(6)∇·∇u1+∇u2=∇u·∇11+∇u2+11+∇u2Δu,we can obtain a new form of the MC model which is(7)minu∫Ω∇u·∇11+∇u2+11+∇u2Δudxdy+12λ∫Ωu−f2dxdy.Because this model is nonconvex, high order, and highly nonlinear, its complete study is nontrivial. However, one can get some insight into its features.


## 2.2. Features of Model (7)

In this subsection, the features of model (7) will be introduced. One can refer to [9] for more details.If∇u≪1 and Δu≪1 for each point in a region, the evolution equation corresponding to the gradient descent flow of model (7) behaves like the biharmonic heat equation in that region. Owing to the property that the biharmonic heat equation can damp high frequency signals quickly, model (7) can remove noise efficiently.In order to illustrate why model (7) can preserve image contrasts, one may consider a simple image f=hχB0,Rx,y defined on Ω=−2R,2R×−2R,2R with χ being the characteristic function, B0,R being an open disk in ℝ2 centered at the origin with radius R, and h>0. Since f is radial symmetric, and one can approximate it by a sequence of smooth radial symmetric function gn. Referring to [9], the properties of preserving image contrasts and sharp edges about model (7) are guaranteed by the following lemma and theorem:Lemma 1.
Assuming thatf=hχB0,Rx,y is an image defined as above, we have(8)∫ΩHfdxdy=4πR,where Hf=∇f·∇1/1+∇f2+1+∇f2Δf.The integral indicates that the regularizer of model (7) does not rely on heights of signals. With the conclusion of this lemma, one can further obtain the following theorem:Theorem 1.
Letf=hχB0,Rx,y be an image defined as above. Then, there exists a constant C>0 such that if λ<C, f is a minimizer of model (7).This theorem demonstrates that model (7) is able to preserve image contrasts, once λ is chosen to be small enough. In contrast, a loss of image contrast may appear while using the ROF model, for whatever small λ. Besides, this theorem illustrates that the model can keep sharp edges as the ROF model.To illustrate why the model can preserve object corners, one needs to consider a different imagef=hχE defined on Ω=−R,R×−R,R with E=0,R×0,R [9]. Similar results are presented in the following lemma and theorem:Lemma 2.
Letf=hχE be an image defined on a rectangle Ω=−R,R×−R,R with E=0,R×0,R. Then,(9)∫ΩHfdxdy=4R.Theorem 2.
Letf=hχE be defined as the above lemma, then there exists a constant C>0 such that if λ<C, f is a minimizer of model (7).Theorem2 indicates another key feature of model (7) on preserving object edges once λ is small enough.


## 2.3. Proposed Relaxation Model and Connections to the Previous Models

Although the new form of the MC model has good properties in preserving image contrasts, sharp edges, and object corners, it is still nonconvex and highly nonlinear. Motivated by the success of relaxation in [8], by relaxing the regularization term, we reformulate (7) into a spatially adaptive TV and Laplacian regularization model, which can be expressed as(10)minu∫Ωαx,y∇u+βx,yΔudxdy+12λ∫Ωu−f2dxdy,where(11)αx,y=∇11+∇u2,βx,y=11+∇u2.It is straightforward to see that model (10) is just the relaxation version of the MC model (1). In [39], an effective model which combines the total variation with the Laplacian regularizer (TVL) was introduced:(12)minuμ1∇u1+μ2Δu1+12KAu−d22,where K and A represent the blurring operator and the Abel transform, respectively. μ1 and μ2 are two constants. The regularization of our model (10) can be seen as an adaptive form of this TVL regularization. Replacing Δ in (10) by ∇2, it becomes the spatially adaptive first- and second-order regularization model (3) proposed in [8]. Compared with the isotropic model, the corresponding anisotropic model is faster when it achieves almost the same denoising effect. Therefore, in order to improve the computational complexity, it is reasonable to use the Laplacian operator Δ to replace the isotropic second-order gradient operator ∇2.


## 3. The Proposed Algorithm

In this section, we propose an efficient numerical algorithm based on the popular ADMM for solving the proposed model (10). Besides, we analyze its convergence and complexity.


### 3.1. Algorithm Description

Firstly, we introduce two auxiliary variablesp and q to reformulate (10) as follows:(13)minu,p,q∫Ωαx,yp+βx,yqdxdy+12λ∫Ωu−f2dxdy,s.t.p=∇u,q=Δu.The augmented Lagrangian functional corresponding to (13) is(14)ℒu,p,q;λ1,λ2=∫Ωαx,yp+βx,yqdxdy+12λ∫Ωu−f2dxdy,+∫Ωλ1·p−∇udxdy+r12∫Ωp−∇u2dxdy,+∫Ωλ2q−Δudxdy+r22∫Ωq−Δu2dxdy,where λ1andλ2 are Lagrangian multipliers and r1andr2 are positive penalty parameters. By applying the ADMM, one needs to minimize three subproblems with other variables fixed and update two Lagrangian multipliers at each iteration which is(15)uk+1=argminu12λ∫Ωu−f2dxdy+r12∫Ωpk−∇u+λ1kr12dxdy+r22∫Ωqk−Δu+λ2kr22dxdy+δ12∫Ωu−uk2dxdy,pk+1=argminp∫Ωαk+1x,ypdxdy+r12∫Ωp−∇uk+1+λ1kr12dxdy+δ22∫Ωp−pk2dxdy,qk+1=argminq∫Ωβk+1x,yqdxdy+r22∫Ωq−Δuk+1+λ2kr22dxdy+δ32∫Ωq−qk2dxdy,λ1k+1=λ1k+r1pk+1−∇uk+1,λ2k+1=λ2k+r2qk+1−Δuk+1,where δ1,δ2,δ3>0 are three penalty parameters and αk+1x,yandβk+1x,y can be explicitly updated as(16)αk+1x,y=∇11+∇uk+12,βk+1x,y=11+∇uk+12.To solve theu-subproblem in (15), we consider its Euler–Lagrange equation, which has the form(17)u−fλ+r1divpk−∇u+λ1kr1−r2Δqk−Δu+λ2kr2+δ1u−uk=0,namely,(18)1λ−r1div∇+r2ΔΔ+δ1u=fλ−divr1pk+λ1k+Δr2qk+λ2k+δ1uk.In order to efficiently solve this nonlinear equation, we use the fast Fourier transform (FFT) with periodic boundary conditions. More concretely,uk+1 is given by(19)uk+1=ℱ−1ℱf/λ−divr1pk+λ1k+Δr2qk+λ2k+δ1uk1/λ−r1ℱdiv∇+r2ℱΔΔ+δ1,where ℱ is the FFT.Bothp and q subproblems have closed form solutions. Specifically, the solutions of p-subproblem and q-subproblem are given by the shrinkage operator:(20)pk+1=shrinkager1∇uk+1−λ1k/r1+δ2pkr1+δ2,αk+1x,yr1+δ2,and(21)qk+1=shrinkager2Δuk+1−λ2k/r2+δ3qkr2+δ3,βk+1x,yr2+δ3,where shrinkages,γ=maxs−γ,0°sgns. Here, ° denotes the component-wise multiplication.We summarize the overall algorithm for solving the spatially adaptive relaxation model (10) in Algorithm 1.Algorithm 1: (ADMM for solving (10)).
Step 1. Input λ,r1,andr2. Initialization: u0=f,p0=0,q0=0, λ10=0,andλ20=0. Let k≔0.Step 2. Compute uk+1 by (19).Step 3. Update αk+1x,y,βk+1x,y using uk+1 according to (16).Step 4. Compute pk+1 by (20).Step 5. Compute qk+1 by (21).Step 6. Update λ1k+1,λ2k+1 by the last two equations of (15).Step 7. Stop if the stopping criterion is satisfied. Otherwise, let k≔k+1, and go to Step 2.


### 3.2. Convergence Analysis

In this subsection, we give a convergence theory for the proposed algorithm (Algorithm1) under the assumption that a solution of the proposed model (10) exists.Theorem 3.
Assume thatuk,pk,qk,λ1k,λ2kk∈N is a sequence generated by Algorithm 1 and u¯,p¯,q¯,λ¯1,andλ¯2k are the points satisfying the following first-order optimality conditions:(22)u−fλ+∇·λ1−Δλ2=0,αx,yg+λ1=0,g∈∂pandαx,y=∇11+∇u2,βx,yh+λ2=0,h∈∂qandβx,y=11+∇u2,p−∇u=0,q−Δu=0.Letuek=uk−u¯,pek=pk−p¯,qek=qk−q¯,λ1ek=λ1k−λ¯1,andλ2ek=λ2k−λ¯2 be the error sequences. If for any gk∈∂pk,g¯∈∂p¯, and hk∈∂qk,h¯∈∂q¯,(23)Δ1k≔αkx,y−α¯x,ygk,pk−p¯,α¯x,y=∇11+∇u¯2,Δ2k≔βkx,y−β¯x,yhk,qk−q¯,β¯x,y=11+∇u¯2,are nonnegative. Then, the generated sequence uk,pk,qk,λ1k,λ2kk∈N converges to a limit point u∗,p∗,q∗,λ1∗,λ2∗ which satisfies the first-order optimality conditions (22).Our proposed model (10) can be seen as a variation of the previous model in [8]. Specifically, the operator ∇2 in (3) is changed to Δ. Furthermore, an ADMM-based algorithm same to [8] is applied to solve the proposed model (10). Therefore, the rigorous convergence proof of Algorithm 1 can refer to the proof presented in [8].


### 3.3. Complexity Analysis

In the following, we will give some simple analysis of complexity for the proposed method.For our proposed method, the cost of updatinguk is 4NlogN+26N+1 at each iteration, where N=mn is the total number of grid points. If a root operation is equal to 12 times addition or subtraction or multiplication or division, at each iteration it will spend 61N (3N times root operations), 11N, 6N, and 3N operations on updating pk,qk,λ1k,andλ2k, respectively. Therefore, the total cost of one step is 4NlogN+107N+1, which gives rise to the ONlogN complexity.


## 3.1. Algorithm Description

Firstly, we introduce two auxiliary variablesp and q to reformulate (10) as follows:(13)minu,p,q∫Ωαx,yp+βx,yqdxdy+12λ∫Ωu−f2dxdy,s.t.p=∇u,q=Δu.The augmented Lagrangian functional corresponding to (13) is(14)ℒu,p,q;λ1,λ2=∫Ωαx,yp+βx,yqdxdy+12λ∫Ωu−f2dxdy,+∫Ωλ1·p−∇udxdy+r12∫Ωp−∇u2dxdy,+∫Ωλ2q−Δudxdy+r22∫Ωq−Δu2dxdy,where λ1andλ2 are Lagrangian multipliers and r1andr2 are positive penalty parameters. By applying the ADMM, one needs to minimize three subproblems with other variables fixed and update two Lagrangian multipliers at each iteration which is(15)uk+1=argminu12λ∫Ωu−f2dxdy+r12∫Ωpk−∇u+λ1kr12dxdy+r22∫Ωqk−Δu+λ2kr22dxdy+δ12∫Ωu−uk2dxdy,pk+1=argminp∫Ωαk+1x,ypdxdy+r12∫Ωp−∇uk+1+λ1kr12dxdy+δ22∫Ωp−pk2dxdy,qk+1=argminq∫Ωβk+1x,yqdxdy+r22∫Ωq−Δuk+1+λ2kr22dxdy+δ32∫Ωq−qk2dxdy,λ1k+1=λ1k+r1pk+1−∇uk+1,λ2k+1=λ2k+r2qk+1−Δuk+1,where δ1,δ2,δ3>0 are three penalty parameters and αk+1x,yandβk+1x,y can be explicitly updated as(16)αk+1x,y=∇11+∇uk+12,βk+1x,y=11+∇uk+12.To solve theu-subproblem in (15), we consider its Euler–Lagrange equation, which has the form(17)u−fλ+r1divpk−∇u+λ1kr1−r2Δqk−Δu+λ2kr2+δ1u−uk=0,namely,(18)1λ−r1div∇+r2ΔΔ+δ1u=fλ−divr1pk+λ1k+Δr2qk+λ2k+δ1uk.In order to efficiently solve this nonlinear equation, we use the fast Fourier transform (FFT) with periodic boundary conditions. More concretely,uk+1 is given by(19)uk+1=ℱ−1ℱf/λ−divr1pk+λ1k+Δr2qk+λ2k+δ1uk1/λ−r1ℱdiv∇+r2ℱΔΔ+δ1,where ℱ is the FFT.Bothp and q subproblems have closed form solutions. Specifically, the solutions of p-subproblem and q-subproblem are given by the shrinkage operator:(20)pk+1=shrinkager1∇uk+1−λ1k/r1+δ2pkr1+δ2,αk+1x,yr1+δ2,and(21)qk+1=shrinkager2Δuk+1−λ2k/r2+δ3qkr2+δ3,βk+1x,yr2+δ3,where shrinkages,γ=maxs−γ,0°sgns. Here, ° denotes the component-wise multiplication.We summarize the overall algorithm for solving the spatially adaptive relaxation model (10) in Algorithm 1.Algorithm 1: (ADMM for solving (10)).
Step 1. Input λ,r1,andr2. Initialization: u0=f,p0=0,q0=0, λ10=0,andλ20=0. Let k≔0.Step 2. Compute uk+1 by (19).Step 3. Update αk+1x,y,βk+1x,y using uk+1 according to (16).Step 4. Compute pk+1 by (20).Step 5. Compute qk+1 by (21).Step 6. Update λ1k+1,λ2k+1 by the last two equations of (15).Step 7. Stop if the stopping criterion is satisfied. Otherwise, let k≔k+1, and go to Step 2.


## 3.2. Convergence Analysis

In this subsection, we give a convergence theory for the proposed algorithm (Algorithm1) under the assumption that a solution of the proposed model (10) exists.Theorem 3.
Assume thatuk,pk,qk,λ1k,λ2kk∈N is a sequence generated by Algorithm 1 and u¯,p¯,q¯,λ¯1,andλ¯2k are the points satisfying the following first-order optimality conditions:(22)u−fλ+∇·λ1−Δλ2=0,αx,yg+λ1=0,g∈∂pandαx,y=∇11+∇u2,βx,yh+λ2=0,h∈∂qandβx,y=11+∇u2,p−∇u=0,q−Δu=0.Letuek=uk−u¯,pek=pk−p¯,qek=qk−q¯,λ1ek=λ1k−λ¯1,andλ2ek=λ2k−λ¯2 be the error sequences. If for any gk∈∂pk,g¯∈∂p¯, and hk∈∂qk,h¯∈∂q¯,(23)Δ1k≔αkx,y−α¯x,ygk,pk−p¯,α¯x,y=∇11+∇u¯2,Δ2k≔βkx,y−β¯x,yhk,qk−q¯,β¯x,y=11+∇u¯2,are nonnegative. Then, the generated sequence uk,pk,qk,λ1k,λ2kk∈N converges to a limit point u∗,p∗,q∗,λ1∗,λ2∗ which satisfies the first-order optimality conditions (22).Our proposed model (10) can be seen as a variation of the previous model in [8]. Specifically, the operator ∇2 in (3) is changed to Δ. Furthermore, an ADMM-based algorithm same to [8] is applied to solve the proposed model (10). Therefore, the rigorous convergence proof of Algorithm 1 can refer to the proof presented in [8].


## 3.3. Complexity Analysis

In the following, we will give some simple analysis of complexity for the proposed method.For our proposed method, the cost of updatinguk is 4NlogN+26N+1 at each iteration, where N=mn is the total number of grid points. If a root operation is equal to 12 times addition or subtraction or multiplication or division, at each iteration it will spend 61N (3N times root operations), 11N, 6N, and 3N operations on updating pk,qk,λ1k,andλ2k, respectively. Therefore, the total cost of one step is 4NlogN+107N+1, which gives rise to the ONlogN complexity.


## 4. Numerical Experiments

In this section, we give some experimental results for the proposed model (10) in terms of our proposed algorithm. For simplicity of presentation, we call it SATVL in the following. We will compare our results with those of the previous algorithm for solving the spatially adaptive TV and high-order TV regularization model (3), which we call SATVTV2. It should be noted that the authors in [8] proposed the same ADMM for solving model (3). All numerical experiments are conducted in MATLAB environment on a PC with a 3.85 GHz CPU processor. To assess the denoising performance qualitatively, the peak signal to noise ratio (PSNR) and the structural similarity (SSIM) [40] are adopted. Their definitions are expressed as(24)PSNR=10log1025521/MN∑i=1M∑j=1Nui,j−Ii,j2,SSIM=2μIμu+c12σIu+c2μI2+μu2+c1σI2+σu2+c2,where I and u represent the original clear image and the restored image, respectively, μI and μu denote the local mean values of images I and u, σI and σu express the standard deviations, σIu signifies the covariance value between I and u, and c1andc2 are two constants which can avoid near-zero denominator values. In general, higher PSNR and SSIM values signify better performance in image denoising.The relative error of the solutionuk is defined by(25)Ruk=uk−uk−11uk−11.It can give some important information about the convergence of iterations. Therefore, we track its value during iteration process. In experiments, the stopping criterion isRuk≤5×10−5 or the iteration reaches 500.The test images are shown in Figure1, which includes one synthetic and four natural images. Since the proximal terms are added for convergence analysis, we set δ1=δ2=δ3=0.Figure 1
Test images of (a) shape, (b) peppers, (c) cameraman, (d) Barbara, and (e) owl.
(a)(b)(c)(d)(e)The values of parameters significantly affect the convergence rates and other characteristics of the methods. In our experiments, we study the influence of parametersλ,r1,andr2 on the SATVTV2. Concretely speaking, after taking a reasonable guess for the parameters: r1andr2, we consider different values of λ. We observe which value of λ could contribute to the biggest value of PSNR, and then we set it to be the value of λ. And in turn, we use the same strategy for the choice of other parameters. Similar strategy is adopted for tuning the parameters of our SATVL. For the parameter tuning, the algorithms are terminated after the stopping criterion is satisfied.Experiment 1.
In this example, we use the shape image as the test image. We first consider to restore the image corrupted by Gaussian white noise withσ=10. For the SATVTV2, we set λ=7.4,r1=0.11,andr2=0.35. For our SATVL, we use λ=12,r1=0.03,andr2=0.25. Besides, the Gaussian white noises at levels σ=20 and σ=30 are added for comparison, respectively. The numerical results are arranged in Table 1, while the corresponding visual results for noise level σ=10 are shown in Figure 2. In Table 1, Iter. and CPU (s) denote the number of iterations and the CPU time in seconds required for the whole denoising process.Table 1
Numerical results.
ImageIndexσ=10σ=20σ=30SATVTV2SATVLSATVTV2SATVLSATVTV2SATVLShape128×128PSNR34.18834.01729.18628.85625.77725.536SSIM0.9880.9860.9760.9710.9590.954Iter.500500500500500500CPU (s)21.56211.56220.23411.34319.42112.500Peppers256×256PSNR33.35533.09829.42829.15227.08126.768SSIM0.9590.9560.9210.9140.8850.874Iter.203329410500500500CPU (s)21.21820.03142.40628.26546.59328.250Cameraman256×256PSNR33.11732.82329.02528.83827.24826.961SSIM0.9540.9490.9030.8990.8730.863Iter.205376283500367500CPU (s)25.07822.37529.23427.50037.59326.875Barbara512×512PSNR31.89831.55927.19726.85824.67224.557SSIM0.9430.9360.8600.8410.7680.761Iter.272349490500500500CPU (s)145.984140.890271.546195.187280.484203.578Owl321×481PSNR33.37533.09629.60029.37827.68527.498SSIM0.9520.9470.8950.8880.8480.841Iter.379461500500500500CPU (s)104.93689.436140.68692.030140.53092.751Figure 2
(a) Noisy shape image withσ=10; (b, c) restored shape images by SATVTV2 and SATVL, respectively; (d, e) the plots of PSNR and relative error versus CPU time, respectively.
(a)(b)(c)(d)(e)Experiment 2.
In our second example, the pepper image is adopted as the test image. The original image is degraded by Gaussian white noise withσ=10. The parameters are set to be λ=9,r1=0.005,andr2=0.7 for the SATVTV2 and λ=14.5,r1=0.001,andr2=0.7 for the SATVL. Then, we consider more Gaussian noises at σ=20 and σ=30, respectively. Here, the values of parameters need to be reset. Numerical results are also listed in Table 1. The visual results about noise level σ=10 are presented in Figure 3.Figure 3
(a) Noisy pepper image withσ=10; (b, c) restored pepper images by SATVTV2 and SATVL, respectively; (d, e) the plots of PSNR and relative error versus CPU time, respectively.
(a)(b)(c)(d)(e)Experiment 3.
In our third example, we test the cameraman image. The original image is corrupted by Gaussian white noise withσ=10. Here, parameters for the SATVTV2 and our SATVL are set to be λ=7.2,r1=0.008,andr2=0.67 and λ=12.4,r1=0.002,andr2=0.71, respectively. In order to show better visual results, we also test the cases of noise levels σ=20 and σ=30. Results of this experiment are given in Table 1. The visual results about the noise level σ=20 are presented in Figure 4.Figure 4
(a) Noisy cameraman image withσ=20; (b, c) restored cameraman images by SATVTV2 and SATVL, respectively; and (d–f) the plots of energy, PSNR, and relative error versus CPU time, respectively.
(a)(b)(c)(d)(e)(f)Experiment 4.
In this example, we choose the Barbara image. Same to the above three experiments, we test the noisy images with different noise levels. Related numerical results are reported in Table1, and the visual results about the noise level σ=30 are displayed in Figure 5.Figure 5
(a) Noisy Barbara image withσ=30; (b, c) restored Barbara images by SATVTV2 and SATVL, respectively; and (d–f) the plots of energy, PSNR, and relative error versus CPU time, respectively.
(a)(b)(c)(d)(e)(f)Experiment 5.
In the last example, the owl image is used. To show the advantage of our proposed method, we also test noisy images with three different noise levels. For each algorithms, related parameters have been tuned to obtain the best restoration performance uniformly.
In terms of PSNR and SSIM, the denoising results by these two methods are similar. It is easy to see from Figures2–5 that the visual results of the SATVL are comparable to the SATVTV2. As shown in Table 1, our SATVL is faster than the SATVTV2. Experiments have shown this situation at different noise levels. From the plots of PSNR and relative error versus CPU time in Figures 2–5, we can also draw this conclusion. Besides, it is shown from these plots that the SATVL has converged, just as the SATVTV2 did. Actually, it is clear that the energy and PSNR have converged to a steady state.


## 5. Conclusion

In this paper, we have proposed a new relaxation model based on mean curvature for adaptive image denoising. We analyzed the connections to the existing models. Furthermore, we proposed an efficient ADMM to solve the proposed relaxation model and studied its convergence and complexity. Compared with the previous method for an adaptive first- and second-order regularization-based image denoising model, experimental results show that our proposed method is fast. Obviously, the proposed regularization term can be applied to other image processing problems, such as image deblurring, and image segmentation.



---

*Source: 3416907-2020-12-01.xml*