A Higher-Order MRF Based Variational Model for Multiplicative Noise Reduction

The Fields of Experts (FoE) image prior model, a filter-based higher-order Markov Random Fields (MRF) model, has been shown to be effective for many image restoration problems. Motivated by the successes of FoE-based approaches, in this letter we propose a novel variational model for multiplicative noise reduction based on the FoE image prior model. The resulting model corresponds to a non-convex minimization problem, which can be efficiently solved by a recently published non-convex optimization algorithm. Experimental results based on synthetic speckle noise and real synthetic aperture radar (SAR) images suggest that the performance of our proposed method is on par with the best published despeckling algorithm. Besides, our proposed model comes along with an additional advantage, that the inference is extremely efficient. Our GPU based implementation takes less than 1s to produce state-of-the-art despeckling performance.


I. INTRODUCTION
Images generated by coherent imaging modalities, e.g., synthetic aperture radar (SAR), ultrasound and laser imaging, inevitably come with multiplicative noise (also known as speckle), due to the coherent nature of the scattering phenomena.The presence of this noise prevents us from interpreting valuable information of images, such as textures, edges and point target, and therefore speckle reduction is often a necessary preprocessing step for successful use of classical image processing algorithms involving image segmentation and classification.The topic of speckle noise reduction (despeckling) has attracted a lot of research attention since early 1980s [17], [16].At present it has been extensively studied.Roughly speaking, the major despeckling techniques fall into four categories: filtering based methods in (1) spatial domain; or (2) a transform domain, e.g., wavelet domain; (3) nonlocal filtering and (4) variational methods.
Early filtering techniques in the spatial domain are developed under the minimum mean square error (MMSE) criterion [17], and then progress to more sophisticated and promising maximum a posterior (MAP) approaches [16].Recently, bilateral filtering has also been modified for despeckling [18].The emergence of wavelet transform in the early of 1990s, opened the way to a new generation of despeckling techniques.Y.J. Chen, R. Ranftl and T. Pock are with the Institute for Computer Graphics and Vision, Graz University of Technology, Inffeldgasse 16, A-8010 Graz, Austria.e-mail: ({cheny, ranftl, pock}@icg.tugraz.at)Wensen Feng is with University of Science and Technology Beijing, Beijing, 100190, China, also with the Institute of Automation, Chinese Academy of Sciences, Beijing, 100190, China.e-mail: (sanmumuren@gmail.com) Hong Qiao is with the Institute of Automation, Chinese Academy of Sciences, Beijing, 100190, China.e-mail: (hong.qiao@mail.ia.ac.cn)There were intensive studies of wavelet based despeckling approaches, see, for instance [5], [3] and references therein.
The nonlocal approaches, which can take the advantage of self-similarity commonly present in natural as well as SAR images, have been already introduced to SAR despeckling [13], [20], [10].By taking into account the peculiar features of multiplicative noise, the so-called SAR-BM3D algorithm [20], which is a SAR-oriented version of the well-known BM3D algorithm [12], exhibits the best published despeckling performance at present.
The last class of methods are variational ones, which minimize some appropriate energy functionals, consisting of a regularizer (also called image prior) and a data fitting term.Up to now, the well-known total variation (TV) has been widely used as a regularizer [22], [24], [4], and the total generalized variation (TGV) regularizer [6] also has been investigated in a recent work [14].
There is a long history of research on regularization technique for image processing problems.The recently proposed Fields of Experts (FoE) [21] image prior model (a higher-order filter-based MRF model), which explicitly characterizes the statistics properties of natural images, defines more effective variational models than hand-crafted regularizers, such as TV and TGV models.The variational model based on a discriminatively trained FoE prior can generate clear stateof-the-art performance for many image restoration problems, such as the additive white Gaussian noise (AWGN) denoising task, see for instance [7], [8].
Motivated by the results of [7], [8], it is interesting to investigate the FoE prior based model for despeckling.In this letter, we propose a novel variational approach for speckle removal, which involves an expressive image prior model -FoE and a highly efficient non-convex optimization algorithm -iPiano, recently proposed in [19].Our proposed method obtains strongly competitive despeckling performance w.r.t. the state-of-the-art method -SAR-BM3D [20], meanwhile, preserve the property of computational efficiency.

II. PROPOSED VARIATIONAL MODELS FOR DESPECKLING
In this section, we propose FoE prior based variational models for speckle noise removal and in particular for SAR images.We introduce an efficient algorithm to solve the corresponding optimization problems.
Given an observation image f corrupted by multiplicative noise, the FoE prior based despeckling model is defined as the following energy minimization problem  where u is the underlying unknown image, the first term is the FoE prior model, and the second part is the data fidelity term, derived from the multiplicative noise model.

A. The FoE prior utilized in our despeckling model
The FoE model is defined by a set of linear filters.According to [8], [7], the student-t distribution based FoE image prior model for an image u is formulated as , N is the number of pixels in image u, N f is the number of linear filters, k i is a set of learned filters with the corresponding weights θ i > 0, k i * u denotes the convolution of the filter k i with a two-dimensional image u, and ρ(•) denotes the non-convex Lorentzian potential function, which is derived from the student-t distribution.Note that the FoE energy is non-convex w.r.t.u.In this work, we make use of the learned filters of a previous work [8], as shown in Figure 1.

B. Modeling of the data term
Let f be the observed SAR image amplitude, which follows a Nakagami distribution depending on the underlying true image amplitude u, the square root of the reflectivity [15] with L the number of looks of the image (i.e., number of independent values averaged) and Γ is the classical Gamma function.
According to the Gibbs function, this likelihood leads to the following energy term via E = −logp(f |u) Combining this data term with the FoE prior model (II.2), we reach the following variational model where •, • denotes the common inner product.Note that the data term is not convex w.r.t.u > 0, which will generally make the corresponding optimization problem harder to solve.
There exists an alternative method to define a convex data term by using the classical Csiszár I-divergence model [11].Although the I-divergence data fitting term typically used in the context of Poisson noise, this seemingly inappropriate data term performs very well in the TV and TGV regularized variational models for despeckling [22], [14].Therefore, we also study this variant for data term modeling in this work.Following previous works of modeling the SAR image intensity [22], [14], the I-divergence based data term for the amplitude model is given by which is strictly convex w.r.t.u for u > 0. Then the variational model involving this convex data term is given by

C. Solving the variational despeckling models
Due to the non-convexity of the prior term, the proposed variational models impose generally very hard optimization problems, especially for the model (II.4) with a non-convex data term.In this work, we resort to a recently published nonconvex optimization algorithm -iPiano [19] to solve them.
The iPiano algorithm is designed for a structured nonsmooth non-convex optimization problem, which is composed of a smooth (possibly non-convex) function F and a convex (possibly non-smooth) function G: (II.7) iPiano is an inertial force enhanced forward-backward splitting algorithm with the following basic update rule where α and β are the step size parameters.
For the model (II.4) with a non-convex data term, we can convert the data term to a convex function via commonly used logarithmic transformation (i.e., w = logu).Therefore, the minimization problem is rewritten as with u = e w .Casting (II.8) in the form of (II.7), we see that , 1 .In order to use the iPiano algorithm, we need to calculate the gradient of F and the proximal map w.r.t.G.It is easy to check that where K i is an N × N highly sparse matrix, implemented as 2D convolution of the image u with filter kernel k i , i.e.,, The proximal map w.r.t.G is given as the following minimization problem (II.9) Instead of using a direct solver for (II.9) in terms of the Lambert W function [9], we utilized Newton's method.We found that this scheme has a quite fast convergence (less than 10 iterations).
For the variational model (II.6), and the proximal map w.r.t.G is given by the following pointwise calculation

III. EXPERIMENTAL RESULTS
In this section, we evaluated the performance of our proposed variational models based on images corrupted by synthetic speckle noise and real noisy SAR images.For synthetic experiments, we calculated the common measurements PSNR and SSIM index [23], and compared the results with the stateof-the-art algorithm -SAR-BM3D [20].

A. The influence of data term
We started with conducting a preliminary despeckling test on a commonly used natural image contaminated by multiplicative noise with L = 8, using our proposed variational models (II.8) and (II.6).We carefully tuned the parameter λ to insure that the corresponding variational models achieve the best performance.The despeckling results are shown in Figure 2(b) and (c).
A first impression from this result is: the variational model with an accurate data fitting term, which is exactly derived from the multiplicative noise model can lead to better result than the model with an "inappropriate" data term.But after having a closer look at the despeckling images, we found an interesting phenomena: these two methods possess complementary strengths and failure modes.For example, for the highly  textured region (highlighted by the white rectangle), (II.8) works much better than (II.6); however, for the homogeneous region (highlighted by the red rectangle), (II.6) generates preferable result.Then an intuitive idea arises to integrate these two data terms so as to leverage their advantages.As a result, a new variational model incorporating these two data terms comes out as follows with u = e w .The data fitting term is still convex, and we can utilize iPiano to solve it.The proximal map w.r.t.G is also solved using Newton's method.
We manually tuned the parameters λ 1 and λ 2 , and conducted the same despeckling experiment.The final result is shown in Figure 2(d).One can see that the combined model indeed leads to significant improvement in terms of both PSNR and SSIM index value.From now on, we will use the combined model (III.1)for despeckling experiments.

B. Results on synthetic speckle noise
We first evaluated the performance of our proposed variational model (III.1)for despeckling based on synthetic noisy images.The results are compared to the best published despeckling algorithm -SAR-BM3D [20].First of all, we considered three standard test images widely used in image processing community, see in Figure 3.The despeckling results at three representative numbers of look L = 1, 3, 8 are summarized in Table I.We made the following empirical choice for the parameters λ 1 and λ 2 for different L: if L = 8, λ 1 = 550, λ 2 = 0.02; if L = 3, λ 1 = 310, λ 2 = 0.008; and if L = 1, λ 1 = 160, λ 2 = 0.004.A practical guideline to tune λ 1 and λ 2 , as well as the implementation of the proposed approach can be found at our homepage [2].
From Table I, one can see that these two competing approaches generate very similar results, with almost identical PSNR and SSIM values.We present three despeckling examples obtained by these two algorithms in Figure 4.
As we know, the despeckling performance of one particular method varies greatly for different image contents, in order to  We present the results in Table II.Again, two competing algorithms behave similarly, which implies that our proposed variational model based on the FoE prior is on par with the best published despeckling algorithm -SAR-BM3D.
In experiments, we found that our method introduces blocktype artifacts for the case of low L, e.g., L = 1 in Figure 4(c).The main reason is that our method is a local model, which becomes less effective to infer the underlying structure solely from the local neighborhoods, if the input image is too noisy.In this case, the advantage of non-locality comes through.

C. Results on a real SAR image
In order to demonstrate the effectiveness of our proposed method on SAR image despeckling task, we conducted a despeckling test on a real noisy SAR image, which is taken by the radar system equipped with the JSTARS aircraft [1] (the number of looks L = 5).For this experiment, we set λ 1 = 50, λ 2 = 0.15 for our model.The results of different algorithms are shown in Figure 5.One can see that the tiny structures in the image, especially the region highlighted by the white rectangle, are well-preserved by our approach after despeckling; however, they are smoothed out to different extents by other algorithms.

D. Run time
Efficiency is also an important aspect for real SAR despeckling task.On the server platform of Intel X5675 3.07GHz, for images of size 481×321 exploited in our experiments, a typical run time of the SAR-BM3D algorithm is about 65.4s, which can be reduced to 5.6s, by its improved version -FANS [10], at the expense of slight performance degradation.Our method typically consumes 27s with a pure Matlab code.
However, our model is a local model, which solely contains convolution of linear filters with an image, in contrast to nonlocal models, e.g., SAR-BM3D and its variant FANS.Therefore, our model is well-suited to GPU parallel computation.Our GPU implementation based on a NVIDIA Geforce GTX 680 accelerates the inference procedure significantly; for a despeckling task with L = 8, it typically takes 0.6s for images of size 512 × 512, 0.41s for 481 × 321 and 0.2s for 256 × 256.

IV. CONCLUSION
In this letter, we have proposed a novel variational model for speckle noise reduction, based on an expressive image prior model -FoE model.Our new variational model poses a generally demanding non-convex optimization problem and we have used the recently proposed algorithm -iPiano to solve it efficiently.Numerical results on synthetic images corrupted by speckle noise and a real SAR image demonstrate that the performance of our method is clearly on par with the best published despeckling algorithm -SAR-BM3D.Furthermore, our model comes along with the additional advantage of simplicity and therefore well-suited to GPU programming.The GPU based implementation can conduct despeckling in less than 1s with state-of-the-art performance.

Fig. 1 .
Fig. 1.48 learned filters of size 7 × 7 exploited in our despeckling model.The first number in the bracket is the norm of the filter and the second one is the weight θ i .
Fig. 2. Despeckling results for a widely used natural image corrupted by multiplicative noise with L = 8, using our proposed variational models with different data terms.The recovery quality is measured by PSNR/SSIM index.

Fig. 4 .
Fig.4.Performance comparison to state-of-the-art algorithm -SAR-BM3D[20] for different L. The results are reported by PSNR/SSIM index.make a comprehensive comparison, we conducted additional despeckling experiments over a standard test dataset -68 Berkeley test images identified by Roth and Black[21], which is widely used for AWGN denoising test.All the results were computed per image and then averaged over the test dataset.We present the results in TableII.Again, two competing algorithms behave similarly, which implies that our proposed variational model based on the FoE prior is on par with the best published despeckling algorithm -SAR-BM3D.In experiments, we found that our method introduces blocktype artifacts for the case of low L, e.g., L = 1 in Figure4(c).The main reason is that our method is a local model, which becomes less effective to infer the underlying structure solely from the local neighborhoods, if the input image is too noisy.In this case, the advantage of non-locality comes through.

Fig. 5 .
Fig. 5. Performance comparison of different algorithms on a real SAR image.

DESPECKLING
RESULTS OF SAR-BM3D [20] AND OUR APPROACH.OUR RESULTS ARE MARKED WITH BLUE COLOR (WITH (III.1)).THE RESULTSARE REPORTED WITH PSNR AND SSIM VALUES.

TABLE II DESPECKLING
RESULTS ON 68 BERKELEY TEST IMAGES.THE RESULTS ARE REPORTED WITH AVERAGE PSNR AND SSIM VALUES.