Adaptive parameter selection for weighted-TV image reconstruction problems

TL;DR

This paper introduces an efficient method for automatically selecting locally adaptive Total Variation regularization parameters in image reconstruction, combining local likelihood estimation with a global discrepancy principle, outperforming existing strategies.

Contribution

It presents a novel hybrid approach that integrates local maximum-likelihood estimation with a global discrepancy principle for adaptive TV regularization parameter selection.

Findings

Outperforms state-of-the-art parameter estimation methods

Effective in various image reconstruction problems

Demonstrates improved image quality and noise handling

Abstract

We propose an efficient estimation technique for the automatic selection of locally-adaptive Total Variation regularisation parameters based on an hybrid strategy which combines a local maximum-likelihood approach estimating space-variant image scales with a global discrepancy principle related to noise statistics. We verify the effectiveness of the proposed approach solving some exemplar image reconstruction problems and show its outperformance in comparison to state-of-the-art parameter estimation strategies, the former weighting locally the fit with the data (Dong et al. '11), the latter relying on a bilevel learning paradigm (Hinterm\"uller et al., '17)

Tables

Table 1. Table 1 : Maximum ISNR/SSIM values achieved by TV-L 2 , SATV and HWTV-L 2 on the skyscraper test image in Fig. 1(a) (top) corrupted by AWGN of two different levels.

		$\sigma =0.02$				$\sigma =0.05$
	TV	SATV	HWTV		TV	SATV	HWTV
ISNR	3.4701	3.6625	4.3331		1.9433	2.0414	2.5408
SSIM	0.8733	0.8966	0.9007		0.7335	0.7797	0.8099

Equations

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

$${\mathrm{L}}_2(u)$$

$$\underset{u\in {\mathbb{R}}^n}{\min }\left\{\alpha \text{TV}(u)+{\text{L}}_2(u)\right\}\text{or, equivalently,}\underset{u\in {\mathbb{R}}^n}{\min }\left\{\text{TV}(u)+\mu {\text{L}}_2(u)\right\}.$$

$$\underset{u\in {\mathbb{R}}^n}{\min }\left\{\mathrm{WTV}(u)+{\mathrm{L}}_2(u)\right\},$$

$$\underset{u\in {\mathbb{R}}^n}{\min }\left\{\mathrm{TV}(u)+{\mathrm{WL}}_2(u)\right\},$$

$$u^{\ast }=\underset{u\in {\mathbb{R}}^n}{\mathrm{argmin}}\left\{\mathrm{WTV}(u)+\mu {\text{L}}_2(u)\right\}.$$

$$u^{\ast }\in \underset{u\in {\mathbb{R}}^n}{\mathrm{argmax}}P(u|g;K)=\underset{u\in {\mathbb{R}}^n}{\mathrm{argmin}}\left\{-\log P(g|u;K)-\log P(u)\right\},$$

$$P(u)=\frac{1}{Z}\text{exp}\left(-\sum _{i=1}^n{\alpha }_i||(\mathrm{D}u{)}_i|{|}_p\right)=\frac{1}{Z}\text{exp}\left(-\mathrm{WTV}(u)\right),p\in \left\{1,2\right\},$$

$$P(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 },w^{\ast },t^{\ast }\}$$

$$\mathcal{L}(u,w,t;{\rho }_w,{\rho }_t;{\alpha }_1,\dots ,{\alpha }_n,\mu )$$

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

$$q_i^{(k)}=(\mathrm{D}{u}^{(k)}{)}_i+\frac{1}{{\beta }_t}({\rho }_t^{(k)}{)}_i\in {\mathbb{R}}^2,$$

$$t_i^{(k+1)}=q_i\max \left(1-\frac{{\alpha }_i^{(k)}}{{\beta }_t\Vert q_i^{(k)}{\Vert }_p},0\right),p\in \{1,2\},i=1,\dots ,n.$$

$$z^{(k)}=K{u}^{(k)}-g+\frac{1}{{\beta }_w}{\rho }_w^{(k)},$$

$$w^{(k+1)}=\frac{{\beta }_w}{{\mu }^{(k)}+{\beta }_w}{z}^{(k)}.$$

$$\left({\mathrm{D}}^T\mathrm{D}+\frac{{\beta }_w}{{\beta }_t}{K}^{T}K\right)u={\mathrm{D}}^T\left(t^{(k+1)}-\frac{1}{{\beta }_t}{\rho }_t^{(k)}\right)+\frac{{\beta }_w}{{\beta }_t}{K}^T\left(w^{(k+1)}-\frac{1}{{\beta }_w}{\rho }_w^{(k)}+g\right),$$

$$\mathcal{L}(\alpha_{i};\mathcal{S}_{i})\;=\;\prod_{x_{i,j}\in\mathcal{S}_{i}}P(x_{i,j};\alpha_{i})=\prod_{j=1}^{N}P(x_{i,j};\alpha_{i})\;=\;\alpha_{i}^{N}\exp\bigg{(}-\sum_{j=1}^{N}\alpha_{i}x_{i,j}\bigg{)}.$$

$$\alpha_{i}=\bigg{(}\frac{1}{N}\sum_{j=1}^{N}x_{i,j}\bigg{)}^{-1},$$

$$\Vert z^{(k)}{\Vert }_2\le \delta$$

$$\Vert z^{(k)}{\Vert }_2>\delta$$

$${\alpha }_i^{(k)}\text{ by }\text{eqref}est\text{ for every }i=1,\dots ,n,$$

$${\mu }^{(k)}\text{ by }\text{eqref}eq:updat{e}_{m}u,$$

$$t^{(k+1)}=\underset{t}{\mathrm{argmin}}\mathcal{L}(u^{(k)},w^{(k)},t;{\rho }_w^{(k)},{\rho }_t^{(k)};{\alpha }_1^{(k)},\dots ,{\alpha }_n^{(k)},{\mu }^{(k)})$$

$$w^{(k+1)}=\underset{r}{\mathrm{argmin}}\mathcal{L}(u^{(k)},w,t^{(k+1)};{\rho }_w^{(k)},{\rho }_t^{(k)};{\alpha }_1^{(k)},\dots ,{\alpha }_n^{(k)},{\mu }^{(k)})$$

$$u^{(k+1)}=\text{nolinebreak}[4]\underset{u}{\mathrm{argmin}}\mathcal{L}(u,w^{(k+1)},t^{(k+1)};{\rho }_w^{(k)},{\rho }_t^{(k)};{\alpha }_1^{(k)},\dots ,{\alpha }_n^{(k)},{\mu }^{(k)})$$

$$\displaystyle\rho_{w}^{(k+1)}=\rho_{w}^{(k)}-\beta_{w}\big{(}w^{(k+1)}-(Ku^{(k+1)}-g)\big{)},$$

$$\displaystyle\rho_{t}^{(k+1)}=\rho_{t}^{(k)}-\beta_{t}\big{(}t^{(k+1)}-\mathrm{D}u^{(k+1)}\big{)},$$

Full text

1 Introduction

In this paper, we are interested in restoring images corrupted by known blur and additive white Gaussian noise (AWGN), i.e. we assume a degradation model of the form $\,g\;{=}\;Ku+\epsilon,$ with $g,u,\epsilon\in\mathbb{R}^{n}$ vectorised forms of the discrete observed image, target image and noise realisation, respectively, and with $K\in\mathbb{R}^{n\times n}$ the blur operator. Such inverse problem is typically ill-posed.

The variational approach to solve ill-posed image restoration problems consists in minimising a composite functional which is the sum of a regularisation term encoding a-priori assumptions on the unknown image $u$ and a data fitting term describing noise statistics. A very popular edge-preserving regulariser firstly proposed in [15] for image denoising is the Total Variation (TV) semi-norm, while the so-called L2 data term is known to be suited for AWGN. They read

[TABLE]

where $(\mathrm{D}_{h}u)_{i},(\mathrm{D}_{v}u)_{i}$ denote the horizontal and vertical discrete gradient components at pixel $i$, respectively. The two instances $p\>{=}1,2$ in (TV) are referred to as anisotropic and isotropic TV, respectively. By taking a weighted average of (TV) and (L2), one gets the TV-L2 model:

[TABLE]

Both the parameters $\alpha$ and $\mu=1/\alpha$ are often referred to as regularisation parameters since their size weights the amount of the regularisation against the trust in the data. Note that the equivalence of the two formulations in (TV-L2) allows in fact to use indifferently either of the two models. To estimate an optimal regularisation parameter, several strategies can be used. When the noise level is known, a classical approach is based on the use of the discrepancy principle [5], while in blind scenarios, optimisation techniques learning the optimal amount of regularisation from training data can be used, see, e.g., [2] and the references therein.

In order to overcome the well-known artefacts of TV-based reconstructions, higher-order [1] and/or locally-adaptive anisotropic regularisers [3, 11, 14] have been proposed in the literature. A simple, though powerful, extension enforcing the TV smoothing to locally adapt to the underlying image structures (such as texture, cartoon…) consists in weighting at any pixel the amount of regularisation [8, 9, 7, 16] or data fit [4]. This reflects in considering two locally-weighted models, referred to as WTV-L2 and TV-WL2, which represent space-variant extensions of the two equivalent formulations of the (TV-L2) model and read as

[TABLE]

where, note, the positive parameters $\alpha_{i}$ and $\mu_{i}$ are now space-variant and both control locally the amount of smoothing. The estimation of parameters $\alpha_{i}$ in (WTV-L2) has been done by inferring local geometries [16] or by means of computationally expensive bilevel-optimisation approaches [8, 9], whereas parameters $\mu_{i}$ in (TV-WL2) have been estimated based on the use of a local discrepancy principle [4]. We remark that despite the aforementioned equivalence of the two approaches in the (TV-L2) scalar parameter case, the two locally-weighted (WTV-L2) and (TV-WL2) models do show significant differences when used for image reconstruction problems as it has been rigorously studied in [7].

Contribution

We propose an image restoration approach based on an hybrid version of the two space-variant (WTV-L2) and (TV-WL2) variational models, with variable regularisation parameters $\alpha_{i}$ and global fidelity parameter $\mu$ referred to as HWTV-L2 model, see (HWTV-L2). We propose a simple yet effective automatic Maximum-Likelihood (ML) estimation procedure of the $\alpha_{i}$ weights in the WTV regulariser as well as with the use of a standard discrepancy principle. The statistical prior assumption motivating our ML estimation approach is that image gradients norms are locally drawn from an half-Laplacian distribution with space-variant scale parameters $\alpha_{i}$. The local closed-form formula obtained by our ML approach is extremely handy and, together with a minimisation algorithm based on the Alternating Directions Method of Multipliers (ADMM), renders our proposal very efficient. The proposed approach outperforms by far both the the classical TV-L2 and the SATV [4] restoration methods in terms of standard image quality indexes (ISNR, SSIM) and, compared to recent bilevel optimisation strategies used in [8, 9] to estimate local TV weights, it slightly improves the restoration quality while at the same time being much more efficient.

2 The proposed hybrid space-variant model

For a given image $g\in\mathbb{R}^{n}$ corrupted by AWGN and blur generated by a known blur operator $K\in\mathbb{R}^{n\times n}$, we propose the following variational model for image restoration

[TABLE]

We note that in addition to the space-variant parameters $\alpha_{i}$ contained in the WTV regulariser, a further scalar data weight $\mu>0$ appears in (HWTV-L2). This makes our model an hybrid version of the two space-variant (WTV-L2) and (TV-WL2) models, where local parameters $\alpha_{i}$ describing local image scales in a statistical sense (see Section 3) are used together with global parameter $\mu$ which codifies the discrepancy w.r.t. the given AWGN level. The redundancy of such parameter is therefore only apparent in (HWTV-L2) as its value is computed depending on the global noise statistics in comparison with the local regularisation strength encoded by the parameters $\alpha_{i}$.

2.1 Statistical derivation

We now justify the choice of the space-variant WTV regulariser in (WTV-L2) by means of statistical arguments.

A common paradigm in image restoration is the Maximum A Posteriori (MAP) approach by which the restored image is obtained as a global minimiser of the negative log-likelihood distribution of the observed image $g$ given the blurring operator $K$ combined with some prior PDF $P(u)$ on the unknown target image $u$. In formulas:

[TABLE]

where equality comes from the Bayes’ formula after dropping the normalisation term $P(g)$. The terms $P(g|u;K)$ and $P(u)$ in (1) are commonly referred to as the likelihood and prior distribution, respectively: they encode available information on the statistics of the noise and on the solution we seek, respectively.

A standard prior model for the gradient magnitude of the unknown image $u$ is TV-Gibbs prior so that $P(u)=\frac{1}{Z}\exp\left(-\alpha\sum_{i=1}^{n}\|(\mathrm{D}u)_{i}\|_{p}\right)$, where $p\in\{1,2\}$. As pointed out in [12], such choice can be equivalently interpreted by saying that each $\|(\mathrm{D}u)_{i}\|_{p}$ distributes according to an half-Laplacian PDF with parameter $\alpha>0$. The use of a one-parameter distribution may be restrictive in modelling images with local properties at different scales (edges, texture…). To allow more flexibility, in [12] a space-variant model where gradient norms distribute according to a half-Laplacian distribution with locally varying scale parameter $\alpha_{i}>0$ has been proposed. The prior associated to such choice is

[TABLE]

where WTV is the regulariser defined in (WTV-L2).

For the sake of completeness, we recall that in the case of AWGN the likelihood term in (1) takes the following form

[TABLE]

where $\sigma>0$ denotes the AWGN standard deviation and $W>0$ is a normalisation constant. By plugging the expression (2) for $P(u)$ and (3) $P(g|u;K)$ in (1), we derive the variational model (WTV-L2).

In our modelling, in order to describe local image features together with global noise discrepancy, we further weight the data fitting term by a global parameter $\mu$, thus obtaining the hybrid reconstruction model (HWTV-L2).

3 ADMM optimisation & automatic parameter selection

In order to solve numerically the image restoration problem (HWTV-L2), we use in the following an ADMM-based algorithm combined with an adaptive estimation procedure of model parameters along the iterations. To do so, we introduce two auxiliary variables $w\in\mathbb{R}^{n}$ and $t\in\mathbb{R}^{2n}$ and rewrite the model in the following constrained form:

[TABLE]

We first 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. The solution $(u^{*},w^{*},t^{*})$ of problem (6) is a saddle point for $\mathcal{L}$ in (5). Hence, we can alternate a minimisation step with respect to the primal variables $t,u,w$ with a maximisation step with respect to the dual variables $\rho_{t},\rho_{w}$, in combination with an iterative update of the space variant parameters $\alpha_{i}$ and $\mu$, which hence will be denoted by $\alpha^{(k)}_{i}$ and $\mu^{(k)}$. In particular, for what concerns $\alpha_{i}^{(k)}$ we use the easy ML estimation strategy described next, whereas for $\mu^{(k)}$ we will rely on a global discrepancy principle

Primal variables update.

The three primal sub-problems can be solved efficiently and in closed-form by simple shrinkage/projection operators (in both cases $p=1,2$) and linear system solvers - see [13, 12]. More in details, the sub-problem with respect to the $t$ primal variable, after some algebraic manipulations, can be written as,

[TABLE]

Denoting by,

[TABLE]

the solution of each one-dimensional separable problem is given by

[TABLE]

Introducing

[TABLE]

we have that the updating formula for $w$ reads:

[TABLE]

Imposing a first order optimality condition with respect to the primal variable $u$, leads to the following linear system,

[TABLE]

that can be solved since the coefficient matrix has full rank - see, e.g, [3].

Parameters update.

In Algorithm 1 both the local space-variant parameters $\alpha_{i}$ and the global parameter $\mu$ are updated along the iterations. This is a standard strategy for this type of optimisation problems (see, e.g., [6, 3]), especially in the case of a cheap update of parameters adapting to the image quality improvement as the one considered in the following. Despite the iterative change in the expression of the cost functional corresponding to such update, we remark that nonetheless we observed empirical convergence of the algorithm. A theoretical proof of such result is left for future research.

For any pixel $i=1,\ldots,n$, we consider the set $\mathcal{S}_{i}:=\{x_{i,j}\}_{j=1}^{N}$ with $x_{i,j}=\|(\mathrm{D}u^{(k)})_{j}\|_{p}$, where $(\mathrm{D}u^{(k)})_{j}$ are gradients in the square neighbourhood $\mathcal{N}_{i}^{r}$ centred in $i$ with side $2r+1$ whose norm is drawn from a half-Laplacian distribution with scale parameter $\alpha_{i}$. The likelihood function of $\alpha_{i}$ thus reads:

[TABLE]

We now look for an $\alpha_{i}>0$ maximising $\mathcal{L}$, or equivalently, minimising $\mathcal{F}(\alpha_{i};\mathcal{S}_{i}):=-\log\leavevmode\nobreak\ \mathcal{L}(\alpha_{i};\mathcal{S}_{i})$. By imposing first order optimality on $\mathcal{F}$ with respect to $\alpha_{i}$, we obtain the closed formula:

[TABLE]

which can be handily updated along the iterations $k\geq 0$ to estimate the local regularisation parameters $\alpha^{(k)}_{i}$ at each pixel $i=1,\ldots,n$ by taking as samples $x^{(k)}_{i,j}=\|(\mathrm{D}u^{(k)})_{j}\|_{p}$, $j=1,\ldots,N$ i.e. the norms of the image gradients in the neighbourhood $\mathcal{N}_{i}^{r}$. We remark that the estimates $\alpha_{i}^{(k)}$ in (13) can be efficiently computed based on 2D convolution (realised by a fast 2D discrete transform) of the map of gradient norms with a square $(2r+1)\times(2r+1)$ averaging kernel.

The parameter $\mu$ is updated along the iterations so as to fulfil the global discrepancy principle as described in [6]: we ask each iterate $u^{(k)}$ to satisfy the condition $||Ku^{(k)}-g||_{2}\leq\delta:=\tau\sigma\sqrt{n}$, where $\sigma$ is the noise standard deviation and the parameter $\tau\approx 1$ is set a priori (see Section 4 for some experiments describing the sensitivity of the model to this parameter). In particular, recalling the definition of $z^{(k)}$ given in (9), the update reads:

[TABLE]

The ADMM pseudo-code is reported in Algorithm 1.

4 Numerical results

In this section we report some numerical results obtained by solving the image reconstruction model (HWTV-L2) via the ADMM Algorithm 1 with fixed penalty parameters $\beta_{t}=20$ and $\beta_{w}=100$. In our experiments we observed that the convergence properties of the algorithm are not affected by this choice, if not in terms of convergence speed. The value $r>0$ denotes the radius of the neighbourhoods $\mathcal{N}_{i}^{r}$ defined in Sect. 3 and used to estimate the space-variant parameters $\alpha_{i}$. Denoting by $u\in[0,1]^{n}$ the ground-truth image, we assess the quality of the reconstruction $u^{*}$ by means of the Improved Signal-to-Noise Ratio $\mathrm{ISNR}(g,u,u^{*}):=10\log_{10}(\|g-u\|_{2}^{2}/\|u^{*}-u\|_{2}^{2})$ and in terms of the Structural Similarity Index (SSIM). We compare our results with the ones obtained by the standard (TV-L2) model, the SATV approach [4] based on coupling the (TV-WL2) model with a local discrepancy-based procedure for automatically selecting the parameters $\mu_{i}$ 111We used the MATLAB code available at: https://www.math.hu-berlin.de/~hp_hint/software/satv.html. and the bilevel learning strategy used in [8, 9] to estimate the parameters $\alpha_{i}$ of (WTV-L2) model via a nested optimisation procedure.

Image deblurring.

We consider the skyscraper test image ($256\times 256$) corrupted by AWGN of two levels $\sigma=0.02,0.05$, and Gaussian blur of band = $5$ and sigma = $1$. The ground-truth image $u$ and the observed image $g$ for $\sigma=0.05$ are shown in Fig.1(a)-1(b), respectively. In this test we highlight the improvements obtained by our (HWTV-L2) Algorithm 1 in comparison to the standard (TV-L2) model and the SATV method, both solved by means of ADMM for comparisons. As mentioned above, for the automatic adaption of the parameter $\mu$ along the iterations, a value for the parameter $\tau$ needs to be chosen. For the three models considered, we observed that the value of $\tau$ maximising the ISNR does not necessarily correspond to the value maximising the SSIM, see Fig.2(a). For the TV-L2 model the maximum SSIM is reached for $\tau\approx 1$, while the ISNR achieves its maximum when $\tau\approx 0.9$, the latter being the case in which texture is better preserved but noise is not completely removed. For SATV, the maximum ISNR and SSIM values are reached approximately for the same $\tau$. As remarked in [4], the SATV method is robust w.r.t. the choice of the radius $r$ of the neighbourhoods used for the estimation. Thus, we set such parameter as the default value $r=5$ in our tests. We performed similar sensitivity tests for our HWTV-L2 model for different $(\tau,r)$ values. Results are shown in Figs.2(c)-2(d). For each method, we then selected the parameter(s) yielding the maximum ISNR/SSIM values and compared the results obtained. In Table 1 we report the achieved ISNR/SSIM values, whereas in Figs.3-4 we show the associated restored images for the case of AWGN with $\sigma=0.05$ - see Fig.1(a)(bottom). We observe that our HWTV-L2 method results in higher quality reconstructions w.r.t. TV-L2 and SATV. Visual inspection confirms the effectiveness of our approach in distinguishing between textured and homogeneous regions, see Figs.3(g)-4(g).

Image denoising.

We now consider the test image turtle222Photo courtesy of K. Papafitsoros. (150$\times$200) corrupted by AWGN of level $\sigma=0.1$ (see Figs.5(a) -5(b)) and focus on the quality and computational improvements of our HWTV-L2 method in the case of anisotropic TV (i.e. $p=1$ in (TV)) in comparison to the alternative bilevel optimisation strategy used [8, 9] for estimating the space-variant parameters $\alpha_{i}$. After optimising the HWTV-L2 method over $\tau$ as discussed above, the maximum achieved value is SSIM = 0.7708 (for $r=40$, $\tau=0.86$), in comparison to SSIM = 0.7602 obtained by using a bilevel optimisation strategy. The reconstructions are shown in Fig.5(c)-5(d). We remark that, in addition to the obtained SSIM and visual improvements, our approach exhibits a very high computational efficiency, whereas bilevel codes are known to be computational expensive and hardly applicable to high-resolution images. For instance, in this experiments the proposed ADMM Algorithm 1 for the HWTV-L2 model required only 40 seconds on a standard laptop, compared to the 1429 seconds required by the bilevel algorithm [9].

5 Conclusions

We proposed an image restoration method based on an hybrid, locally-weighted TV-L2 variational model where local regularisation parameters are combined with a global data fidelity weight. Numerically, we solve the model by means of an ADMM-type algorithm combined with an effective and efficient automatic update of the local parameters via ML estimation and a global discrepancy constraint. Compared to standard as well as state-of-the-art competing models, the proposed approach outperforms in terms of standard image quality measures (ISNR, SSIM) as well as computational efficiency.

Acknowledgements

LC acknowledges the support of the Fondation Mathématique Jacques Hadamard (FMJH). LC and MP are thankful to the organisers of the special trimester The mathematics of imaging held at the IHP (Paris, France) where part of this research was carried out and to G. Peyré for the financial support provided within the ERC project NORIA. Research of AL, MP and FS was supported by the “National Group for Scientific Computation (GNCS-INDAM)” and by the ex60 project “Funds for selected research topics”. The authors are deeply grateful to K. Papafitsoros for the computation of the bilevel WTV reconstruction used for comparison in Fig. 5.