Space-variant TV regularization for image restoration

TL;DR

This paper introduces space-variant TV regularization models that adaptively estimate local shape parameters to improve image restoration quality for various noise types and blurs.

Contribution

The paper presents a novel space-variant generalization of TV regularization with automatic local parameter estimation, enhancing restoration performance over traditional TV models.

Findings

Effective for images with diverse gradient distributions

Outperforms traditional TV in noisy and blurred conditions

Compatible with L2 and L1 fidelity terms

Abstract

We propose two new variational models aimed to outperform the popular total variation (TV) model for image restoration with L$_2$ and L$_1$ fidelity terms. In particular, we introduce a space-variant generalization of the TV regularizer, referred to as TV$_p^{SV}$, where the so-called shape parameter $p\,$ is automatically and locally estimated by applying a statistical inference technique based on the generalized Gaussian distribution. The restored image is efficiently computed by using an alternating direction method of multipliers procedure. We validated our models on images corrupted by Gaussian blur and two important types of noise, namely the additive white Gaussian noise and the impulsive salt and pepper noise. Numerical examples show that the proposed approach is particularly effective and well suited for images characterized by a wide range of gradient distributions.

Figures

Tables

Table 1. Table 1 : Example 1: ISNR values for different models with different noise level (BSNR) on test images.

		geometric			mandrill
BSNR	TV-L2	TVp-L2	TV${}_{}{}^{\mathrm{sv}}{}_p{}^{}$-L2	TV-L2	TVp-L2	TV${}_{}{}^{\mathrm{sv}}{}_p{}^{}$-L2
20	7.77	7.92	8.36	1.38	1.64	1.78
30	9.01	9.87	10.30	2.90	3.04	3.31
40	11.58	12.98	13.47	5.32	5.56	6.09

Equations

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

$$\begin{array}[]{ccc}\mathrm{AWGN:}&&\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}$$

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

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

$$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^{\mathrm{sv}}(u)+\mu {\mathrm{L}}_q(u;g)\right\},q\in \{1,2\},$$

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

$$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{)}\,,$$

$$\{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)$$

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

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

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

$${\lambda }_r^{(k+1)}$$

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

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

$$v^{(k)}=K{u}^{(k)}-g+{\lambda }_r^{(k)}/{\beta }_r.$$

$$\mathrm{case}q=1:r^{(k+1)}$$

$$\mathrm{case}q=2:r^{(k+1)}$$

$$\delta^{(k+1)}(\mu)\,\;{:=}\;\,\|r^{(k+1)}\|_{2}\;{=}\;\>\big{(}\beta_{r}/(\beta_{r}+\mu)\big{)}\,\|\,v^{(k)}\|_{2}\;,$$

$$\begin{array}[]{llll}\|\,v^{(k)}\|_{2}\;{\leq}\;\bar{\delta}&\;\Longrightarrow&\mu^{(k+1)}\;{=}\;0,&\;\;r^{(k+1)}\;{=}\;v^{(k)}\vspace{0.2cm}\\ \|\,v^{(k)}\|_{2}\;{>}\;\bar{\delta}&\;\Longrightarrow&\mu^{(k+1)}\;{=}\;\beta_{r}\big{(}\|\,v^{(k)}\|_{2}/\bar{\delta}-1\big{)},&\;\;r^{(k+1)}\;{=}\;\bar{\delta}\,\,v^{(k)}/\|\,v^{(k)}\|_{2}\end{array}$$

$$\mathrm{BSNR}(g,u)=10{\log }_{10}\Vert Ku-E[Ku]{\Vert }_2^2/\Vert g-Ku{\Vert }_2^2$$

$$\mathrm{ISNR}(g,u,u^{\ast })=10{\log }_{10}\Vert g-u{\Vert }_2^2/\Vert u^{\ast }-u{\Vert }_2^2,$$

Full text

1 Introduction

During the image acquisition and transmission processes, degradation effects such as those due to blur and noise always occur. The goal of image restoration is to eliminate these unwanted effects and to recover clean images from the acquired blurred and noisy ones. We consider grayscale images with rectangular $\,d_{1}\!{\times}\,d_{2}$ domain, such that $n\,{:=}\,d_{1}d_{2}$ is the total number of pixels in the images. The general discrete 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$.

In this paper, we are interested in two important types of noise, namely the additive (zero-mean) white Gaussian noise (AWGN) which typically appears, e.g., in Magnetic Resonance Tomography, 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 images, for these two kinds of noise the general degradation model in (1) reads as

[TABLE]

For what concerns AWGN, the additive corruptions $n_{i}\in\mathbb{R}$, $i\in\Omega$, represent independent realizations from the same univariate Gaussian distribution with zero mean and standard deviation $\sigma$. 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 [10] or it could be estimated [4]. Just like AWGN is 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. 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 $J:\mathbb{R}^{n}\to\mathbb{R}$ such that, typically, restoration is casted as an optimization problem of the form

[TABLE]

where the functionals $R(u)$ and $F(u;g)$, commonly referred to 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. In particular, the functional form of the fidelity term is strictly connected to the characteristics of the noise corruption.

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

[TABLE]

For what regards the regularization term in (3), a very popular choice is represented by the Total Variation semi-norm, 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 (a) it is convex and (b) it allows for restored images with sharp, neat edges. By substituting the TV regularizer (5) and the L2 or L1 fidelity terms (4) for $R$ and $F$ in (3), respectively, one obtains the so-called TV-L2 [12] - or ROF - and TV-L1 [11] restoration models; in formulas:

[TABLE]

The TV-L2 and TV-L1 models in (6) are non-smooth convex and allows to obtain good quality restorations of images corrupted by AWGN and SPN, respectively, such that they are regarded as sort of baseline models. The goal of this paper is to devise two new variational models which are able to outperform the TV-L2 and TV-L1 models, in particular by designing a new, better performing regularizer, and also to propose an efficient minimization algorithm for the solution of these models based on the Alternating Direction Method of Multipliers (ADMM) strategy [3].

The two proposed models are as follows:

[TABLE]

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

[TABLE]

A different value $p_{i}$ for each pixel $i\,$ is thus allowed by the proposed regularizer (8), such that local, space-variant properties of the clean image $u$ can be potentially addressed. The usefulness of this great flexibility is however conditioned to the existence of effective procedures for the automatic estimation of the $p_{i}$ values. As it will be discussed in the paper, the algorithm used in [8] for estimating a unique, global $p$ value is not sufficiently robust to be used for inferring our local $p_{i}$ values. Hence, in the paper we also propose a new suitable estimation procedure of the $p_{i}$ values based on the statistical inference technique described in [7]. The regularization term in (8) is a space-variant version of the TVp regularizer proposed in  [8] where the estimation of a global fixed $p$-value relied on the gradient magnitudes of the image and such a distribution is in general too rigid for effectively modeling the actual gradient magnitudes distribution of real images. In the proposed model (7)–(8) with $q=2$, we also set automatically the regularization parameter $\mu$ based on the well-known discrepancy principle [13].

The paper is organized as follows. In Sect. 2 we briefly outline the procedure proposed for the automatic estimation of the $p_{i}$ parameters. The ADDM-based minimization algorithm is illustrated in Sect. 3 and numerical results are reported in Sect. 4.

2 Estimation of the space-variant parameters

The method proposed in [8] 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 patches for the estimation of local $p$ values. In the following we briefly outline our proposal based on the statistical inference procedure illustrated in [7].

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 [7] to the local data set consisting of the computed gradient magnitudes in a neighborhood of the pixel $i$. In particular, we use square neighborhoods $N_{\,i}^{\,s}$ of size $s\in\{3,5,\ldots\}$ centered at pixel $i\in\Omega$. Following [7], the values $p_{i}$, $i\in\Omega$, parameters of the Generalized Gaussian Distributions, are estimated as follows:

[TABLE]

where $\mathrm{card}(A)$ denotes the cardinality of set $A$ and where the function $h:\mathbb{R}_{+}^{*}\to\mathbb{R}_{+}^{*}$, referred to as the generalized Gaussian ratio function in [7], is defined by

[TABLE]

with $\Gamma(\,\cdot\,)$ indicating the Gamma function [1]. The function $h$ in (11) 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 (10) can be carried out very efficiently. In Fig. 1 the maps of local $p$ values, obtained with neighborhoods of size $s=3$ (b) and $s=11$ (c) starting from the original test image geometric (a) are shown. Both maps are scaled in the same range for visual comparison. As the size $s$ increases, we acquire different kind of details, 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. 4 numerical experiments have been carried out by computing the $p$-map starting from the corrupted images.

3 Applying ADMM to the proposed model

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 V$ and $t\in Q$, with $V:=\mathbb{R}^{n}$, $Q:=\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 (12)–(13) by ADMM [3], we define the augmented Lagrangian functional

[TABLE]

where $\beta_{r},\beta_{t}>0$ are scalar penalty parameters and $\lambda_{r}\in V$, $\lambda_{t}\in Q$ are the vectors of Lagrange multipliers associated with the linear constraints $r=Ku-g$ and $t=Du$ in (13), respectively. 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 associated with the augmented Lagrangian in (14) - minimization for the primal variables $u,r,t$, maximization for the dual variables $\lambda_{r},\lambda_{t}$ - reads as follows:

[TABLE]

In the following we describe how to solve the minimization sub-problem (15) - in both cases $q\in\{1,2\}$ - for the primal variable $r$ only. In fact, thanks to the preliminary ADMM variable splitting procedure, sub-problems (16)–(17) for the variables $t$ and $u$ are identical in the two cases $q\in\{1,2\}$ and, more importantly, their solution can be obtained based on formulas given in [8] for the same sub-problems.

Solving the sub-problem for $\mathbf{r}$

Recalling definition (14) and carrying out some simple algebraic manipulations, the minimization sub-problem (15) reads as

[TABLE]

with the constant (w.r.t. the optimization variable $r$) vector $v^{(k)}\in V$ given by

[TABLE]

Since $\mu\geq 0$, $\beta_{r}>0$, in both cases $q\in\{1,2\}$ the cost function in (20) is strictly convex and its (unique) global minimizer - that is, the solution $r^{(k+1)}$ of (20) - can be computed, depending on $q$, by means of the following closed-form formulas:

[TABLE]

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

In case that the regularization parameter $\mu$ is regarded as a constant - that is, it is fixed a priori - then formulas (22)–(23) 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 [13]. To this aim, in the following we propose a procedure which builds upon those presented in [5, 9] 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 (23) as a function $\delta^{(k+1)}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ of the regularization parameter $\mu$:

[TABLE]

where the second equality comes from (23). The discrepancy function in (24) is continuous, non-negative and monotonically decreasing over its entire domain $\mu\in\mathbb{R}_{+}$ 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 norm of the vector $v^{(k)}$ in (21).

In case that $\,\|\,v^{(k)}\|_{2}\leq\bar{\delta}$, with $\bar{\delta}$ denoting the noise level, then from (24) 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\mathbb{R}_{+}$, that is the discrepancy principle is satisfied for any $\mu\geq 0$. We set $\mu^{(k+1)}=0$, such that, according to (23), 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)}$ in (24) guarantee that there exists a unique value $\mu^{(k+1)}$ of $\,\mu$ such that $\delta^{(k+1)}(\mu^{(k+1)})=\bar{\delta}$. Recalling (24), 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 (23), 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 (22) for the case $q=1$ whereas for the case $q=2$ it is determined as follows:

[TABLE]

4 Numerical results

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

• •

TV-Lq, $q=1,2$, defined in (6) with $p=1$ fixed,

• •

TVp-Lq, $q=1,2$, with $p\in(0,2]$ fixed.

We remark that the TVp-L2 model has been introduced in [8], whereas the TVp-L1 model has not been proposed before and can be regarded as a further contribution of this paper, together with the automatic selection procedure for the space-variant $p$ parameters. For what concerns the TV${}_{p}^{\mathrm{sv}}$-L1 model in order to have a robust evaluation of the $p$-map, the image is preliminarily processed by an adaptive filter. 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 [4]. We replace the corrupted pixels with the mean of non-corrupted pixels of its neighborhood. The size of the neighborhood is variable and depends on the percentage of non-corrupted pixels in it. The image obtained is then used to compute the $p$-map. The described strategy has been introduced instead of the simple median filter, whose smoothing effects is quite high. 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 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 parameters $\beta_{t}$ and $\beta_{r}$ are suitably set. Usually good choices are $(\beta_{t},\beta_{r})=(1,1),(10,5)$. The iterations of the algorithms are stopped as soon as two successive iterates satisfy $\;\|u^{(k)}-u^{(k-1)}\|_{2}/\|u^{(k-1)}\|_{2}\,\;{<}\;\,10^{-4}$. 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 for that test. In the following, we report numerical results concerning the restoration of images corrupted by AWGN (Example 1) and SPN (Example 2).

Example 1: restoration of images corrupted by AWGN

In this subsection we are testing the performance of TV${}_{p}^{\mathrm{sv}}$-L2 on piecewise constant (geometric ($256\times 256$) Fig. 2(a)) and textured images (mandrill ($512\times 512$) Fig. 2(c)) with different noise levels. In Table 1 the results are compared in terms of ISNR with the ones obtained by TV-L2 and TVp-L2. Both geometric and mandrill images have been corrupted by a Gaussian blur of band=5 and standard deviation sigma=1.0, and by an AWGN, with BSNR=20,30,40. The $p$-maps have been computed by setting the size of the neighborhoods $s=3$. The good quality of the reconstructed images can be appreciated by a visual inspection of Fig. 2(b),(d) and by comparing the ISNR values reported in Table 1.

Example 2: restoration of images corrupted by SPN

In this subsection we report the performance of TV${}_{p}^{\mathrm{sv}}$-L1 on a $200\times 200$ medical image representing a particular of a CT scan of an abdomen - see Fig. 3(a). It has been corrupted by a SPN of level $\gamma=0.35$ and by a Gaussian blur of band=9 and sigma=2.5 (Fig. 3 (b)). The $p$-map in Fig.3(c), computed by setting the size of the neighborhood s=25, presents higher values in the textured regions. A comparison of the methods TV-L1, TVp-L1, TV${}_{p}^{\mathrm{sv}}$-L1 leads to ISNR=$11.81,12.97,13.60$, respectively. The quality of the reconstructed images can be visually appreciated in Fig. 3(d),(e),(f).

5 Conclusions

We have proposed two new variational models which are able to outperform the popular TV model for image restoration with L2 and L1 fidelity terms. In particular, we introduced the TV${}_{p}^{SV}$ regularizer, a space-variant generalization of the popular TV prior, where the shape parameter $p$ is automatically and locally estimated by an effective procedure based on the statistical inference technique in [7]. The restored image is efficiently computed by using an ADMM-based algorithm. Numerical examples show that the proposed approach is particularly effective and well suited for images corrupted by Gaussian blur and two important types of noise, the AWGN and SPN. As future work, we plan to extensively test our models on a new immunofluorescence portable diagnostic systems where low-cost complementary metal-oxide semiconductor (CMOS) sensors are used. In this system different noise sources affect a noise-free image acquired by the CMOS-based imaging system:the Photo-Response Non-Uniformity is usually modeled as an AWGN while the signal dependent Photon Shot Noise is more properly modeled as a Poisson noise and the Analog-to-Digital Converter noise as an SPN with known positions.