An Automatic Pixel-Wise Multi-Penalty Approach to Image Restoration

This work tackles the problem of image restoration, a crucial task in many fields of applied sciences, focusing on removing degradation caused by blur and noise during the acquisition process. Drawing inspiration from the multi-penalty approach based on the Uniform Penalty principle, discussed in previous work, here we develop a new image restoration model and an iterative algorithm for its effective solution. The model incorporates pixel-wise regularization terms and establishes a rule for parameter selection, aiming to restore images through the solution of a sequence of constrained optimization problems. To achieve this, we present a modified version of the Newton Projection method, adapted to multi-penalty scenarios, and prove its convergence. Numerical experiments demonstrate the efficacy of the method in eliminating noise and blur while preserving the image edges.


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] Au + e = b (1) where b ∈ R M and u ∈ R N , respectively, are vectorized forms of the observed m x × m y image and the exact n x × n y image to be restored, A ∈ R M×N is the linear operator modeling the imaging system, and e represents Gaussian white noise with mean zero and standard deviation σ.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 u from b.
In a variational framework, image restoration can be reformulated as a constrained optimization problem of the form min u 1 2 Au − b 2 + R(u) whose objective function contains a L 2 -based term, imposing consistency of the model with the data, and a regularization term R(u), forcing the solution to satisfy some a priori properties.Here and henceforth, the symbol • 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 Ω 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 λ i are regularization parameters and the ψ i (u) are regularization functions for i = 1, . . ., 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 λ i , i = 1, . . ., 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: R(u) = λψ(u). ( 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 ψ(u) 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 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] proposed a modified discrepancy principle for multi-penalty regularization and provided 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 ≤ q < 2 and 2 ≤ p < ∞.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 extend the methodology developed in [18] to image restoration problems and we perform a comparative analysis with state-of-the-art regularization methods for this application.We propose to find an estimate u * of u satisfying where is a positive scalar and L ∈ R N×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.We formulate an iterative algorithm for computing the solution (u * , λ * ) of ( 6), where λ * = (λ * 1 , . . ., λ * N ) T .Once the regularization parameters are set in every inner iteration, 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 multipenalty regularization.Consequently, the convergence of the modified NP method can be established.
The principal contributions of this work are summarized as follows: 1.
We propose a variational pixel-wise regularization model tailored for image restoration and derived from the theoretical model developed in [18].2.
We devise an algorithm capable of effectively and efficiently solving the proposed model.

3.
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.

Materials and Methods
In this section, we present the iterative algorithm that generates the sequence (u (k) , λ (k) ) converging to the solution (u * , λ * ) in (6).
Starting from an initial guess u (0) taken as the observed image b, the correspondent initial guess of the regularization parameters is computed as: where and N i is a neighborhood of size R × R, (with R odd and R ≥ 1) of the i-th pixel with coordinates (x i , y i ).
The successive terms (u (k+1) , λ (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: λ (0) ∈ R N , Tol, A, b, , Output: (u * , λ * ) Algorithm 1 is well defined, and we experimentally observe a converging behaviour.Its formal convergence proof is obtained in [18] (theorem 3.4) for the case R = 1, since in this case, Algorithm 1 corresponds to UPenMM.Otherwise, to preserve convergence, we should introduce a correction as proposed in the Generalized Uniform Penalty method (GUPenMM) [18].However, even without this correction, we obtained good-quality results and we prefer here to investigate Algorithm 1 because, in the case of large scale image restoration problems, it is much more convenient from a computational point of view.Moreover, we verified that the results obtained with such a correction are qualitatively comparable with those given by Algorithm 1, as the human eye cannot distinguish differences 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 J (k) (u) the function to be minimized at step 3 in Algorithm 1: and by g its gradient ∇J (k) (u) where the iteration index k has been omitted for easier notation.Moreover, let g I denote the reduced gradient: where I(u) is the set of indices [20]: and ε is a small positive parameter.The Hessian matrix ∇ 2 J (u) has the form where Λ (k) is the diagonal matrix with diagonal elements λ N .A general iteration of the proposed NP-like method has the form: where p ( ) is the search direction, α ( ) is the steplength and [•] + denotes the projection on the positive orthant.At each iteration , the computation of p ( ) requires the solution of the linear system where H (k) is the following approximation to ∇ 2 J (k) (u) Under periodic boundary conditions, H (k) is a BCCB matrix and system (11) can be efficiently solved in the Fourier space by using Fast Fourier Transforms.Therefore, despite its simplicity, the BCCB approximation H (k) is efficient, since it allows to solve the linear system in O(N log 2 N) operations, and effective, as is shown by the numerical results.Finally, given the solution d ( ) of ( 11), the search direction p ( ) is obtained as The step length α ( ) is computed with the variation of the Armijo rule discussed in [20] as the first number of the sequence where u ( ) (β m ) = [u ( ) − β m u ( ) ] + and η ∈ (0, 1  2 ).We observe that the approximated Hessian H (k) is constant for each inner iteration and it is positive definite, then it satisfies Following that, 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.

Numerical Experiments
All the experiments were performed under Windows 10 and MATLAB R2021a running on a desktop (Intel(R) Core(TM) i5-8250CPU@1.60GHz).Quantitatively, we evaluated 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 û is the restored image, u is the reference image, and b is the blurred, noisy image.Four reference images were used in the experiments: galaxy, mri, leopard, and elaine, shown in Figures 1-4.The first three images have size 256 × 256, while the elaine image is 512 × 512.In order to define the test problems, each reference image was convolved with two PSFs corresponding 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     We compared the proposed pixel-wise multi-penalty regularization model (MULTI) with some commonly used state-of-the-art methods based on a variational approach.In particular, we considered the Tikhonov method (TIKH) [22], the Total Variation (TV) [2], and Total Generalized Variation (TGV) [3] regularization with nonnegative constraints.Tikhonov and TV regularization are quite classic regularization terms.It is well known that Tikhonov regularization tends to make images overly smooth and it fails to preserve sharp edges.On the contrary, TV regularization better preserves sharp edges but often produces staircase effects.TGV has been recently proposed to overcome the drawbacks of Tikhonov and TV regularization such as blurring and the staircasing effect.Therefore, we compared MULTI with TGV in order to demonstrate the capacity of MULTI to preserve sharp features as well as smooth transition variations.
In our numerical experiments, the regularization parameter values for TIKH, TV, and TGV were 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 [23].Regarding the TGV regularization, the RESPOND method [24] was 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 at the url https://github.com/diserafi/respond(accessed on 18 September 20).
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 N i in (8) was 5 × 5 pixels for all tests except for galaxy, where a 3 × 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 S i .We note that at the cost of adjusting a single parameter > 0, it is possible to achieve point-wise optimal regularization.Table 1.Values for the parameter in (7).
Tables 2-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 (λ 1 , . . ., λ 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 Tables 2-5, we note that for the case of Gaussian blur, MULTI consistently achieves the best results, as highlighted in bold.This is evident from its higher MSSIM and ISNR values and lower RE values.However, for the Out-of-focus case, there are three instances (representing 25%) where TGV exhibits superior error parameters.Furthermore, our observations indicate that TGV consistently outperforms both TIKH and TV in terms of accuracy.Therefore, in Figures 1-4 we only represent the images obtained by MULTI and TGV in the out-of-focus case, with δ = 10 −2 , as this is a very challenging case.
The strength of MULTI is evident by observing some details of the reconstructed images.In Figures 5-8 we show some cropped details of the original images and compare it with MULTI and TGV reconstructions.Figure 5 shows a detail of the galaxy with a few stars over a dark background.In this case the image sparsity is better preserved by MULTI. Figure 6 shows the galaxy centre: it is a smooth area which is well recovered by MULTI while TGV shows staircasing.In Figure 7 we observe that the leopard's whiskers and fur spots are better reproduced by MULTI.Moreover, from the images provided in Figure 8, it can be observed that MULTI method better preserves the local characteristics of the image, avoiding flattening the smooth areas and optimally preserving the sharp contours.We observe that a smooth area such as the cheek is better represented by MULTI, avoiding the staircase effect.Moreover, an area with strong contours, such as the teeth and the eyes, is better depicted.In summary, these examples show the good capacity of MULTI to preserve the different image structures, narrow peaks, and smooth areas by using local regularization parameters that are inversely proportional to the local curvature approximated by the discrete Laplacian.Finally, we show in Figure 10 an example of the algorithm behavior, reporting the history of the regularization parameter norm, relative error, and residual norm (top row), in the case of leopard test with out-of-focus blur and noise level δ = 10 −2 .In the bottom row we show the decreasing behavior of the objective function and projected gradient norm.The relative error flattens after a few iterations, and the same behavior can be observed in all the other tests.Therefore, we used a large tolerance value (Tol = 10 −1 ) in the outer loop of Algorithm 1, making it possible to obtain good regularization parameters and accurate restorations in a few outer iterations.We observe that even in the most difficult case, Table 4 row 12, the total computation time is 285 s, proving the algorithm efficiency.

Conclusions
Despite the interest of recent literature on multi-penalty regularization, its drawback lies in the difficult computation of the regularization parameters.Our work proposes the pixel-wise regularization model to tackle the significant task of image restoration, concentrating on eliminating degradation originating from blur and noise.We show that multi-penalty regularization can be realized by an algorithm that is able to compute efficiently and automatically a large number of regularization parameters.The numerical results confirm the algorithm's proficiency in eliminating noise and blur while concurrently preserving the edges of the image.Such an approach can be exploited in different real-world imaging applications, such as computed tomography, super-resolution, and biomedical imaging in general.Finally, further analyses of the properties of the proposed algorithm will be the subject of future works.
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 δ = e / b .The values δ = 2.5 × 10 −2 , 10 −2 , 5 × 10 −3 were used.

Figure 8 .
Figure 8. Elaine test problem: out-of-focus blur; δ = 10 −2 .A detail of the original image (left), MULTI restoration (centre), and TGV restoration (right).Red arrows highlight the different image features.The regularization parameters computed by MULTI are represented in Figure 9; 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 range of the parameters automatically adjusts according to the different test types.

Table 2 .
Numerical results for the galaxy test problem.Column Iters shows the number of RESPOND, ADMM and SGP iterations for TGV, TV and TIKH, respectively.For the MULTI algorithm, it shows the number of outer iterations and NP iterations in parenthesis.The best results are highlighted in bold.Column λ 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 (λ 1 , . . ., λ N ) computed by Algorithm 1.

Table 3 .
Numerical results for the mri test problem.Column Iters shows the number of RESPOND, ADMM, and SGP iterations for TGV, TV, and TIKH, respectively.For the MULTI algorithm, it shows the number of outer iterations and NP iterations in parenthesis.Column λ 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 (λ 1 , . . ., λ N ) computed by Algorithm 1.The best results are highlighted in bold.

Table 4 .
Numerical results for the leopard test problem.Column Iters shows the number of RE-SPOND, ADMM, and SGP iterations for TGV, TV, and TIKH, respectively.For the MULTI algorithm, it shows the number of outer iterations and NP iterations in parenthesis.Column λ 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 (λ 1 , . . ., λ N ) computed by Algorithm 1.The best results are highlighted in bold.

Table 5 .
Numerical results for the elaine test problem.Column Iters shows the number of RESPOND, ADMM, and SGP iterations for TGV, TV, and TIKH, respectively.For the MULTI algorithm, it shows the number of outer iterations and NP iterations in parenthesis.Column λ 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 (λ 1 , . . ., λ N ) computed by Algorithm 1.The best results are highlighted in bold.