A new difference of anisotropic and isotropic total variation regularization method for image restoration

1.   Introduction

Total variation ^[1] regularization method is widely used in image processing, such as image reconstruction and denoising. Generally, there are two types: isotropic and anisotropic. TV has a strong edge preserving ability and can process piecewise constant images. However, the image obtained by the TV often suffers from staircasing artifacts. Therefore, in order to suppress artifacts and improve image quality, TV variants have been studied by many scholars ^[2,3,4,5,6].

TV consists of the $L_1$ norm and gradient operator. From the factor of statistics, the $L_1$ norm may have deviation estimations for larger coefficients, where it will lose the oracle property ^[7]. Instead of the $L_1$ norm, the $L_0$ norm is the appropriate choice to ensure sparseness. However, the minimization of $L_0$ is a combination optimization problem and NP hard. Then, several penalty functions that integrate the benefits of $L_1$ and $L_0$ norms while circumventing their limitations have been suggested, such as smoothly clipped absolute deviation ^[8], capped-$L_1$ ^[9], transformed-$L_1$ ^[10], and $L_1-L_2$ ^[11]. Numerous algorithms employing convex norms on gradients have also proposed, including the alternating minimization algorithm ^[12] and the augmented Lagrange multiplier algorithm ^[13]. The alternating direction method of multipliers (ADMM) ^[14] and primal-dual algorithm ^[15] are two benchmark methods. In fact, the primal-dual hybrid gradient method introduced by Chambolle and Pock in ^[15] is equivalent to the preconditioned ADMM. Additionally, the widely used split Bregman method ^[16] can be considered as an instance of ADMM with two variables. Additionally, there are also many algorithms for solving nonconvex models. As variational models are inherently nonconvex, obtaining the solution for global optimization is difficult. Some approaches, such as the proximal alternating linearized minimization method ^[17] and the difference-of-convex algorithm (DCA) ^[18], typically converge to local minimizers. In ^[19], DCA was firstly created by Dinh et al. for solving the DC programming. Because of its simplicity and efficiency, it has been widely used and extensively studied.

DC programming and DCA are the natural extension from modern convex analysis to nonconvex analysis, broad enough to include most real world nonconvex optimizations. During the last decades, these theoretical and algorithmic tools have been greatly enriched. He et al. ^[20] proposed a unified Douglas-Rachford algorithm for generalized DC programming. Niu et al. ^[21] developed a refined inertial DC algorithm for DC programming. Artacho et al. ^[22,23] investigated the line search idea in a DC program with linear constraint. Wen et al. ^[24] used the Nesterov acceleration technique for a special class of DC program. Moreover, DC programming and DCA have always been popular for nonconvex models in image processing, such as image restoration, image reconstruction, and compressive sensing. In ^[25,26], Lou et al. proved that $L_1-L_2$ is closer to $L_0$ than $L_1$ and developed a weighted difference of anisotropic and isotropic TV as a regularization for image denoising, deblurring, segmentation and magnetic resonance imaging (MRI). Its solution was obtained by DCA and the split Bregman. In ^[27], Li et al. introduced a multiplicative noise removal model, whose objective function can be written as the difference of two convex functions. DCA with a primal-dual hybrid gradient method was used for solving nonconvex programming. Sun et al. ^[28] studied the DCA to address the Log-norm TV image restoration problem and ADMM was used to solve the resulting subproblem. In ^[29], for impulsive noise removal, DCA with adaptive proximal parameters was proposed. DCA is a simple descent algorithm. Interested readers are encouraged to consult these state of the art developments in ^[30].

The contributions of this study are as follows. First, motivated by the idea that $L_1-L_2$ can be considered as a Fischer-Burmeister (FB) function, we propose a new nonconvex TV regularization term based on a general FB function. Second, two efficient DCAs are designed for solving nonconvex optimization problems. The proposed DCA needs to solve the TV type of subproblem, which can be efficiently handled by ADMM. Last, experiments on image denoising and MRI reconstruction show that the proposed methods improve on the classical TV model and are compared with state-of-the-art methods. Especially, for MRI, the new model can successfully reconstruct the Shepp-Logan image using merely 7 radial lines.

The rest of this paper is organized as follows. We show our new models in Section 2. Section 3 presents the DCA with ADMM for solving the new nonconvex models. Section 4 shows the effectiveness of the new models and algorithms through the presentation of experimental results. The conclusion is provided in the last section.

2.   New models

The problem of restoring image with additive noise is the following:

where $f\in R^m$ is noise image, $A\in R^{m\times n}$ is a linear operator, $u\in R^n$ is original image and $e\in R^m$ is Gaussian noise. To simplify the presentation, we express both the original and observation images or data as vectors.

Our aim is to restore the initial image $u$ from the noisy data $f$. Image restoration is a typical inverse problem, and it's impossible to obtain the original image $u$ from $f$ by just the direct inversion of (2.1). A classic method is variational regularization for addressing the ill-posed problem. This approach can make the solutions stable. Then, for the Gaussian noise removal, the famous TV model is as follows:

where $\|\cdot\|$ represents the Euclidean norm. As discussed above, TV can be classified into two categories: isotropic and anisotropic:

where $D$ is the gradient operator, $D_x$ and $D_y$ are the horizontal and vertical derivative operators, respectively, and $\|\cdot\|_1$ means the $L_1$ norm.

Illuminated by $L_1-L_2$ minimization in compression sensing ^[11], Lou et al. ^[26] proposed the well-known $L_1-\gamma L_2$ regularization function, which is the coupling of anisotropic and isotropic TV:

where $\gamma$ is an element of the interval(0, 1]. In ^[31], the authors find that $J(u)$ can be considered as the application of the FB function to the gradient of TV (FBTV), where the scalar-valued FB function is defined as

FBTV considers the non-sparse gradient vectors along the boundaries in the image. The effectiveness of FBTV was demonstrated to surpass that of either anisotropic or isotropic TV.

Motivated by $J(u)$ can be described as a FB function and $L_1-L_2$ is closer to $L_0$ than $L_1$, we propose a new difference regularization method by using the following generalized FB (GFB) function:

Now, we present the nonconvex variational regularization models, named GFBTV-C, as

and GFBTV as

where $\alpha > 0$ is the regularization parameter. It is observed that $\psi(v_1, v_2)$ takes into account other gradient directions, not just the horizontal or vertical. This allows it to regulate the impact of anisotropic TV for detecting sharp edges while reducing the blocky artifacts. To better understand the different metrics, we plot the level curves corresponding to $L_0$, $L_1$, $L_1-0.5L_2$, $L_1-L_2$ and GFB with $\theta$ = 0.1 and 0.01 in Figure 1. From the figure, it can observe that the contour lines of GFB are bending more inward and closer to $L_0$. This illustrates that GFB can enhance sparsity, and $\theta$ acts like a parameter controlling to what extent.

3.   Numerical algorithms

Though $\|\sqrt{(1-\theta)(|D_xu|^2+|D_yu|^2)+\theta(|D_xu|-|D_yu|)^2}\|_1$ is not convex, DCA can also be used to easily solve the GFBTV-C (2.2) and GFBTV (2.3) models. In the following, we apply the DCA with ADMM to handle the two nonsmooth and nonconvex models.

3.1.   DCA for solving the model GFBTV-C

First, in order to use the DCA, we separate into (2.2) as follows:

where

DCA optimizes (3.1) via linearizing $\Psi(u)$:

where

If the denominator is zero, the corresponding value is set as zero. Therefore, by the form of ADMM, the convex subproblem (3.2) could be rewritten as the unconstrained problems:

where $z$ is a Lagrange multiplier which enforces $Au = f$ and the parameter $\mu$ is the penalty parameter.

The subproblem (3.3) equals to solve the TV problem. Therefore, ADMM is also used to find the solution. In order to utilize ADMM, two auxiliary variables, $d_x$ and $d_y$, are introduced to rephrase the original minimization problem (3.3) as a constrained problem:

As per the classical optimization theory, the augmented Lagrangian function for Eq (3.4) can be expressed as follows:

where $b_x, b_y$ are the Lagrange multipliers and $\lambda > 0$ is the penalty parameter.

Subsequently, our attention turns to obtain the solution for subproblems in (3.5). The $u$-subproblem can be written as follows:

Using the first-order optimality condition of Eq (3.6), we have

In the problem of reconstructing undersampled magnetic resonance imaging (MRI) using compressed sensing ^[32], the composite $A = RF$ is formed by combining a sampling operator $R$ and a Fourier transform operator $F$. When periodic boundary conditions are applied, both matrices $D^TD$ and $A^TA$ become block circulant with circulant blocks. As a result, the coefficient matrix $\mu A^TA - \lambda D^TD$ can be diagonalized using a discrete Fourier transform. This means that Eq (3.7) can be efficiently solved by a fast Fourier transform (FFT):

Here, $\mathcal{F}(\cdot)$ and $\mathcal{F}^{-1}(\cdot)$ refer to the FFT and inverse FFT operations, respectively. The symbol $"*"$ represents a complex conjugation, while $" \odot "$ denotes an elementwise multiplication.

For $d_x, d_y$-subproblem in (3.5), the solution can be updated via soft shrinkage

where

and $\mathrm{sign}$ as the signum function. Additionally, the Lagrange multipliers are updated as follows:

At last, we present the DCA with ADMM for solving the proposed model GFBTV-C in Algorithm 1.

3.2.   DCA for solving the model GFBTV

For the corresponding unconstrained problem (2.3), the DCA subproblem is expressed as

Similar to (3.5), we have

The augmented Lagrangian function for Eq (3.12) can be written as

The next step is to express the $u$-subproblem as

Then, using the first-order optimality condition of Eq (3.14), we get

Moreover, if $D^TD$ and $A^TA$ have special structure, Eq (3.15) can be solved via FFT as follows:

For $d_x, d_y$-subproblem in (3.13), the solution can be updated via soft shrinkage

and the Lagrange multipliers are updated as follows:

Algorithm 2 is presented to demonstrate the DCA for resolving the proposed model GFBTV. To tackle the unconstrained problem (2.3), Algorithm 2 is nearly identical to Algorithm 1, but differs by including an extra update on $z$.

4.   Experiments

We present the experimental results of model GFBTV-C (2.2) utilizing Algorithm 1 for MRI and model GFBTV (2.3) employing Algorithm 2 for image denoising. From Figure 1, we find that $\theta$ = 0.01 and 0.1 can approximate the $L_0$ norm perfectly. Hence, we choose the control parameter $\theta$ = 0.01 for image denoising and $\theta$ = 0.1 for MRI naturally. All experiments were executed on a desktop computer with a 2.9 GHz processor and 16GB RAM. The peak signal-to-noise ratio (PSNR) is used to evaluate the restored quality, which is defined as

where MSE represents the mean square error between $u_k$ and $u_{ori}$, $u_k$ is the recovered image, and $u_{ori}$ represents the original image. Furthermore, it is worth noting that the structure similarity (SSIM) ^[33] has also been reported. It knows that as the PSNR and SSIM values increase, the quality of the reconstructed image improves.

4.1.   Image denoising

In this subsection, we will discuss the issue of image denoising, in which the operator $A$ is set to an identity matrix for the nonconvex model GFBTV. We compare it with several well-known variational models and a deep learning-based approach named DnCNN ^[34]. The first variational model we consider is the TV model, which is also known as the classic ROF model ^[1]. The second model is $L_1-\gamma L_2$ with $\gamma$ equaling 0.5 or 1, as discussed in ^[26]. The last model is the transform total variation (TTV) in ^[35], and DCA is formulated for the TTV minimization problem. We can efficiently solve the new model GFBTV using Algorithm 2, where $u_{k+1}$ can be calculated by (3.16) with $A = I$.

For this experiment, 100 images were randomly selected from the Pascal dataset as the test images. Both of these images were subject to degradation caused by additive Gaussian noise with a mean of zero and a standard deviation of $\sigma$. The values of $\sigma$ chosen for the experiment are 0.05 and 0.1. The penalty parameter $\mu$ was set to 5 for $\sigma$ = 0.1 and 15 for $\sigma$ = 0.05 during testing. For the proposed model, we set $\theta = 0.01$. The default parameter settings are the same as the references ^[26,35] for a fair comparison. We initialize $u_0 = f$, and the maximum number of DCA is 2.

The average comparison results are reported in Table 1. The PSNR and SSIM values of the different algorithms for various levels of $\sigma$ are provided. It is shown from Table 1 that DnCNN is superior to the gradient models in terms of PSNR and SSIM values. Although our new approach has slightly lower PSNR and SSIM values compared to DnCNN, it outperforms TV based models for various noise levels about SSIM. The denoising results of the six models are presented in Figures 2–5. All of the methods effectively remove noise as shown in the four figures and some fine structures and details are easy to see from the zoomed regions. DnCNN demonstrates the highest performance regarding PSNR among the compared approaches. However, when evaluating the four test images, our proposed method exhibits higher SSIM values than both DnCNN and other TV based models, and can reserve more details of the original images. It is also obvious that artifacts appear in the DnCNN, TV and $L_1-L_2$ methods. Hence, the experimental results show that the new model can effectively reduce the blocky artifacts and obtain better visual effects.

4.2.   Magnetic resonance imaging

In this subsection, we consider the problem of MR reconstruction. Compressive sensing based magnetic resonance imaging technology is widely used in images solving inverse problems. The sparse image reconstruction model recovers images from under-sampled k-space data, which effectively solves the problem of data acquisition efficiency in MRI domain. Consequently, we adopt the constrained formulation GFBTV-C (2.2), and Algorithm 1 is used for solving it. Additionally, we compare with TV, $L_1-\gamma L_2$ with $\gamma$ = 0.5 or 1 in ^[26], and TTV ^[35]. In this test, $\theta$ is set as 0.1. The related parameters in the proposed method and competing methods are left to their default settings in ^[26,35].

First, we consider the issue of MRI reconstruction utilizing a Shepp-Logan phantom from seven and eight radial projections. The results of various models in reconstructing the Shepp-Logan phantom are shown in Figure 6. Among these methods, only GFBTV-C can attain precise recovery using 7 lines, and the result is excellent. When it comes to 8 lines, GFBTV-C, TTV and $L_1-0.5L_2$ can achieve perfect reconstruction. TV and $L_1-L_2$ methods have some artifacts in the reconstructed images.

Second, we use two brain images test as the test images: Brain-1 and Brain-2 with dimensions of 256 $\times$ 256. These are generated from the literatures ^[36,37]. Then, we perform undersampling reconstruction under different sampling patterns. The k-space data for Brain-1 were randomly sampled using a Cartesian mask with 87 lines. For Brain-2, a variable density mask with a sampling rate of 30% was used to undersample the k-space data. The brain images reconstructed the residue images between the reconstruction images and original images by the five models, and are shown in Figures 7 and 8. All methods have better recovery effects. Visually, there is no significant difference. The residue images illustrate that the proposed method can visually gain better reconstructions. It appears that the TV model has the most residual images. From the residual images, it can be observed that TTV, $L_1-0.5L_2$, $L_1-L_2$ methods have similar effects on the image details. Table 2 presents the numerical performance of the five compared models. From the table, we see that the TV model needs less CPU time (in second). TTV, $L_1-0.5L_2$, $L_1-L_2$ and our GFBTV-C have similar time costs. Meanwhile, we know that GFBTV-C provides an improved reconstruction quality than TV, $L_1-0.5L_2$, and $L_1-L_2$ as it has the highest PSNR and SSIM values. The PSNR values of the new method are improved by at least 0.82 dB, and the SSIM values around 0.02. As the sampling rate becomes smaller, the performance of the proposed model improves more significantly.

5.   Conclusions

This paper introduces a novel nonconvex TV term that utilizes the GFB function to restore images corrupted by Gaussian noise. Additionally, we use DCA to solve the optimization problem. The results of our experiments on MRI and image denoising exhibit that our proposed models are capable of achieving better solutions in comparison to both convex models and other nonconvex models, and produce more preferable results.

Besides, we observe that the new regularization function has an improved to encourage sparsity compared to $L_1$, $L_1-\gamma L_2$ norm, and approximates $L_0$ norm more closely. This property leads to a better solution. Furthermore, the proposed model succeeds in reconstructing the Shepp-Logan phantom image from merely 7 radial lines. Hence, the superiority of our approach is more apparent when the sampling rate is relatively low or the available measurements are limited. The proposed model in this paper is mainly developed for image denoise and MRI, in the future, we will intend to extend the new method to computed tomography, Parallel MRI, electrical impedance tomography and so on. It is also essential to highlight that the convergence of our algorithm remains an open challenge and will be addressed in the future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgement

This work was supported in part by the National Natural Science Foundation of China under Grant 11901137 and Grant 61967004, in part by the Guangxi Key Laboratory of Automatic Detecting Technology and Instruments under Grant YQ22108, and in part by the Innovation Project of GUET Graduate Education under Grant 2022YCXS161.

Conflict of interest

The authors declare no conflict of interest.