IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

1831

Speckle Reduction via Higher Order Total Variation Approach Wensen Feng, Hong Lei, and Yang Gao

Abstract— Multiplicative noise (also known as speckle) reduction is a prerequisite for many image-processing tasks in coherent imaging systems, such as the synthetic aperture radar. One approach extensively used in this area is based on total variation (TV) regularization, which can recover significantly sharp edges of an image, but suffers from the staircase-like artifacts. In order to overcome the undesirable deficiency, we propose two novel models for removing multiplicative noise based on total generalized variation (TGV) penalty. The TGV regularization has been mathematically proven to be able to eliminate the staircasing artifacts by being aware of higher order smoothness. Furthermore, an efficient algorithm is developed for solving the TGV-based optimization problems. Numerical experiments demonstrate that our proposed methods achieve state-of-theart results, both visually and quantitatively. In particular, when the image has some higher order smoothness, our methods outperform the TV-based algorithms. Index Terms— Speckle, total generalized variation, total variation, inexact split Uzawa.

I. I NTRODUCTION

S

PECKLE noise in many coherent imaging systems (e.g., SAR [1], ultrasound and laser imaging) keeps us from interpreting valuable information such as texture, edges, and point targets. Therefore, the reduction of speckle is critical for many image processing tasks involving image segmentation, target detection and classification. Without loss of generality, we focus on SAR images with single frequency and single polarization, where the speckle is assumed to be fully developed, implying that the number of point scatterers in a resolution cell is large [1]. Up to now, the major despeckling techniques can be loosely grouped in four categories: 1) spatial; 2) wavelet-based; 3) nonlocal filtering; and 4) variational. A typical spatial approach is the multi-look technique which amounts to incoherently averaging independent observations of the same resolution cell, thus reducing the noise intensity. However, for the averaging approaches, the reduction of noise is always at the cost of a clear loss in image resolution. To overcome this deficiency, a great deal of research has been conducted to develop suitable filters which can reduce the noise, yet preserve details and edges [35]. Early filters of this

Manuscript received February 26, 2013; revised October 22, 2013 and January 29, 2014; accepted February 6, 2014. Date of publication February 26, 2014; date of current version March 12, 2014. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. A. N. Rajagopalan. The authors are with the Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China (e-mail: [email protected]; [email protected]; [email protected]). Digital Object Identifier 10.1109/TIP.2014.2308432

kind [39] were developed under the minimum-mean-squareerror (MMSE) criterion. Afterwards, the more sophisticated and promising maximum a posteriori (MAP) approach was proposed [30]. Nonetheless, these filters may cause the loss of details if the statistical description of the SAR image is inaccurate. The despeckling techniques in the wavelet domain [2], [8], [9] resort to statistical wavelet shrinkage with MAP Bayesian approach, and focus on the selection of the most suitable prior for the signal coefficients [2]–[7]. In general, the wavelet-based methods guarantee a superior ability to preserve signal resolution. However, they often suffer from the ringing effects near the edges of the images, or isolated patterns in flat areas, resulting in visually annoying. Recently, the modern denoising methods, e.g., nonlocal mean (NLM) [10], block-matching 3D (BM3D) [11] and K-SVD [12], have been readily extended to SAR despeckling [13]–[16]. As an extension of the NLM algorithm, the probabilistic patch-based (PPB) filter [13] provides promising results by developing an efficient similarity measure well suited to SAR images in a data-driven way. A drawback of this filter is the suppression of thin and dark details in the regularized images. By taking into account the peculiar features of SAR images, Parrilli et al. [15] derived a SAR-oriented version of BM3D. It exhibits an objective performance comparable or superior to other techniques on simulated speckled images, and guarantees a very good subjective quality on real SAR images. However, the nonlocal filtering methods also present their typical artifacts in the form of structured signal-like patches in flat areas, originated from the randomness of speckle and reinforced through the patch selection process. In addition, the huge computational burden is another disadvantage encountered with these methods. Meanwhile, the variational approaches have attracted much interest in removing the speckle by directly estimating the reflectance of the underlying scene. There are two terms in the variational model: a data fidelity term and a regularizer. Due to the edge sharpness preserving and noise smoothing properties [17], TV has been widely used as a regularizer in the speckle removal task [18]–[24], [26], [27]. A. Variational Models for Speckle Reduction and Corresponding Algorithms In this subsection, the major TV-based variational models and the corresponding algorithms are briefly reviewed.

1057-7149 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

1832

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Let 0 ⊂ R2 and Y be the observed SAR image intensity. It is known that Y can be assumed to be the product of the underlying true image intensity u : 0 → R+ and the multiplicative noise η: Y = uη.

(1)

The probability density function (pdf) of η for the L-look SAR image is given by the following Gamma distribution: p(η) =

where H is Heaviside function and  is the classical Gamma function. Using the MAP criterion, Aubert and Aujol [18] derived a new model (AA-model), and aimed to solve the energy minimization problem as below:    ∗ u = arg min λ logu + Yu −1 d x + (u), (2) 0

u∗

where denotes the despeckled result, λ > 0 denotes the regularization parameter and  (u) is related to the prior probability distribution of u. The TV regularizer TV(u) =  |∇u| d x is adopted as φ(u) in the AA-model, leading to 0 the final problem:     |∇u| d x. (3) u ∗ = arg min λ logu + Yu −1 d x + u

0

0

As to (3), the global solution is hard to find because of its nonconvexity. To resolve this problem, logarithmic transformation (i.e., u l = log(u)) is commonly used [19], [22]. Then, the transformed variational model becomes     |∇u l | d x, u l + Ye−ul d x + u l∗ = arg min λ ul

0

0

u l∗

u∗ = e .

(4)

However, since the logarithmic transformation is nonlinear, the values in dark area of the image are expanded while those in the bright part are compressed. To address this deficiency, Steidl and Teuber [24] showed that the classical I-divergence model given in (5) does not require nonlinear transformation, and shares the same solution of model (4) theoretically.   ∗ |∇u| d x. (5) u = arg min λ (u − Ylogu)d x + u

0

0

In what follows, the discrete formulations of problems (4) and (5), and the solving algorithms are reviewed. For simplicity, we assume that 1 corresponds to the discretized image domain 1 = {(i, j) |i, j ∈ Z, 1 ≤ i ≤ M, 1 ≤ j ≤ N }, where Z denotes the set of positive integers, and (i, j ) is the pixel location. Then, the standard isotropic discrete TV regularizer can be represented as M N 2 ∇u1 = ((∇u)h(i, j ) ) + ((∇u)v(i, j ) )2 , i=1

(∇u)v(i, j )

ul

u l∗

u∗ = e ,

j =1

and (∇u)h(i, j ) denote the vertical and horiwhere zontal first order difference at pixel (i, j ) ∈ 1 , respectively.

(6)

where , denotes the inner product, and u ∗ = arg min λ u − Y log u, 1 + ∇u1 . u>0

1 LL ηL−1 e−Lη H (η) , (L)

u

Therefore, the exponential model (4) and the I-divergence model (5) are expressed as the following discrete formulations:

 u l∗ = arg min λ u l + Ye−ul , 1 + ∇u l  1 ,

(7)

To cope with the exponential model of (6), Bioucas-Dias and Figueiredo [19] employed the alternating direction method of multipliers (ADMM) [25]. This algorithm is named as the multiplicative image denoising by augmented Lagrangian (MIDAL) which requires inner iterations. Huang et al. [22] used a variable splitting to address the problem of (6). Nevertheless, if the penalty parameter is unsuitable, their model becomes severely ill-conditioned. For computing the minimizer of the I-divergence model (7), Steidl and Teuber [24] applied the technique named as alternating direction method of multipliers (ADMM). It needs evaluating the inverse operation, which consumes rather expensive computational cost. Recently, Hyenkyun Woo and Sangwoon Yun [26] proposed a general framework combining an alternating minimization algorithm (AMA) [28] with a shifting technique (AMAST). Their proposed method has a simple closed-form solution at each iteration, thus avoiding any inner loops or inverse operation. B. Proposed Approach Indeed, solutions of variational problems with TV regularization admit many desirable properties, most notably the appearance of sharp edges. However, the regularization with TV also has some shortcomings, most notably the so-called staircasing artifact, i.e., the undesired appearance of edges. This is originated from the model assumption that an image is consist of piecewise-constant areas, even in regions with smooth transitions of pixel values. To overcome the staircasing artifact, several modified models have been suggested to combine first- and higherorder TV functionals, in particular the infimal convolution (ICTV) [29] and total generalized variation (TGV) [31]. As detailed in [32], near the part where the jump of the image is not constant across the boundary, plugging the functional ICTV into variational problems gives solutions which still admit staircasing artifacts. However, these undesirable effects are eliminated by the TGV-based methods. Therefore, the TGV regularizer has been widely used in the variational models for reducing the additive noise, deblurring, and some other applications [36]. In this study, we incorporate the TGV penalty into the existing data fidelity term for speckle removal, and develop two novel variation despeckling models. As demonstrated in our numerical experiments, the TGV-based despeckling method outperforms the traditional TV methods by reducing the staircasing artifact, since TGV involves and balances higher-order derivatives of the image. Meanwhile, to solve the

FENG et al.: SPECKLE REDUCTION VIA HIGHER ORDER TOTAL VARIATION APPROACH

proposed models efficiently, we propose an algorithm based on the primal-dual algorithm [38]. The remainder of the paper is organized as follows. In Section II, we present a brief review of TGV. Then, the new variational models for speckle removal are proposed along with the analysis of the existence of the solutions. Afterwards, we present, in Section III, the explicit numerical schemes to solve the proposed models. Subsequently, Section IV describes the experiment results. In the final Section V, the concluding remarks are drawn. II. P ROBLEM F ORMULATION In this section, we start with a brief review of TGV, and then propose two new variational models based on it, followed by the analysis of the existence of a solution to the models. A. Total Generalized Variation For notational convenience, we assume that  ⊂ Rd is a non-empty, open and connected set. Note that, d ∈ Z, d ≥ 1 is a fixed space dimension. For the images used in this paper, we have d = 2. Then, the TGV of order k [31] and positive weights α = (α0 , . . . , αk−1 ) is defined as   k vd x|v ∈ C k (, Sym k (Rd )), c  udiv

(8) TGVkα (u) = sup

divl v ≤ αl , l = 0, . . . , k − 1, ∞ where Symk (Rd ) denotes the space of symmetric tensors of order k with arguments in Rd , Cck (, Symk (Rd )) denotes the vector space of compactly supported symmetric tensor field, and sup(x) or supremum(x) denotes the least upper bound of a set x. In the case of k = 2, Sym2 (Rd ) corresponds to the space of all symmetric d × d matrices. When k = 3, Sym3 (Rd ) is equivalent to the space S d×d×d of all symmetric d × d × d tensors. Note that TGVkα (u) is a semi-norm that is 0 for all polynomials of degree less than or equal to k − 1. As a consequence, the image reconstruction with TGV regularization leads to piecewise polynomial intensities. More details on TGV can be found in [31]. Specifically, we use the second order TGV as follows,   udiv2 wd x|w ∈ Cc2 (, S d×d ), 2  (9) TGVα (u) = sup w∞ ≤ α0 , divw∞ ≤ α1 , where the infinity norm of w and div(w) are defined as ⎛ ⎞1 2 d   2   ⎝ ⎠ w∞ = sup wl j (x) , x∈

divw∞ = sup

x∈

l, j =1

 d 

 21 |(divw)l |2

and the divergences are given by (divw)l = l ≤ d, and div2 w =

d  l, j =1

,

l=1

∂ 2 wl j ∂ xl ∂ x j

d  j =1

∂wl j ∂x j

,1 ≤

.

To analyze the properties of TGV2α (u), we need its discrete version following essentially the presentation in [36]. For this

1833

purpose, the basic related operators should be reviewed beforehand [36]. First, we denote the space Cc2 (1 , R) by U , the space Cc2 (1 , R2 ) by V , and Cc2 (1 , S 2×2 ) by W . In other words, the components of v ∈ V are v 1 ∈ U and v 2 ∈ U . Likewise, the components of w ∈ W are w11 ∈ U , w12 ∈ U and w22 ∈ U respectively. Second, the first-order forward and backward difference operators for u ∈ U are defined as ∂x+ u, ∂ y+ u, ∂x− u and ∂ y− u, whose concrete forms can be referred to [36]. Third, the gradient, the symmetrized gradient as well as the divergence operator are expressed as follows [36]  +  ∂ u ∇ :U → V, ∇u = x+ , ∂y u  1 − 1 − 2  ∂ −v 1 2 (∂ y v + ∂ x v ) , E : V → W, E(v) = 1 − 1x − 2 ∂ y− v 2 2 (∂ y v + ∂ x v ) div1 : V →U, div1 v = ∂x− v 1 + ∂ y− v 2,   + 11 ∂x w + ∂ y+ w12 . div2 : W → V, div2 w = ∂x+ w12 + ∂ y+ w22 Note that, the aforementioned operators are adjoint to each other, i.e., (div1 )∗ = −∇ and (div2 )∗ = −E. With the constraint div2 w = v in (9), the discretized TGV2α (u) of u ∈ U can now be given by  u, div1 v |(v, w) ∈ V × W, div2 w = v, 2 TGVα (u) = max . w∞ ≤ α0 , v∞ ≤ α1 , Based on [33], the discretized TGV2α (u) can be rewritten as TGV2α (u) = min α1 ∇u − v1 + α0 E(v)1 . v∈V

(10)

Here, the definition in (10) provides a way of balancing between the first and second derivatives of u. An intuitive interpretation of (10) is given as

follows. In the neighborhood

of edges, it holds that ∇ 2 u 1 is considerably larger than ∇u1 , so it is beneficial for choosing v = 0 to minimize (10), and hence the function TGV2α (u) resembles TV. On the other hand, ∇ 2 u is locally small in smooth regions of the image, it is therefore favorable to choose v = ∇u in these regions, resulting locally in the penalization of ∇ 2 u. B. New Variational Model In this subsection, we propose two novel variational models based on TGV2α (u) for despeckling. Before describing the models, a brief review of the derivation process for the functional (2) is provided, as shown in [18]. If the radar cross section (RCS) u of a target is given, the conditional probability density function [34] of the observed pixel intensity Y can be represented as p(Y|u) = 1 L L L−1 − LY e u . Aubert and Aujol [18] assumed that u (L) ( u ) Y follows a Gibbs prior: p(u) = Z1 exp(−γ φ(u)) where Z is a normalizing constant, γ denotes a related parameter, and φ is a non negative given function. Then, with the premise that p(Y) is constant [18], maximizing the log-likelihood log( p(u|Y)) = log( p(Y|u)) + log( p(u)) − log( p(Y)) amounts to minimizing the formulation as below  L (logu + Yu −1 )d x + (u), γ 0

1834

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

 where (u) = 0 φ(u)d x. Aubert and Aujol [18] supposed (u) to be TV(u), which is necessarily suboptimal, especially when the underlying images should be assumed piecewise linear instead of constant. Therefore, we incorporate the TGV penalty into the existing data fidelity term for speckle removal, and develop two novel variation despeckling models. Due to the MAP criterion with TGV regularization, our despeckling model can be formulated as  Y (log u + )d x + TGV2α (u) (11) u ∗ = arg min λ u u>0 0 where λ = γL and TGV2α (u) is defined as (9). As is well known, TV results are significantly corrupted by staircasing artifacts in regions where the assumption that image consists of piecewise-constant areas is violated. By the analysis of (8) and (10), TGV works equally well in these regions. Moreover, the global solution of (11) is hard to find because the term of logu + Yu −1 is non-convex [18]. To resolve this problem, we apply the logarithmic transformation [19], [22] which yields the exponential model:  ∗ ∗ (u l +Ye−ul )d x +TGV2α (u l ), u ∗f inal = eul . u l = arg min λ ul

0

(12) For notational convenience, we replace u l by u in (12), which means that the final solution can be simply expressed as ∗ u ∗f inal = eu . Accordingly, the proposed exponential model can be finally expressed as  ∗ ∗ u = arg min λ (u + Ye−u )d x + TGV2α (u), u ∗f inal = eu .

Likewise, G(u) for the model (13) can be represented as  λ 0 (u + Ye−u )d x. Then, the general problem [36] reads as u ∗ = arg min G(u) + TGV2α (u) u∈L p (0 )

(15)

where p ∈ (1, ∞), and p ≤ d/(d − 1). Note that, we have d = 2 in our application, and therefore p ≤ 2. Next, the kernel of TGV2α (u) [36] which is abbreviated as ker(TGV2α (u)), should be defined as the set of u such that TGV2α (u) = 0. The definition is given by   u : 0 → R|u(x) = Ax + b a.e. for some A ∈ R1×d , b ∈ R . Then, a linear mapping P [36] is defined as a projection operator onto the kernel of TGV2α (u):   P : L p (0 ) → ker TGV2α (u) . Besides, it holds that P|ker(TGV2 (u)) = I and P 2 = P where I α denotes the identity matrix. To prove the existence of solutions for (15), the following coercivity assumption on G should be made [36]: For any sequence {u n } in L p (0 ), it follows that 



n  

Pu → ∞ and (I− P) u n bounded⇒ G u n → ∞. p p (16)

If (16) is satisfied, there exists at least one minimizer for (15) [36]. Under these prerequisites, we have the following proposition. Proposition 1: Under the assumptions of Y ≥ 0 and Y ∈ L ∞ (0 ) that Y generally satisfies, the proposed model (13) satisfies the requirement (16). As a consequence, there must u 0 (13) exist a minimizer for (13). Proof: To verify the existence of a minimizer of model p with p ≤ 2, we first Furthermore, due to the inherent shortcomings of the logarith- (13) in L (0 ) for some  p ∈ (1, ∞) −u mic transformation, we also propose the classical I-divergence observe that G (u) = λ 0 (u + Ye )d x is proper, convex and lower semi-continuous [37]. Moreover, the functional model with TGV regularization: G (u) is strictly convex if Y > 0. We can also note that  ∞ u ∗ = arg min λ (u − Y log u)d x + TGV2α (u). (14) Y ∈ L (0 ) holds, which is suitable since the noisy image Y always satisfies this constraint. Next, to establish property u>0 0 (16), choose P : L 1 (0 ) → ker(TGV2α (u)) as an arbitrary Note that, due to the scale invariant property of TV [24], linear continuous projection. If, for a sequence {u n } ⊂ L p (0 ) the solution of (5) is equal to that of (4). Nevertheless, this it holds that Pu n  p → ∞, then also Pu n 1 → ∞ kind of special property is not satisfied by TGV, which means since all norms are equivalent on the finite-dimensional space that the I-divergence model (14) with TGV reguarization does ker(TGV2α (u)). Consequently, the following inequality holds not share the same solution of (13). As a matter of fact, the as P is continuous on L 1 (0 ), model (14) is more appropriate for Poisson noise reduction     n than despeckling. However, numerical experiments show that (u n + Ye−u )d x ≥ λ un d x G un = λ the model (14) obtains satisfactory results, and is as efficient 0

0

as (13). Therefore, we study both of the models (13) and (14) = λ u n 1 ≥ λc Pu n 1 for better observing the performance of TGV. To obtain existence of solutions for (13) and (14), we for some constant c > 0. Accordingly, G (u n ) → ∞ as follow the approach proposed in [36]. In order to describe n → ∞, and therefore (16) is satisfied. Together with the variational problems with TGV2α regularization, a proper, the observation that G is proper, convex, and lower semiconvex and lower semi-continuous functional G : L p (0 ) → continuous, this implies that there is a solution to (13). Proposition 2: Under the assumptions of Y ≥ 0 and Y ∈ (−∞, ∞] should be defined in a general form. Note that, the functional G has different formulations for different L ∞ (0 ) that Y generally satisfies, the proposed model (14) applications. While removing the additive Gaussian noise, satisfies the requirement (16). As a consequence, there must G(u) is given by u − f 2 where f is the noisy image. exist a minimizer for (14).

FENG et al.: SPECKLE REDUCTION VIA HIGHER ORDER TOTAL VARIATION APPROACH

Proof: To verify the existence of a minimizer of model (14) in L p (0 ) for some  p ∈ (1, ∞) with p ≤ 2, we first observe that G (u) = λ 0 (u − Y log u)d x is proper, convex and lower semi-continuous if u > 0 and Y ≥ 0 [37]. To establish (16), Y ∈ L ∞ (0 ) holds. It is suitable as the noisy image Y always satisfies this constraint. Next, a linear continuous projection P : L 1 (0 ) → ker(TGV2α (u)) is defined. Assume a sequence {u n } ⊂ L p (0 ) satisfying Pu n  p → ∞, then Pu n 1 → ∞ holds since all norms are equivalent on the finite-dimensional space ker(TGV2α (u)). For the coercivity, one may assume that u n > 0 for each n without loss of generality. Please note that for u, μ > 0, due to the convexity of − log, the subgradient inequality gives log u ≤ log μ + (u − μ)/μ. Choosing μ such that Yμ∞ ≤ 1/2 and log μ − 1 ≥ 0,   ≤ μ2 logμ + 12 u − 12 μ. we get Ylog(u) ≤ Y log(μ) + u−μ μ For simplicity, we abbreviate μ2 logμ + 12 u − 12 μ as 21 u + M where M = μ2 logμ − 12 μ. This implies that,      un − M dx un − G un ≥ λ 2 0  

λ λc n n

= Md x ≥ dx u 1−λ Pu 1 − λM 2 2 0 0

for some constant c > 0. Accordingly, G (u n ) → ∞ as n → ∞ and (16) is satisfied. Together with the observation that G is proper, convex, and lower semi-continuous, this implies that there is always a solution to (14). For a numerical realization, we discretize the problems in (13) and (14) in a straightforward way:

 ∗ u ∗ = arg min λ u +Ye−u , 1 +TGV2α (u), u ∗f inal = eu , (17) u

u ∗ = arg min λ u − Y log u, 1 + TGV2α (u),

(18)

u>0

where TGV2α (u) is defined as (10). Through the above analysis together with the experimental results, we can safely say that the proposed models are stable and can produce anticipated effects on staircasing artifacts removal. III. E XPLICIT N UMERICAL S CHEME In the following, to compute the minimizers of (17) and (18), we employ the primal-dual ascent-descent algorithm [38] because of its fast convergence and applicability for a wide range of problems. Certainly, other algorithms might be suited and efficient as well. In its general form, the primal-dual method finds a saddle point for the problem min max Kx, y + F (x) − H(y)

x∈X y∈Y

(19)

where X and Y are Hilbert spaces, K : X → Y is a linear and continuous mapping, and F : X → (−∞, ∞], H : Y → (−∞, ∞] are proper, convex and lower semi-continuous functionals. Then, the primal-dual algorithm [38] for the solution of (19) is defined through the iterations as follows. Choose σ, τ > 0   such that σ τ K2 < 1, initial values x 0 , y 0 ∈ X × Y and

1835

Algorithm 1 Prima-Dual Method (PD) for Speckle Removal With TGV Penalty (PDTGV) Input: The observed noisy image Y Initialization: 2√ Choose σ > 0, τ > 0 such that σ τ ≤ ; 17+ 33 Choose λ > 0, k Outer > 0 and a stopping criterion c;     Choose u 0 , v 0 ∈ U × V , p 0 , q 0 ∈ V × W , and set u¯ 0 = u 0 , u˜ 0 = u 0 , v¯ 0 = v 0 ;

u k −u k−1 While k < k Outer and u k−1  2 > c do 2    1) pk+1 ← proj P  pk + σ ∇ u¯ k− v¯ k ; 2) q k+1 ← proj Q q k + σ E(v¯ n ) ; 3) u˜ k+1 ← u k +τ div1 pk+1 ; 4) if F (u) = λ u − Y log u, 1

2 u˜ k+1 −λτ + (λτ −u˜ k+1 ) +4λτ Y k+1 u ← , 2

 if F (u) = λ u + Ye−u , 1 Apply Newton method to solve for u k+1 ; k; 5) u¯ k+1 ← 2u k+1 −   uk+1 k+1 k 6) v ←v +τ p + div2 q k+1 ; 7) v¯ k+1 ← 2v k+1 − v k ; 8) k ← k + 1; end while k Final Solution:u ∗f inal = eu for (17) and u k for (18).

x¯ 0 = x 0 . The iterations read as ⎧ k+1 ⎪ = (I + σ ∂H)−1 (y k + σ K x¯ k ), ⎨y x k+1 = (I + τ ∂F)−1 (x k − τ K∗ y k+1 ), ⎪ ⎩ k+1 x¯ = 2x k+1 − x k .

(20)

Here, the operators (I + τ ∂H)−1 and (I + σ ∂F)−1 can be characterized as the solutions of x − x ¯ 22 + F (x), 2τ x∈X y − y¯ 22 y ∗ = (I + σ ∂H)−1 ( y¯ ) ⇔ y ∗ = arg min + H (y). 2σ y∈Y

x ∗ = (I + τ ∂F )−1 (x) ¯ ⇔ x ∗ = arg min

According to [38], if X and Y are finite dimensional, the primal-dual algorithm converges to a saddle point (x ∗ , y ∗ ) of (19) provided that a saddle point exists. For addressing the problem (17) and (18), we rewrite them as min

(u,v)∈U ×V

F (u) + α1 ∇u − v1 + α0 E(v)1

(21)



where F(u) denotes λ u + Ye−u , 1 for (17) and λ u − Y log u, 1 for (18). In order to solve (21) using the primal-dual algorithm, we reformulate it as a convexconcave saddle-point problem which is derived by the duality principles: min

max

(u,v)∈U ×V ( p,q)∈ P×Q

F (u) + ∇u − v, p + E(v), q

(22)

where p, q are the dual variables associated with the sets given by     P = p ∈ V |  p∞ ≤ α1 , Q = q ∈ W | q∞ ≤ α0 .

1836

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Introducing indicator functionals, i.e., ! 0 if x ∈ Q, IQ (x) = ∞ else,

One can apply Newton method to address this sub-problem since the Hessian matrix of (27) only has values along the diagonal, which yields great convenience for computing. Second, if F(u) = λ u − Y log u, 1 for (18), we should solve

the discrete functional in (22) can be reformulated as max

min

(u,v)∈U ×V ( p,q)∈V ×W

u k+1 = arg min λ u − Y log u, 1 +

F(u)+ ∇u −v, p + E(v), q

−I{·∞ ≤α1 } ( p)−I{·∞ ≤α0 } (q) .

u

(23)

Then, the problem (23) admits the structure (19) if one chooses " # ∇ −I X = U × V, Y = V × W, K = 0 E F (x) = F (u, v) = F(u), H(y) = H( p, q) = I{·∞ ≤α1 } ( p) + I{·∞ ≤α0 } (q). Note that, the functionals F and H can be performed componentwise since they only depend on one component of x and y, respectively. Moreover, the functional H corresponds to projection operators on the sets P and Q. Consequently, the iterations of p and q can be treated as solving for y in (20). Thus, the explicit expressions for pk+1 and q k+1 are given by    p − p ¯ 22 , p¯ = pk +σ ∇ u¯ k − v¯ k p∈P 2σ    = proj P ( p) ¯ = proj P pk + σ ∇ u¯ k − v¯ k ,   q − q ¯ 22 , q¯ = q k +σ E v¯ k = arg min q∈Q 2σ    = projQ (q) ¯ = proj Q q k + σ E v¯ k ,

p k+1 = arg min ⇒ pk+1 q

k+1

⇒ q k+1

where the Euclidean projectors proj P ( p) , proj Q (q) denote the respective projection onto the aforementioned convex sets P and Q. They can be easily computed via component-wise operations [36], ¯ = proj P ( p)

p¯ q¯    , proj Q (q)  . (24) ¯ = | p| ¯ ¯ max 1, α1 max 1, |αq| 0

Likewise, the subproblems of u and v can be treated as solving for x in (20), and can be formulated as

u − u˜ k+1 2 k+1 2 u , u˜ k+1 = u k +τ div1 pk+1 = arg min F (u) + 2τ u∈U (25) as well as   v − v ˜ 22 , v˜ = v k + τ pk+1 + div2 q k+1 v∈V 2τ   (26) = v˜ = v k + τ pk+1 + div2 q k+1 .

v k+1 = arg min ⇒ v k+1

There are two situations to consider for solving the subproblem 

in (25). First, if F(u) = λ u + Ye−u , 1 for (17), we should solve

2

 1

u k+1 = arg min λ u + Ye−u , 1 +

u − u˜ k+1 . (27) 2 2τ u∈U

(28)

This problem is easy to solve by setting its first-order derivative as zero, leading to the following equation: ˜ − λτ Y = 0, u 2 + (λτ − u)u whose solution is u˜ − λτ +

(λτ − u) ˜ 2 + 4λτ Y

. (29) 2 It is worth noting that the computations of (26), (27) and (29) are computed component-wisely. Finally, in order to choose the step sizes σ and τ , the norm of K is √ estimated according to [36] and given by 17+ 33 2 K < < 12. Therefore, the proposed algorithm 2 summarized in Algorithm √ 1 is known to converge for step33), which means that the choice of sizes σ τ ≤ 2/(17 + √ σ = τ = 1/ 12 leads to convergence. Additionally, the PDTGV algorithm for the exponential model (13) is abbreviated as PDTGV_EXP, and PDTGV_Idiv for the I-divergence model (14). In the following, the notation PDTGVs involves PDTGV_EXP and PDTGV_Idiv. u k+1 =

as well as

u − u ˜ 22 . 2τ

IV. E XPERIMENTAL R ESULT In this section, we report numerical experiments on speckle reduction with 10 test images in Fig. 1, including five synthesized images contaminated by speckle noise and five real SAR images. The real SAR image SAR3 is a TerraSAR-X image whose number of looks L is 1. Although the values of L of the other real SAR images from http://www.sandia.gov/radar/sardata.html are unknown, we can estimate the equivalent number of looks (ENL) according to [34] from the homogeneous regions within the numbered white rectangles in Fig. 1(f), (g), (i), and (j). The results are ENL = 5 for SAR1 and SAR2, and ENL = 3 for SAR4 and SAR5. Therefore, we set the values of L just as presented in Fig. 1. Additionally, we also compare the proposed algorithms with two kinds of methods: 1) the nonlocal methods PPBit [13] and SAR-BM3D [11], which usually provide state-of-the-art results; 2) the variational methods with TV regularization, including MIDAL [25], ADMM-III [24] and AMAST [26]. For the AMAST methods, we only provide the despeckled images by AMAST-a-EXP for the sake of simplicity. In the experiment, the number of looks L is set to be 1 or 3 for synthesized images, fulfilling the practical application. For the PPBit filter, a similarity window of size 7 × 7 and a search window of size 21 ×21 are used, and the parameters are set as α = 0.92 and T = 0.2. Note that, four iterations are employed while implementing the PPBit filter. Two quality measures are taken to evaluate the despeckled results. The first one is the peak signal-to-noise ratio (PSNR) 2552 n which is calculated by PSNR(dB) = 10log10 u−u ∗ 2 , where 2

FENG et al.: SPECKLE REDUCTION VIA HIGHER ORDER TOTAL VARIATION APPROACH

1837

Fig. 1. Ten different test images. (a) Pic1(256 × 256). (b) Pepper(256 × 256). (c) Pic2(256 × 256). (d) Rem1(350 × 350). (e) Rem2(420 × 420). (f) SAR1(256 × 256), L = 5. (g) SAR2(256 × 256), L = 5. (h) SAR3(512 × 512), L = 1. (i) SAR4(1180 × 1180), L = 3. (j) SAR5(1195 × 1195), L = 3.

Fig. 2. Performance comparison of different algorithms on Pic1. (a) Noisy image (L = 3). (b) PPBit. (c) SAR-BM3D. (d) MIDAL. (e) PDTGV_EXP. (f) PDTGV_Idiv. TABLE I C OMPARISON OF THE P ERFORMANCE OF THE A LGORITHMS : PPB- IT, SAR-BM3D, MIDAL, AMAST- A -I-DIV, AMAST- A -EXP, ADMM-III, PDTGV S IN T ERMS OF PSNR AND SSIM

n is the size of image, u ∗ and u are the estimated and original images. The second one is the structural similarity index (SSIM) [41], which measures the structural detail similarity between u ∗ and u. The SSIM is given by SSIM(u ∗ , u) =

(2μu ∗ μu )(2σuu ∗ + c2 ) (μ2u ∗ + μ2u + c1 )(σu2 + σu2∗ + c2 )

where μu and μu ∗ are the mean values of image u and u ∗ , σu and σu ∗ denote their standard deviations, respectively, and σuu ∗ is the covariance of u and u ∗ . Besides, it holds that c1 = (0.01D)2 and c2 = (0.03D)2 , where D is the dynamic range (255 for 8-bit grayscale images). The range of SSIM value lies in [0, 1] with 1 for the perfect quality.

1838

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Fig. 3. Performance comparison of different algorithms on the image Pic2. (a) Noisy image (L = 3). (b) PPBit. (c) SAR-BM3D. (d) ADMM-III. (e) MIDAL. (f) AMAST-a-EXP. (g) PDTGV_EXP. (h) PDTGV_Idiv. (i) Noisy image (L = 1). (j) PPBit. (k) SAR-BM3D. (l) ADMM-III. (m) MIDAL. (n) AMAST-a-EXP. (o) PDTGV_EXP. (p) PDTGV_Idiv.

We terminate all the algorithms except the PPBit and

SAR-BM3D by the condition u k − u k−1 2 ≥ c u k−1 2 , where c = 3 × 10−5 . The number of Newton iterations is 15 for the MIDAL and PDTGV_EXP. All simulations listed here are implemented in C laguage with interface to Matlab 7.12 (R2011a) through the mex function on a laptop equipped with 2.40GHz CPU and 4G RAM memory. A. Choice of the Regularization Parameters The proposed TGV-based models, like TV-based models, need a proper balance between the data fidelity term and the regularization term. Underweighting of the TGV term leads to residual noise, while overweighting removes image features. For the sake of simplicity, we generally set the empirical value

of λ. Although finer tuning of the parameters may lead to better performance, the results under the following parameter settings are consistently promising already. For synthesized images, the choices of λ = 1.5 for L = 3 and λ = 0.9 for L = 1 yield good results empirically. For real SAR images, the following weights are used: λ = 1.9 for L = 1, λ = 3 for L = 3, and λ = 5 for L = 5. For compatibility with the conventional TV semi-norm, α1 = 1 is adopted in this paper. Moreover, the default value α0 = 2 produces satisfactory results for most applications and does not need to be tuned [33]. B. Results on Synthetic Images Fig. 2 presents the comparison between different methods on the piecewise affine image Pic1 contaminated by the

FENG et al.: SPECKLE REDUCTION VIA HIGHER ORDER TOTAL VARIATION APPROACH

1839

Fig. 4. Performance comparison of different algorithms on the image Rem1. (a) Noisy image (L = 3). (b) PPBit. (c) SAR-BM3D. (d) ADMM-III. (e) MIDAL. (f) AMAST-a-EXP. (g) PDTGV_EXP. (h) PDTGV_Idiv. (i)-(p) Zoom of the images in (a)-(h). (i) Noisy image. (j) PPBit. (k) SAR-BM3D. (l) ADMM-III. (m) MIDAL. (n) AMAST-a-EXP. (o) PDTGV_EXP. (p) PDTGV_Idiv.

speckle noise with L = 3. All of the methods remove noise while maintaining sharp edges efficiently. But since the original image consists of regions with vanishing second derivative, it is not surprising that TV results are significantly corrupted by staircasing artifacts. The proposed TGV-based models, on the other hand, are capable of recovering almost perfectly by taking into account the second derivative as well as edges. Note that, the nonlocal technique SAR-BM3D are affected by ghost artifacts, while PPBit exhibits some smooth brushstrokes, because of the attempt to recognize structures even when they are absent. This is more evident for homogeneous areas. Such problems encountered with the

nonlocal methods are also confirmed by the results presented in Figs. 3–5. In this experiment, the PSNR index for SAR-BM3D and PDTGVs is evidently better than that of others. In terms of SSIM, PDTGVs exceeds all the others. Therefore, it is reasonable to conclude that TGV result is best in the sense that it has piecewise smooth intensities and sharp edges without staircasing artifacts. Fig. 3 presents the despeckled results of the image Pic2 corrupted by speckle noise with an equivalent number of look L = 1 and 3. From observation, in regions where the piecewise constant assumption is violated, the results of the TV methods

1840

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Fig. 5. Performance comparison of different algorithms on the image Rem2. (a) Noisy image (L = 3). (b) PPBit. (c) SAR-BM3D. (d) ADMM-III. (e) MIDAL. (f) AMAST-a-EXP. (g) PDTGV_EXP. (h) PDTGV_Idiv. (i)-(p) Zoom of the images in (a)-(h). (i) Noisy image. (j) PPBit. (k) SAR-BM3D. (l) ADMM-III. (m) MIDAL. (n) AMAST-a-EXP. (o) PDTGV_EXP. (p) PDTGV_Idiv.

are significantly degraded by staircasing artifacts. However, PDTGVs works equally well in these regions, as pointed out in Fig. 2. Moreover, by contrasting Fig. 3(f) and (g), it can be noticed that the edges and piecewise constant regions introduced in the TV results are more emphasized than the TGV results. In this test, the SAR-BM3D performs best on detail preservation, although it brings in some annoying artifacts. In Figs. 4 and 5, one can observe that PDTGVs is able to recover smooth regions as well as discontinuities at object boundaries. It is worth taking a closer look at the zoom of the despeckled images shown in Figs. 4(i)–(p) and 5(i)–(p). In Figs. 4(o) and (p) and 5(o) and (p), PDTGVs overcomes

the staircasing artifacts around the areas with higher-order smoothness. Meanwhile, the tendency of TV to produce staircasing artifacts is again presented in Figs. 4(l)–(n) and 5(l)–(n). Additionally, if the noisy image has the periodic patterns, the nonlocal methods perform better than the other ones, just as presented in the bottom right corner of Fig. 4(j) and (k). Nevertheless, the nonlocal methods bring in some tiny distortions on the surfaces of the buildings which yet recover better in the TGV results. Note that, the centers of the domes in image Rem2 should be assumed to be piecewise-constant while the other parts should be assumed to be piecewise linear. As described in Fig. 5, PDTGVs restores the dome to the piecewise-linear results since the

FENG et al.: SPECKLE REDUCTION VIA HIGHER ORDER TOTAL VARIATION APPROACH

1841

Fig. 6. Performance comparison of different algorithms on the image SAR2. (a) Noisy image (L = 5). (b) PPBit. (c) SAR-BM3D. (d) MIDAL. (e) PDTGV_EXP. (f)–(j) Zoom of the images in (a)–(e). (f) Noisy image. (g) PPBit. (h) SAR-BM3D. (i) MIDAL. (j) PDTGV_EXP.

Fig. 7. Performance comparison on the parts of SAR1, SAR4 and SAR5. Left column: Noisy image. Middle column: MIDAL. Right column: PDTGV_EXP.

piecewise-constant parts are totally hidden by the severe noise. The TV-based methods recover piecewise-constant areas which seem unfaithful. Therefore, both of the results produced by PDTGVs and TV-based methods can not accord with the underlying true image, leading to their similar performance in PSNR/SSIM. However, the PDTGVs removes the staircasing artifacts produced by TV. Table I lists the PSNR (in dB) and SSIM of the different algorithms on five synthetic images in Fig. 1. Comparing the indexes in Table I and the results in Figs. 2–5, we can easily draw the following conclusions. 1) Overall, the best PSNR is guaranteed by SAR-BM3D, followed by PDTGVs and PPBit. If the image contains

some piecewise linear features, PDTGVs becomes more competitive. Extremely, while testing with the images Pic1 and Pic2 consisting of regions with vanishing second derivative, the PSNRs of PDTGVs are even higher than those of SAR-BM3D. 2) In most cases, PDTGVs generates a slightly higher PSNR than the TV-based methods. Since TGV does not incorporate neighborhood information, PDTGVs does not lead to a significantly higher PSNR, although it provides more accurate recovery of smooth regions. Actually, there are still some cases in which the performances of TV-based methods and PDTGVs are comparable, just as shown in Fig. 5. Note that, in order to preserve more details, the automatic spatially dependent selection procedure of λ [40] can be adopted in our model. 3) In terms of SSIM, the best performance is provided by PDTGVs for recovering the piecewise smooth images, indicating that our method is the most powerful in geometry-preserving in this case, which can noticeably be visually perceived in Figs. 2–5. 4) In the visual quality, the staircasing artifacts, which are typical for TV-based methods, do not appear when the proposed PDTGVs are used. Besides, the nonlocal methods produce the so-called ghosts, structured signal-like patches that appear in homogeneous areas of the image. This phenomenon is more evident for SAR-BM3D, as shown in Figs. 2(c) and 3(c) and (k). C. Results on Real SAR Images In this subsection, we take five tests with real SAR images to observe the performance of the proposed method. Note that, we only provide the despeckled images by PDTGV_EXP since the TGV results are visually similar. In Fig. 6, all of the methods remove the speckle in homogeneous regions and preserve edges well. By observation,

1842

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 4, APRIL 2014

Fig. 8.

Performance comparison of different algorithms on the image SAR3(L = 1). (a) SAR-BM3D. (b) MIDAL. (c) PDTGV_EXP. TABLE II T HE C OMPUTATION T IME (S ECONDS ) OF D IFFERENT M ETHODS

the TV-based result in Fig. 6(i) introduces additional staircasing artifacts, leading to a more blocky, less authentic result compared to the original image SAR2. Besides, it is notable that in the PPB filtering result as shown in Fig. 6(b), some weak objects especially in the lower left corner are neglected or attenuated, which are yet recovered better via the other approaches. Fig. 7 presents the close-up results of the areas in the white squares as shown in Fig. 1(f), (g), (i), and (j), which again demonstrate the advantages of PDTGV_EXP over TV methods. While dealing with the one-look scene SAR3 in Fig. 8, some obscure differences can be observed that TV-based algorithms brings in staircasing artifacts which are almost absent when using PDTGV_EXP. Here, the underlying image intensity of SAR3 should be assumed piecewise linear instead of constant. Therefore, it is not surprising that TV-based result recovers piecewise-constant cartoon-like image. The PDTGV_EXP outputs better results as it takes into consideration the higher order smoothness. D. Comparison of the Computation Time Table II reports the CPU time of different methods. The AMAST algorithms outperform all others in CPU time consuming, thanks to their simple but effective iterative steps. In addition, the speed of PDTGVs is remarkably higher than that of SAR-BM3D. Particularly, the method PDTGV_Idiv avoids any inner iterations or inverse operations, whereby it achieves huge speed advantage. V. C ONCLUSION In this study, we propose two novel despeckling models based on TGV regularization, and an efficient solving algorithm originated from the primal-dual method. Equipped with

options to adjust the high degrees of smoothness, the proposed approaches PDTGV_EXP and PDTGV_Idiv reduce the staircasing artifacts commonly seen in traditional TV-based algorithms. In comparison with the state-of-the-art method SAR-BM3D, PDTGVs especially the method PDTGV_Idiv, achieve much computational advantages, and with comparable or even better results for the images containing piecewise linear features. Moreover, PDTGVs can be effectively accelerated through parallelization, and hence can be efficiently applied on hugesize SAR images. In addition, in order to preserve more details, we suggest two subjects for future investigations. First, some sparse transforms can be brought into the model, e.g., the Shearlet transform [43]. Second, the automatic spatially dependent selection procedure of λ [40] can be introduced.

R EFERENCES [1] J. Goodman, “Some fundamental properties of speckle,” J. Opt. Soc. Amer., vol. 66, no. 11, pp. 1145–1150, 1976. [2] J. J. Ranjani and S. Thiruvengadan, “Dual-tree complex wavelet transform based SAR despeckling using interscale dependence,” IEEE Trans. Geosci. Remote Sensing, vol. 48, no. 6, pp. 2723–2731, Jun. 2010. [3] A. Achim, E. E. Kuruoglu, and J. Zerubia, “SAR image filtering based on the heavy-tailed Rayleigh model,” IEEE Trans. Image Process., vol. 15, no. 9, pp. 2686–2693, Sep. 2006. [4] M. I. H. Bhuiyan, M. O. Ahmad, and M. N. S. Swamy, “Spatially adaptive wavelet-based method using the Cauchy prior for denoising the SAR images,” IEEE Trans. Circuits Syst. Video Technol., vol. 17, no. 4, pp. 500–507, Apr. 2007. [5] T. Bianchi, F. Argenti, and L. Alparone, “Segmentation-based MAP despeckling of SAR images in the undecimated wavelet domain,” IEEE Trans. Geosci. Remote Sensing, vol. 46, no. 9, pp. 2728–2742, Sep. 2008. [6] J. J. Ranjani and S. J. Thiruvengadam, “Generalized SAR despeckling based on DTCWT exploiting interscale and intrascale dependences,” IEEE Geosci. Remote Sensing Lett., vol. 8, no. 3, pp. 552–556, May 2011.

FENG et al.: SPECKLE REDUCTION VIA HIGHER ORDER TOTAL VARIATION APPROACH

[7] F. Argenti, T. Bianchi, A. Lapini, and L. Alparone, “Fast MAP despeckling based on Laplacian–Gaussian modeling of wavelet coefficients,” IEEE Geosci. Remote Sensing Lett., vol. 9, no. 1, pp. 13–17, Jan. 2012. [8] H. Guo, J. Odegard, M. Lang, R. Gopinath, I. Selesnick, and C. Burrus, “Wavelet based speckle reduction with application to SAR based ATD/R,” in Proc. IEEE Int. Conf. Image Process., vol. 1. Jan. 1994, pp. 75–79. [9] L. Gagnon, “Speckle filtering of SAR images: A comparative study between complex wavelet-based and standard filters,” Proc. SPIE, vol. 3169, pp. 80–91, Oct. 1997. [10] A. Buades, B. Coll, and J. M. Morel, “A review of image denoising algorithms, with a new one,” Multiscale Model. Simul., vol. 4, no. 2, pp. 490–530, Jul. 2005. [11] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3D transform-domain collaborative filtering,” IEEE Trans. Image Process., vol. 16, no. 8, pp. 2080–2095, Aug. 2007. [12] M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” IEEE Trans. Image Process., vol. 15, no. 12, pp. 3736–3745, Dec. 2006. [13] C. A. Deledalle, L. Denis, and F. Tupin, “Iterative weighted maximum likelihood denoising with probabilistic patch-based weights,” IEEE Trans. Image Process., vol. 18, no. 12, pp. 2661–2672, Dec. 2009. [14] T. Teuber and A. Lang, “Nonlocal filters for removing multiplicative noise,” in Scale Space and Variational Methods. New York, NY, USA: Springer-Verlag, 2012, pp. 50–61. [15] S. Parrilli, M. Poderico, C. V. Angelino, and L. Verdoliva, “A nonlocal SAR image denoising algorithm based on LLMMSE wavelet shrinkage,” IEEE Trans. Geosci. Remote Sensing, vol. 50, no. 2, pp. 606–616, Feb. 2012. [16] Y. M. Huang, L. Moisan, M. K. Ng, and T. Y. Zeng, “Multiplicative noise removal via a learned dictionary,” IEEE Trans. Image Process., vol. 21, no. 11, pp. 4534–4543, Dec. 2012. [17] L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D, Nonlinear Phenomena, vol. 60, pp. 259–268, Nov. 1992. [18] G. Aubert and J. F. Aujol, “A variational approach to remove multiplicative noise,” SIAM J. Appl. Math., vol. 68, no. 4, pp. 925–946, 2008. [19] J. M. Bioucas-Dias and M. A. T. Figueiredo, “Multiplicative noise removal using variable splitting and constrained optimization,” IEEE Trans. Image Process., vol. 19, no. 7, pp. 1720–1730, Jul. 2010. [20] L. Denis, F. Tupin, J. Darbon, and M. Sigelle, “SAR image regularization with fast approximation discrete minimization,” IEEE Trans. Image Process., vol. 18, no. 7, pp. 1588–1600, Jul. 2009. [21] S. Durand, J. Fadili, and M. Nikolova, “Multiplicative noise removal using L1 fidelity on frame coefficients,” J. Math Imag. Vis., vol. 36, no. 3, pp. 201–226, Mar. 2010. [22] Y. M. Huang, M.K. Ng, and Y. W. Wen, “A new total variation method for multiplicative noise removal,” SIAM J. Imag. Sci., vol. 2, no. 1, pp. 20–40, 2009. [23] J. Shi and S. Osher, “A nonlinear inverse scale space method for a convex multiplicative noise model,” SIAM J. Imag. Sci., vol. 1, no. 3, pp. 294–321, 2008. [24] G. Steidl and T. Teuber, “Removing multiplicative noise by DouglasRachford splitting,” J. Math Imag. Vis., vol. 36, no. 2, pp. 168–184, Feb. 2010. [25] J. Eckstein and D. Bertsekas, “On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators,” Math. Program., Ser. A, B, vol. 55, no. 3, pp. 293–318, Jun. 1992. [26] H. Woo and S. Yun, “Alternating minimization algorithm for speckle reduction with a shifting technique,” IEEE Trans. Image Process., vol. 21, no. 4, pp. 1701–1714, Apr. 2012. [27] H. Woo and S. Yun, “A new multiplicative denoising variational model based on mth root transformation,” IEEE Trans. Image Process., vol. 21, no. 5, pp. 2523–2533, May 2012. [28] P. Tseng, “Applications of a splitting algorithm to decomposition in convex programming and variational inequalities,” SIAM J. Control Optim., vol. 29, no. 1, pp. 119–138, Jan. 1991. [29] A. Chambolle and P.-L. Lions, “Image recovery via total variation minimization and related problems,” Numer. Math., vol. 76, no. 2, pp. 167–188, 1997. [30] D. T. Kuan, A. A. Sawchuk, T. C. Strand, and P. Chavel, “Adaptive restoration of images with speckle,” IEEE Trans. Acoust., Speech, Signal Process., vol. 35, no. 3, pp. 373–383, Jan. 1987. [31] K. Bredies, K. Kunisch, and T. Pock, “Total generalized variation,” SIAM J. Imag. Sci., vol. 3, no. 3, pp. 492–526, 2010.

1843

[32] M. Benning, C. Brune, M. Burger, and J. Muller, “Higher-order TV methods–enhancement via Bregman iteration,” J. Sci. Comput., vol. 54, nos. 2–3, pp. 269–310, 2012. [33] F. Knoll, K. Bredies, T. Pock, and R. Stollberger, “Second order total generalized variation (TGV) for MRI,” Magn. Res. Med., vol. 65, no. 22, pp. 480–491, 2011. [34] C. Oliver and S. Quegan, Understanding Synthetic Aperture Radar Imaging. Raleigh, NC, USA: SciTech Publ., 2004. [35] V. Frost, K. Shanmugan, and J. Holtzman, “A model for radar images and its applications to adaptive digital filtering of multiplicative noise,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 2, no. 2, pp. 157–166, Mar. 1982. [36] K. Bredies, “Recovering piecewise smooth multichannel images by minimization of convex functionals with total generalized variation penalty,” Inst. Math. Sci. Comput., Graz Univ. Technol., Graz, Austria, Tech. Rep. 12-006, 2012. [37] R. T. Rockafellar and R. J.-B. Wets, Variational Analysis. New York, NY, USA: Springer-Verlag, 1998. [38] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” J. Math. Imag. Vis., vol. 40, no. 1, pp. 120–145, 2011. [39] J. S. Lee, “Digital image enhancement and noise filtering by use of local statistics,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 2, no. 2, pp. 165–168, Mar. 1980. [40] K. Bredies, Y. Dong, and M. Hintermüller, “Spatially dependent regularization parameter selection in total generalized variation models for image restoration,” Inst. Math. Sci. Comput., Graz Univ. Technol., Graz, Austria, Tech. Rep. 12-002, 2012. [41] Z. Wang, C. B. Alan, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Image Process., vol. 13, no. 4, pp. 600–612, Apr. 2004. [42] A. Mittal, A. K. Moorthy, and A. C. Bovik, “No-reference image quality assessment in the spatial domain,” IEEE Trans. Image Process., vol. 21, no. 12, pp. 4695–4708, Dec. 2012. [43] D. Labate, W.-Q. Lim, G. Kutyniok, and G. Weiss, “Sparse multidimensional representation using shearlets,” Proc. SPIE, vol. 5914, pp. 254–262, Sep. 2005.

Wensen Feng received the B.S. degree from Xidian University, Xi’an, China, in 2005. He is currently pursuing the Ph.D. degree with the Institute of Electronics, Chinese Academy of Sciences, Beijing, China. His research interests include image denoising, segmentation and restoration, and radar image processing.

Hong Lei received the B.S. degree from Xidian University, Xi’an, China, and the M.S. degree from the Beijing Institute of Technology, Beijing, China, in 1984 and 1991, respectively. He has been with the Institute of Electronic, Chinese Academy of Science, since 1991. His current research interests include microwave communication technology, microwave imaging, airborne and spaceborne SAR system design and their signal processing, and antenna theory and its engineering application.

Yang Gao received the B.S. degree in electronic engineering from the Nanjing University of Aeronautics and Astronautics, Nanjing, China, in 2009. He is currently pursuing the Ph.D. degree in highresolution synthetic aperture radar (SAR) signal processing technology with the University of Chinese Academy of Science, Beijing. His research interests include the high-resolution SAR imaging, SAR autofocus techniques, and image processing.