## Space-adaptive anisotropic bivariate Laplacian regularization for image restoration

Luca Calatroni, Alessandro Lanza, Monica Pragliola and Fiorella Sgallari

Abstract In this paper we present a new regularization term for variational image restoration which can be regarded as a space-variant anistropic extension of the classical Total Variation (TV) regularizer. The proposed regularizer comes from the statistical assumption that the gradients of the unknown target image distribute locally according to space-variant bivariate Laplacian distributions. The high flexibility of the proposed regularizer holds the potential for effectively modelling local image properties, in particular driving in an adaptive manner the strength and direction of non-linear TV-diffusion. The free parameters of the regularizer are automatically set -and, eventually, updated - based on a robust Maximumum Likelihood estimation procedure. A minimization algorithm based on the Alternating Direction Method of Multipliers is presented for the efficient numerical solution of the proposed variational model. Some experimental results are reported which demonstrate the high-quality of restorations achievable by the proposed model, in particular with respect to classical TV-regularized models.

Keywords : Image restoration, variational methods, ADMM, TV regularization.

Luca Calatroni

CMAP, CNRS, École Polytechnique, Institut Polytechnique de Paris, Palaiseau, 91128, Route de Saclay, France e-mail: luca.calatroni@polytechnique.edu

Alessandro Lanza

University of Bologna, Piazza di Porta San Donato 5, Bologna, Italy e-mail: alessandro.lanza2@ unibo.it

Monica Pragliola

University of Bologna, Piazza di Porta San Donato 5, Bologna, Italy e-mail: monica.pragliola2@ unibo.it

Fiorella Sgallari

University of Bologna, Piazza di Porta San Donato 5, Bologna, Italy e-mail: fiorella.sgallari@ unibo.it

## 1 Introduction

Image restoration is the task of recovering a sharp image u starting from a blurred and noisy observation g . In this work, we consider a degradation model of the form

$$g = K u + e,$$

where g ; u 2 R n are vectorized images, K 2 R n GLYPH&lt;2&gt; n is the linear blur operator and e 2 R n 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 u GLYPH&lt;3&gt; , estimate of the original u , which solves a well-posed problem. In the variational approach, u GLYPH&lt;3&gt; is the minimizer of a cost functional J ' u K ; ; g ' : R n ! R . In formula,

$$u ^ { * } \in \arg \min _ { u \in \mathbb { R } ^ { n } } \left \{ \mathfrak { D } ( u ; K, g ) \coloneqq \mathcal { R } ( u ) + \mu \mathcal { F } ( u ; K, g ) \right \}.$$

The functionals R and F are commonly referred to as the regularization and the data fidelity term, respectively. While R encodes prior information on the desired image u , F is a data term which measures the 'distance' between the given image g and u after the action of the operator K with respect to some norm corresponding to the noise statistics in the data, cf., e.g., [10]. The regularization parameter GLYPH&lt;22&gt; &gt; 0 controls the trade-off between the two terms.

In this paper, we consider an Additive White Gaussian Noise (AWGN) corrupting the blurred image Ku , i.e. e GLYPH&lt;24&gt; N ' 0 ; GLYPH&lt;27&gt; 2 I n ' , where I n is the n -dimensional identity matrix. It is well known that, in presence of AWGN, a suitable choice for F ' u K ; ; g ' is the so-called L2 fidelity term, reading as,

$$\mathcal { F } ( u ; K, g ) = L _ { 2 } ( u ; K, g ) = \frac { 1 } { 2 } \| K u - g \| _ { 2 } ^ { 2 }.$$

Apopular choice for the regularization term R ' u ' is given by the TV semi-norm [9],

$$\mathcal { R } ( u ) = T V ( u ) = \sum _ { i = 1 } ^ { n } \| ( \nabla u ) _ { i } \| _ { 2 } \,,$$

where 'r u ' i : = GLYPH&lt;0&gt; ' D u h ' i ; ' D u v ' i GLYPH&lt;1&gt; T 2 R 2 denotes the discrete gradient of image u at pixel i , with D h ; D v 2 R n GLYPH&lt;2&gt; n linear operators representing finite difference discretizations of the first-order horizontal and vertical partial derivatives, respectively. Couplingthe L2 data term with the TV regularizer leads to one of the most widespread variational models for image restoration problem, the TV-L2 (or ROF) model,

$$u ^ { * } \in \arg \min _ { u \in \mathbb { R } ^ { n } } \left \{ \sum _ { i = 1 } ^ { n } \| ( \nabla u ) _ { i } \| _ { 2 } + \frac { \mu } { 2 } \| K u - g \| _ { 2 } ^ { 2 } \right \}.$$

The global perspective of the TV-L2 model does not allow to diversify the action of the regularizer on regions of the image presenting different properties. In [2, 3, 5, 6],

the authors have proposed space-variant regularization term based on statistical assumptions on the distribution of the ' 2-norm of the gradients and on the gradients themselves.

In this paper, we propose a space-variant anisotropic extension of the TV regularizer in (2) which, as it will be illustrated in Sect. 2, comes from the a priori assumption that the gradients of the target image u distribute locally according to a Bivariate Laplace Distribution (BLD). The proposed BLTV regularizer takes the form

$$BL T V ( u ; \lambda _ { 1 }, \lambda _ { 2 }, \theta ) = \sum _ { j } ^ { n } \| \Lambda _ { i } R _ { \theta _ { i } } ( \nabla u ) _ { i } \| _ { 1 }$$

$$\text{BLTV} ( u ; \lambda _ { 1 }, \lambda _ { 2 }, \theta ) & = \sum _ { i = 1 } ^ { n } \| \Lambda _ { i } R _ { \theta _ { i } } ( \nabla u ) _ { i } \| _ { 1 } & ( 3 ) \\ & = \sum _ { i = 1 } ^ { n } \left [ \lambda _ { i } ^ { ( 1 ) } | \langle r _ { i } ^ { ( 1 ) }, ( \nabla u ) _ { i } \rangle | + \lambda _ { i } ^ { ( 2 ) } | \langle r _ { i } ^ { ( 2 ) }, ( \nabla u ) _ { i } \rangle | \right ], & ( 4 )$$

where GLYPH&lt;3&gt; i is a 2 GLYPH&lt;2&gt; 2 positive definite diagonal matrix and R GLYPH&lt;18&gt; i is the rotation matrix corresponding to the angle GLYPH&lt;0&gt; GLYPH&lt;18&gt; i . Mathematically,

$$\Lambda _ { i } = \begin{pmatrix} \lambda _ { i } ^ { ( 1 ) } & 0 \\ 0 & \lambda _ { i } ^ { ( 2 ) } \end{pmatrix}, \quad R _ { \theta _ { i } } = \begin{pmatrix} \cos \theta _ { i } & \sin \theta _ { i } \\ - \sin \theta _ { i } & \cos \theta _ { i } \end{pmatrix} = \begin{pmatrix} r _ { i } ^ { ( 1 ) } \\ r _ { i } ^ { ( 2 ) } \end{pmatrix}.$$

We denote by GLYPH&lt;21&gt; 1 ; GLYPH&lt;21&gt; 2 ; GLYPH&lt;18&gt; 2 R n the maps of the parameters defining the local distributions. Hence, the proposed BLTV-L2 variational restoration model reads as

$$u ^ { * } \in \arg \min _ { u \in \mathbb { R } ^ { n } } \left \{ \text{BLTV} ( u ; \lambda _ { 1 }, \lambda _ { 2 }, \theta ) + \frac { \mu } { 2 } \| K u - g \| _ { 2 } ^ { 2 } \right \}.$$

The 3 n free parameters defining the BLTV regularizer in (3)-(5) 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 (6), we also propose an efficient and robust procedure for automatically estimating the parameter maps GLYPH&lt;21&gt; ' 1 ' , GLYPH&lt;21&gt; ' 2 ' , GLYPH&lt;18&gt; 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

Fig. 1: Space-invariant isotropic TV-diffusion (a), space-variant anisotropic BLTVdiffusion (b), level curves of TV (c) and BLTV (d) regularization terms.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00801v1-figure-c59c85e1-4acb-46a2-a268-18bd0aa4d2ab.jpg)

is the algorithm adopted here to solve the minimization problem. Notice that the convexity of the BLTV-L2 model ensures the convergence of the ADMM.

## 2 Deriving the model via MAP

Applying the MaximumAPosteriori (MAP)estimation approach to image restoration consists in computing the restored image as a global maximizer of the posterior probability Pr ' u j g ; K ' of the unknown target image u given the observation g , namely:

$$u ^ { * } \, \in \, \arg \max _ { u \in \mathbb { R } ^ { n } } \, \Pr ( u | g ; K ) = \arg \min _ { u \in \mathbb { R } ^ { n } } \, \left \{ \, - \ln \Pr ( g | u ; K ) - \ln \Pr ( u ) \, \right \} \,, \quad ( 7 )$$

where, after applying the Bayes' rule, we drop the evidence term Pr ' g ' and extracted GLYPH&lt;0&gt; ln of the objective function. The two terms Pr ' u ' and Pr ' g j u K ; ' in (7) are referred to as the prior and the likelihood . The likelihood term associated with AWGN corruption takes the form

$$\Pr ( g | u ; K ) = \prod _ { i = 1 } ^ { n } \, \frac { 1 } { \sqrt { 2 \pi } \sigma } \exp \left ( - \, \frac { ( K u - g ) _ { i } ^ { 2 } } { 2 \sigma ^ { 2 } } \, \right ) = W \exp \left ( - \, \frac { \| K u - g \| _ { 2 } ^ { 2 } } { 2 \sigma ^ { 2 } } \, \right ), \quad ( 8 )$$

where GLYPH&lt;27&gt; &gt; 0 denotes the AWGN standard deviation and W &gt; 0 is a normalization constant. For what concerns the prior, a common choice is to model the unknown image u as a Markov Random Field (MRF) such that the image can be characterized by its Gibbs prior distribution, whose general form is:

$$\Pr ( u ) = \prod _ { i = 1 } ^ { n } \Pr _ { i } ( u ) = \prod _ { i = 1 } ^ { n } z _ { i } \exp \left ( \, - \alpha \, V _ { c _ { i } } ( u ) \, \right ) = Z \exp \left ( \, - \alpha \, \sum _ { i = 1 } ^ { n } V _ { c _ { i } } ( u ) \, \right ), \quad ( 9 )$$

where GLYPH&lt;11&gt; &gt; 0 is the MRF parameter, f c i g n i = 1 is the set of all cliques (a clique is a set of neighboring pixels) for the MRF, V c i is the potential function defined on the clique c i and Z = ˛ n i = 1 z i is a normalization constant. Choosing as potential function at the generic i -th pixel the magnitude of the discrete gradient at the same pixel, i.e.

V c i = k'r u ' i k 2 for any i 2 f 1 ; : : : ; n g , the Gibbs prior in (9) reduces to the TV prior which, plugged into (7), yields the popular TV regularizer. We remark that the TV prior corresponds to choosing

$$\Pr _ { i } ( u ) = z _ { i } \exp \left ( - \alpha \| ( \nabla u ) _ { i } \| _ { 2 } \right ) \,.$$

Here, we propose a space-variant anisotropic generalization of the local prior Pr i ' u ' above. More in detail, we assume that the discrete gradient at any pixel i distributes according to a space-variant BLD, such that our prior reads as

$$\kappa _ { 1 } \quad \cdot & \quad \cdot & \quad \cdot \\ & \Pr ( u ) & = \prod _ { i = 1 } ^ { n } \Pr _ { i } ( u ) = \prod _ { i = 1 } ^ { n } z _ { i } \, \exp \left ( - \, \| \Lambda _ { i } R _ { \theta _ { i } } ( \nabla u ) _ { i } \| _ { 1 } \right ) \\ & = Z \, \exp \left ( - \sum _ { i = 1 } ^ { n } \| \Lambda _ { i } R _ { \theta _ { i } } ( \nabla u ) _ { i } \| _ { 1 } \right ), & & ( 1 0 ) \\ & \quad \cdot & \quad \cdot$$

where Z is the normalization constant and GLYPH&lt;3&gt; i ; R GLYPH&lt;18&gt; i 2 R 2 GLYPH&lt;2&gt; 2 are defined in (5).

Plugging the prior (10) and the likelihood (8) into (7) and neglecting the constant terms, the proposed BLTV-L2 variational model (6) is obtained.

## 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 GLYPH&lt;21&gt; ' 1 ' ; GLYPH&lt;21&gt; ' 2 ' ; GLYPH&lt;18&gt; 2 R n identifying all the local BLDs has to be proposed as well. To this aim, we resort to the ML approach. Consider a set of N 2-dimensional samples S : = f s 1 ; : : : ; s N g drawn from a BLD with parameters ' GLYPH&lt;21&gt; 1 ; GLYPH&lt;21&gt; 2 ; GLYPH&lt;18&gt; ' . Here, the samples play the role of image gradients at pixels of a neighborhood of radius r centered at a generic pixel i . Assuming independence of the samples, according to the definition of BLD, the likelihood function is defined by

$$\Pr ( \mathcal { S } | \lambda _ { 1 }, \lambda _ { 2 }, \theta ) = \left ( \frac { \lambda _ { 1 } \lambda _ { 2 } } { 4 } \right ) ^ { N } \exp \left ( - \sum _ { i = 1 } ^ { N } \left ( \lambda _ { 1 } | \langle r _ { 1 }, s _ { i } \rangle | + \lambda _ { 2 } | \langle r _ { 2 }, s _ { i } \rangle | \right ) \right ),$$

where, clearly, r 1 and r 2 depend on GLYPH&lt;18&gt; - see (5). The goal here is to find

$$( \lambda _ { 1 } ^ { * }, \lambda _ { 2 } ^ { * }, \theta ^ { * } ) & \in \arg \max _ { \lambda _ { 1 }, \lambda _ { 2 }, \theta } \Pr ( \mathcal { S } | \lambda _ { 1 }, \lambda _ { 2 }, \theta ) = \arg \min _ { \lambda _ { 1 }, \lambda _ { 2 }, \theta } - \ln \Pr ( \mathcal { S } | \lambda _ { 1 }, \lambda _ { 2 }, \theta ) \\ & = - N \ln \frac { \lambda _ { 1 } \, \lambda _ { 2 } } { 4 } + \sum _ { i = 1 } ^ { N } \left ( \lambda _ { 1 } | \langle r _ { 1 }, s _ { i } \rangle | + \lambda _ { 2 } | \langle r _ { 2 }, s _ { i } \rangle | \right ).$$

Imposing a first-order optimality condition with respect to GLYPH&lt;21&gt; 1 ; GLYPH&lt;21&gt; 2 leads to the following closed-form estimation formulas:

$$\lambda _ { 1 } = \left ( \frac { 1 } { N } \sum _ { i = 1 } ^ { N } | \langle r _ { 1 }, s _ { i } \rangle | \right ) ^ { - 1 }, \quad \lambda _ { 2 } = \left ( \frac { 1 } { N } \sum _ { i = 1 } ^ { N } | \langle r _ { 2 }, s _ { i } \rangle | \right ) ^ { - 1 }. \quad \quad ( 1 2 )$$

Substituting the expressions in (12) into the objective function (11), we thus obtain the simplified minimization problem in the only variable GLYPH&lt;18&gt; :

$$\theta ^ { * } \in \arg \min _ { \theta } \left \{ \ln \left ( \sum _ { i = 1 } ^ { N } | \langle r _ { 1 }, s _ { i } \rangle | \right ) + \ln \left ( \sum _ { i = 1 } ^ { N } | \langle r _ { 2 }, s _ { i } \rangle | \right ) \right \}. \quad \quad ( 1 3 )$$

## 4 ADMM

In order to solve numerically the proposed image restoration model (6), we use an ADMM-basedalgorithm - see [1]. We first introduce two auxiliary variables w 2 R n , t 2 R 2 n and rewrite the model in the equivalent linearly constrained form:

$$\Phi _ { 1 } \quad \Phi _ { 2 } \quad \Phi _ { 3 } \quad \Phi _ { 4 } \quad \Phi _ { 5 } \quad \Phi _ { 6 } \quad \Phi _ { 7 } \quad \Phi _ { 8 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \\ \Phi _ { 1 } \quad \Phi _ { 1 } \quad \Phi _ { 2 } \quad \Phi _ { 2 } \quad \Phi _ { 3 } \quad \Phi _ { 4 } \quad \Phi _ { 5 } \quad \Phi _ { 6 } \quad \Phi _ { 6 } \quad \Phi _ { 7 } \quad \Phi _ { 7 } \quad \Phi _ { 8 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 } \quad \Phi _ { 9 }$$

$$\text{subject to} \, \colon \, \ w = K u - g, \ t = D u, \quad \quad ( 1 5 )$$

where the space-variant matrices GLYPH&lt;3&gt; i , R GLYPH&lt;18&gt; i can be estimated via the ML procedure described in Sect. 3 based only on the observed image g (i.e. as a preliminary pre-processing step) or also updated along the ADMM iterations. We define the augmented Lagrangian functional:

$$\mathcal { L } ( u, w, t ; \rho _ { w }, \rho _ { t } ) & \coloneqq \sum _ { i = 1 } ^ { n } \| \Lambda _ { i } R _ { \theta _ { t } } t _ { i } \| _ { 1 } + \frac { \mu } { 2 } \| w \| _ { 2 } ^ { 2 } - \rho _ { t } ^ { T } ( t - D u ) + \frac { \beta _ { t } } { 2 } \| t - D u \| _ { 2 } ^ { 2 } \\ & - \rho _ { w } ^ { T } ( w - ( K u - g ) ) + \frac { \beta _ { w } } { 2 } \| w - ( K u - g ) \| _ { 2 } ^ { 2 }, \\ & \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \,$$

where GLYPH&lt;12&gt; w ; GLYPH&lt;12&gt; t &gt; 0arescalar penalty parameters and GLYPH&lt;26&gt; w 2 R n , GLYPH&lt;26&gt; t 2 R 2 n are the vectors of Lagrange multipliers associated with the given linear constraints. The solution f u GLYPH&lt;3&gt; ; w GLYPH&lt;3&gt; ; t GLYPH&lt;3&gt; g of problem (14) is a saddle point for L in (16), see, e.g., [1]. Hence, we can alternate a minimization step with respect to t ; u ; w with a maximization step with respect to GLYPH&lt;26&gt; t ; GLYPH&lt;26&gt; w . Mathematically,

$$u ^ { ( k + 1 ) } \, \leftarrow \, \arg \min _ { u \in \mathbb { R } ^ { n } } \, \mathcal { L } ( u, w ^ { ( k ) }, t ^ { ( k ) } ; \rho _ { w } ^ { ( k ) }, \rho _ { t } ^ { ( k ) } ),$$

$$w ^ { ( k + 1 ) } \, \leftarrow \, \arg \min _ { r \in \mathbb { R } ^ { n } } \, \mathcal { L } ( u ^ { ( k + 1 ) }, w, t ^ { ( k ) } ; \rho _ { w } ^ { ( k ) }, \rho _ { t } ^ { ( k ) } ),$$

$$t ^ { ( k + 1 ) } \, \leftarrow \, \arg \min _ { t \in \mathbb { R } ^ { 2 n } } \, \mathcal { L } ( u ^ { ( k + 1 ) }, w ^ { ( k + 1 ) }, t ; \rho _ { w } ^ { ( k ) }, \rho _ { t } ^ { ( k ) } ), \, \quad \quad ( 1 9 )$$

$$\rho _ { w } ^ { ( k + 1 ) } \, \leftarrow \, \rho _ { w } ^ { ( k ) } - \beta _ { r } \left ( \, w ^ { ( k + 1 ) } - ( K u ^ { ( k + 1 ) } - g ) \, \right ),$$

$$\rho _ { t } ^ { ( k + 1 ) } \, \leftarrow \, \rho _ { t } ^ { ( k ) } - \beta _ { t } \left ( t ^ { ( k + 1 ) } - D u ^ { ( k + 1 ) } \right ).$$

The solution of the primal sub-problem (17) can be efficiently computed by means of standard linear Fast Fourier Transform (FFT) solvers. The sub-problem (18) can be solved in closed-form by following [3, Section 3]. Finally, the sub-problem (19) can be solved by computing efficiently the proximal operator of the anisotropic 1 norm, for which the proof of [2, Proposition 6.3] can be easily adapted. The regularization parameter GLYPH&lt;22&gt; is updated along the iterations so as to fulfill the global discrepancy principle as described in [4]. We refer the reader also to [5, 7] for more details on the numerical solution of the algorithm.

## 5 Experimental results

In this section, we evaluate the performance of the proposed BLTV-L2 restoration model compared with the baseline TV-L2 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

$$k \geq 1 5 0 0 \text{ or } \delta ^ { ( k ) } \coloneqq \frac { \| u ^ { ( k ) } - u ^ { ( k - 1 ) } \| } { \| u ^ { ( k - 1 ) } \| } \geq 1 0 ^ { - 6 }.$$

The quality of the restored images u GLYPH&lt;3&gt; is measured by means of the Improved Signalto-Noise Ratio ISNR ' u GLYPH&lt;3&gt; ; g ; u ' = 10 log 10 k g GLYPH&lt;0&gt; u k 2 2 k u GLYPH&lt;3&gt; GLYPH&lt;0&gt; u k 2 2 , with u denoting the original uncorrupted image, and of the Structural-Similarity-Index (SSIM) [11].

We consider the test image brain in Fig.2 (a) (570 GLYPH&lt;2&gt; 430) and the test image abdomen in Fig.4 (a) (350 GLYPH&lt;2&gt; 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 GLYPH&lt;27&gt; 2 f 10 20 ; g -see, e.g., Fig.2 (b) and Fig.4 (b). The parameter maps are computed at the beginning starting from the observed image g 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 GLYPH&lt;27&gt; = 20 and of abdomen for GLYPH&lt;27&gt; = 10 via TV-L2 and BLTV-L2 models are shown in Figs.2(c)-(d) and Figs.4(c)-(d), respectively.

From a visual inspection, the restoration via BLTV-L2 seems to be more neat and less cartooned than the TV reconstructions. As reported in Tab.1, the ISNR and

Fig. 2: First row . Original test image brain (a), observed image corrupted by Gaussian blur and AWGN with GLYPH&lt;27&gt; = 20 (b), TV reconstruction (c), BLTV reconstruction (d). Second row . Close-up(s) of the first row.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00801v1-figure-59bf13b7-a77e-477a-aa0c-fb0374d4cc85.jpg)

Table 1: Maximum ISNR/SSIM values achieved by TV-L2 and BLTV-L2 on brain and abdomen test images corrupted by AWGN of two different levels.

|      | brain          | brain          | brain          | brain          | abdomen        | abdomen        | abdomen        | abdomen        |
|------|----------------|----------------|----------------|----------------|----------------|----------------|----------------|----------------|
|      | GLYPH<27> = 10 | GLYPH<27> = 10 | GLYPH<27> = 20 | GLYPH<27> = 20 | GLYPH<27> = 10 | GLYPH<27> = 10 | GLYPH<27> = 20 | GLYPH<27> = 20 |
|      | TV-L 2         | BLTV-L         | TV-L 2         | BLTV-L 2       | TV-L 2         | BLTV-L         | TV-L 2         | BLTV-L 2       |
| 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           |

SSIM values for the two test images and for different noise levels obtained by the BLTV-L2 model outperform the ones reached by the TV-L2 model.

The final parameter maps computed by BLTV-L2 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 (12) is explicit, thus very cheap, whereas the computation of GLYPH&lt;18&gt; GLYPH&lt;3&gt; in (13) 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00801v1-figure-61abbb6f-a5e0-4a8d-95e4-c625a610c58e.jpg)

GLYPH&lt;21&gt;

GLYPH&lt;21&gt;

GLYPH&lt;18&gt;

Fig. 3: Final parameter maps for brain test image.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00801v1-figure-875524ea-d101-42a8-868b-979c73037ccd.jpg)

Fig. 4: First row . Original test image abdomen (a), observed image corrupted by Gaussian blur and AWGN with GLYPH&lt;27&gt; = 10 (b), TV reconstruction (c), BLTV reconstruction (d). Second row . Close-up(s) of the first row.

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 distribu-

Fig. 5: Final parameter maps for abdomen test image.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00801v1-figure-ecf7831c-0672-46a7-aa37-fda0dded0b44.jpg)

tions. 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 [8] - is a matter being studied.

## References

- 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 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 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.: 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 p -' 2 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 &amp;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).
Luca Calatroni Alessandro Lanza Monica Pragliola Fiorella Sgallari
Space-adaptive anisotropic bivariate Laplacian regularization for image restoration
In this paper we present a new regularization term for variational image
restoration which can be regarded as a space-variant anisotropic extension of
the classical isotropic Total Variation (TV) regularizer. The proposed
regularizer comes from the statistical assumption that the gradients of the
target image distribute locally according to space-variant bivariate Laplacian
distributions. The highly flexible variational structure of the corresponding
regularizer encodes several free parameters which hold the potential for
faithfully modelling the local geometry in the image and describing local
orientation preferences. For an automatic estimation of such parameters, we
design a robust maximum likelihood approach and report results on its
reliability on synthetic data and natural images. A minimization algorithm
based on the Alternating Direction Method of Multipliers (ADMM) is presented
for the efficient numerical solution of the proposed variational model. Some
experimental results are reported which demonstrate the high-quality of
restorations achievable by the proposed model, in particular with respect to
classical Total Variation regularization.