Space-variant Generalized Gaussian Regularization for Image Restoration

Space-v arian t Generalized Gaussian Regularization for Image Restoration A Lanza a , S. Morigi a , M. Pragliola a and F. Sgallari a a Dep artment of Mathematics, University of Bolo gna, Piazza di Porta San Donato 5, Bolo gna, IT Abstract W e propose a new space-v ariant regularization term for v ariational image restoration based on the assumption that the gradient magnitudes of the target image distribute lo cally according to a half-Generalized Gaussian distribution. This leads to a highly flexible regularizer characterized b y t wo p er-pixel free parameters, whic h are automatically estimated from the observ ed image. The prop osed regularizer is coupled with either the L 2 or the L 1 fidelit y terms, in order to effectively deal with additiv e white Gaussian noise or impulsive noises suc h as, e.g, additive white Laplace and salt and p epp er noise. The restored image is efficiently computed b y means of an iterative n umerical algorithm based on the alternating direction metho d of multipliers. Numerical examples indicate that the prop osed regularizer holds the p otential for achieving high quality restorations for a wide range of target images characterized by different gradien t distributions and for the differen t types of noise considered. 1 In tro duction Image restoration refers to the reco very of a clean sharp image from a noisy , and p otentially blurred, observ ation. In this pap er, we consider the problem of restoring images corrupted by known blur and differen t types of noise. W e consider gra y lev el images of size d 1 × d 2 , suc h that n := d 1 d 2 is the total num b er of pixels in the images. The general mo del of the image degradation pro cess under blur and noise corruptions can b e written as g = N ( K u ) , (1) where u, g ∈ R n represen t vectorized forms of the unknown clean image and of the observed corrupted image, resp ectively , K ∈ R n × n is a kno wn linear blurring op erator and N ( · ) denotes the noise corruption op erator, which in most cases is of random nature. Giv en K and g , the goal of image restoration is to solv e the ill-conditioned - or ev en singular, dep ending on K - in verse problem of recov ering an as accurate as possible estimate u ∗ of the unkno wn clean image u . The class of variational metho ds for image restoration relies on determining restored images u ∗ ∈ R n as the minimizers of suitable cost functionals J ( u ) : R n → R such that, t ypically , restoration is casted as an optimization problem of the form u ∗ ← arg min u ∈ R n J ( u ) , J ( u ) := R ( u ) + µ F ( u ; g ) , (2) where the functionals R ( u ) and F ( u ; g ) , commonly referred as the r e gularization and the fidelity term, enco de prior information on the clean image u and on the observ ation mo del (1), resp ectively , with the so-called regularization parameter µ > 0 con trolling the trade-off b etw een the t w o terms. The functional form of the fidelity term is strictly connected with the c haracteristics of the noise corruption. In this pap er, we are interested in three imp ortant types of noise, namely the additive (zero- mean) white Gaussian noise (A W GN) which typically app ears, e.g., in Magnetic Resonance T omography , the additive (zero-mean) white Laplace noise (A WLN) and the impulsive salt and p epp er noise (SPN) usually due to transmission errors or malfunctioning pixel elemen ts in camera sensors. Denoting by 1 Ω := { 1 , . . . , n } the set of all pixel p ositions in the vectorized images, for these three kinds of noise the general degradation model in (1) reads as A W GN and A WLN : SPN : g i = ( K u ) i + n i ∀ i ∈ Ω , g i =  ( K u ) i for i ∈ Ω 0 ⊆ Ω n i ∈ { 0 , 1 } for i ∈ Ω 1 := Ω \ Ω 0 . (3) F or what concerns A WGN and A WLN, the additive corruptions n i ∈ R , i ∈ Ω , represent indep endent realizations from the same univ ariate Gaussian and Laplace distribution with zero mean and standard deviation σ , resp ectiv ely . In the case of SPN, only a subset Ω 1 of the pixels is corrupted by noise, whereas the complementary subset Ω 0 is noise-free. In particular, the corrupted pixels can take only the t wo p ossible extreme v alues { 0 , 1 } (we assume that images ha ve range [0 , 1] ), with the same probability . The subset Ω 1 is known in some applications [11] or it could be estimated [5]. As the zero-mean A W GN and A WLN are fully c haracterized from a probabilistic p oin t of view by the unique scalar parameter σ , SPN is characterized by the parameter γ ∈ [0 , 1] which represents the probability for a pixel to be noise-corrupted. It is well known that A WGN and the impulsive A WLN and SPN are suitably dealt with by the so- called L 2 and L 1 fidelit y terms, which are related to the ` 2 and ` 1 norm of the residue image, resp ectively [15]; in form ulas: F ( u ; g ) = L q ( u ; g ) := 1 q k K u − g k q q , q ∈ { 1 , 2 } . (4) F or what regards the regularization term in (2), a v ery popular choice is represen ted by the T otal V ariation (TV) semi-norm [13], that is R ( u ) = TV( u ) := n X i =1 k ( ∇ u ) i k 2 , (5) where ( ∇ u ) i :=  ( D h u ) i , ( D v u ) i  T ∈ R 2 denotes the discrete gradien t of image u at pixel i , with D h , D v ∈ R n × n linear op erators representing finite difference discretizations of the first-order horizontal and v ertical partial deriv atives, resp ectively . Popularit y of TV regularizer for image restoration is mainly due to tw o facts, namely it is con vex and allows for restored images with sharp, neat edges. By substituting the TV regularizer (5) and the L 2 or L 1 fidelit y terms (4) for R and F in (2), resp ectively , one obtains the so-called TV-L 2 [13] - or R OF - and TV-L 1 [12] restoration models, which reads as u ∗ ← arg min u ∈ R n { TV( u ) + µ L q ( u ; g ) } , q ∈ { 1 , 2 } . (6) The TV-L 2 and TV-L 1 mo dels in (6) are non-smo oth conv ex and allo w to obtain go o d qualit y restorations of images corrupted by A WGN and A WLN/SPN, resp ectively , such that they are regarded as sort of baseline models. The con tribution of this pap er consists of a new space-v arian t regularization term which, coupled with the L 2 or L 1 fidelit y terms, gives rise to a generalization of the models in (6) of the form u ∗ ← arg min u ∈ R n  TV sv p,α ( u ) + µ L q ( u ; g )  , q ∈ { 1 , 2 } , (7) where the new space-v ariant TV sv p,α regularizer is defined b y TV sv p,α ( u ) := n X i =1 α i k ( ∇ u ) i k p i 2 , α i ∈ ]0 , + ∞ [ , p i ∈ ]0 , 2] ∀ i ∈ Ω . (8) The prop osed regularizer in (8) is highly flexible as it is characterized b y t wo p er-pixel free parameters p i , α i , such that local, space-v ariant prop erties of the target clean image u can be p otentially addressed. The usefulness of suc h a great flexibilit y in the prop osed regularizer is ho wev er conditioned to the existence of effective pro cedures for the automatic estimation of the p i and α i parameters. Hence, in this paper we also prop ose a suitable method for the automatic estimation of such parameters from the observ ed image partially based on the statistical inference technique described in [8]. As outlined in Section 2, the rationale of our prop osal is that the distribution of the gradient mag- nitudes of the unkno wn clean image is space-v arian t and it is well modeled lo cally by a t w o-parameters Generalized Gaussian distribution. As highlighted in [9], the TV regularizer in (5) comes from implic- itly assuming a space-inv ariant (that is, frame-based), one-parameter half-Laplacian distribution for the 2 gradien t magnitudes. Based on the observ ation that such a distribution is not sufficiently flexible, the authors in [9] proposed a generalization of the TV regularizer, referred to as TV p , which relies on a space-in v ariant, tw o-parameters half-Generalized Gaussian distribution with p denoting the additional parameter, called the shap e parameter. Finally , the v ery recently prop osed TV sv p regularizer [16] further generalizes the TV p b y assuming a space-v ariant shap e parameter in the tw o-parameters half-Generalized Gaussian distribution. The authors in [16] demonstrated exp erimentally ho w using a lo cal, space-v arian t mo del holds the p otential for ac hieving higher quality restorations than using a global, space-inv arian t mo del. The numerical solution of the tw o prop osed v ariational mo dels (7)–(8) are obtained by means of an efficien t iterative minimization algorithm based on the Alternating Direction Metho d of Multipliers (ADMM) strategy [4]. The pap er is organized as follo ws. The deriv ation of the mo dels (7) and a more detailed motiv ation of their in tro duction are prop osed in Section 2. In Section 3 we briefly outline the statistical inference pro cedure used for automatically estimating the p i and α i parameters. The ADDM-based minimization algorithm is describ ed in detail in Section 4. Some meaningful numerical results are rep orted in Section 5 and, finally , conclusions are drawn in Section 6. 2 Motiv ation via MAP estimator The Maxim um A P osteriori (MAP) Estimation approach [3] relies on the maximization of the p osterior probabilit y density function Pr ( u | g ; K ) asso ciated to the clean unkno wn image u : u ∗ ← arg max u ∈ R n Pr ( u | g ; K ) . (9) Relying on the Bay es’ formula, and dropping the evidenc e term Pr ( g ) , this is equiv alen t to maximize the pro duct of the prior Pr ( u ) and the likeliho o d Pr ( g | u ; K ) probability densit y functions. By taking the negativ e logarithm of this product, problem (9) can b e reformulated as follows: u ∗ ← arg min u ∈ R n { − log Pr ( g | u ; K ) − log Pr ( u ) } . (10) A t first, we fo cus on the setting of the prior. A common choice is to mo del the unknown image u as a Mark ov Random Field (MRF) such that the image can be c haracterized b y its Gibbs prior distribution, whose general form is: Pr ( u ) = 1 Z n Y i =1 exp ( − α V c i ( u ) ) = 1 Z exp  − α n X i =1 V c i ( u )  , (11) where α > 0 is the MRF parameter, { c i } n i =1 is the set of all cliques (a clique is a set of neighboring pixels) for the MRF, V c i is the p otential function defined on the clique c i and Z is the partition function, that is a function not dep ending on u which allows for the normalization of the prior. Cho osing as potential function at the generic i -th pixel the magnitude of the discrete gradien t at the same pixel, i.e. V c i = k ( ∇ u ) i k 2 , the Gibbs prior in (11) reduces to the p opular TV prior: Pr ( u ) = 1 Z exp  − α n X i =1 k ( ∇ u ) i k 2  = 1 Z exp  − α TV ( u )  , where Z is the normalization constant not dep ending on u . The adoption of a TV prior can b e further in terpreted as assuming that the ` 2 norm of the gradien t at any pixel of the unknown clean image, k ( ∇ u ) i k 2 , follo ws a space-in v ariant half-Laplacian (or exp onential) distribution: Pr ( x ; α ) =    α exp ( − α x ) for x ≥ 0 0 for x < 0 . In [9], a deep in vestigation ab out the effect of replacing the half-Laplacian distribution with the more flexible half-Generalized Gaussian distribution Pr ( x ; p, α ) =    αp Γ(1 /p ) exp( − ( αx ) p ) for x ≥ 0 0 for x < 0 (12) 3 has b een carried out. The presence of a second parameter p allo ws for a b etter approximation of the ` 2 norm gradien t distribution and leads to the introduction of the TV p prior: Pr ( u ) = 1 Z exp  − α n X i =1 k ( ∇ u ) i k p 2  = 1 Z exp  − α TV p ( u )  . In this paper, we prop ose a prior, consisting in a non-stationary (space-v arian t) Marko v Random Field. The parameters α, p of the half-Generalized Gaussian distribution of the magnitude of the discrete gradien ts change as the clique c i c hanges. Therefore, the prior takes the follo wing form: Pr ( u ) = 1 Z exp  − n X i =1 α i k ( ∇ u ) i k p i 2  = 1 Z exp  − TV sv p,α ( u )  . (13) The adoption of a space-v arian t approach is exp ected to b e more flexible for the restoration of images presen ting cliques with differen t prop erties, i.e. images in which texture, smo oth, piecewise constant regions, and edges co-exist. In order to justify the reason why a space-v arian t mo del should b e adopted in general, we consider the test image skyscraper illustrated in Fig. 1(a). W e selected tw o regions c haracterized by smooth and texture structures - see the cyan-bordered and the yello w-bordered b oxes, respectively , in Figs. 1(d),1(g). The histogram of the gradient magnitudes of the whole image is shown in Fig 1(b) and zo omed in Fig. 1(c). The sup erimp osed green dashed lines, whic h hav e b een repro duced in eac h sub-figure for comparison, represent the half-Generalized Gaussian distributions that b est fit the histograms and whose parameters ha ve b een computed considering all the pixels of the test image. The histogram of the gradient magnitudes in the t wo b ordered regions are sho wn in Figs. 1(e), 1(h) and zo omed in Figs. 1(f ), 1(i). The red solid lines represent the half-Generalized Gaussian distributions that b est fit the histograms and whose parameters ha ve b een computed considering only the pixels in the b o xes. It is worth noticing how the histograms of the gradient magnitudes in the t wo selected regions are v ery different from each other and also differ from the one of the whole test image. As a result, the red lines fit the histogram shap es in Figs. 1(e), 1(h) b etter than the green ones - see also the zooms in Figs. 1(f ),1(i). This is the b enefit of the space-v ariant strategy , whic h is able to mo del space-v ariant image features. Going back to the MAP inference formula (10), in particular to the likelihoo d term Pr ( g | u ; K ) , assuming the noise to be additive indep endent iden tically distributed, w e hav e: Pr ( g | u ; K ) = n Y i =1 Pr ( g i | u ; K ) . (14) The likelihoo d term (14) clearly takes different form according to the distribution of the noise. In the follo wing we sp ecify the lik eliho o d asso ciated to the noises considered in this pap er. A dditive White Gaussian Noise . If the noise is known to b e A WG with standard deviation σ , the lik eliho o d term in (14) is as follo ws: Pr ( g | u ; K ) = n Y i =1 1 √ 2 π σ exp  − ( K u − g ) 2 i 2 σ 2  = 1 W exp  − k K u − g k 2 2 2 σ 2  , (15) where W is the normalization constan t not depending on u . Therefore, after replacing our prior (13) and the Gaussian lik eliho o d (15) in the MAP inference formula (10), and dropping the constan t terms, w e obtain our TV sv p,α -L 2 mo del in (7) with q = 2 , that is in extended form: u ∗ ← arg min u ∈ R n ( n X i =1 α i k ( ∇ u ) i k p i 2 + µ 2 k K u − g k 2 2 ) , where w e set µ = 1 /σ 2 . 4 (a) test image (b) global histogram (c) zoom of (b) (d) smooth region (e) local histogram for (d) (f ) zo om of (e) (g) texture region (h) local histogram for (g) (i) zoom of (h) Figure 1: Gradien t magnitudes histograms on the whole test image, on a smo oth region and on a texture region. A dditive White Laplace Noise . If the noise is kno wn to b e A WL with scale parameter β , the lik eliho o d term in (14) takes the following form: Pr ( g | u ; K ) = n Y i =1 1 2 β exp  − | K u − g | i β  = 1 W exp  − k K u − g k 1 β  , (16) Therefore, after replacing our prior (13) and the Laplace likelihoo d (16) in the MAP inference formula (10), and dropping the constan t terms, w e obtain our TV sv p,α -L 1 mo del in (7) with q = 1 , that is in extended form: u ∗ ← arg min u ∈ R n ( n X i =1 α i k ( ∇ u ) i k p i 2 + µ k K u − g k 1 ) , (17) where w e set µ = 1 /β . Salt and Pepper noise . The SPN can b e classified as a sp arse 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 p opular c hoice is to adopt the ` 0 pseudo-norm of the residual K u − g as the fidelity term. Nevertheless, it is v ery common to substitute the ` 0 pseudo-norm with the ` 1 norm, which is easier to deal with - since it is conv ex - and still allows a go o d sparsification effect. Hence, also in this case, the problem to whic h we are referring is (17). 5 3 Estimation of the space-v arian t parameters The prop osed regularization term (8) is derived by assuming that, for eac h pixel position i = 1 , ..., n , the magnitude - that is, the ` 2 norm - of the gradients of the target image in a surrounding neighbor- ho o d distributes according to a half-Generalized Gaussian (hGG) distribution, whose probabilit y densit y function is giv en in (12). This means that the distribution of the ` 2 norm of the gradien ts in the target image is defined pixel-wise as follows: Pr  k ( ∇ u ) i k ; p i , α i  = α i p i Γ(1 /p i ) exp  −  α i k ( ∇ u ) i k  p i  . In order to use the prop osed regularization term, we thus need to generate the p -map and the α -map. The metho d prop osed in [9] for estimating a global, image-based p v alue requires a v ery large n um b er of samples in order to provide statistically reliable estimates, therefore it could not b e generalized to our prop osal since we use small size image neighborho o ds for the estimation of lo cal p v alues. In [16] the authors prop osed a new metho d based on the statistical inference pro cedure illustrated in [8] which is sufficien tly robust to our purposes. F or completeness, in the follo wing w e briefly outline the method. Let u ∈ R n b e the v ectorized form of an image for whic h w e wan t to estimate the asso ciated vector of space-v arian t parameters p i , i ∈ Ω . First, we compute the v ector m ∈ R n con taining the magnitudes of the gradien ts of the image u ; in formulas: m i := k ( ∇ u ) i k 2 , i ∈ Ω . Then, we estimate each parameter p i b y applying the statistical inference technique in [8] to the lo cal data set consisting of the computed gradient magnitudes in a neighborho o d of the pixel i . In particular, w e use symmetric square neighborho o ds N s i of size s ∈ { 3 , 5 , . . . } centered at pixel i ∈ Ω . F ollo wing [8], the v alues p i , i ∈ Ω , shap e parameters of the hGG distributions, are estimated as follows: p i = h − 1 ( ρ i ) , ρ i = card  N s i   X j ∈ N s i m 2 j  /  X j ∈ N s i | m j |  2 , i ∈ Ω , (18) where card( A ) denotes the cardinality of set A and where the function h :]0 , + ∞ [ → ]0 , + ∞ [ , referred to as the gener alize d Gaussian r atio function in [8], is defined b y h ( z ) =  Γ(1 /z ) Γ(3 /z )  /  Γ 2 (2 /z )  , (19) with Γ( · ) indicating the Gamma function [1]. The function h in (19) is con tinuous, monotonically decreasing and surjective, hence inv ertible. Moreov er, since h is not data-dependent, its inv erse h − 1 , represen ting the v alues p i , can b e pre-computed off-line and stored as a lo okup-table, restricted to (0 , 2] , suc h that at run-time the final step of the estimation in (18) can be carried out very efficiently . The k ey no velt y of our prop osal relies on exploiting all the adv an tages of using a space-v arian t hGG distribution mo del, hence w e compute also the map of lo cal scale parameters α i , i = 1 , . . . , n . W e prop ose to estimate such scale parameters by means of a Maximum Lik eliho o d approach. Once p i for a pixel is estimated, the local likelihoo d function is giv en by L ( α, p i ; x 1 , ..., x n ) = n Y i =1  αp i Γ(1 /p i )  exp( − ( αx i ) p i ) =  αp i Γ(1 /p i )  n exp  − n X i =1 ( αx i ) p i  , (20) suc h that the v alue of the lo cal scale parameter is obtained by maximizing (20), that is by solving the follo wing optimization problem: α i = arg max α log L ( α, p i ; x 1 , ..., x n ) = arg max α  n log α − n X i =1 ( αx i ) p i  . (21) 6 By imp osing the first order optimality condition for problem (21), we obtain the following closed form estimation form ula: α i =  p i n n X i =1 x p i i  − 1 p i . (22) In Fig. 2 the maps of lo cal p v alues, obtained with neigh b orho o ds 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 an y case the metho d asso ciates very low p v alues with flat regions and higher v alues with edges. It is w orth remarking that in Sect. 5 numerical experiments ha ve been carried out by computing the p -map starting from the corrupted images. (a) (b) (c) Figure 2: Original test image skyscraper (a), p -map for s = 3 (b) and s = 10 (c). In Fig. 3 the maps of lo cal α v alues, obtained with neighborho o ds 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. (a) (b) (c) Figure 3: Original test image skyscraper (a), α -map for s = 3 (b) and s = 10 (c). 4 Numerical solution b y ADMM In this section, we illustrate the ADMM-based iterativ e algorithm used to numerically solve the prop osed mo del (7)–(8) for b oth cases q = 2 and q = 1 . T o this purp ose, first we resort to the v ariable splitting tec hnique [2] and introduce tw o auxiliary v ariables r ∈ R n and t ∈ R 2 n , suc h that mo del (7)–(8) is rewritten in the follo wing equiv alen t constrained form: { u ∗ , r ∗ , t ∗ } ← arg min u,r,t  n X i =1 α i k t i k p i 2 + µ q k r k q q  , q ∈ { 1 , 2 } , (23) sub ject to : r = K u − g , t = D u , (24) where D := ( D T h , D T v ) T ∈ R 2 n × n and t i :=  ( D h u ) i , ( D v u ) i  T ∈ R 2 represen ts the discrete gradient of image u at pixel i . T o solv e problem (23)–(24) by ADMM [4], we define the augmented Lagrangian functional L ( u, r , t ; λ r , λ t ) = n X i =1 α i k t i k p i 2 + µ q k r k q q − h λ t , t − D u i + β t 2 k t − D u k 2 2 − h λ r , r − ( K u − g ) i + β r 2 k r − ( K u − g ) k 2 2 , (25) 7 where β r , β t > 0 are scalar p enalt y parameters and λ r ∈ R n , λ t ∈ R 2 n are the vectors of Lagrange m ultipliers asso ciated with the linear constrain ts r = K u − g and t = D u in (24), resp ectively . Solving (23)–(24) is th us equiv alen t to seek for the solutions of the following saddle p oint problem: Find ( x ∗ ; λ ∗ ) ∈ X × Λ suc h that L ( x ∗ ; λ ) ≤ L ( x ∗ ; λ ∗ ) ≤ L ( x ; λ ∗ ) ∀ ( x ; λ ) ∈ X × Λ , (26) with the augmen ted lagrangian functional L defined in (25) and where, for simplicity of notations, w e set x := ( u, r, t ) , λ := ( λ r , λ t ) , X := R n × R n × R 2 n and Λ := R n × R 2 n . Giv en the previously computed (or initialized for k = 0 ) vectors u ( k ) , λ ( k ) r and λ ( k ) t , the k -th iteration of the prop osed ADMM-based iterative scheme applied to the solution of the saddle-p oint problem (26) - minimization for the primal v ariables u, r, t , maximization for the dual v ariables λ r , λ t - reads as follows: r ( k +1) ← arg min r ∈ R n L ( u ( k ) , r , t ( k ) ; λ ( k ) r , λ ( k ) t ) , (27) t ( k +1) ← arg min t ∈ R 2 n L ( u ( k ) , r ( k +1) , t ; λ ( k ) r , λ ( k ) t ) , (28) u ( k +1) ← arg min u ∈ R n L ( u, r ( k +1) , t ( k +1) ; λ ( k ) r , λ ( k ) t ) , (29) λ ( k +1) r ← λ ( k ) r − β r  r ( k +1) − ( K u ( k +1) − g )  , (30) λ ( k +1) t ← λ ( k ) t − β t  t ( k +1) − D u ( k +1)  . (31) In the following three subsections we describe how to solve the minimization sub-problems (27), (28) and (29) for the primal v ariables r , t and u , resp ectively , in b oth cases q ∈ { 1 , 2 } . In particular, we remark that thanks to the preliminary ADMM v ariable splitting pro cedure, sub-problems (28) and (29) for the v ariables t and u are identical in the tw o cases q ∈ { 1 , 2 } and their solution can b e obtained based on form ulas given in [9] for the same sub-problems. 4.1 Minimization sub-problem for the primal v ariable r Recalling the definition of the augmented Lagrangian functional in (25) and carrying out some simple algebraic manipulations, the minimization sub-problem (27) for the primal v ariable r can b e written as r ( k +1) ← arg min r ∈ R n  µ q k r k q q + β r 2   r − v ( k )   2 2  , q ∈ { 1 , 2 } , (32) with the constan t (with respect to the optimization v ariable r ) vector v ( k ) ∈ R n giv en by v ( k ) = K u ( k ) − g + 1 β r λ ( k ) r . (33) Since µ ≥ 0 , β r > 0 , in both cases q ∈ { 1 , 2 } the cost function in (32) is strongly conv ex, hence it admits a unique global minimizer. In particular, the unique solution r ( k +1) of (32) can b e computed, dep ending on q , b y means of the follo wing closed-form formulas: case q = 1 : r ( k +1) = sign  v ( k )   max  | v ( k ) | − µ/β r , 0  , (34) case q = 2 : r ( k +1) =  β r / ( β r + µ )  v ( k ) , (35) where sign ( · ) and | · | in (34) denote the comp onent-wise signum and absolute v alue functions and  indicates the component-wise vectors pro duct. W e remark that form ula (34) represen ts a well-kno wn comp onen t-wise soft-thresholding op erator - see e.g. [12] - whereas (35) comes easily from first-order optimalit y conditions of (32). In case that the regularization parameter µ is regarded as a constant - that is, it is fixed a priori - then form ulas (34)–(35) allow to determine very efficien tly the solution r ( k +1) of this sub-problem. How ever, as previously stated, in the case q = 2 w e aim also at automatically adjusting µ along iterations - that is, µ b ecomes µ ( k ) - such that the final solution u ∗ of our mo del (7)–(8) satisfies the discrepancy principle [14]. T o this aim, in the following w e prop ose a pro cedure which builds up on those presented in [6, 10] but, due to a differen t ADMM initial v ariable splitting, needs to be adapted and is worth to b e outlined in detail. 8 W e consider the discrepancy asso ciated with the solution r ( k +1) in (35) as a function δ ( k +1) : [0 , + ∞ [ → [0 , + ∞ [ of the regularization parameter µ : δ ( k +1) ( µ ) :=   r ( k +1) ( µ )   2 = β r β r + µ   v ( k )   2 , (36) where the second equality comes from (35). The discrepancy function in (36) is clearly contin uous, non- negativ e and monotonically decreasing ov er its entire domain µ ∈ [0 , + ∞ [ and at the extremes w e hav e δ ( k +1) ( µ = 0) = k v ( k ) k 2 , δ ( k +1) ( µ → + ∞ ) = 0 . In order to set a v alue µ ( k +1) suc h that the discrepancy principle is satisfied here for the auxiliary v ariable r (recall that r = K u − g represents the residue of the restoration), w e consider t wo complemen tary cases based on the v alue of the norm of the vector v ( k ) defined in (33). In case that k v ( k ) k 2 ≤ ¯ δ , with ¯ δ denoting the noise lev el, then from (36) and from the fact that 0 < β r / ( β r + µ ) ≤ 1 , it follows that δ ( k +1) ( µ ) ≤ ¯ δ ∀ µ ∈ [0 , + ∞ [ , that is the discrepancy principle is satisfied for an y µ . In this case we th us set µ ( k +1) = 0 , such that, according to (35), the sub-problem solution is r ( k +1) = v ( k ) . In case that k v ( k ) k 2 > ¯ δ , the prop erties of the discrepancy function δ ( k +1) ( µ ) in (36) guaran tee that there exists a unique v alue µ ( k +1) of µ suc h that δ ( k +1) ( µ ( k +1) ) = ¯ δ . Recalling (36), we hav e  β r / ( β r + µ ( k +1) )  k v ( k ) k 2 = ¯ δ ⇐ ⇒ µ ( k +1) = β r  k v ( k ) k 2 / ¯ δ − 1  . Replacin g this expression for µ in (35), the sub-problem solution is r ( k +1) = ¯ δ v ( k ) / k v ( k ) k 2 . T o summarize, the solution of this sub-problem at any iteration k is computed by (34) for the case q = 1 whereas for the case q = 2 it is determined as follo ws: k v ( k ) k 2 ≤ ¯ δ = ⇒ µ ( k +1) = 0 , r ( k +1) = v ( k ) k v ( k ) k 2 > ¯ δ = ⇒ µ ( k +1) = β r  k v ( k ) k 2 / ¯ δ − 1  , r ( k +1) = ¯ δ v ( k ) / k v ( k ) k 2 (37) 4.2 Minimization sub-problem for the primal v ariable t Giv en the definition of the augmen ted Lagrangian functional in (25), the minimization sub-problem for the primal v ariable t in (28) can b e written as follo ws: t ( k +1) ← arg min t ∈ R 2 n ( n X i =1 α i k t i k p i 2 − h λ ( k ) t , t − D u ( k ) i + β t 2    t − D u ( k )    2 2 ) ← arg min t ∈ R 2 n ( n X i =1 α i k t i k p i 2 + β t 2     t −  D u ( k ) + 1 β t λ ( k ) t      2 2 ) ← arg min t ∈ R 2 n n X i =1 ( α i k t i k p i 2 + β t 2     t i −   D u ( k )  i + 1 β t  λ ( k ) t  i      2 2 ) . (38) Note that in (38) the minimized functional is written in explicit comp onent-wise (or pixel-wise) form, with  D u ( k )  i ,  λ ( k ) t  i ∈ R 2 denoting the discrete gradien t and the Lagrange m ultipliers at pixel i , resp ectiv ely . Solving the 2 n -dimensional minimization problem in (38) is thus equiv alen t to solve the n follo wing indep endent 2 -dimensional problems: t ( k +1) i ← arg min t i ∈ R 2  k t i k p i 2 + ( β t /α i ) 2    t i − q ( k ) i    2 2  , i = 1 , . . . , n , (39) with the constan t v ectors q ( k ) i ∈ R 2 defined b y q ( k ) i :=  D u ( k )  i + 1 β t  λ ( k ) t  i , i = 1 , . . . , n . (40) The solutions of the n optimization problems in (39) can b e obtained based on the results rep orted in Proposition 1 of [9], that is: t ( k +1) i = ξ ( k ) i q ( k ) i , i = 1 , . . . , n , (41) where, in particular, the shrinkage co efficients ξ ( k ) i ∈ [0 , 1] , i = 1 , . . . , n, are given b y formulas (50)–(52) in [9]. The o verall computational cost of this subproblem is linear in the num b er of pixels n . 9 4.3 Minimization sub-problem for the primal v ariable u As illustrated in [9], the minimization sub-problem (29) for the primal v ariable u reduces to the solution of the follo wing n × n system of linear equations  D T D + β r β t K T K  u = D T  t ( k +1) − 1 β t λ ( k ) t  + β r β t K T  r ( k +1) − 1 β r λ ( k ) r + g  , (42) whic h is solv able if the co efficien t matrix has full-rank, that is if the following condition holds: Ker  D T D  ∩ Ker  K T K  = { 0 } , (43) where Ker( M ) denotes the n ull space of matrix M and 0 is the n -dimensional n ull v ector. In our case, condition (43) is satisfied. In fact, K represents a blurring op erator, which is a lo w-pass filter, whereas the regularization matrix D is a first-order difference op erator and, hence, is a high-pass filter. Moreo ver, since β t , β r > 0 , the co efficient matrix in (42) is symmetric p ositive definite and typically highly sparse. Hence, the linear system in (42) can b e solved quite efficiently by the iterative (even tually preconditioned) conjugate gradient metho d. Moreov er, under appropriate assumptions ab out the solution u near the image b oundary , the linear system can b e solved even more efficiently . W e assume p erio dic b oundary conditions for u , so that b oth D T D and K T K are blo ck circulan t matrices with circulant blo c ks and, hence, the co efficient matrix in (42) can b e diagonalized by the 2D discrete F ourier transform (FFT implementation). Pro vided that the penalty parameters β t , β r are kept fixed during the ADMM iterations, the co efficient matrix in (42) do es not change and it can b e diagonalized once for all at the b eginning. Therefore, at any ADMM iteration the linear system (42) can b e solved by one forward 2D FFT and one in v erse 2D FFT, each at a cost of O ( n log n ) . 4.4 ADMM-based iterativ e sc heme T o summarize previous results, in Algorithm 1 we rep ort the main steps of the prop osed ADMM-based iterativ e scheme used to solve the saddle-p oint problem (25)–(26) and, hence, to compute solutions of the proposed mo del (7)–(8). In the field of image and sign al pro cessing the ADMM has b een one of the most pow erful and successful metho ds for solving v arious conv ex or nonconv ex optimization problems. In con vex settings, n umerous con vergence results hav e b een established for ADMM as well as its v arieties, see for example [18] and references therein . In particular, con vergence results co ver the proposed TV p,α -L q mo dels, q ∈ { 1 , 2 } , in the sp ecial case of p i ≥ 1 ∀ i . Ho wev er, in case that one or more p i < 1 , the ADMM is under nonconv ex settings, where there hav e b een a few studies on the conv ergence prop erties. T o the b est of our knowledge, existing con v ergence results of ADMM for nonconv ex 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 v arious applications inv olving nonconv ex optimization problems, and this is a practical justification of its wide use. 5 Numerical results In this section, we ev aluate exp erimentally the p erformance of the t wo prop osed mo dels TV sv p,α -L q , q ∈ { 1 , 2 } , defined in (7)–(8), when applied to the restoration of gray level images synthetically corrupted b y kno wn blur and b y A W GN - in the case of TV sv p,α -L 2 mo del - and SPN or A WLN - in the case of TV sv p,α -L 1 mo del. In particular, the t wo prop osed models are compared with: • TV-L q , q ∈ { 1 , 2 } , defined in (6); see [13], [15]; • TV p -L q , q ∈ { 1 , 2 } , with global p ∈ (0 , 2] ; see [9], [16]; • TV sv p -L q , q ∈ { 1 , 2 } , with lo cal p i ∈ (0 , 2] , i ∈ { 1 , . . . , n } ; see [16]. F or what concerns the preliminary estimation of the p i and α i parameters, we directly apply the pro cedure outlined in Section 3 to the observ ed corrupted image g for the A W GN and A WLN cases. Instead, for the SPN case, in order to ha ve a robust estimation of the parameters, a preliminary pro cessing b y an adaptive filter is required. In particular, we assume that the p osition of the pixels corrupted by the SPN is known a priori, otherwise it can b e easily detected as suggested in [5]. W e replace the corrupted pixels with the mean of the non-corrupted pixels of its neigh b orho o d. The size of the neighborho o d is 10 Algorithm 1 ADMM-based scheme for mo dels (7)–(8) input : observ ed image g ∈ R n output : appro ximate solution u ∗ ∈ R n of (7)–(8) 1. initialize: 2. · estimate parameters p i and α i , i = 1 , . . . , n , b y (18) and (22), resp ectiv ely 3. · set u (0) = g , λ (0) r = λ (0) t = 0 2. for k = 1, 2, 3, . . . until c onver genc e do : 3. · up date primal v ariables: 4. · compute r ( k +1) b y (33) and (34) for q = 1 , (37) for q = 2 5. · compute t ( k +1) b y (40), (41) and form ulas (50)–(52) in [9] 5. · compute u ( k +1) b y solving (42) 6. · up date dual v ariables: 7. · compute λ ( k +1) r , λ ( k +1) t b y (30), (31) 8. end for 9. u ∗ = u ( k +1) v ariable and dep ends on the p ercentage P of non-corrupted pixels in it. If P is b elow a fixed threshold P (usually P = 0 . 4 ), then the neigh b orho o d is enlarged, in order to incorp orate a greater num b er of uncorrupted pixels. The obtained image is then used to compute the p -map and the α -map. The describ ed strategy has b een introduced instead of the classic median filter, whose smo othing effects is quite high. Clearly , the same approac h is adopted for the TV sv p -L 1 mo del to estimate the p -map only . The quality of the observed corrupted images g and of the restored images u ∗ is measured - in dB - b y means of the Blurred Signal-to-Noise Ratio BSNR( g , u ) = 10 log 10 k K u − E [ K u ] k 2 2 k g − K u k 2 2 and the Impro ved Signal-to-Noise Ratio ISNR( g , u, u ∗ ) = 10 log 10 k g − u k 2 2 k u ∗ − u k 2 2 , resp ectiv ely , with u denoting the original uncorrupted image and E [ K u ] the a verage intensit y of the blurred image K u . In general, the larger the ISNR v alue, the higher the qualit y of restoration. F or all the ADMM-based minimization algorithms and for all the tests, the p enalty parameters β t and β r are suitably set. Moreo ver, for all tests, the ADMM iterations of all the compared algorithms are stopp ed as soon as tw o successiv e iterates satisfy   u ( k ) − u ( k − 1)   2   u ( k − 1)   2 < 10 − 4 . F or the mo dels with the L 2 fidelit y term, the regularization parameter µ has b een automatically set based on the discrepancy principle. F or the mo dels with the L 1 fidelit y term, µ has been hand-tuned indep enden tly in eac h test so as to provide the highest p ossible ISNR v alue. In the follo wing, w e rep ort n umerical results concerning the restoration of blurred images corrupted by A WGN (Example 1) and by impulsiv e SPN and A WLN (Example 2). Example 1: restoration of images corrupted by blur and A W GN. In this example, we ev aluate experimentally the p erformance of the prop osed TV sv p,α -L 2 mo del on a purely piecewise constant test image - geometric ( 256 × 256 ), Fig. 4(a) - and a partially textured test image - skyscraper ( 256 × 256 ), Fig. 4(d). Both images ha v e been synthetically corrupted b y a Gaussian blur of parameters band=5 and sigma=1.0 and b y A WGN characterized b y different noise levels. The p, α -maps hav e b een computed b y using neigh b orho o ds of size s = 3 . 11 T able 1: Example 1: achiev ed ISNR v alues. geometric skyscrap er BSNR TV-L 2 TV p -L 2 TV sv p -L 2 TV sv p,α -L 2 TV-L 2 TV p -L 2 TV sv p -L 2 TV sv p,α -L 2 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 In T able 1 the performance of our mo del are compared in terms of achiev ed ISNR v alues with those of the TV-L 2 , TV p -L 2 and TV sv p -L 2 mo dels. The go o d quality of the restored image b y our mo del can b e appreciated b y a visual insp ection of Figs. 4(c),(f ) and by comparing the ISNR v alues rep orted in T able 1. (a) original (b) corrupted (BSNR = 20) (c) TV sv p,α -L 2 (ISNR = 8.60) (d) original (e) corrupted (BSNR = 20) (f ) TV sv p,α -L 2 (ISNR = 3.31) Figure 4: Example 1: restoration of the test images geometric and skyscraper corrupted by blur and A WGN. Example 2: restoration of images corrupted by blur and SPN or A WLN. In this example w e ev aluate the performance of the prop osed TV sv p,α -L 1 mo del on three medical test images lungs ( 468 × 591 ), Fig. 5 (a), ecography ( 401 × 511 ), Fig. 6 (a), and aneurism ( 701 × 766 ), Fig. 7 (a), synthetically corrupted b y Gaussian blur of parameters band=5 and sigma=1 and b y tw o t yp es of impulsive noise, namely SPN and A WLN. The images are provided in the rep ository at https://medpix.nlm.nih.go v. First, for what concerns corruptions by SPN, in Figs. 5, 6, 7 we rep ort for the three considered test images the original and corrupted image together with the estimated p, α -maps in the first column (with the size s of the neighborho o ds used for the p, α -maps estimation rep orted in the captions), the restoration results, obtained b y the four compared metho ds, in the second column (with the ac hiev ed ISNR v alues in the captions) and a zo omed detail of the restored images - green- b ordered in Figs. 5 (a), 6 (a), 7 (a) - in the last column. The rep orted ISNR v alues as w ell as the visual insp ection of the restored images and of the zoomed details strongly indicate how the prop osed space-v arian t regularizer allows for higher quality restorations. In particular, it is w orth remarking how, with resp ect to the space-v arian t TV sv p mo del, the additional de- grees of freedom represen ted by the scale parameters α i used in our prop osal, yield a sufficient additional 12 T able 2: Example 2 (A WLN): ac hieved ISNR v alues. TV-L 1 TV p -L 1 TV sv p -L 1 TV sv p,α -L 1 lungs 6.20 6.80 7.30 7.85 ecograph y 5.93 6.40 7.88 8.32 aneurism 9.10 9.44 10.13 10.70 flexibilit y for av oiding un wan ted spurious effects - see, e.g., spikes in Figs. 5 (i), 6 (i), 7 (i). In the second part of this example, w e consider the restoration of the same three medical test images corrupted b y the same blur of parameters band=5 , sigma=1 and b y a different impulsiv e noise, namely A WLN of lev el yielding BSNR=10. In T able 2 w e report the ISNR v alues ac hieved by the compared metho ds and in Fig. 8 we show the original images, the corrupted images and the restored images by our mo del. The results in T able 2 confirm that, also in case of images corrupted b y A WLN, the prop osed TV sv p,α -L 1 mo del outp erforms its comp etitors in terms of ISNR. Moreo ver, the restored images depicted in the last column of Fig. 8 provide further evidence of the go o d quality restorations achiev able by our prop osal. 6 Conclusions W e presen ted a new space-v arian t regularizer for v ariational image restoration based on the assumption that the gradient magnitudes of the target image distribute lo cally according to a half-Generalized Gaus- sian distribution. Thanks to the high num b er of free parameters in volv ed in the regularizer and to the fact that such parameters can b e automatically and robustly estimated from the observed image, our prop osal exhibits a very high flexibilit y which p otentially allows for an effective mo deling of space-v arian t prop erties of images. By coupling the prop osed regularizer with either the L 1 or L 2 fidelit y terms, we tested our prop osal on images corrupted by blur and different types of noise, namely A WGN, A WLN and SPN. The restored images, obtained by means of an efficient ADMM-based n umerical algorithm, strongly indicate that the proposed regularizer holds the p otential for achieving high qualit y restoration results for a wide range of target images c haracterized b y different gradien t distributions and for the differen t types of noise considered. A ckno wledgmen ts : Researc h w as supported by the “National Group for Scien tific Computation (GNCS-IND AM)” and by ex60 pro ject by the Universit y of Bologna “F unds for selected research topics” . References [1] Abramo witz M, and Stegun IA. 1970. Handbo ok of Mathematical F unctions. New Y ork. [2] Bioucas-Dias J and Figueredo M. 2010. F ast image Recov ery Using V ariable Splitting and Con- strained Optimization. IEEE T rans Image Pro c 19: 2345–2356. [3] Bo vik A.C. 2010. Handb o ok of Image and Video Pro cessing. A cademic press. [4] Bo yd S, Parikh N, Ch u E, Peleato B and Ec kstein J. 2011. Distributed Optimization and Statistical Learning via the Alternating Direction Metho d of Multipliers F oundations and T rends in Machine Learning 3: 1–122. [5] Cai JF, Chan RH and Nikolo v a M. 2010. F ast t wo-phase image deblurring under impulse noise. Journ Math Imaging Vision 36: 46–53. [6] He C, Hu C, Zhang W and Shi B. 2014. A F ast Adaptiv e P arameter Estimation for T otal V ariation Image Restoration. IEEE T rans Image Pro c 23: 4954–4967. [7] Kai-Sheng S(2006) A globally conv ergen t and consistent metho d for estimating the shap e parameter of a generalized Gaussian distribution. IEEE T rans Infor 52: 510–527. 13 [8] Karnran S and Leon-Garcia A. 1995. Estimation of shap e parameter for generalized Gaussian distri- butions in subband decomp ositions of video. IEEE T rans Circuits and Systems for Video T ec hnology 5: 52–56. [9] Lanza A, Morigi S and Sgallari F. 2016. Constrained T V p - ` 2 Mo del for Image Restoration, Journ. Scien tific Computing, 68: 64–91. [10] Lanza A, Morigi S and Sgallari F. 2016. Conv ex Image Denoising via Non-con v ex Regularization with P arameter Selection. Journ Math Imag Vision 56: 195–220. [11] Lazzaro D, Morigi S, Melpignano P , Loli Piccolomini E and Benini L. 2017. Image enhancemen t v ariational metho ds for Enabling Strong Cost Reduction in OLED-based Poin t-of-Care Immunofluo- rescen t Diagnostic Systems. In ternational Journal for Numerical Metho ds in Biomedical Engineering 34(3), e2932. [12] Min T, Y ang J and He B. 2009. Alternating direction algorithms for total v ariation deconv olution in image reconstruction. TR0918, Dept Math, Nanjing Universit y . [13] R udin LI, Osher S and F atemi E. 1992: Nonlinear total v ariation based noise remo v al algorithms. Ph ysics D, 60: 259–268 [14] W en Y and Chan RH. 2012. P arameter Selection for T otal V ariation Based Image Restoration Using Discrepancy Principle. IEEE T rans Image Pro c. 21: 1770–1781. [15] T ao M, Y ang J. 2009. Alternating Direction Algorithm for T otal V ariation Deconv olution in Image Reconstruction, Departmen t of Mathematics. Nanjing Univ ersity , T ech. Rep. TR0918. [16] Lanza A, Morigi S, Pragliola M, Sgallari F. 2018. Space-v arian t TV regularization for image restora- tion. Lecture Notes in Computational Vision and Biomechanics. [17] Hong M, Luo Z, Raza viya yn M. 2014. Con vergence Analysis of Alternating Direction Metho d of Multipliers for a F amily of Nonconv ex Problems. Preprin t, [18] He B, Y uan X. 2012. On the O(1/n) conv ergence rate of the Douglas-Rachford alternating direction metho d. SIAM Journal on Numerical Analysis, vol. 50, no. 2, pp. 700–709. 14 (a) original (b) TV-L 1 (ISNR = 11.04) (c) zoom of (b) (d) corrupted (e) TV p -L 1 (ISNR = 12.48) (f ) zoom of (e) (g) p -map ( s = 3 ) (h) TV sv p -L 1 (ISNR = 15.30) (i) zoom of (h) (l) α -map ( s = 3 ) (m) TV sv p,α -L 1 (ISNR = 16.56) (n) zoom of (m) Figure 5: Example 2 (SPN): visual restoration results for the test image lungs corrupted b y a γ = 0 . 1 lev el noise. 15 (a) original (b) TV-L 1 (ISNR = 22.13) (c) zo om of (b) (d) corrupted (e) TV p -L 1 (ISNR = 23.15) (f ) zoom of (e) (g) p -map ( s = 10 ) (h) TV sv p -L 1 (ISNR = 25.46) (i) zo om of (h) (l) α -map ( s = 10 ) (m) TV sv p,α -L 1 (ISNR = 28.01) (n) zo om of (m) Figure 6: Example 2 (SPN): visual restoration results for the test image ecography corrupted b y a γ = 0 . 35 lev el noise. 16 (a) original (b) TV-L 1 (ISNR = 18.55) (c) zoom of (b) (d) corrupted (e) TV p -L 1 (ISNR = 19.10) (f) zo om of (e) (g) p -map ( s = 3 ) (h) TV sv p -L 1 (ISNR = 21.14) (i) zoom of (h) (l) α -map ( s = 3 ) (m) TV sv p,α -L 1 (ISNR = 24.47) (n) zoom of (m) Figure 7: Example 2 (SPN): visual restoration results for the test image aneurism corrupted by a γ = 0 . 1 lev el noise. 17 (a) original (b) original (c) original (d) corrupted (BSNR=10) (e) corrupted (BSNR=10) (f ) corrupted (BSNR=10) (g) TV sv p,α -L 1 (h) TV sv p,α -L 1 (i) TV sv p,α -L 1 Figure 8: Example 2 (A WLN): visual restoration results. 18