Space-variant Generalized Gaussian Regularization for Image Restoration

TL;DR

This paper introduces a space-variant regularizer based on a local half-Generalized Gaussian distribution for image restoration, adaptable to various noise types and image gradient distributions, with automatic parameter estimation.

Contribution

It proposes a novel, flexible regularizer with per-pixel parameters, coupled with efficient algorithms for high-quality image restoration under different noise conditions.

Findings

Effective in restoring images with diverse gradient distributions.

Performs well with Gaussian, Laplace, and salt-and-pepper noise.

Achieves high-quality results in numerical experiments.

Abstract

We propose a new space-variant regularization term for variational image restoration based on the assumption that the gradient magnitudes of the target image distribute locally according to a half-Generalized Gaussian distribution. This leads to a highly flexible regularizer characterized by two per-pixel free parameters, which are automatically estimated from the observed image. The proposed regularizer is coupled with either the $L_2$ or the $L_1$ fidelity terms, in order to effectively deal with additive white Gaussian noise or impulsive noises such as, e.g, additive white Laplace and salt and pepper noise. The restored image is efficiently computed by means of an iterative numerical algorithm based on the alternating direction method of multipliers. Numerical examples indicate that the proposed regularizer holds the potential for achieving high quality restorations for a wide range of target images characterized by different gradient distributions and for the different types of noise considered.

Tables

Table 1. Table 1 : Example 1: achieved ISNR values.

			geometric				skyscraper
BSNR	TV-L2	TVp-L2	TV${}_{}{}^{\mathrm{sv}}{}_p{}^{}$-L2	TV${}_{}{}^{\mathrm{sv}}{}_{p,\alpha }{}^{}$-L2	TV-L2	TVp-L2	TV${}_{}{}^{\mathrm{sv}}{}_p{}^{}$-L2	TV${}_{}{}^{\mathrm{sv}}{}_{p,\alpha }{}^{}$-L2
20	7.77	7.92	8.36	8.60	2.76	3.00	3.07	3.31
30	9.01	9.87	10.30	10.57	5.12	5.52	5.94	6.40

Table 2. Table 2 : Example 2 (AWLN): achieved ISNR values.

	TV-L1	TVp-L1	TV${}_{}{}^{\mathrm{sv}}{}_p{}^{}$-L1	TV${}_{}{}^{\mathrm{sv}}{}_{p,\alpha }{}^{}$-L1
lungs	6.20	6.80	7.30	7.85
ecography	5.93	6.40	7.88	8.32
aneurism	9.10	9.44	10.13	10.70

Equations

$$g=\mathcal{N}\left(Ku\right),$$

$$u^{\ast }\leftarrow \arg \underset{u\in {\mathbb{R}}^n}{\min }\mathcal{J}(u),\mathcal{J}(u):=\mathcal{R}(u)+\mu \mathcal{F}(u;g),$$

$$\begin{array}[]{ccc}\mathrm{AWGN}$\;$\mathrm{and}$\;$\mathrm{AWLN:}&&\mathrm{SPN:}\vspace{0.1cm}\\ g_{i}\>\;{=}\;\>(Ku)_{i}\;{+}\;n_{i}\;\;\>\forall\,i\in\Omega\,,&&g_{i}\>\;{=}\;\>\left\{\begin{array}[]{ll}(Ku)_{i}&\mathrm{for}\;\;\>i\in\Omega_{0}\subseteq\Omega\\ n_{i}\in\{0,1\}&\mathrm{for}\;\;\>i\in{\Omega_{1}}:=\Omega\setminus\Omega_{0}\end{array}\right.\,.\end{array}$$

$$\mathcal{F}(u;g)={\mathrm{L}}_q(u;g):=\frac{1}{q}\Vert Ku-g{\Vert }_q^q,q\in \{1,2\}.$$

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

$$u^{\ast }\leftarrow \arg \underset{u\in {\mathbb{R}}^n}{\min }\left\{\mathrm{TV}(u)+\mu {\mathrm{L}}_q(u;g)\right\},q\in \{1,2\}.$$

$$u^{\ast }\leftarrow \arg \underset{u\in {\mathbb{R}}^n}{\min }\left\{{\mathrm{TV}}_{p,\alpha }^{\mathrm{sv}}(u)+\mu {\mathrm{L}}_q(u;g)\right\},q\in \{1,2\},$$

$${\mathrm{TV}}_{p,\alpha }^{\mathrm{sv}}(u):=\sum _{i=1}^n{\alpha }_i\Vert (\nabla u{)}_i{\Vert }_2^{p_i},{\alpha }_i\in ]0,+\infty [,p_i\in ]0,2]\forall i\in \Omega .$$

$$u^{\ast }\leftarrow \arg \underset{u\in {\mathbb{R}}^n}{\max }\text{Pr}(u|g;K).$$

$$u^{\ast }\leftarrow \arg \underset{u\in {\mathbb{R}}^n}{\min }\left\{-\log \text{Pr}(g|u;K)-\log \text{Pr}(u)\right\}.$$

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

$$\text{Pr}(u)=\frac{1}{Z}\,\text{exp}\,\bigg{(}\,-\alpha\,\sum_{i=1}^{n}\lVert(\nabla u)_{i}\rVert_{2}\,\bigg{)}=\frac{1}{Z}\,\text{exp}\bigg{(}-\alpha\,\text{TV}(u)\,\bigg{)},$$

$$\text{Pr}(x;\alpha)\>\;{=}\;\>\left\{\begin{array}[]{ll}\alpha\,\exp\,(\,-\alpha\,x\,)&\mathrm{for}\;\;\>x\geq 0\\ &\\ 0&\mathrm{for}\;\;\>x<0\end{array}\right.\,.$$

$$\text{Pr}(x;p,\alpha)\>\;{=}\;\>\left\{\begin{array}[]{ll}\frac{\alpha p}{\Gamma(1/p)}\exp(-(\alpha x)^{p})&\mathrm{for}\;\;\>x\geq 0\\ &\\ 0&\mathrm{for}\;\;\>x<0\end{array}\right.\,$$

$$\text{Pr}(u)=\frac{1}{Z}\,\text{exp}\,\bigg{(}\,-\alpha\,\sum_{i=1}^{n}\lVert(\nabla u)_{i}\rVert_{2}^{p}\,\bigg{)}=\frac{1}{Z}\,\text{exp}\bigg{(}-\alpha\,\text{TV}_{p}(u)\,\bigg{)}.$$

$$\text{Pr}(u)=\frac{1}{Z}\,\text{exp}\,\bigg{(}\,-\,\sum_{i=1}^{n}\alpha_{i}\lVert(\nabla u)_{i}\rVert_{2}^{p_{i}}\,\bigg{)}=\frac{1}{Z}\,\text{exp}\bigg{(}-\,\text{TV}_{p,\alpha}^{\mathrm{sv}}(u)\,\bigg{)}.$$

$$\text{Pr}(g|u;K)=\prod _{i=1}^n\text{Pr}(g_i|u;K).$$

$$\text{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{)}=\frac{1}{W}\,\text{exp}\bigg{(}-\frac{\lVert Ku-g\rVert_{2}^{2}}{2\sigma^{2}}\,\bigg{)},$$

$$u^{\ast }\leftarrow \arg \underset{u\in {\mathbb{R}}^n}{\min }\left\{\sum _{i=1}^n{\alpha }_i\Vert (\nabla u{)}_i{\Vert }_2^{p_i}+\frac{\mu }{2}\Vert Ku-g{\Vert }_2^2\right\},$$

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

$$u^{\ast }\leftarrow \arg \underset{u\in {\mathbb{R}}^n}{\min }\left\{\sum _{i=1}^n{\alpha }_i\Vert (\nabla u{)}_i{\Vert }_2^{p_i}+\mu \Vert Ku-g{\Vert }_1\right\},$$

$$\text{Pr}\big{(}\,\lVert(\nabla u)_{i}\rVert\,;p_{i},\alpha_{i}\,\big{)}\>\;{=}\;\>\frac{\alpha_{i}p_{i}}{\Gamma(1/p_{i})}\,\exp\big{(}-\big{(}\,\alpha_{i}\,\lVert(\nabla u)_{i}\rVert\,\big{)}^{p_{i}}\,\big{)}.$$

$$m_i:={\Vert (\nabla u{)}_i\Vert }_2,i\in \Omega .$$

$$p_{i}\,\;{=}\;\,h^{-1}(\rho_{i}),\quad\;\rho_{i}\,\;{=}\;\;\mathrm{card}\big{(}N_{\,i}^{\,s}\big{)}\,\bigg{(}\sum_{j\in N_{\,i}^{\,s}}\!m_{j}^{2}\bigg{)}\,/\,\bigg{(}\sum_{j\in N_{\,i}^{\,s}}\!|m_{j}|\bigg{)}^{\!\!2},\quad i\in\Omega\,,$$

$$h(z)\,\;{=}\;\,\big{(}\Gamma(1/z)\,\,\Gamma(3/z)\big{)}\,/\,\big{(}\Gamma^{2}(2/z)\big{)}\,,$$

$$\mathcal{L}(\alpha ,p_i;x_1,...,x_n)$$

$${\alpha }_i$$

$$\alpha_{i}=\bigg{(}\frac{p_{i}}{n}\sum_{i=1}^{n}x_{i}^{p_{i}}\bigg{)}^{-\frac{1}{p_{i}}}.$$

$$\{u^{\ast },r^{\ast },t^{\ast }\}\leftarrow \arg \underset{u,r,t}{\min }$$

$$\mathrm{subject}\mathrm{to}:$$

$$\mathcal{L}(u,r,t;{\lambda }_r,{\lambda }_t)$$

Full text

1 Introduction

Image restoration refers to the recovery of a clean sharp image from a noisy, and potentially blurred, observation. In this paper, we consider the problem of restoring images corrupted by known blur and different types of noise.

We consider gray level images of size $\,d_{1}\!{\times}\,d_{2}$, such that $n\,{:=}\,d_{1}d_{2}$ is the total number of pixels in the images. The general model of the image degradation process under blur and noise corruptions can be written as

[TABLE]

where $u,g\in\mathbb{R}^{n}$ represent vectorized forms of the unknown clean image and of the observed corrupted image, respectively, $K\in\mathbb{R}^{n\times n}$ is a known linear blurring operator and $\mathcal{N}(\,\cdot\,)$ denotes the noise corruption operator, which in most cases is of random nature. Given $K$ and $g$, the goal of image restoration is to solve the ill-conditioned - or even singular, depending on $K$ - inverse problem of recovering an as accurate as possible estimate $u^{*}$ of the unknown clean image $u$.

The class of variational methods for image restoration relies on determining restored images $u^{*}\in\mathbb{R}^{n}$ as the minimizers of suitable cost functionals $\mathcal{J}(u):\mathbb{R}^{n}\to\mathbb{R}$ such that, typically, restoration is casted as an optimization problem of the form

[TABLE]

where the functionals $\mathcal{R}(u)$ and $\mathcal{F}(u;g)$, commonly referred as the regularization and the fidelity term, encode prior information on the clean image $u$ and on the observation model (1), respectively, with the so-called regularization parameter $\mu>0$ controlling the trade-off between the two terms.

The functional form of the fidelity term is strictly connected with the characteristics of the noise corruption. In this paper, we are interested in three important types of noise, namely the additive (zero-mean) white Gaussian noise (AWGN) which typically appears, e.g., in Magnetic Resonance Tomography, the additive (zero-mean) white Laplace noise (AWLN) and the impulsive salt and pepper noise (SPN) usually due to transmission errors or malfunctioning pixel elements in camera sensors. Denoting by $\Omega:=\{1,\ldots,n\}$ the set of all pixel positions in the vectorized images, for these three kinds of noise the general degradation model in (1) reads as

[TABLE]

For what concerns AWGN and AWLN, the additive corruptions $n_{i}\in\mathbb{R}$, $i\in\Omega$, represent independent realizations from the same univariate Gaussian and Laplace distribution with zero mean and standard deviation $\sigma$, respectively. In the case of SPN, only a subset $\Omega_{1}$ of the pixels is corrupted by noise, whereas the complementary subset $\Omega_{0}$ is noise-free. In particular, the corrupted pixels can take only the two possible extreme values $\{0,1\}$ (we assume that images have range $[0,1]$), with the same probability. The subset $\Omega_{1}$ is known in some applications [11] or it could be estimated [5]. As the zero-mean AWGN and AWLN are fully characterized from a probabilistic point of view by the unique scalar parameter $\,\sigma$, SPN is characterized by the parameter $\gamma\in[0,1]$ which represents the probability for a pixel to be noise-corrupted.

It is well known that AWGN and the impulsive AWLN and SPN are suitably dealt with by the so-called L2 and L1 fidelity terms, which are related to the $\ell_{2}$ and $\ell_{1}$ norm of the residue image, respectively [15]; in formulas:

[TABLE]

For what regards the regularization term in (2), a very popular choice is represented by the Total Variation (TV) semi-norm [13], that is

[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. Popularity of TV regularizer for image restoration is mainly due to two facts, namely it is convex and allows for restored images with sharp, neat edges. By substituting the TV regularizer (5) and the L2 or L1 fidelity terms (4) for $\mathcal{R}$ and $\mathcal{F}$ in (2), respectively, one obtains the so-called TV-L2 [13] - or ROF - and TV-L1 [12] restoration models, which reads as

[TABLE]

The TV-L2 and TV-L1 models in (6) are non-smooth convex and allow to obtain good quality restorations of images corrupted by AWGN and AWLN/SPN, respectively, such that they are regarded as sort of baseline models.

The contribution of this paper consists of a new space-variant regularization term which, coupled with the $L_{2}$ or $L_{1}$ fidelity terms, gives rise to a generalization of the models in (6) of the form

[TABLE]

where the new space-variant TV${}_{p,\alpha}^{\mathrm{sv}}$ regularizer is defined by

[TABLE]

The proposed regularizer in (8) is highly flexible as it is characterized by two per-pixel free parameters $p_{i}$, $\alpha_{i}$, such that local, space-variant properties of the target clean image $u$ can be potentially addressed. The usefulness of such a great flexibility in the proposed regularizer is however conditioned to the existence of effective procedures for the automatic estimation of the $p_{i}$ and $\alpha_{i}$ parameters. Hence, in this paper we also propose a suitable method for the automatic estimation of such parameters from the observed image partially based on the statistical inference technique described in [8].

As outlined in Section 2, the rationale of our proposal is that the distribution of the gradient magnitudes of the unknown clean image is space-variant and it is well modeled locally by a two-parameters Generalized Gaussian distribution. As highlighted in [9], the TV regularizer in (5) comes from implicitly assuming a space-invariant (that is, frame-based), one-parameter half-Laplacian distribution for the gradient magnitudes. Based on the observation that such a distribution is not sufficiently flexible, the authors in [9] proposed a generalization of the TV regularizer, referred to as TVp, which relies on a space-invariant, two-parameters half-Generalized Gaussian distribution with $p$ denoting the additional parameter, called the shape parameter. Finally, the very recently proposed TV${}_{p}^{\mathrm{sv}}$ regularizer [16] further generalizes the TVp by assuming a space-variant shape parameter in the two-parameters half-Generalized Gaussian distribution. The authors in [16] demonstrated experimentally how using a local, space-variant model holds the potential for achieving higher quality restorations than using a global, space-invariant model.

The numerical solution of the two proposed variational models (7)–(8) are obtained by means of an efficient iterative minimization algorithm based on the Alternating Direction Method of Multipliers (ADMM) strategy [4].

The paper is organized as follows. The derivation of the models (7) and a more detailed motivation of their introduction are proposed in Section 2. In Section 3 we briefly outline the statistical inference procedure used for automatically estimating the $p_{i}$ and $\alpha_{i}$ parameters. The ADDM-based minimization algorithm is described in detail in Section 4. Some meaningful numerical results are reported in Section 5 and, finally, conclusions are drawn in Section 6.

2 Motivation via MAP estimator

The Maximum A Posteriori (MAP) Estimation approach [3] relies on the maximization of the posterior probability density function Pr$(u|g;K)$ associated to the clean unknown image $u$:

[TABLE]

Relying on the Bayes’ formula, and dropping the evidence term Pr$(g)$, this is equivalent to maximize the product of the prior Pr$(u)$ and the likelihood Pr$(g|u;K)$ probability density functions. By taking the negative logarithm of this product, problem (9) can be reformulated as follows:

[TABLE]

At first, we focus on the setting of 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$ is the partition function, that is a function not depending on $u$ which allows for the normalization of the prior.

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}$, the Gibbs prior in (11) reduces to the popular TV prior:

[TABLE]

where $Z$ is the normalization constant not depending on $u$. The adoption of a TV prior can be further interpreted as assuming that the $\ell_{2}$ norm of the gradient at any pixel of the unknown clean image, $\lVert(\nabla u)_{i}\rVert_{2}$, follows a space-invariant half-Laplacian (or exponential) distribution:

[TABLE]

In [9], a deep investigation about the effect of replacing the half-Laplacian distribution with the more flexible half-Generalized Gaussian distribution

[TABLE]

has been carried out. The presence of a second parameter $p$ allows for a better approximation of the $\ell_{2}$ norm gradient distribution and leads to the introduction of the TVp prior:

[TABLE]

In this paper, we propose a prior, consisting in a non-stationary (space-variant) Markov Random Field. The parameters $\alpha,p$ of the half-Generalized Gaussian distribution of the magnitude of the discrete gradients change as the clique $c_{i}$ changes. Therefore, the prior takes the following form:

[TABLE]

The adoption of a space-variant approach is expected to be more flexible for the restoration of images presenting cliques with different properties, i.e. images in which texture, smooth, piecewise constant regions, and edges co-exist.

In order to justify the reason why a space-variant model should be adopted in general, we consider the test image skyscraper illustrated in Fig. 1(a). We selected two regions characterized by smooth and texture structures - see the cyan-bordered and the yellow-bordered boxes, respectively, in Figs. 1(d),1(g).

The histogram of the gradient magnitudes of the whole image is shown in Fig 1(b) and zoomed in Fig. 1(c). The superimposed green dashed lines, which have been reproduced in each sub-figure for comparison, represent the half-Generalized Gaussian distributions that best fit the histograms and whose parameters have been computed considering all the pixels of the test image.

The histogram of the gradient magnitudes in the two bordered regions are shown in Figs. 1(e), 1(h) and zoomed in Figs. 1(f), 1(i). The red solid lines represent the half-Generalized Gaussian distributions that best fit the histograms and whose parameters have been computed considering only the pixels in the boxes.

It is worth noticing how the histograms of the gradient magnitudes in the two selected regions are very different from each other and also differ from the one of the whole test image. As a result, the red lines fit the histogram shapes in Figs. 1(e), 1(h) better than the green ones - see also the zooms in Figs. 1(f),1(i). This is the benefit of the space-variant strategy, which is able to model space-variant image features.

Going back to the MAP inference formula (10), in particular to the likelihood term Pr$(g|u;K)$, assuming the noise to be additive independent identically distributed, we have:

[TABLE]

The likelihood term (16) clearly takes different form according to the distribution of the noise. In the following we specify the likelihood associated to the noises considered in this paper.

Additive White Gaussian Noise. If the noise is known to be AWG with standard deviation $\sigma$, the likelihood term in (16) is as follows:

[TABLE]

where $W$ is the normalization constant not depending on $u$.

Therefore, after replacing our prior (15) and the Gaussian likelihood (17) in the MAP inference formula (10), and dropping the constant terms, we obtain our TV${}_{p,\alpha}^{\mathrm{sv}}$-L2 model in (7) with $q=2$, that is in extended form:

[TABLE]

where we set $\mu=1/\sigma^{2}$.

Additive White Laplace Noise. If the noise is known to be AWL with scale parameter $\beta$, the likelihood term in (16) takes the following form:

[TABLE]

Therefore, after replacing our prior (15) and the Laplace likelihood (19) in the MAP inference formula (10), and dropping the constant terms, we obtain our TV${}_{p,\alpha}^{\mathrm{sv}}$-L1 model in (7) with $q=1$, that is in extended form:

[TABLE]

where we set $\mu=1/\beta$.

Salt and Pepper noise. The SPN can be classified as a sparse noise, since it corrupts only a subset of pixels according to (3) . In this case, in order to strongly promote the sparsity of the noise, a popular choice is to adopt the $\ell_{0}$ pseudo-norm of the residual $Ku-g$ as the fidelity term. Nevertheless, it is very common to substitute the $\ell_{0}$ pseudo-norm with the $\ell_{1}$ norm, which is easier to deal with - since it is convex - and still allows a good sparsification effect. Hence, also in this case, the problem to which we are referring is (20).

3 Estimation of the space-variant parameters

The proposed regularization term (8) is derived by assuming that, for each pixel position $i=1,...,n$, the magnitude - that is, the $\ell_{2}$ norm - of the gradients of the target image in a surrounding neighborhood distributes according to a half-Generalized Gaussian (hGG) distribution, whose probability density function is given in (13). This means that the distribution of the $\ell_{2}$ norm of the gradients in the target image is defined pixel-wise as follows:

[TABLE]

In order to use the proposed regularization term, we thus need to generate the $p$-map and the $\alpha$-map.

The method proposed in [9] for estimating a global, image-based $p$ value requires a very large number of samples in order to provide statistically reliable estimates, therefore it could not be generalized to our proposal since we use small size image neighborhoods for the estimation of local $p$ values. In [16] the authors proposed a new method based on the statistical inference procedure illustrated in [8] which is sufficiently robust to our purposes. For completeness, in the following we briefly outline the method.

Let $u\in\mathbb{R}^{n}$ be the vectorized form of an image for which we want to estimate the associated vector of space-variant parameters $p_{i}$, $i\in\Omega$. First, we compute the vector $m\in\mathbb{R}^{n}$ containing the magnitudes of the gradients of the image $u$; in formulas:

[TABLE]

Then, we estimate each parameter $p_{i}$ by applying the statistical inference technique in [8] to the local data set consisting of the computed gradient magnitudes in a neighborhood of the pixel $i$. In particular, we use symmetric square neighborhoods $N_{\,i}^{\,s}$ of size $s\in\{3,5,\ldots\}$ centered at pixel $i\in\Omega$. Following [8], the values $p_{i}$, $i\in\Omega$, shape parameters of the hGG distributions, are estimated as follows:

[TABLE]

where $\mathrm{card}(A)$ denotes the cardinality of set $A$ and where the function $h:]0,+\infty[\to]0,+\infty[$, referred to as the generalized Gaussian ratio function in [8], is defined by

[TABLE]

with $\Gamma(\,\cdot\,)$ indicating the Gamma function [1]. The function $h$ in (24) is continuous, monotonically decreasing and surjective, hence invertible. Moreover, since $h$ is not data-dependent, its inverse $h^{-1}$, representing the values $p_{i}$, can be pre-computed off-line and stored as a lookup-table, restricted to $(0,2]$, such that at run-time the final step of the estimation in (23) can be carried out very efficiently.

The key novelty of our proposal relies on exploiting all the advantages of using a space-variant hGG distribution model, hence we compute also the map of local scale parameters $\alpha_{i}$, $i=1,\ldots,n$. We propose to estimate such scale parameters by means of a Maximum Likelihood approach. Once $p_{i}$ for a pixel is estimated, the local likelihood function is given by

[TABLE]

such that the value of the local scale parameter is obtained by maximizing (25), that is by solving the following optimization problem:

[TABLE]

By imposing the first order optimality condition for problem (26), we obtain the following closed form estimation formula:

[TABLE]

In Fig. 2 the maps of local $p$ values, obtained with neighborhoods of size $s=3$ (b) and $s=11$ (c) starting from the original test image skyscraper (a) are shown. Both maps are scaled in the same range for visual comparison. As the size $s$ increases, image features of increasing scale are highlighted, but in any case the method associates very low $p$ values with flat regions and higher values with edges. It is worth remarking that in Sect. 5 numerical experiments have been carried out by computing the $p$-map starting from the corrupted images.

In Fig. 3 the maps of local $\alpha$ values, obtained with neighborhoods of size $s=3$ (b) and $s=11$ (c) starting from the original test image skyscraper (a) are shown. Both maps are scaled in the range $[0,1]$ for visual comparison.

4 Numerical solution by ADMM

In this section, we illustrate the ADMM-based iterative algorithm used to numerically solve the proposed model (7)–(8) for both cases $q=2$ and $q=1$. To this purpose, first we resort to the variable splitting technique [2] and introduce two auxiliary variables $r\in\mathbb{R}^{n}$ and $t\in\mathbb{R}^{2n}$, such that model (7)–(8) is rewritten in the following equivalent constrained form:

[TABLE]

where $D:=(D_{h}^{T},D_{v}^{T})^{T}\in\mathbb{R}^{2n\times n}$ and t_{i}\>{:=}\>\big{(}(D_{h}u)_{i}\,,\,(D_{v}u)_{i}\big{)}^{T}\in\mathbb{R}^{2} represents the discrete gradient of image $u$ at pixel $i$.

To solve problem (28)–(29) by ADMM [4], we define the augmented Lagrangian functional

[TABLE]

where $\beta_{r},\beta_{t}>0$ are scalar penalty parameters and $\lambda_{r}\in\mathbb{R}^{n}$, $\lambda_{t}\in\mathbb{R}^{2n}$ are the vectors of Lagrange multipliers associated with the linear constraints $r=Ku-g$ and $t=Du$ in (29), respectively.

Solving (28)–(29) is thus equivalent to seek for the solutions of the following saddle point problem:

[TABLE]

with the augmented lagrangian functional $\mathcal{L}$ defined in (30) and where, for simplicity of notations, we set $x:=(u,r,t)$, $\lambda:=(\lambda_{r},\lambda_{t})$, $X:=\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{2n}$ and $\Lambda:=\mathbb{R}^{n}\times\mathbb{R}^{2n}$.

Given the previously computed (or initialized for $k=0$) vectors $u^{(k)}$, $\lambda_{r}^{(k)}$ and $\lambda_{t}^{(k)}$, the $k$-th iteration of the proposed ADMM-based iterative scheme applied to the solution of the saddle-point problem (31) - minimization for the primal variables $u,r,t$, maximization for the dual variables $\lambda_{r},\lambda_{t}$ - reads as follows:

[TABLE]

In the following three subsections we describe how to solve the minimization sub-problems (32), (33) and (34) for the primal variables $r$, $t$ and $u$, respectively, in both cases $q\in\{1,2\}$. In particular, we remark that thanks to the preliminary ADMM variable splitting procedure, sub-problems (33) and (34) for the variables $t$ and $u$ are identical in the two cases $q\in\{1,2\}$ and their solution can be obtained based on formulas given in [9] for the same sub-problems.

4.1 Minimization sub-problem for the primal variable $r$

Recalling the definition of the augmented Lagrangian functional in (30) and carrying out some simple algebraic manipulations, the minimization sub-problem (32) for the primal variable $r$ can be written as

[TABLE]

with the constant (with respect to the optimization variable $r$) vector $v^{(k)}\in\mathbb{R}^{n}$ given by

[TABLE]

Since $\mu\geq 0$, $\beta_{r}>0$, in both cases $q\in\{1,2\}$ the cost function in (37) is strongly convex, hence it admits a unique global minimizer. In particular, the unique solution $r^{(k+1)}$ of (37) can be computed, depending on $q$, by means of the following closed-form formulas:

[TABLE]

where $\mathrm{sign}(\,\cdot\,)$ and $|\,\cdot\,|$ in (39) denote the component-wise signum and absolute value functions and $\,\odot$ indicates the component-wise vectors product. We remark that formula (39) represents a well-known component-wise soft-thresholding operator - see e.g. [12] - whereas (40) comes easily from first-order optimality conditions of (37).

In case that the regularization parameter $\mu$ is regarded as a constant - that is, it is fixed a priori - then formulas (39)–(40) allow to determine very efficiently the solution $r^{(k+1)}$ of this sub-problem. However, as previously stated, in the case $q=2$ we aim also at automatically adjusting $\mu$ along iterations - that is, $\mu$ becomes $\mu^{(k)}$ - such that the final solution $u^{*}$ of our model (7)–(8) satisfies the discrepancy principle [14]. To this aim, in the following we propose a procedure which builds upon those presented in [6, 10] but, due to a different ADMM initial variable splitting, needs to be adapted and is worth to be outlined in detail.

We consider the discrepancy associated with the solution $r^{(k+1)}$ in (40) as a function $\delta^{(k+1)}:[0,+\infty[\rightarrow[0,+\infty[$ of the regularization parameter $\mu$:

[TABLE]

where the second equality comes from (40). The discrepancy function in (41) is clearly continuous, non-negative and monotonically decreasing over its entire domain $\mu\in[0,+\infty[$ and at the extremes we have $\,\delta^{(k+1)}(\mu=0)=\|v^{(k)}\|_{2}$, $\,\delta^{(k+1)}(\mu\to+\infty)=0$. In order to set a value $\mu^{(k+1)}$ such that the discrepancy principle is satisfied here for the auxiliary variable $r$ (recall that $r=Ku-g$ represents the residue of the restoration), we consider two complementary cases based on the value of the norm of the vector $v^{(k)}$ defined in (38).

In case that $\,\|\,v^{(k)}\|_{2}\leq\bar{\delta}$, with $\bar{\delta}$ denoting the noise level, then from (41) and from the fact that $\,0<\beta_{r}/(\beta_{r}+\mu)\leq 1$, it follows that $\,\delta^{(k+1)}(\mu)\;{\leq}\;\bar{\delta}\;\>\forall\,\mu\in[0,+\infty[$, that is the discrepancy principle is satisfied for any $\mu$. In this case we thus set $\mu^{(k+1)}=0$, such that, according to (40), the sub-problem solution is $r^{(k+1)}=v^{(k)}$.

In case that $\,\|v^{(k)}\|_{2}>\bar{\delta}$, the properties of the discrepancy function $\delta^{(k+1)}(\mu)$ in (41) guarantee that there exists a unique value $\mu^{(k+1)}$ of $\,\mu$ such that $\delta^{(k+1)}(\mu^{(k+1)})=\bar{\delta}$. Recalling (41), we have \big{(}\beta_{r}/(\beta_{r}+\mu^{(k+1)})\big{)}\|\,v^{(k)}\|_{2}\;{=}\;\bar{\delta}\>\;\;\;{\Longleftrightarrow}\;\;\;\mu^{(k+1)}=\beta_{r}\big{(}\,\|\,v^{(k)}\|_{2}/\bar{\delta}\>\;{-}\;1\,\big{)}. Replacing this expression for $\mu$ in (40), the sub-problem solution is $\,r^{(k+1)}\;{=}\;\>\bar{\delta}\,v^{(k)}/\|\,v^{(k)}\|_{2}$.

To summarize, the solution of this sub-problem at any iteration $k$ is computed by (39) for the case $q=1$ whereas for the case $q=2$ it is determined as follows:

[TABLE]

4.2 Minimization sub-problem for the primal variable $t$

Given the definition of the augmented Lagrangian functional in (30), the minimization sub-problem for the primal variable $t$ in (33) can be written as follows:

[TABLE]

Note that in (43) the minimized functional is written in explicit component-wise (or pixel-wise) form, with $\left(Du^{(k)}\right)_{i},\left(\lambda_{t}^{(k)}\right)_{i}\in\mathbb{R}^{2}$ denoting the discrete gradient and the Lagrange multipliers at pixel $i$, respectively. Solving the $2n$-dimensional minimization problem in (43) is thus equivalent to solve the $n$ following independent $2$-dimensional problems:

[TABLE]

with the constant vectors $q^{(k)}_{i}\in\mathbb{R}^{2}$ defined by

[TABLE]

The solutions of the $n$ optimization problems in (44) can be obtained based on the results reported in Proposition 1 of [9], that is:

[TABLE]

where, in particular, the shrinkage coefficients $\,\xi_{i}^{(k)}\in[0,1],\;\,i=1,\ldots,n,\,$ are given by formulas (50)–(52) in [9].

The overall computational cost of this subproblem is linear in the number of pixels $n$.

4.3 Minimization sub-problem for the primal variable $u$

As illustrated in [9], the minimization sub-problem (34) for the primal variable $u$ reduces to the solution of the following $n\times n$ system of linear equations

[TABLE]

which is solvable if the coefficient matrix has full-rank, that is if the following condition holds:

[TABLE]

where $\mathrm{Ker}(M)$ denotes the null space of matrix $M$ and [math] is the $n$-dimensional null vector. In our case, condition (48) is satisfied. In fact, $K$ represents a blurring operator, which is a low-pass filter, whereas the regularization matrix $D$ is a first-order difference operator and, hence, is a high-pass filter. Moreover, since $\beta_{t},\beta_{r}>0$, the coefficient matrix in (47) is symmetric positive definite and typically highly sparse. Hence, the linear system in (47) can be solved quite efficiently by the iterative (eventually preconditioned) conjugate gradient method. Moreover, under appropriate assumptions about the solution $u$ near the image boundary, the linear system can be solved even more efficiently. We assume periodic boundary conditions for $u$, so that both $D^{T}D$ and $K^{T}K$ are block circulant matrices with circulant blocks and, hence, the coefficient matrix in (47) can be diagonalized by the 2D discrete Fourier transform (FFT implementation). Provided that the penalty parameters $\beta_{t}$, $\beta_{r}$ are kept fixed during the ADMM iterations, the coefficient matrix in (47) does not change and it can be diagonalized once for all at the beginning. Therefore, at any ADMM iteration the linear system (47) can be solved by one forward 2D FFT and one inverse 2D FFT, each at a cost of $O(n\log n)$.

4.4 ADMM-based iterative scheme

To summarize previous results, in Algorithm 1 we report the main steps of the proposed ADMM-based iterative scheme used to solve the saddle-point problem (30)–(31) and, hence, to compute solutions of the proposed model (7)–(8).

In the field of image and signal processing the ADMM has been one of the most powerful and successful methods for solving various convex or nonconvex optimization problems. In convex settings, numerous convergence results have been established for ADMM as well as its varieties, see for example [18] and references therein. In particular, convergence results cover the proposed TVp,α-Lq models, $q\in\{1,2\}$, in the special case of $p_{i}\geq 1\;\forall\,i$. However, in case that one or more $p_{i}<1$, the ADMM is under nonconvex settings, where there have been a few studies on the convergence properties. To the best of our knowledge, existing convergence results of ADMM for nonconvex problems is very limited to particular classes of problems and under certain conditions of the dual step size [17]. Nevertheless, the ADMM works extremely well for various applications involving nonconvex optimization problems, and this is a practical justification of its wide use.

5 Numerical results

In this section, we evaluate experimentally the performance of the two proposed models TV${}_{p,\alpha}^{\mathrm{sv}}$-Lq, $q\in\{1,2\}$, defined in (7)–(8), when applied to the restoration of gray level images synthetically corrupted by known blur and by AWGN - in the case of TV${}_{p,\alpha}^{\mathrm{sv}}$-L2 model - and SPN or AWLN - in the case of TV${}_{p,\alpha}^{\mathrm{sv}}$-L1 model. In particular, the two proposed models are compared with:

• •

TV-Lq, $\;\;\,q\in\{1,2\}$, defined in (6); see [13], [15];

• •

TVp-Lq, $\;q\in\{1,2\}$, with global $p\in(0,2]$; see [9], [16];

• •

TV${}_{p}^{\mathrm{sv}}$-Lq, $q\in\{1,2\}$, with local $p_{i}\in(0,2]$, $i\in\{1,\ldots,n\}$; see [16].

For what concerns the preliminary estimation of the $p_{i}$ and $\alpha_{i}$ parameters, we directly apply the procedure outlined in Section 3 to the observed corrupted image $g$ for the AWGN and AWLN cases. Instead, for the SPN case, in order to have a robust estimation of the parameters, a preliminary processing by an adaptive filter is required. In particular, we assume that the position of the pixels corrupted by the SPN is known a priori, otherwise it can be easily detected as suggested in [5]. We replace the corrupted pixels with the mean of the non-corrupted pixels of its neighborhood. The size of the neighborhood is variable and depends on the percentage $\mathcal{P}$ of non-corrupted pixels in it. If $\mathcal{P}$ is below a fixed threshold $\overline{\mathcal{P}}$ (usually $\overline{\mathcal{P}}=0.4$), then the neighborhood is enlarged, in order to incorporate a greater number of uncorrupted pixels. The obtained image is then used to compute the $p$-map and the $\alpha$-map. The described strategy has been introduced instead of the classic median filter, whose smoothing effects is quite high. Clearly, the same approach is adopted for the TV${}_{p}^{\mathrm{sv}}$-L1 model to estimate the $p$-map only.

The quality of the observed corrupted images $g$ and of the restored images $u^{*}$ is measured - in dB - by means of the Blurred Signal-to-Noise Ratio

[TABLE]

and the Improved Signal-to-Noise Ratio

[TABLE]

respectively, with $u$ denoting the original uncorrupted image and $E\,[Ku]$ the average intensity of the blurred image $Ku$. In general, the larger the ISNR value, the higher the quality of restoration.

For all the ADMM-based minimization algorithms and for all the tests, the penalty parameters $\beta_{t}$ and $\beta_{r}$ are suitably set. Moreover, for all tests, the ADMM iterations of all the compared algorithms are stopped as soon as two successive iterates satisfy

[TABLE]

For the models with the L2 fidelity term, the regularization parameter $\mu$ has been automatically set based on the discrepancy principle. For the models with the L1 fidelity term, $\mu$ has been hand-tuned independently in each test so as to provide the highest possible ISNR value. In the following, we report numerical results concerning the restoration of blurred images corrupted by AWGN (Example 1) and by impulsive SPN and AWLN (Example 2).

Example 1: restoration of images corrupted by blur and AWGN. In this example, we evaluate experimentally the performance of the proposed TV${}_{p,\alpha}^{\mathrm{sv}}$-L2 model on a purely piecewise constant test image - geometric ($256\times 256$), Fig. 4(a) - and a partially textured test image - skyscraper ($256\times 256$), Fig. 4(d). Both images have been synthetically corrupted by a Gaussian blur of parameters band=5 and sigma=1.0 and by AWGN characterized by different noise levels. The $p,\alpha$-maps have been computed by using neighborhoods of size $s=3$.

In Table 1 the performance of our model are compared in terms of achieved ISNR values with those of the TV-L2, TVp-L2 and TV${}_{p}^{\mathrm{sv}}$-L2 models. The good quality of the restored image by our model can be appreciated by a visual inspection of Figs. 4(c),(f) and by comparing the ISNR values reported in Table 1.

Example 2: restoration of images corrupted by blur and SPN or AWLN. In this example we evaluate the performance of the proposed TV${}_{p,\alpha}^{\mathrm{sv}}$-L1 model on three medical test images lungs ($468\times 591$), Fig. 5 (a), ecography ($401\times 511$), Fig. 6 (a), and aneurism ($701\times 766$), Fig. 7 (a), synthetically corrupted by Gaussian blur of parameters band=5 and sigma=1 and by two types of impulsive noise, namely SPN and AWLN. The images are provided in the repository at https://medpix.nlm.nih.gov.

First, for what concerns corruptions by SPN, in Figs. 5, 6, 7 we report for the three considered test images the original and corrupted image together with the estimated $p,\alpha$-maps in the first column (with the size $s$ of the neighborhoods used for the $p,\alpha$-maps estimation reported in the captions), the restoration results, obtained by the four compared methods, in the second column (with the achieved ISNR values in the captions) and a zoomed detail of the restored images - green- bordered in Figs. 5 (a), 6 (a), 7 (a) - in the last column.

The reported ISNR values as well as the visual inspection of the restored images and of the zoomed details strongly indicate how the proposed space-variant regularizer allows for higher quality restorations. In particular, it is worth remarking how, with respect to the space-variant TV${}_{p}^{\mathrm{sv}}$ model, the additional degrees of freedom represented by the scale parameters $\alpha_{i}$ used in our proposal, yield a sufficient additional flexibility for avoiding unwanted spurious effects - see, e.g., spikes in Figs. 5 (i), 6 (i), 7 (i).

In the second part of this example, we consider the restoration of the same three medical test images corrupted by the same blur of parameters band=5 , sigma=1 and by a different impulsive noise, namely AWLN of level yielding BSNR=10. In Table 2 we report the ISNR values achieved by the compared methods and in Fig. 8 we show the original images, the corrupted images and the restored images by our model. The results in Table 2 confirm that, also in case of images corrupted by AWLN, the proposed TV${}_{p,\alpha}^{\mathrm{sv}}$-L1 model outperforms its competitors in terms of ISNR. Moreover, the restored images depicted in the last column of Fig. 8 provide further evidence of the good quality restorations achievable by our proposal.

6 Conclusions

We presented a new space-variant regularizer for variational image restoration based on the assumption that the gradient magnitudes of the target image distribute locally according to a half-Generalized Gaussian distribution. Thanks to the high number of free parameters involved in the regularizer and to the fact that such parameters can be automatically and robustly estimated from the observed image, our proposal exhibits a very high flexibility which potentially allows for an effective modeling of space-variant properties of images. By coupling the proposed regularizer with either the L1 or L2 fidelity terms, we tested our proposal on images corrupted by blur and different types of noise, namely AWGN, AWLN and SPN. The restored images, obtained by means of an efficient ADMM-based numerical algorithm, strongly indicate that the proposed regularizer holds the potential for achieving high quality restoration results for a wide range of target images characterized by different gradient distributions and for the different types of noise considered.

Acknowledgments: Research was supported by the “National Group for Scientific Computation (GNCS-INDAM)” and by ex60 project by the University of Bologna “Funds for selected research topics”.