Space-adaptive anisotropic bivariate Laplacian regularization for image restoration

TL;DR

This paper introduces a space-adaptive anisotropic bivariate Laplacian regularizer for image restoration, improving upon classical Total Variation by modeling local image geometry and orientation with a flexible, statistically grounded variational approach.

Contribution

It proposes a novel space-variant regularizer based on bivariate Laplacian distributions, with an automatic parameter estimation method and an efficient ADMM-based minimization algorithm.

Findings

Achieves higher quality image restorations compared to classical TV.

Demonstrates reliable parameter estimation on synthetic and natural images.

Shows the effectiveness of the proposed model through experimental results.

Abstract

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.

Tables

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

	brain				abdomen
	$\sigma =10$		$\sigma =20$		$\sigma =10$		$\sigma =20$
	TV-L2	BLTV-L2	TV-L2	BLTV-L2	TV-L2	BLTV-L2	TV-L2	BLTV-L2
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

Equations

$$g=Ku+e,$$

$$u^{*}\in\arg\min_{u\in\mathbb{R}^{n}}\;\big{\{}\mathcal{J}(u;K,g):=\mathcal{R}(u)+\mu\mathcal{F}(u;K,g)\big{\}}.$$

$$\mathcal{F}(u;K,g)={\mathrm{L}}_2(u;K,g)=\frac{1}{2}\Vert Ku-g{\Vert }_2^2.$$

$$\mathcal{R}(u)=\mathrm{TV}(u)=\sum _{i=1}^n\Vert (\nabla u{)}_i{\Vert }_2,$$

$$u^{*}\in\arg\min_{u\in\mathbb{R}^{n}}\bigg{\{}\sum_{i=1}^{n}\lVert(\nabla u)_{i}\rVert_{2}+\frac{\mu}{2}\lVert Ku-g\rVert_{2}^{2}\bigg{\}}.$$

$$\mathrm{BLTV}(u;{\lambda }_1,{\lambda }_2,\theta )$$

$$=$$

$${\Lambda }_i=\left(\begin{array}{cc}{\lambda }_i^{(1)}& 0\\ 0& {\lambda }_i^{(2)}\end{array}\right),R_{{\theta }_i}=\left(\begin{array}{cc}\cos {\theta }_i& \sin {\theta }_i\\ -\sin {\theta }_i& \cos {\theta }_i\end{array}\right)=\left(\begin{array}{c}{r}_i^{(1)}\\ r_i^{(2)}\end{array}\right).$$

$$u^{*}\in\arg\min_{u\in\mathbb{R}^{n}}\bigg{\{}\mathrm{BLTV}(u;\lambda_{1},\lambda_{2},\theta)+\frac{\mu}{2}\lVert Ku-g\rVert_{2}^{2}\bigg{\}}.$$

$$u^{\ast }\in \arg \underset{u\in {\mathbb{R}}^n}{\max }\text{Pr}(u|g;K)=\arg \underset{u\in {\mathbb{R}}^n}{\min }\left\{-\ln \text{Pr}(g|u;K)-\ln \text{Pr}(u)\right\},$$

$$\mathrm{Pr}(g|u;K)=\prod_{i=1}^{n}\,\frac{1}{\sqrt{2\pi}\sigma}\,\text{exp}\bigg{(}-\frac{(Ku-g)_{i}^{2}}{2\sigma^{2}}\,\bigg{)}=W\,\text{exp}\bigg{(}-\frac{\lVert Ku-g\rVert_{2}^{2}}{2\sigma^{2}}\,\bigg{)},$$

$$\text{Pr}(u)=\prod_{i=1}^{n}\text{Pr}_{i}(u)=\prod_{i=1}^{n}z_{i}\,\text{exp}\,(\,-\alpha\,V_{c_{i}}(u)\,)=Z\,\text{exp}\,\bigg{(}\,-\alpha\,\sum_{i=1}^{n}V_{c_{i}}(u)\,\bigg{)},$$

$${\text{Pr}}_i(u)=z_i\text{exp}\left(-\alpha \Vert (\nabla u{)}_i{\Vert }_2\right).$$

$$\mathrm{Pr}(u)$$

$$\mathrm{Pr}(\mathcal{S}|\lambda_{1},\lambda_{2},\theta)=\bigg{(}\frac{\lambda_{1}\,\lambda_{2}}{4}\bigg{)}^{N}\exp\bigg{(}-\sum_{i=1}^{N}\bigg{(}\lambda_{1}|\langle r_{1},s_{i}\rangle|+\lambda_{2}|\langle r_{2},s_{i}\rangle|\bigg{)}\bigg{)}\,,$$

$$({\lambda }_1^{\ast },{\lambda }_2^{\ast },{\theta }^{\ast })$$

$$\lambda_{1}=\bigg{(}\frac{1}{N}\sum_{i=1}^{N}|\langle r_{1},s_{i}\rangle|\bigg{)}^{-1}\,,\quad\lambda_{2}=\bigg{(}\frac{1}{N}\sum_{i=1}^{N}|\langle r_{2},s_{i}\rangle|\bigg{)}^{-1}.$$

$$\theta^{*}\in\arg\min_{\theta}\bigg{\{}\mathrm{ln}\bigg{(}\sum_{i=1}^{N}|\langle r_{1},s_{i}\rangle|\bigg{)}+\mathrm{ln}\bigg{(}\sum_{i=1}^{N}|\langle r_{2},s_{i}\rangle|\bigg{)}\bigg{\}}\,.$$

$$\{u^{\ast },w^{\ast },t^{\ast }\}$$

$$\mathcal{L}(u,w,t;{\rho }_w,{\rho }_t)$$

$$u^{(k+1)}$$

$$w^{(k+1)}$$

$$t^{(k+1)}$$

$${\rho }_w^{(k+1)}$$

$${\rho }_t^{(k+1)}$$

$$k\ge 1500\mathrm{or}{\delta }^{(k)}:=\frac{\Vert u^{(k)}-u^{(k-1)}\Vert }{\Vert u^{(k-1)}\Vert }\ge 10^{-6}.$$

Full text

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

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

[TABLE]

where $g,u\in\mathbb{R}^{n}$ are vectorized images, $K\in\mathbb{R}^{n\times n}$ is the linear blur operator and $e\in\mathbb{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^{*}$, estimate of the original $u$, which solves a well-posed problem. In the variational approach, $u^{*}$ is the minimizer of a cost functional $\mathcal{J}(u;K,g):\mathbb{R}^{n}\to\mathbb{R}$. In formula,

[TABLE]

The functionals $\mathcal{R}$ and $\mathcal{F}$ are commonly referred to as the regularization and the data fidelity term, respectively. While $\mathcal{R}$ encodes prior information on the desired image $u$, $\mathcal{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., stuart . The regularization parameter $\mu>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\sim\mathcal{N}(0,\sigma^{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 $\mathcal{F}(u;K,g)$ is the so-called L2 fidelity term, reading as,

[TABLE]

A popular choice for the regularization term $\mathcal{R}(u)$ is given by the TV semi-norm ROF ,

[TABLE]

where (\nabla u)_{i}:=\big{(}(D_{h}u)_{i},(D_{v}u)_{i}\big{)}^{T}\in\mathbb{R}^{2} denotes the discrete gradient of image $u$ at pixel $i$, with $D_{h},D_{v}\in\mathbb{R}^{n\times n}$ linear operators representing finite difference discretizations of the first-order horizontal and vertical partial derivatives, respectively.

Coupling the 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,

[TABLE]

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 siam ; ncmip ; cmbbe ; vip , the authors have proposed space-variant regularization term based on statistical assumptions on the distribution of the $\ell_{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 (4) 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

[TABLE]

where $\Lambda_{i}$ is a $2\times 2$ positive definite diagonal matrix and $R_{\theta_{i}}$ is the rotation matrix corresponding to the angle $-\theta_{i}$. Mathematically,

[TABLE]

We denote by $\lambda_{1},\lambda_{2},\theta\in\mathbb{R}^{n}$ the maps of the parameters defining the local distributions. Hence, the proposed BLTV-L2 variational restoration model reads as

[TABLE]

The $3n$ 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 $\lambda^{(1)}$, $\lambda^{(2)}$, $\theta$ 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-L2 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$(u|g;K)$ of the unknown target image $u$ given the observation $g$, namely:

[TABLE]

where, after applying the Bayes’ rule, we drop the evidence term Pr$(g)$ and extracted $-\ln$ of the objective function. The two terms Pr$(u)$ and Pr$(g|u;K)$ in (10) are referred to as the prior and the likelihood. The likelihood term associated with AWGN corruption takes the form

[TABLE]

where $\sigma>0$ denotes the AWGN standard deviation and $W>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:

[TABLE]

where $\alpha>0$ is the MRF parameter, $\{c_{i}\}_{i=1}^{n}$ 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=\prod_{i=1}^{n}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}}=\lVert(\nabla u)_{i}\rVert_{2}$ for any $i\in\{1,\ldots,n\}$, 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

[TABLE]

Here, we propose a space-variant anisotropic generalization of the local prior $\text{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

[TABLE]

where $Z$ is the normalization constant and $\Lambda_{i},R_{{\theta}_{i}}\in\mathbb{R}^{2\times 2}$ are defined in (8).

Plugging the prior (14) and the likelihood (11) into (10) and neglecting the constant terms, the proposed BLTV-L2 variational model (9) 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 $\lambda^{(1)},\lambda^{(2)},\theta\in\mathbb{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 $\mathcal{S}:=\left\{s_{1},\ldots,s_{N}\right\}$ drawn from a BLD with parameters $(\lambda_{1},\lambda_{2},\theta)$. 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

[TABLE]

where, clearly, $r_{1}$ and $r_{2}$ depend on $\theta$ - see (8). The goal here is to find

[TABLE]

Imposing a first-order optimality condition with respect to $\lambda_{1},\lambda_{2}$ leads to the following closed-form estimation formulas:

[TABLE]

Substituting the expressions in (17) into the objective function (16), we thus obtain the simplified minimization problem in the only variable $\theta$:

[TABLE]

4 ADMM

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 $w\in\mathbb{R}^{n}$, $t\in\mathbb{R}^{2n}$ and rewrite the model in the equivalent linearly constrained form:

[TABLE]

where the space-variant matrices $\Lambda_{i}$, $R_{\theta_{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:

[TABLE]

where $\beta_{w},\beta_{t}>0$ are scalar penalty parameters and $\rho_{w}\in\mathbb{R}^{n}$, $\rho_{t}\in\mathbb{R}^{2n}$ are the vectors of Lagrange multipliers associated with the given linear constraints. The solution $\{u^{*},w^{*},t^{*}\}$ of problem (20) is a saddle point for $\mathcal{L}$ in (21), see, e.g., boyd . Hence, we can alternate a minimization step with respect to $t,u,w$ with a maximization step with respect to $\rho_{t},\rho_{w}$. Mathematically,

[TABLE]

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 $1$ norm, for which the proof of (siam, , Proposition 6.3) can be easily adapted. The regularization parameter $\mu$ 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-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

[TABLE]

The quality of the restored images $u^{*}$ is measured by means of the Improved Signal-to-Noise Ratio $\mathrm{ISNR}(u^{*},g,u)=10\log_{10}\frac{\|g-u\|_{2}^{2}}{\|u^{*}-u\|_{2}^{2}}$, with $u$ denoting the original uncorrupted image, and of the Structural-Similarity-Index (SSIM) ssim .

We consider the test image brain in Fig.2 (a) (570 $\times$ 430) and the test image abdomen in Fig.4 (a) (350 $\times$ 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 $\sigma\in\{10,20\}$ - 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 $\sigma=20$ and of abdomen for $\sigma=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 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 (17) is explicit, thus very cheap, whereas the computation of $\theta^{*}$ 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.