Edge adaptive hybrid regularization model for image deblurring

As mentioned above, the harmonic model is effective in suppressing the noise in smooth areas of the image, however, it also brings about the blur to the edges inevitably. On the other hand, the total variation-based methods protect the image edges efficiently but still have a staircase artifacts in the smooth regions. In this paper, we combine these two models to make up for their respective shortcomings. Note that the simple combination is not distinctive enough for image restoration. The fixed parameter in the regularization term operates uniformly on the whole image, while ideal edges and the smooth areas have opposite demand of regularization. Specifically, the larger parameter of regularization term is conducive to restrain noise in the smooth area, however has a side effect of blurring the edges. On the other hand, a smaller regularization parameter retains edge information but brings about a staircase effect to the smooth areas. Accordingly, we introduce two edge adaptive parameters for TV and Tikhonov regularization terms respectively, which treat the edges and the smooth areas flexibly and hence improve the denoising performance in a reasonable way. On this basis, a new edge adaptive hybrid regularization (EAHR) model is proposed. The algorithm first detects the edges of the image to build an edge information matrix. It then adjusts parameters for TV and Tikhonov terms according to the local information of each pixel. Overall, the Edge Adaptive Hybrid Regularization (EAHR) model for image restoration we consider is given as follows

where $\bm{\alpha}_{1}$ and $\bm{\alpha}_{2}\in\mathbb{R}^{M\times N}$ are spatially varying non-negative parameter matrices, $\|\nabla\mathbf{u}\|_{1,\bm{\alpha}_{1}}$ is given by (1.9) with $\bm{\alpha}_{1}$ instead of $\bm{\alpha}$, and $\|\nabla\mathbf{u}\|_{2,\bm{\alpha}_{2}}$ is given by

In the proposed model, the edge adaptive parameters $\alpha_{1}$ and $\alpha_{2}$ for TV and Tikhonov regularization terms play key roles in the final recovered images. As explained in  [30], the edge detection operator is a function of $|\nabla\mathbf{u}|$ for detecting the image edges. Following it, we define the edge information matrix $\mathbf{E}$ as

where $\mathbf{G}$ is Gaussian kernel, and $\mathbf{\mathbf{G}\otimes\nabla\mathbf{u}}$ is the convolution of $\mathbf{G}$ with $\nabla\mathbf{u}$ as follows

In order to improve the robustness of edge adaptive algorithm, we binarize the edge adaptive parameters according to the edge information matrix $\mathbf{E}$. The spatially varying non-negative parameter matrices $\bm{\alpha}_{1}$ and $\bm{\alpha}_{2}$ in the model (2.1) are defined by

and

where $\alpha_{1}$ and $\alpha_{2}$ are positive parameters, $\theta_{1},\theta_{2}\in(0,1)$ are scaling parameters, and $\tau>0$ is a threshold. With $\bm{\alpha}_{1}$ and $\bm{\alpha}_{2}\in\mathbb{R}^{M\times N}$ fixed, the proposed model is convex. Therefore, it can be solved efficiently by sPADMM.

The proposed model (2.1) balances the relationship between the edge and smooth areas well by adopting the two parameters based on the edge information. By the proposed weights, it achieves the goal of retaining the edge and eliminating the staircase effect simultaneously. Moreover, while the TV-based models sometimes erroneously treats the noise in the smooth area as the image edge, the proposed model can effectively eliminate these misjudged pixels by edge detection process.

As mentioned, our proposed model can be efficiently solved by sPADMM. To apply this algorithm, we first introduce an auxiliary variable $\mathbf{k}$ to take the place of $\nabla\mathbf{u}$ in the both terms of $\|\nabla\mathbf{u}\|_{1,\bm{\alpha}_{1}}$ and $\|\nabla\mathbf{u}\|_{2,\bm{\alpha}_{2}}^{2}$, then the model (2.1) can be reformulated as the following constrained optimization problem

The augmented Lagrangian function of the problem (2.6)  is defined as

where $\bm{\lambda}=[\bm{\lambda}_{1},\bm{\lambda}_{2}]\in\mathbb{R}^{2M\times N}$ is the Lagrange multiplier parameter matrix with $\bm{\lambda}(i,j)=(\bm{\lambda}_{1}(i,j),\bm{\lambda}_{2}(i,j))$, and $\beta>0$ is the penalty parameter. In each iteration of sPADMM, we minimize $L$ with respect to $\mathbf{u}$ for fixed $\mathbf{k}$ and then minimize $L$ with respect to $\mathbf{k}$ for fixed $\mathbf{u}$. After these two steps, we update the multiplier $\lambda$. Hence, the solution to the problem (2.7) is approached by iterating the following three equations:

Here, for any $\mathbf{x}\in\mathbb{R}^{M\times N}$, $\mathbf{S}\in\mathbb{R}^{M\times N}$ is the self-adjoint positive semidefinte matrix and $\mathbf{S}\mathbf{x}$ is the matrix multiplication of $\mathbf{S}$ and $\mathbf{x}$, $\|\mathbf{x}\|_{\mathbf{S}}=\sqrt{\langle\mathbf{x},\mathbf{S}\mathbf{x}\rangle}$ denotes the matrix norm. If we take $\mathbf{S}_{1}=0$ and $\mathbf{S}_{2}=0$, the sPADMM will be the alternating direction method of multipliers (ADMM) [5]. In our algorithm, $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ are positive definite matrices. The sPADMM for our EAHR model is given by Algorithm 1.

In each iteration of sPADMM, the subproblems (2.8) and (2.9) need to be solved. The $\mathbf{u}-$subproblem is written as

In our experiment, we choose the self-adjoint positive definite linear operators $\mathbf{S}_{1}=\rho_{1}\mathbf{I}$, where $\mathbf{I}$ is the unit matrix and $\rho_{1}>0$. Note that the subproblem (2.11) is a quadratic optimization problem subject to the optimality condition

which is solved by Fourier transform  [40] as follows

Next, we analyze the nonsmooth subproblem (2.9) in detail to find its global minimizer.

In order to solve subproblem (2.9), we first introduce a proposition as follows.

Proposition 2.1.  For any  $\alpha,\beta>0$, and $z,t\in\mathbb{R}^{2}$ , the minimizer of

is given by

where $\left(a\right)_{+}=\max\left(a,0\right)$, $\frac{\mathbf{0}}{0}=\mathbf{0}$ and $\mathbf{0}=(0,0)$.

Proof: If $t\neq\mathbf{0}$, we have

Let  $\frac{\partial f}{\partial z}=0$, we get

The equation (2.15) implies that

It follows from (2.16)  and  (2.15) that

Then the equation (2.14) holds when $t\neq\mathbf{0}$. The equation (2.14) also holds when $t=\mathbf{0}$. Therefore, Proposition 4.1 is proved.

The subproblem (2.9) is written as

The problem above is equivalently expressed as

In our experiment, we use the self-adjoint positive definite linear operators $\mathbf{S}_{2}=\rho_{2}\mathbf{I}$, where $\mathbf{I}$ is the unit matrix and $\rho_{2}>0$. Letting $z=\mathbf{k}(i,j)$, $q=\nabla\mathbf{u}(i,j)$, $p_{1}=\bm{\alpha}_{1}(i,j)$, $p_{2}=\bm{\alpha}_{2}(i,j)$ and $\widetilde{\lambda}=\bm{\lambda}(i,j)$, we rewrite the problem above as follows

where  $|z|=\sqrt{z_{1}^{2}+z_{2}^{2}}$. Let

the problem (2.18) is equivalent to

From Proposition  2.1, the optimal solution of this problem is

By the definition of $z^{*}$ and $t$, we have

where

Then $\mathbf{k}^{n}$ with $\mathbf{k}^{n}(i,j)$ given by (2.22) is the solution of (2.9).

The pseudo-code of the proposed algorithm is given in Algorithm 1. The linear-rate convergence of the algorithm is characterized in Theorem 2.1.

Theorem 2.1.  Let $\{(\mathbf{u}^{n},\mathbf{k}^{n},\mathbf{\lambda}^{n})\}$ be the sequence generated by the algorithm sPADMM. If $\eta\in(0,\frac{1+\sqrt{5}}{2})$, then the sequence $\{(\mathbf{u}^{n},\mathbf{k}^{n},\mathbf{\lambda}^{n})\}$ converges to the KKT point of (2.6).

The detailed proof of the theorem can be found in the Appendix.

In this section, we introduce all the parameters used in our algorithm. From (2.5), $\theta_{1},\theta_{2}\in(0,1)$ are scaling parameters, and $\tau>0$ is a threshold which is used to adjust the value of $\alpha_{1}$ and $\alpha_{2}$. There are three parameters $(\mu,\beta,\tau)$ in our Algorithm 1 that need to be adjusted. We choose the self-adjoint positive definite linear operators $\mathbf{S}_{1}=\mathbf{S}_{2}=\rho\mathbf{I}$ with $\rho=0.1$ and the maximum iterative number (denoted by $MaxIter$ ) is set as $500$. Through a multitude of experiments, we choose better parameters to improve the quality of restoring images damaged by blur and noise.

We have done a lot of simulations to find the best value of parameters $\mu$, $\tau$ and $\beta$ for each blurring kernel with different noise levels. As a result, we get the parameters $\mu$, $\tau$ and $\beta$ as functions of $\sigma\in(0,100]$ for each blurring kernel. Taking Gaussian blur as a example, the parameters $\mu$, $\tau$ and $\beta$ are given as:

where $\sigma\in(0,100]$. When $\sigma=3$, the corresponding values of $\mu$, $\tau$ and $\beta$ are 5000, 0.9 and 0.1 respectively. In the case of Motion blur and Average blur, we have similar functions for the parameters $\mu$, $\tau$ and $\beta$ which are given in the code.

The parameter settings for the other compared algorithms are given below. The TV model [38] has two parameters $\mu$ and $\lambda$, where $\mu$ is the parameter to balance the data fidelity and regularity terms and $\lambda$ plays an important role in analyzing the convergence. In our tests, we let $\mu\in\{200,300,400,700,1500\}$ and $\lambda\in\{10,15,20,50,80,100,200,300\}$. For the DCA model [31], we fix the parameter $\lambda=1$, and then choose the parameter $\mu$ from the sequence $\{50,100,150,200,300,380,700,800,1500\}$. The parameters used by the TRL2 model [42] are set as $\alpha\in\{5,10,15,20,30,50,100,150,200\}$, $\tau\in\{0.1,0.2,0.3,0.5\}$ and $\beta\in\{0.1,0.5,1,2,3,5,10\}$. In the SOCF model [28], we take fixed values $\mu_{1}=0.2,\mu_{2}=1e-1$ and $\eta=1.6$ represented the step-length. We choose $\mu\in\{1e-4,1e-3\}$ and $\lambda\in\{2,3,4,7,11,13\}$. The BM3D model [15] has two parameters: regularized inversion (RI) which is the collaborative hard-thresholding for BM3D filtering and regularized Wiener inversion (RWI) for collaborative Wiener filtering. We set these two parameters as $RI\in\{4e-4,4e-3,1e-3\}$ and $RWI\in\{5e-3,3e-2\}$.

In order to verify the performance of our proposed model, we have tested the recovery results of gray images and color images under different blur kernels and different Gaussian noise levels. Tables 1,  2 and  3 show the values of PSNR, SSIM of the classic TV  [38], DCA with $L_{1}-0.5L_{2}$  [31], TRL2  [42], SOCF  [28] and BM3D  [15] on Gaussian blur $\mathrm{GB(9,5)}/\sigma=5$ and Motion blur $\mathrm{MB(20,60)}/\sigma=5$ and Average blur $\mathrm{AB(20,60)}/\sigma=3$. The number in bold means the best result. Obviously, our algorithm gives the best results in most of the cases. We also show the average of PSNR/SSIM value of the restored image in the Table 4. Comprehensive evaluations on the image set with three kinds of kernels and four different noise levels demonstrate the superiority of our proposed method when compared with other state-of-the-art methods.

In order to further demonstrate the superiority of our algorithm on image restoration more intuitively, we compare the visual quality of restored images of our proposed method and other state-of-the-art methods. We also zoom in key parts of the image for better illustrating. Figures $\ref{fig:GBg}-\ref{fig:ABc}$ show that our method yields the best image quality in terms of removing noise, preserving edges, and maintaining image sharpness. Looking carefully at the details in the figures, we can see that TV  [38] oversmoothes images, TRL2  [42] introduces staircase effect, DCA  [31] sharpens edges but loses some detail information, there exits the ringing phenomenon by SOCF  [28] and BM3D model  [15] introduces some artifacts.

Table 5 lists the computation time on a desktop (Intel(R) Core(TM) i5-8250 CPU @1.60 GHz), which reveal that our method is faster than DCA  [31], but is slower than TV  [38], BM3D  [15], TRL2 [42] and SOCF [28]. Acceleration will be considered in the future work.

In brief, our algorithm (EAHR) has competitive performance in sharpening the edges and removing the noise, which outperforms other mainstream methods both in PSNR/SSIM values and visual quality.