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]
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 illposed 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
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 \xb7\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
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 multipenalty 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 singlepenalty regularization, also known as Tikhonovlike regularization, which corresponds to the choice
$p=1$:
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 multipenalty regularization, with a significant number of researchers focusing on scenarios involving two penalty terms. Notably, the widelyused elastic regression in Statistics serves as an example of a multipenalty 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 multiparameter version of the Lcurve 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
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 multipenalty regularization and provide theoretical background for this a posteriori rule.
Works such as [
10,
11,
12,
13,
14] also explore multipenalty 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 twopenalty regularization, incorporating
${L}_{0}$ and
${L}_{2}$ penalty terms, to tackle nonlinear illposed 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 spacevariant generalized Gaussian regularization approach for image restoration, emphasizing its applicative potential.
In [
18], a multipenalty pointwise 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
where
$\u03f5$ 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 pixelwise 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 singlepenalty image restoration under Poisson noise and it is adapted here to the context of multipenalty 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 pixelwise 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:
where
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 35 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},\u03f5$, Output: $({\mathbf{u}}^{*},{\mathsf{\lambda}}^{*})$

 1:
Set $k=0$.
 2:
repeat
 3:
$\mathbf{u}}^{(k+1)}=arg\underset{\mathbf{u}\in \Omega}{min}\phantom{\rule{0.277778em}{0ex}}\left(\right)open="\{"\; close="\}">\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$
 4:
${\mathcal{S}}_{i}^{(k+1)}=\underset{i\in {\mathcal{N}}_{i}}{max}{\left(L{\mathbf{u}}^{(k+1)}\right)}_{i}^{2}),\phantom{\rule{2.em}{0ex}}i=1,\dots N$
 5:
${\lambda}_{i}^{(k+1)}=\frac{\parallel A{\mathbf{u}}^{(k+1)}{\mathbf{b}\parallel}^{2}}{N({\mathcal{S}}_{i}^{(k+1)}+\u03f5)},\phantom{\rule{0.277778em}{0ex}}i=1,\dots ,N$
 6:
$k=k+1$
 7:
until$\parallel {\mathsf{\lambda}}^{\left(k\right)}{\mathsf{\lambda}}^{(k1)}\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:
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:
where
$\mathcal{I}\left(\mathbf{u}\right)$ is the set of indices [
20]:
with
and
$\overline{\epsilon}$ is a small positive parameter.
The Hessian matrix
${\nabla}^{2}\mathcal{J}\left(\mathbf{u}\right)$ has the form
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 NPlike method has the form:
where
${\mathbf{p}}^{(\ell )}$ is the search direction,
${\alpha}^{(\ell )}$ is the steplength and
${[\xb7]}^{+}$ denotes the projection on the positive orthant.
At each iteration
ℓ, the computation of
${\mathbf{p}}^{(\ell )}$ requires the solution of the linear system
where
${H}^{\left(k\right)}$ is the following approximation to
${\nabla}^{2}{\mathcal{J}}^{\left(k\right)}\left(\mathbf{u}\right)$
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
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
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
Then, the results given in [
19] for singlepenalty image restoration under Poisson noise, can be applied here to prove the convergence of the NLlike iterations to critical points.
The stopping criteria for the NPlike 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) i58250CPU@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:
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 RestoreTool [
1], and an outoffocus 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\xb7{10}^{2},{10}^{2},5\xb7{10}^{3}$ were used.
We compare the proposed pixelwise multipenalty 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 TVbased 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 TGVbased 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
$\u03f5$ 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
$\u03f5$ used in each test are reported in
Table 1, where we observe that the value of
$\u03f5$ is proportional to the noise level. The parameter
$\u03f5$ 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
$\u03f5>0$, it is possible to achieve pointwise 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 outoffocus 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 outoffocus 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.