An automatic pixel-wise multi-penalty approach to image restoration

1. Introduction

Image restoration is an important task in many areas of applied sciences since digital images are frequently degraded by blur and noise during the acquisition process. Image restoration can be mathematically formulated as the linear inverse problem [1]

$$A\mathbf{u}+\mathbf{e}=\mathbf{b}$$

(1)

where $\mathbf{b}\in {\mathbb{R}}^M$ and $\mathbf{u}\in {\mathbb{R}}^N$ respectively are vectorized forms of the observed $m_x\times m_y$ image and the exact $n_x\times n_y$ image to be restored, $A\in {\mathbb{R}}^{M\times N}$ is the linear operator modeling the imaging system and $\mathbf{e}$ represents Gaussian white noise with mean zero and standard deviation $\sigma$. The image restoration problem (1) is inherently ill-posed and regularization strategies, based on the prior information on the unknown image, are usually employed in order to effectively restore the image $\mathbf{u}$ from $\mathbf{b}$.

In a variational framework, image restoration can be reformulated as a constrained optimization problem of the form

$$\begin{array}{cc}\underset{\mathbf{u}}{min}\hfill & {\displaystyle \frac{1}{2}{\parallel A\mathbf{u}-\mathbf{b}\parallel }^2+\mathcal{R}\left(\mathbf{u}\right)}\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & \mathbf{u}\in \Omega \hfill \end{array}$$

(2)

whose objective function contains a $L_2$-based term, imposing consistency of the model with the data, and a regularization term $\mathcal{R}\left(\mathbf{u}\right)$, forcing the solution to satisfy some a priori properties. Here and henceforth, the symbol $\parallel ·\parallel$ denotes the Euclidean norm. The constraint imposes some characteristics on the solution which are often given by the physics underlying the data acquisition process. Since image pixels are known to be nonnegative, a typical choice for $\Omega$ is the positive orthant.

The quality of the restored images strongly depends on the choice of the regularization term which, in a very general framework, can be expressed as

$$\mathcal{R}\left(\mathbf{u}\right)=\sum _{i=1}^p{\lambda }_i{\psi }_i\left(\mathbf{u}\right)$$

(3)

where the positive scalars ${\lambda }_i$ are regularization parameters and the ${\psi }_i\left(\mathbf{u}\right)$ are regularization functions for $i=1,\dots ,p$. The multi-penalty approach (3) allows to impose several regularity properties on the desired solution, however a crucial issue with its realization is the need to define reliable strategies for the choice of the regularization parameters ${\lambda }_i$, $i=1,\dots ,p$.

Therefore, in the literature, the most common and famous regularization approach is single-penalty regularization, also known as Tikhonov-like regularization, which corresponds to the choice $p=1$:

$$\mathcal{R}\left(\mathbf{u}\right)=\lambda \psi \left(\mathbf{u}\right).$$

(4)

In image restoration, smooth functions based on the $L_2$-norm or convex nonsmooth functions like the Total Variation, the $L_1$ norm or the Total Generalized Variation are usually used for $\psi \left(\mathbf{u}\right)$ in (4) [2,3]. Even in the case $p=1$, the development of suitable parameter choice criteria is still an open question. The recent literature has demonstrated a growing interest in multi-penalty regularization, with a significant number of researchers focusing on scenarios involving two penalty terms. Notably, the widely-used elastic regression in Statistics serves as an example of a multi-penalty regularization technique, integrating the $L_1$ and $L_2$ penalties from the Lasso and Ridge methods. However, the majority of the literature primarily addresses the development of suitable rules for parameter selection.

Lu, Pereverzev et al. [4,5] have extensively investigated two $L_2$-based terms, introducing a refined discrepancy principle to compute dual regularization parameters, along with its numerical implementation.

The issue of parameter selection is further discussed in [6], where a generalized multi-parameter version of the L-curve criterion is proposed, and in [7], which suggests a methodology based on the GCV method.

Reichel and Gazzola [8] propose regularization terms of the form

$${\varphi }_i\left(\mathbf{u}\right)={\parallel D_i\mathbf{u}\parallel }^2,i=1,\dots ,p$$

(5)

where $D_i$ are suitable regularization matrices. They present a method to determine the regularization parameters utilizing the discrepancy principle, with a special emphasis on the case $p=2$. Fornasier et al. [9] propose a modified discrepancy principle for multi-penalty regularization and provide theoretical background for this a posteriori rule.

Works such as [10,11,12,13,14] also explore multi-penalty regularization for unmixing problems, employing two penalty terms based on $L_q$ and $L_p$ norms, $0\le q<2$ and $2\le p<\infty$. The latter specifically concentrates on the $L_1$ and $L_2$ norms.

The study [15] assesses two-penalty regularization, incorporating $L_0$ and $L_2$ penalty terms, to tackle nonlinear ill-posed problems and analyzes its regularizing characteristics.

In [16], an automated spatially adaptive regularization model combining harmonic and Total Variation (TV) terms is introduced. This model is dependent on two regularization parameters and two edge information matrices. Despite the dynamic update of the edge information matrix during iterations, the model necessitates fixed values for the regularization parameters.

Calatroni et al. [17] present a space-variant generalized Gaussian regularization approach for image restoration, emphasizing its applicative potential.

In [18], a multipenalty point-wise approach based on the Uniform Penalty principle is considered and analyzed for general linear inverse problems, introducing two iterative methods, UpenMM and GUpenMM, and analyzing their convergence.

Here we focus on image restoration and, following [18], we propose to find an estimate ${\mathbf{u}}^{*}$ of $\mathbf{u}$ satisfying

$$\begin{array}{cc}\hfill {\mathbf{u}}^{*}=& {\displaystyle \arg \underset{\mathbf{u}\ge 0}{min}\frac{1}{2}{\parallel A\mathbf{u}-\mathbf{b}\parallel }^2+\frac{1}{2}\sum _{i=1}^N{\lambda }_i^{*}{\left(L\mathbf{u}\right)}_i^2,}\hfill \\ \hfill {\lambda }_i^{*}=& {\displaystyle \frac{\parallel A{\mathbf{u}}^{*}{-\mathbf{b}\parallel }^2}{N({\left(L{\mathbf{u}}^{*}\right)}_i^2+ϵ)}.}\hfill \end{array}$$

(6)

where $ϵ$ is a positive scalar and $L\in {\mathbb{R}}^{N\times N}$ is the discrete Laplacian operator. This model, named MULTI, is specifically tailored for the image restoration problem. Observe that MULTI incorporates a pixel-wise regularization term and includes a rule for choosing the parameters. Building upon UPenMM and GUpenMM methods analyzed in [18], we formulate an iterative algorithm for computing the solution $({\mathbf{u}}^{*},{\mathsf{\lambda }}^{*})$ of (6), where ${\mathsf{\lambda }}^{*}={({\lambda }_1^{*},\dots ,{\lambda }_N^{*})}^T$. In every inner iteration, once the regularization parameters are set, the constrained minimization subproblem is efficiently solved by a customized version of the Newton Projection (NP) method. Here, the Hessian matrix is approximated by a Block Circulant with Circulant Blocks (BCCB) matrix, which is easily invertible in the Fourier space. This modified version of NP was designed in [19] for single-penalty image restoration under Poisson noise and it is adapted here to the context of multi-penalty regularization. Consequently, the convergence of the modified NP method can be established.

The principal contributions of this work are summarized as follows:

• We propose a novel variational pixel-wise regularization model for image restoration.
• We devise an algorithm capable of effectively and efficiently solving the proposed model.
• Through numerical experiments, we demonstrate that the proposed approach can proficiently eliminate noise and blur in smooth areas of an image while preserving its edges.

The structure of this paper is as follows: Section 2 introduces the proposed algorithm. The results of numerical experiments are presented in Section 3, and finally, the conclusions are drawn in Section 4.

2. Materials and Methods

In this section we present the iterative algorithm that generates the sequence $({\mathbf{u}}^{\left(k\right)},{\mathsf{\lambda }}^{\left(k\right)})$ converging to the solution $({\mathbf{u}}^{*},{\mathsf{\lambda }}^{*})$ in (6).

Starting from an initial guess ${\mathbf{u}}^{\left(0\right)}$ taken as the observed image $\mathbf{b}$, the correspondent initial guess of the regularization parameters is computed as:

$${\lambda }_i^{\left(0\right)}=\frac{\parallel A{\mathbf{u}}^{\left(0\right)}{-\mathbf{b}\parallel }^2}{N({\mathcal{S}}_i^{\left(0\right)}+ϵ)},i=1,\dots ,N$$

(7)

where

$${\mathcal{S}}_i^{\left(0\right)}=\underset{i\in {\mathcal{N}}_i}{max}{\left(L{\mathbf{u}}^{\left(0\right)}\right)}_i^2,i=1,\dots ,N$$

(8)

and ${\mathcal{N}}_i$ is a neighborhood of size $R\times R$, (with R odd and $R\ge 1$) of the $i-$th pixel with coordinates $(x_i,y_i)$.

The successive terms $({\mathbf{u}}^{(k+1)},{\mathsf{\lambda }}^{(k+1)})$ are obtained by the update formulas reported in steps 3-5 of Algorithm 1. The iterations are stopped when the relative distance between two successive regularization vectors is smaller than a fixed tolerance $Tol>0$.

Algorithm 1 Input: ${\mathsf{\lambda }}^{\left(0\right)}\in {\mathbb{R}}^N$, $Tol$, $A,\mathbf{b},ϵ$, Output: $({\mathbf{u}}^{*},{\mathsf{\lambda }}^{*})$

1: Set $k=0$.   2: repeat 3:    ${\displaystyle {\mathbf{u}}^{(k+1)}=arg\underset{\mathbf{u}\in \Omega }{min}\left\{\frac{1}{2}{\parallel A\mathbf{u}-\mathbf{b}\parallel }^2+\sum _{i=1}^N{\lambda }_i^{\left(k\right)}{\left(L\mathbf{u}\right)}_i^2\right\}}$    4:    ${\displaystyle {\mathcal{S}}_i^{(k+1)}=\underset{i\in {\mathcal{N}}_i}{max}{\left(L{\mathbf{u}}^{(k+1)}\right)}_i^2),i=1,\dots N}$    5:     ${\displaystyle {\lambda }_i^{(k+1)}=\frac{\parallel A{\mathbf{u}}^{(k+1)}{-\mathbf{b}\parallel }^2}{N({\mathcal{S}}_i^{(k+1)}+ϵ)},i=1,\dots ,N}$    6:    $k=k+1$    7: until${\displaystyle \parallel {\mathsf{\lambda }}^{\left(k\right)}-{\mathsf{\lambda }}^{(k-1)}\parallel \le Tol\parallel {\mathsf{\lambda }}^{\left(k\right)}\parallel }$    8: ${\mathsf{\lambda }}^{*}={\mathsf{\lambda }}^{\left(k\right)}$, ${\mathbf{u}}^{*}={\mathbf{u}}^{\left(k\right)}$

Algorithm 1 is well defined, and it is possible to prove its convergence. Let’s assume for now that the maximum value of $\left(L{\mathbf{u}}^{\left(k\right)}\right)$ in the neighborhood ${\mathcal{N}}_i$ of the $i-$th pixel is exactly in the center of ${\mathcal{N}}_i$ with coordinates $(x_i,y_i)$. In this case, Algorithm 1 corresponds to UPenMM and, using Theorem 3.4 in [18] we can prove that the sequence $({\mathbf{u}}^{\left(k\right)},{\mathsf{\lambda }}^{\left(k\right)})$ converges to $({\mathbf{u}}^{*},{\mathsf{\lambda }}^{*})$ in (6). In the general case, when ${\left(L{\mathbf{u}}^{\left(k\right)}\right)}_i<{\mathcal{S}}_i^{\left(k\right)}$, it is possible to preserve convergence through a correction obtained from a convex combination between ${\left(L{\mathbf{u}}^{\left(k\right)}\right)}_i$ and ${\mathcal{S}}_i^{\left(k\right)}$. This approach is presented as Generalized Uniform Penalty method (GUPenMM) [18].

However, such algorithmic correction is not necessary in image restoration problems where it is not essential to compute solutions with a high level of accuracy. Indeed the human eye cannot distinguish the difference smaller than a few gray levels.

At each inner iteration, the constrained minimization subproblem (step 3 in Algorithm 1) is solved efficiently by a tailored version of the NP method where the Hessian matrix is approximated by a BCCB matrix easily invertible in the Fourier space.

Let us denote by ${\mathcal{J}}^{\left(k\right)}\left(\mathbf{u}\right)$ the function to be minimized at step 3 in Algorithm 1:

$${\mathcal{J}}^{\left(k\right)}\left(\mathbf{u}\right)=\frac{1}{2}{\parallel A\mathbf{u}-\mathbf{b}\parallel }^2+\sum _{i=1}^N{\lambda }_i^{\left(k\right)}{\left(L\mathbf{u}\right)}_i^2$$

and by $\mathbf{g}$ its gradient $\nabla {\mathcal{J}}^{\left(k\right)}\left(\mathbf{u}\right)$ where the iteration index k has been omitted for easier notation. Moreover, let ${\mathbf{g}}_{\mathcal{I}}$ denote the reduced gradient:

$${\left\{{\mathbf{g}}_{\mathcal{I}}\right\}}_i=\left\{\begin{array}{cc}{g}_i,\hfill & \mathrm{if}i\notin \mathcal{I}\left(\mathbf{u}\right);\hfill \\ 0,\hfill & \mathrm{otherwise};\hfill \end{array}\right.$$

(9)

where $\mathcal{I}\left(\mathbf{u}\right)$ is the set of indices [20]:

$$\mathcal{I}\left(\mathbf{u}\right)=\left\{i\right|0\le u_i\le \epsilon \mathrm{and}{g}_i>0\}$$

with

$$\epsilon =min\{\overline{\epsilon },w\},w=\parallel \mathbf{u}-{[\mathbf{u}-\nabla \mathcal{J}\left(\mathbf{u}\right)]}^{+}\parallel$$

and $\overline{\epsilon }$ is a small positive parameter.

The Hessian matrix ${\nabla }^2\mathcal{J}\left(\mathbf{u}\right)$ has the form

$${\nabla }^2{\mathcal{J}}^{\left(k\right)}\left(\mathbf{u}\right)=A^{T}A+L^T{\Lambda }^{\left(k\right)}L$$

where ${\Lambda }^{\left(k\right)}$ is the diagonal matrix with diagonal elements ${\lambda }_1^{\left(k\right)},\dots ,{\lambda }_N^{\left(k\right)}$.

A general iteration of the proposed NP-like method has the form:

$${\mathbf{u}}^{(\ell +1)}={\left[{\mathbf{u}}^{(\ell )}-{\alpha }^{(\ell )}{\mathbf{p}}^{(\ell )}\right]}^{+}$$

(10)

where ${\mathbf{p}}^{(\ell )}$ is the search direction, ${\alpha }^{(\ell )}$ is the steplength and ${[·]}^{+}$ denotes the projection on the positive orthant.

At each iteration ℓ, the computation of ${\mathbf{p}}^{(\ell )}$ requires the solution of the linear system

$$H^{\left(k\right)}{\mathbf{d}}^{(\ell )}={\mathbf{g}}_{\mathcal{I}}^{(\ell )},$$

(11)

where $H^{\left(k\right)}$ is the following approximation to ${\nabla }^2{\mathcal{J}}^{\left(k\right)}\left(\mathbf{u}\right)$

$$H^{\left(k\right)}=A^{T}A+{\mu }^{\left(k\right)}{L}^{T}L,{\mu }^{\left(k\right)}=\mathrm{mean}\left({\lambda }_1^{\left(k\right)},\dots ,{\lambda }_N^{\left(k\right)}\right).$$

(12)

Under periodic boundary conditions, $H^{\left(k\right)}$ is a BCCB matrix and system (11) can be efficiently solved in the Fourier space by using Fast Fourier Transforms. Hence, despite its simplicity, the BCCB approximation $H^{\left(k\right)}$ is efficient, since it allows to solve the linear system in $O\left(N{log}_{2}N\right)$ operations, and effective as it is shown by the numerical results. Finally, given the solution ${\mathbf{d}}^{(\ell )}$ of (11), the search direction ${\mathbf{p}}^{(\ell )}$ is obtained as

$$p_i^{(\ell )}=\left\{\begin{array}{cc}{d}_i^{(\ell )},\hfill & \mathrm{if}i\notin \mathcal{I}\left({\mathbf{u}}^{(\ell )}\right);\hfill \\ g_i^{(\ell )},\hfill & \mathrm{otherwise};\hfill \end{array}\right.i=1,\dots ,N.$$

(13)

The step length ${\alpha }^{(\ell )}$ is computed with the variation of the Armijo rule discussed in [20] as the first number of the sequence ${\left\{{\beta }^m\right\}}_{m\in \mathbb{N}}$, $0<\beta <1$, such that

$$\mathcal{J}\left({\mathbf{u}}^{(\ell )}\right)-\mathcal{J}\left({\mathbf{u}}^{(\ell )}\left({\beta }^m\right)\right)\ge \eta \left({\beta }^m\sum _{i\notin \mathcal{I}\left({\mathbf{u}}^{(\ell )}\right)}{g}_i^{(\ell )}{p}_i^{(\ell )}+\sum _{i\in \mathcal{I}\left({\mathbf{u}}^{(\ell )}\right)}{g}_i^{(\ell )}\left(u_i^{(\ell )}-u_i^{(\ell )}\left({\beta }^m\right)\right)\right)$$

(14)

where ${\mathbf{u}}^{(\ell )}\left({\beta }^m\right)={[{\mathbf{u}}^{(\ell )}-{\beta }^m{\mathbf{u}}^{(\ell )}]}^{+}$ and $\eta \in (0,\frac{1}{2})$.

We observe that the approximated Hessian $H^{\left(k\right)}$ is constant for each inner iteration ℓ and it is positive definite, then it satisfies

$$c_1{\parallel y\parallel }^2\le y^T{\left(H^{\left(k\right)}\right)}^{-1}y\le c_2{\parallel y\parallel }^2\forall y\in R^N,\forall \ell >0.$$

Then, the results given in [19] for single-penalty image restoration under Poisson noise, can be applied here to prove the convergence of the NL-like iterations to critical points.

The stopping criteria for the NP-like method are based on the relative distance between two successive iterates and the relative projected gradient norm. In addition, a maximum number of NP iterations have been fixed.

3. Numerical Experiments

All the experiments are performed under Windows 10 and MATLAB R2021a running on a desktop (Intel(R) Core(TM) i5-8250CPU@1.60 GHz). Quantitatively, we evaluate the quality of image restoration by the relative error (RE), improved signal to noise ratio (ISNR) and mean structural similarity index (MSSIM) measures. The MSSIM is defined by Wang et al. [21] and ISNR is calculated as:

$$\mathrm{ISNR}=20{log}_{10}\left(\frac{\parallel \mathbf{b}-\mathbf{u}\parallel }{\parallel \widehat{\mathbf{u}}-\mathbf{u}\parallel }\right)$$

where $\widehat{\mathbf{u}}$ is the restored image, $\mathbf{u}$ is the reference image and $\mathbf{b}$ is the blurred, noisy image.

Four reference images were used in the experiments: galaxy, mri, leopard and elaine, shown in Figure 1, Figure 2, Figure 3 and Figure 4. The first three images have size $256\times 256$, while the elaine image is $512\times 512$. In order to define the test problems, each reference image was convolved with two PSFs corrsponding to a Gaussian blur with variance 2, generated by the psfGauss function from the MATLAB toolbox Restore-Tool [1], and an out-of-focus blur with radius 5, obtained with the function fspecial from the MATLAB Image Processing Toolbox. The resulting blurred image was then corrupted by Gaussian noise with different values of the noise level $\delta =\parallel \mathbf{e}\parallel /\parallel \mathbf{b}\parallel$. The values $\delta =2.5·{10}^{-2},{10}^{-2},5·{10}^{-3}$ were used.

We compare the proposed pixel-wise multi-penalty regularization model (MULTI) with Tikhonov (TIKH), Total Variation (TV) and Total Generalized Variation (TGV) regularization with nonnegative constraints. The regularization parameter values for TIKH, TV and TGV have been chosen heuristically by minimizing the relative error values. The Alternating Direction Method of Multipliers (ADMM) was used for the solution of the TV-based minimization problem, while for TIKH, we used the Scaled Gradient Projection (SGP) method with Barzilai and Borwein rules for the step length selection. Regarding the TGV regularization, the RESPOND method [22] has been used. We remark that RESPOND has been originally proposed for the restoration of images corrupted by Poisson noise, by using Directional Total Generalized Variation regularization. It has been adapted here to deal with TGV-based restoration of images under Gaussian noise. The MATLAB implementation for Poisson noise is available on GitHub.

The tolerance $Tol$ in the outer loop of MULTI in Algorithm 1 step 7 was ${10}^{-1}$, while the maximum number of iterations was 20. Regarding the NP method, a tolerance of ${10}^{-5}$ was used and the maximum number of iterations was 1000.

The size of the neighborhood ${\mathcal{N}}_i$ in (8) was $5\times 5$ pixels for all tests except for Galaxy, where a $3\times 3$ neighborhood was used.

The values of the parameter $ϵ$ in (6), used in the various tests, are in the range $[{10}^{-4},{10}^{-3}]$. In order to compare all the algorithms at their best performance, the values $ϵ$ used in each test are reported in Table 1, where we observe that the value of $ϵ$ is proportional to the noise level. The parameter $ϵ$ represents a threshold and, in general, should have a small value when compared to the non null values of ${\mathcal{S}}_i$. We note that at the cost of adjusting a single parameter $ϵ>0$, it is possible to achieve point-wise optimal regularization.

Table 2, Table 3, Table 4 and Table 5 report the numerical results for all the test problems. The last column of the tables shows the used values of the regularization parameter for TIKH, TV and TGV while, for MULTI, it reports the norm of the regularization parameters vector $({\lambda }_1,\dots ,{\lambda }_N)$ computed by Algorithm 1. Column 7 shows the number of RESPOND, ADMM and SGP iterations for TGV, TV and TIKH, respectively. For the MULTI algorithm, Column 7 shows the number of outer iterations and NP iterations in parenthesis.

In Table 2, Table 3, Table 4 and Table 5, we observe that TGV always outperforms TIKH and TV in terms of accuracy, since it has greater MSSIM and ISNR values and smaller values of RE. Therefore, in Figure 1, Figure 2, Figure 3 and Figure 4 we represent the images obtained by MULTI and TGV in the out-of-focus case, with $\delta ={10}^{-2}$ as this is a very challenging case.

Moreover, from the images provided in Figure 5, it can be observed that MULTI method better preserves the local characteristics of the image, avoiding flattening the shooth areas and optimally preserving the sharp contours. We observe that a smooth area like the cheek is better represented by MULTI, avoiding the presence of the staircasing effect. Moreover, an area with strong contours, such as the teeth and the eyes, is better depicted.

The regularization parameters computed by MULTI are represented in Figure 6, the adherence of the regularization parameters’ values to the image content is clear, showing larger values in areas where the image is flat and smaller values where there are pronounced gradients (edges). The parameters range automatically adjusts according to the different test types.

Finally we show in Figure 7 an example of the algorithm behaviour, reporting the history of regularization parameter norm, relative error and residual norm in the case of leopard test with out-of-focus blur and noise level $\delta ={10}^{-2}$. The relative error flattens after a few iterations, and the same behaviour can be observed in all the other test. Therefore we used a large tolerance value ($Tol={10}^{-1}$) in the outer loop of Algorithm 1.

4. Conclusions

Our work is inspired by the uniform multi-penalty approach developed in [18] and proposes the pixel-wise regularization model to tackle the significant task of image restoration, concentrating on the elimination of degradation originating from blur and noise. The iterative method analized in [18] is adapted an made efficient for image restoration problems, by introducing a modified version of the Newton Projection method, tailored for multi-penalty scenarios.

The numerical results confirm the algorithm’s proficiency in eliminating noise and blur while concurrently preserving the edges of the image. These outcomes encourage us to extend the application of this algorithm to various imaging problems.