1 Introduction
Image restoration is the task of recovering a sharp image starting from a blurred and noisy observation . In this work, we consider a degradation model of the form
(1) 
where
are vectorized images,
is the linear blur operator and is an additive noise vector. A possible strategy to overcome the illposedeness of the linear system in (1) is to reformulate the problem. We thus look for , estimate of the original , which solves a wellposed problem. In the variational approach, is the minimizer of a cost functional . In formula,(2) 
The functionals and are commonly referred to as the regularization and the data fidelity term, respectively. While encodes prior information on the desired image , is a data term which measures the ‘distance’ between the given image and after the action of the operator with respect to some norm corresponding to the noise statistics in the data, cf., e.g., stuart . The regularization parameter controls the tradeoff between the two terms.
In this paper, we consider an Additive White Gaussian Noise (AWGN) corrupting the blurred image , i.e. , where is the
dimensional identity matrix. It is well known that, in presence of AWGN, a suitable choice for
is the socalled L fidelity term, reading as,(3) 
A popular choice for the regularization term is given by the TV seminorm ROF ,
(4) 
where denotes the discrete gradient of image at pixel ,
with linear operators representing finite difference discretizations of the firstorder horizontal and vertical partial derivatives, respectively.
Coupling the L data term with the TV regularizer leads to one of the most widespread variational models for image restoration problem, the TVL (or ROF) model,
(5) 
The global perspective of the TVL model does not allow to diversify the action of the regularizer on regions of the image presenting different properties. In siam ; ncmip ; cmbbe ; vip , the authors have proposed spacevariant regularization term based on statistical assumptions on the distribution of the norm of the gradients and on the gradients themselves.
In this paper, we propose a spacevariant anisotropic extension of the TV regularizer in (4) which, as it will be illustrated in Sect. 2, comes from the a priori assumption that the gradients of the target image distribute locally according to a Bivariate Laplace Distribution (BLD). The proposed BLTV regularizer takes the form
(6)  
(7) 
where is a positive definite diagonal matrix and is the rotation matrix corresponding to the angle . Mathematically,
(8) 
We denote by the maps of the parameters defining the local distributions. Hence, the proposed BLTVL variational restoration model reads as
(9) 
The free parameters defining the BLTV regularizer in (6)–(8) hold the potential for effectively modelling local image properties, in particular driving in a suitable adaptive manner the strength and direction of nonlinear TVdiffusion. In Figs.1
(a)(b) the red ellipses represent TVdiffusion strengths along all possible directions at few sample pixel locations for TV and BLTV regularizers. It is evident how for classical TV such ellipses turn out to be circles (isotropy) of constant radius (spaceinvariance), whereas our BLTV regularizer allows for ellipses (anistropy) of different size (spacevariance). In practice, such flexibility of BLTV can be (and will be) exploited to diffuse in different ways in regions exhibiting different properties: for instance, strong isotropic diffusion in homogeneous regions, strongly anisotropic diffusion in regions characterized by a dominant edge direction. In Figs.
1 (c)(d) we show the level curves of the TV regularizer and one among the infinity of possible configuration of the level curves of the BLTV regularizer, respectively, revealing once again the flexibility of the proposed regularizer.(a)  (b)  (c)  (d) 
Together with the variational model in (9), we also propose an efficient and robust procedure for automatically estimating the parameter maps , , based on a Maximum Likelihood (ML) approach. The parameter maps can be updated along the iterations of the Alternating Direction Methods of Multipliers (ADMM), which is the algorithm adopted here to solve the minimization problem. Notice that the convexity of the BLTVL model ensures the convergence of the ADMM.
2 Deriving the model via MAP
Applying the Maximum A Posteriori
(MAP) estimation approach to image restoration consists in computing the restored image as a global maximizer of the posterior probability Pr
of the unknown target image given the observation , namely:(10) 
where, after applying the Bayes’ rule, we drop the evidence term Pr and extracted of the objective function. The two terms Pr and Pr in (10) are referred to as the prior and the likelihood. The likelihood term associated with AWGN corruption takes the form
(11) 
where
denotes the AWGN standard deviation and
is a normalization constant. For what concerns the prior, a common choice is to model the unknown image as a Markov Random Field (MRF) such that the image can be characterized by its Gibbs prior distribution, whose general form is:(12) 
where is the MRF parameter, is the set of all cliques (a clique is a set of neighboring pixels) for the MRF, is the potential function defined on the clique and is a normalization constant. Choosing as potential function at the generic th pixel the magnitude of the discrete gradient at the same pixel, i.e. for any , the Gibbs prior in (12) reduces to the TV prior which, plugged into (10), yields the popular TV regularizer. We remark that the TV prior corresponds to choosing
(13) 
Here, we propose a spacevariant anisotropic generalization of the local prior above. More in detail, we assume that the discrete gradient at any pixel distributes according to a spacevariant BLD, such that our prior reads as
(14)  
where is the normalization constant and are defined in (8).
3 Parameter estimation via ML
In order to make the introduction of the proposed regularizer actually useful, an efficient, robust, and automatic procedure for the estimation of the parameter maps identifying all the local BLDs has to be proposed as well. To this aim, we resort to the ML approach. Consider a set of 2dimensional samples drawn from a BLD with parameters . Here, the samples play the role of image gradients at pixels of a neighborhood of radius centered at a generic pixel . Assuming independence of the samples, according to the definition of BLD, the likelihood function is defined by
(15) 
where, clearly, and depend on  see (8). The goal here is to find
(16)  
Imposing a firstorder optimality condition with respect to leads to the following closedform estimation formulas:
(17) 
Substituting the expressions in (17) into the objective function (16), we thus obtain the simplified minimization problem in the only variable :
(18) 
4 Admm
In order to solve numerically the proposed image restoration model (9), we use an ADMMbased algorithm  see boyd . We first introduce two auxiliary variables , and rewrite the model in the equivalent linearly constrained form:
(20)  
where the spacevariant matrices , can be estimated via the ML procedure described in Sect. 3 based only on the observed image (i.e. as a preliminary preprocessing step) or also updated along the ADMM iterations. We define the augmented Lagrangian functional:
(21)  
where are scalar penalty parameters and , are the vectors of Lagrange multipliers associated with the given linear constraints. The solution of problem (20) is a saddle point for in (21), see, e.g., boyd . Hence, we can alternate a minimization step with respect to with a maximization step with respect to . Mathematically,
(22)  
(23)  
(24)  
(25)  
(26) 
The solution of the primal subproblem (22
) can be efficiently computed by means of standard linear Fast Fourier Transform (FFT) solvers. The subproblem (
23) can be solved in closedform by following (ncmip, , Section 3). Finally, the subproblem (24) can be solved by computing efficiently the proximal operator of the anisotropic norm, for which the proof of (siam, , Proposition 6.3) can be easily adapted. The regularization parameter is updated along the iterations so as to fulfill the global discrepancy principle as described in ape . We refer the reader also to cmbbe ; tvp for more details on the numerical solution of the algorithm.5 Experimental results
In this section, we evaluate the performance of the proposed BLTVL restoration model compared with the baseline TVL model, also solved by ADMM. The stopping criteria of the ADMM for both models are defined based on the number of iterations as well as on the iterates relative change, i.e. we stop iterating as soon as
(27) 
The quality of the restored images is measured by means of the Improved SignaltoNoise Ratio , with denoting the original uncorrupted image, and of the StructuralSimilarityIndex (SSIM) ssim .
We consider the test image brain in Fig.2 (a) (570 430) and the test image abdomen in Fig.4 (a) (350 480) with pixel values between 0 and 255, synthetically corrupted by spaceinvariant blur with Gaussian kernel of parameters band
= 9, sigma
= 2, and by AWGN of different levels  see, e.g., Fig.2 (b) and Fig.4 (b). The parameter maps are computed at the beginning starting from the observed image and then updated every 300 iterations based on the current iterate. The radius of the neighborhoods used for the local parameter estimation has been set equal to 8 and 5 for the brain and abdomen test images, respectively. The reconstructions of brain for and of abdomen for via TVL and BLTVL models are shown in Figs.2(c)(d) and Figs.4(c)(d), respectively.
(a)  (b)  (c)  (d) 
(e)  (f)  (g)  (h) 
From a visual inspection, the restoration via BLTVL seems to be more neat and less cartooned than the TV reconstructions. As reported in Tab.1, the ISNR and SSIM values for the two test images and for different noise levels obtained by the BLTVL model outperform the ones reached by the TVL model.
The final parameter maps computed by BLTVL are shown in Figs. 35.
brain  abdomen  
TVL  BLTVL  TVL  BLTVL  TVL  BLTVL  TVL  BLTVL  
ISNR  4.27  5.52  5.00  6.67  3.76  4.66  6.10  6.77 
SSIM  0.87  0.88  0.83  0.85  0.78  0.80  0.74  0.76 
(a)  (b)  (c)  (d) 
(e)  (f)  (g)  (h) 
Computational times. We tested the joint parameter estimation + reconstruction model on a standard laptop with inbuilt MATLAB software, version 2016b. As far as the ML parameter estimation of the parameter maps procedure is concerned, we notice that the update (17) is explicit, thus very cheap, whereas the computation of in (18) requires the solution of an optimisation problem. We solve the problem by linesearching upon a suitable discretization of the parameter space. For Fig. 2, the ML parameter estimation procedure took 8.45 secs.
The ADMM algorithmic substeps with automatic parameter update every 300 iterations computes the numerical solution in 381 secs for the highresolution image in Fig. 2. A possible way to accelerate the speed of the algorithm would be the computation the parameter maps only in terms of the given image and not along the iterations, although of course that would render a less accurate result.
6 Conclusions and future works
We presented a new spacevariant anisotropic regularization term for image restoration based on the a priori statistical assumption that the gradients of the unknown target image distribute locally according to spacevariant bivariate Laplace distributions. The high flexibility of the proposed regularizer together with the presented ML parameters estimation procedure and ADMMbased minimization algorithm yield a very effective and efficient approach. Preliminary experiments on images corrupted by blur and AWGN strongly indicate that the proposed variational model achieves highquality restorations and, in particular, outperforms by far the results obtained by classical TVregularized restoration models. Coupling the proposed regularizer with other fidelity terms suitable for dealing with noises other than Gaussian  such as, e.g., Laplace, Poisson and mixed PoissonGaussian MIXED  is a matter being studied.
Bibliography

(1)
Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J.: Distributed optimization and statistical learning via the alternating direction method of multipliers: In: Foundations and Trends in Machine Learning, 3 (2011).
 (2) Calatroni, L., Lanza, A., Pragliola, M., Sgallari, F.: A flexible spacevariant anisotropic regularisation for image restoration with automated parameter selection. In: SIAM Journal of Imaging Sciences, 12 (2019).
 (3) Calatroni, L., Lanza, A., Pragliola, M., Sgallari, F.: Adaptive parameter selection for weightedTV image reconstruction problems. In: Journal of Physics: Conference Series, NCMIP 2019, to appear (2019).
 (4) He, C., Hu, C., Zhang, W. and Shi, B.: A Fast Adaptive Parameter Estimation for Total Variation Image Restoration. In: IEEE Transactions on Image Processing, 23 (2014).
 (5) Lanza, A., Morigi, S., Pragliola, M., Sgallari, F.: Spacevariant generalised gaussian regularisation for image restoration. In: Computational Methods in Biomechanics and Biomedical Engineering: Imaging and Visualization, 13 (2018).
 (6) Lanza, A., Morigi, S., Pragliola, M., Sgallari, F.: Spacevariant TV regularization for image restoration. In: VipIMAGE 2017, J. M. R. Tavares and R. Natal Jorge, eds., Cham, Springer International Publishing, (2018).
 (7) Lanza, A., Morigi, S., Sgallari, F.: Constrained TV model for image restoration. In: Journal of Scientific Computing, 68 (2016).
 (8) Lanza, A., Morigi, S., Sgallari, F., Wen, Y.W.: Image restoration with Poisson  Gaussian mixed noise. In: Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization, 2 (2014).
 (9) Rudin, L.I., Osher, S., Fatemi, E.: Nonlinear total variation based noise removal algorithms. In: Physica D: Nonlinear Phenomena, 60 (1992).
 (10) Stuart, A.M.: Inverse problems: a Bayesian perspective. In: Acta Numerica, 19 (2010).
 (11) Zhou, W., Bovik, A., Sheikh, H., Simoncelli, E.: Image quality assessment: From error visibility to structural similarity. In: IEEE Transactions on Image Processing, 13 (2004).
Comments
There are no comments yet.