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
are vectorized images,is the linear blur operator and is an additive noise vector. A possible strategy to overcome the ill-posedeness of the linear system in (1) is to reformulate the problem. We thus look for , estimate of the original , which solves a well-posed problem. In the variational approach, is the minimizer of a cost functional . In formula,
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 trade-off 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 foris the so-called L fidelity term, reading as,
A popular choice for the regularization term is given by the TV semi-norm ROF ,
where denotes the discrete gradient of image at pixel ,
with linear operators representing finite difference discretizations of the first-order 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 TV-L (or ROF) model,
The global perspective of the TV-L 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 space-variant 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 space-variant 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
where is a positive definite diagonal matrix and is the rotation matrix corresponding to the angle . Mathematically,
We denote by the maps of the parameters defining the local distributions. Hence, the proposed BLTV-L variational restoration model reads as
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 non-linear TV-diffusion. In Figs.1
(a)-(b) the red ellipses represent TV-diffusion 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 (space-invariance), whereas our BLTV regularizer allows for ellipses (anistropy) of different size (space-variance). 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.
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 BLTV-L 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 Prof the unknown target image given the observation , namely:
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
denotes the AWGN standard deviation andis 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:
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
Here, we propose a space-variant anisotropic generalization of the local prior above. More in detail, we assume that the discrete gradient at any pixel distributes according to a space-variant BLD, such that our prior reads as
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 2-dimensional 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
where, clearly, and depend on - see (8). The goal here is to find
Imposing a first-order optimality condition with respect to leads to the following closed-form estimation formulas:
In order to solve numerically the proposed image restoration model (9), we use an ADMM-based algorithm - see boyd . We first introduce two auxiliary variables , and rewrite the model in the equivalent linearly constrained form:
where the space-variant matrices , can be estimated via the ML procedure described in Sect. 3 based only on the observed image (i.e. as a preliminary pre-processing step) or also updated along the ADMM iterations. We define the augmented Lagrangian functional:
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,
The solution of the primal sub-problem (22
) can be efficiently computed by means of standard linear Fast Fourier Transform (FFT) solvers. The sub-problem (23) can be solved in closed-form by following (ncmip, , Section 3). Finally, the sub-problem (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 BLTV-L restoration model compared with the baseline TV-L 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
The quality of the restored images is measured by means of the Improved Signal-to-Noise Ratio , with denoting the original uncorrupted image, and of the Structural-Similarity-Index (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 space-invariant 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 TV-L and BLTV-L models are shown in Figs.2(c)-(d) and Figs.4(c)-(d), respectively.
From a visual inspection, the restoration via BLTV-L 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 BLTV-L model outperform the ones reached by the TV-L model.
The final parameter maps computed by BLTV-L are shown in Figs. 3-5.
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 line-searching upon a suitable discretization of the parameter space. For Fig. 2, the ML parameter estimation procedure took 8.45 secs.
The ADMM algorithmic sub-steps with automatic parameter update every 300 iterations computes the numerical solution in 381 secs for the high-resolution 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 space-variant 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 space-variant bivariate Laplace distributions. The high flexibility of the proposed regularizer together with the presented ML parameters estimation procedure and ADMM-based 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 high-quality restorations and, in particular, outperforms by far the results obtained by classical TV-regularized 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 Poisson-Gaussian MIXED - is a matter being studied.
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 space-variant 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 weighted-TV 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.: Space-variant 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.: Space-variant 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).