Fast adaptive nonlocal SAR despeckling

524

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. 11, NO. 2, FEBRUARY 2014

Fast Adaptive Nonlocal SAR Despeckling Davide Cozzolino, Sara Parrilli, Giuseppe Scarpa, Giovanni Poggi, and Luisa Verdoliva

Index Terms—Despeckling, nonlocal filtering, synthetic aperture radar (SAR).

statistics of SAR images. In particular, several ad hoc solutions are used to improve the efficiency of block matching: 1) an activity-driven variable-size search area; 2) a probabilistic early termination (PET) algorithm; 3) computation of distances by look-up tables. The new despeckling algorithm is conceptually similar to SAR-BM3D and has similar performance, but is much faster. The next section discusses SAR-BM3D bottlenecks. Sections III and IV describe, respectively, fast collaborative filtering and fast block matching. Then, Section V provides experimental results on both simulated and real-world SAR images.

I. I NTRODUCTION

II. A NALYSIS OF SAR-BM3D

XTRACTING information from synthetic aperture radar (SAR) images is complicated by the presence of speckle which reduces the readability of data, and the reliability of data processing algorithms. Therefore, some form of despeckling is routinely used by both photo-interpreters and computer programs before proceeding with specific tasks. Research on this topic has been very intense, and many techniques have been proposed based, for example, on adaptive linear filtering [1], or wavelet shrinkage [2] and, recently, on nonlocal filtering [3]–[6]. By exploiting image self similarity, nonlocal filtering mimics a true statistical averaging of pixels, thus allowing strong speckle reduction and accurate preservation of features. On the down side, it requires the computation of a large number of block-similarity measures and is therefore computationally demanding. Complexity is a serious problem in the SAR field, since systems of the last generation produce high-resolution images of hundreds of megapixels, and several recent papers focus on fast SAR despeckling techniques [7], [8]. It is certainly an issue for SAR-BM3D [6], a nonlocal technique based on the Block-Matching 3-D (BM3D) algorithm for AWGN denoising [9]. SAR-BM3D combines the nonlocal approach with several powerful tools, such as wavelet transforms, transform-domain shrinkage, Wiener filtering, providing a very good performance at the price of a significant processing time. To speed it up we introduce here several innovations, tailored to the peculiar

SAR-BM3D comprises two passes, in each one the whole image is scanned block-wise, with partial overlapping, and at each block three processing steps are performed: block matching, collaborative filtering and aggregation. Block matching is used to locate, in a given search area, the K blocks which best match the target according to a suitable measure. The selected blocks are then stacked in a 3-D structure and filtered all together (collaboratively), in a transform domain, so as to exploit their similarity. Finally, block estimates are aggregated to produce the output image. The first pass of the algorithm filters the noisy input to produce a basic estimate of the clean image. The second pass filters the noisy input anew using, however, the more reliable statistics computed on the basic estimate to improve the filter performance.

Abstract—Despeckling techniques based on the nonlocal approach provide an excellent performance, but exhibit also a remarkable complexity, unsuited to time-critical applications. In this letter, we propose a fast nonlocal despeckling filter. Starting from the recent SAR-BM3D algorithm, we propose to use a variable-size search area driven by the activity level of each patch, and a probabilistic early termination approach that exploits speckle statistics in order to speed up block matching. Finally, the use of look-up tables helps in further reducing the processing costs. The technique proposed conjugates excellent performance and low complexity, as demonstrated on both simulated and real-world SAR images and on a dedicated SAR despeckling benchmark.

E

Manuscript received February 5, 2013; revised April 10, 2013 and May 28, 2013; accepted June 17, 2013. D. Cozzolino, G. Scarpa, G. Poggi, and L. Verdoliva are with the DIETI, Universitá Federico II di Napoli, 80138 Naples, Italy (e-mail: Davide. [email protected]; [email protected]; [email protected]; [email protected]). S. Parrilli is with the Centro Italiano Ricerche Aerospaziali, 81043 Capua, Italy (e-mail: [email protected]). Digital Object Identifier 10.1109/LGRS.2013.2271650

A. Complexity of Collaborative Filtering This step comprises direct and inverse transform of each 3-D group and Wiener filtering. We indicate with N the block size, S the search area size. Transform complexity depends, obviously, on group size, N K, and type of transform. As an example, using FFT for both the 2-D transform of each block of the group and the 1-D transform along the group calls for N K log2 (N K) multiplications/block in each pass, much less than block matching. To improve the accuracy of subsequent estimates, however, SAR-BM3D uses undecimated wavelet transform (UDWT) in the first-pass, with a much larger complexity, proportional to N KH log2 (N ) log2 (K) multiplications/block, with H the filter length, comparable now to that of block matching. Moreover, since the data volume of the transformed group grows significantly, Wiener filtering itself calls for heavy computation. Therefore, to reduce complexity, it is mandatory to replace UDWT with some faster transforms.

1545-598X © 2013 IEEE

COZZOLINO et al.: FAST ADAPTIVE NONLOCAL SAR DESPECKLING

B. Complexity of Block Matching The selection of K nearest neighbors (NN) requires computing the distance between the target block1 x, and all the S blocks of the search area, y, according to a suitable vector  distance dN (x, y) = N i=1 d(xi , yi ), with d(·, ·) the singleletter distance. Therefore, for each analyzed block, SN singleletter distances must be computed. BM3D uses the squared error measure, which is optimal [10] for AWGN denoising under most criteria of interest. However, to take into account the statistics of speckle, SAR-BM3D uses a more complex distance, derived in [3] under a probabilistic criterion for Nakagami distributed speckle. In the first pass the single-letter distance is     2  1 x y x + y2 (1) + d (x, y) = log = log . (1) 2 y x 2xy In [10] this statistic has been proved to be equivalent to the generalized likelihood ratio, and to provide the best performance in nonlocal filtering. In the second pass, under the same approach, the distance takes into account both the observed values x and y of the pixels under analysis, and the basic estimates u  and v output by the first pass, becoming  2    2 x + y2 ( u − v2 )2 γL , v) = log d(2) (x, y, u + 2xy 2L − 1 u 2 v2 (2) with L the number of looks and the weight γ set experimentally to 1 [6]. By switching the sum of logs with the log of products, and neglecting additions and tests, the cost of block matching amounts to 3SN [5SN ] multiplications/block in the first [second] pass, responsible for much of the complexity. III. FAST C OLLABORATIVE F ILTERING To replace UDWT we tested a number of alternatives choosing eventually the Daubechies-4 wavelet in the spatial domain and the Haar wavelet along the group. Even this optimal choice, however, causes some performance impairment, due to the Wiener filtering, which is optimal when the necessary statistics are perfectly known, but very sensitive to estimation errors. The UDWT provided plenty of data to carry out reliable estimates, which are not available anymore with non-redundant transforms, causing the observed performance impairment. To recover this loss, we replaced Wiener filtering with a more robust wavelet thresholding, which also guarantees good results, provided we can estimate, for each group, the local noise 2 variance σG used to set the optimal threshold. Indeed, if we assume all signals to have constant statistics over a group, which is reasonable given its limited size, statistical moments can be estimated through sample averages as shown in [6]. Wavelet thresholding, however, is well-known to produce some artifacts. With hard thresholding, isolated basis functions stand out in homogeneous areas because of random abovethreshold coefficients; on the contrary, soft thresholding oversmooths important signal features in heterogeneous regions. 1 Blocks are stacked as N -vectors, we use boldface for vectors and capital letters for random variables.

525

We therefore devised an adaptive strategy, using soft thresholding for flat homogeneous blocks and hard thresholding for active ones. As decision statistic we use the ratio of arithmetic to geometric means of the observed block intensities originally proposed in [11] for flat/active classification. Given N , L, and the the desired false alarm probability PF A , that a flat block be classified as active, the decision threshold λT is set using the approximation proposed in [12], valid for N  1. In order to preserve relevant image features the statistic is computed on a larger 16 × 16 block centered on the 8 × 8 target block, so as to reduce the probability of error. With these modifications, collaborative filtering provides satisfactory results and its complexity is negligible with respect to block matching. IV. FAST B LOCK M ATCHING A. Variable-Size Search Area Computing distances only for some selected blocks is the most immediate speed-up strategy. Usually, candidate blocks are classified based on some synthetic features, and those not belonging to the target’s class are discarded. With intense speckle, however, classification becomes quite unreliable, so we use this strategy very conservatively. In particular, we rely on the flat/active classification described before and shrink significantly the search area for flat target blocks, likely to be drawn from an homogeneous area of the image. This straightforward action has motivations that go beyond efficiency, as observed in [13]. Nonlocal techniques, in fact, compare very well with local filtering when dealing with edged/textured blocks, since they are able to find blocks that share an underlying signal similar to the target’s, and exploit such similarities. In homogeneous areas, however, all blocks have similar signal content, and the nearest neighbors turn out to be the blocks that have a noise pattern most similar to the target’s. Therefore, we keep using a large search area for active blocks, but reduce it for flat ones, cutting complexity and improving performance. B. Probabilistic Early Termination Let dN (x, y) be an additive distance measure between N -vectors and assume we have already found K candidate NNs to x, with T the distance of the farthest. We will update the list of NNs, inserting a new vector y, only if dN (x, y) < T , which requires computing and accumulating N single-letter distances. However, we can stop the computation as soon as the partial sum dn (x, y) =

n 

d(xi , yi ),

n n] N

 = Pr d(Φi , Ψi ) < r (5)

, v), we In the second pass, involving distance d(2) (x, y, u proceed in a similar way, but take into account the different behavior of terms depending on (x, y), the original noisy signal, and on ( u, v), the first-pass estimate. In fact, since the latter terms present a higher dynamics, lacking the log, we sum them first, in order to increase the chance of early rejection. Since these terms are almost noise-free, due to the first-pass filtering, the decision threshold remains fixed at t (0) until terms of the first type are taken into account. C. Computing Distances Through Look-up Tables Variable-search area and PET largely reduce the number of letter-distances to compute, but these are still computationally demanding, so we replace computation with memory accesses. This requires a quantization of the whole image by a suitable NQ -level quantizer Q. If qn = Q(x) and qm = Q(y), then   qn qm (1) d (x, y)  log + − log 2 = LUT(n, m). (9) qm qn The NQ × NQ look-up table LUT is computed in advance. In order to keep a good accuracy, the quantization must be reasonably dense, especially for small values of the input since the distance diverges when one of them approaches zero. Therefore, we carry out a uniform quantization in the log domain, in the range dictated by the extreme values, xmin > 0 and xmax , of the input image. Once computed the quantization step as Δ = log2 (xmax /xmin )/(NQ − 1) we set qn = xmin 2nΔ ,

n = 0, . . . , NQ − 1.

(10)

With such a quantizer, however, it is not really necessary to compute and store a matrix of distances since  (11) LUT(n, m) = log2 2(n−m)Δ + 2(m−n)Δ − log 2

i=n+1

where the inequality stems from the fact that u = v maximizes Pr[d(1) (uΦ, vΨ) < t]. The random variables Di = d(1) (Φi , Ψi ) are i.i.d., with finite mean μ and variance σ 2 , hence their sum is well approximated (but for n  N , when early stopping makes no sense) by a Gaussian random variable with mean (N − n)μ and variance (N − n)σ 2 , allowing us to compute the upper bound in (5). In conclusion, knowing T , μ and σ, we can compute a threshold value t (n) such that  (1) d(1) (X, Y) ≥ t (n) ⇒ Pr d (X, Y) < T <  (6)  n N that is, the candidate can be discarded with probability of error less than . Note that (5) is a very conservative bound since, by assuming ui = vi for all components to come, it neglects the contribution, often significant, given to the overall distance by the underlying signals. For√our √ first-pass distance (1) the mean and variance of D = d(1) ( X, Y ), with X,Y i.i.d. Gamma RVs, are [12] μ = E[D] = ψ 0 (2L) − ψ 0 (L) − log 2 ψ 1 (L) σ 2 = VAR[D] = − ψ 1 (2L) 2 with ψ m (x) the polygamma function of order m.

(7) (8)

depends only on the difference of the indexes. Hence, only a vector of size 2NQ − 1 is necessary, thus reducing both memory occupation and memory access time. Similar considerations hold for the look-up table computation of the second-pass distance implemented in a similar way. V. E XPERIMENTAL R ESULTS We analyze performance in terms of both complexity and image quality, also in comparison with some relevant reference algorithms, the Enhanced Lee filter [1] with 5 × 5-pixel window, and the nonlocal PPB with 7 × 7-pixel patches, 21 × 21-pixel search area, and 25 iterations, as suggested in [3]. Complexity is measured by the algorithms’ CPU time on a 3.40 GHz 64 bit desktop computer with 8 GB memory. To measure image quality we use three complementary approaches: statistical simulation based on optical images, the benchmarking tool proposed in [16], [17] based on physical-level simulation of SAR canonical scenes, and real-world SAR images. We use 8 × 8-pixel blocks, a 39 × 39-pixel search area, and 16 [32] NNs in the first [second] pass. The major parameters have been tuned through simulation on different images: λT such that PF A  10−3 for each L,  = 0.18 [0.36] for first [second] pass PET, and NQ = 512. Shrinkage threshold is set

COZZOLINO et al.: FAST ADAPTIVE NONLOCAL SAR DESPECKLING

527

TABLE I P ERFORMANCE IN T ERMS OF PSNR ( D B), SSIM-I NDEX , β -I NDEX , AND CPU T IME ( S ) FOR L ENA (512 × 512) OVER 100 R EALIZATIONS

Fig. 1.

Images used in the experiments: Lena, Rosen-2, Rosen-3.

to 2.7 σG as in [9]. All parameters have been set once and for all, obtaining always satisfactory results in the experiments, so the user can forget about them and keep the default values. In Table I we show objective performance indices measured on 1-look, 4-look, and 16-look SAR images simulated based on the Lena image of Fig. 1 (left). We report PSNR and structural similarity (SSIM) index, which measure how similar the despeckled image is to the clean reference, the β index [18], which measures preservation of edges in the filtered image, and the CPU time. We use Monte Carlo simulation on 100 realizations, as suggested in [19], and report both mean value and standard deviation (smaller font). FANS turns out to be over 10 times faster than SAR-BM3D, comparable to it in terms of both PSNR and SSIM, and superior in terms of β-index, due to the classification step. Notice also that the first pass of FANS guarantees a very good quality/complexity tradeoff, with negligible performance losses and a CPU time of just 2–3 s. Moreover, FANS largely outperforms PPB and Enhanced Lee although the latter is much faster than all other algorithms, and preferable if time is the dominant issue. These considerations hold both for the single-look and multilook cases, but the performance gaps narrow somewhat when starting from a better image, as for L = 16. PSNR figures for simulated SAR images help understanding major trends, but must be complemented by visual inspection and, possibly, numerical measures specific of SAR imagery. Therefore, in Fig. 2 we show (top row) the output of the various despeckling algorithms on the Rosen-2 image of Fig. 1 (center), a TerraSAR-X image (Infoterra GmbH) taken over Rosenheim, in Germany. PPB provides the best speckle rejection followed closely by FANS, both by visual inspection and in terms of equivalent number of looks (ENL), computed in the whiteframed area and reported in the caption. On the other hand, PPB, and in lesser measure Enhanced Lee, oversmooth the

Fig. 2. Results for Rosen-2. (Top) Images filtered by Enhanced Lee (ENL = 9.9), PPB (ENL = 73.0), SAR-BM3D (ENL = 7.2), FANS (ENL = 50.6). Middle: ratio images. (Bottom) Edge maps provided by the ratio edge detector.

image, leaving traces of signal structures in the ratio images (middle row), absent for SAR-BM3D and FANS. The FANS image also provides the best edge map (bottom row), obtained by the ratio edge detector [20], catching most actual edges (unlike PPB) with little or no false alarms (unlike Enhanced Lee and SAR-BM3D). Similar considerations hold for the Rosen-3 image of in Fig. 1 (right), a section of which is shown enlarged in Fig. 3 together with the filter outputs to facilitate visual inspection. Notice that, by reinforcing pixel patterns already present in the noisy image, nonlocal techniques find structures in mostly homogeneous areas which might correspond to signal details and textures but also to random speckle patterns, something impossible to decide lacking a ground truth. To gain better insight about this and other phenomena relevant for SAR images, Table II reports the objective indicators proposed in [17] computed on the single-look SAR images of Fig. 4 generated by physical-level simulation. Each indicator is computed w.r.t. a 512-look reference image of the same scene. ENL accounts for speckle suppression in homogeneous regions, the coefficient of variation (Cx ) measures texture preservation, Despeckling Gain (DG) measures the SNR improvement granted by despeckling, Edge smearing (ES) accounts for edge profile degradation, Figure of Merit (FOM) for edge

528

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. 11, NO. 2, FEBRUARY 2014

between despeckling and classification. All techniques have good radiometric accuracy except PPB which has low CBG . In general, FANS seems to keep the good performance already exhibited by SAR-BM3D with a much smaller running time. Given the promising results observed in this work, future research will focus on classification-based nonlocal despeckling and on the extension to the case of correlated speckle. Software is available at http://www.grip.unina.it for full reproducibility of results. R EFERENCES

Fig. 3. 512 × 100 strip of Rosen-3: original, Enhanced Lee (ENL = 12.3), PPB (ENL = 79.8), SAR-BM3D (ENL = 21.8), FANS (ENL = 55.3). TABLE II P ERFORMANCE C OMPARISON ON THE D ESPECKLING B ENCHMARK

Fig. 4. Images of the benchmark: Homogeneous, DEM, Squares, Corner.

detectability after despeckling, finally, Contrast to Background (CBG ) measures radiometric fidelity. Nonlocal techniques guarantee by far the strongest speckle suppression (large ENL), with FANS on top, thanks to the stronger filtering carried out in regions classified as flat. SARBM3D and FANS guarantee a very good preservation of features in textured areas (Cx  2.40, large DG). SAR-BM3D has also an excellent performance on edges, both in terms of smearing (but the best is not to filter at all) and edge detectability (large FOM), while FANS is penalized by the classifier that misses a whole edge, declared as a flat region and hence strongly smoothed, ringing a bell on the delicate interplay

[1] A. Lopes, R. Touzi, and E. Nezry, “Adaptive speckle filters and scene heterogeneity,” IEEE Trans. Geosci. Remote Sens., vol. 28, no. 6, pp. 992– 1000, Nov. 1990. [2] T. Bianchi, F. Argenti, and L. Alparone, “Segmentation-based MAP despeckling of SAR Images in the undecimated wavelet domain,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 9, pp. 2728–2742, Sep. 2008. [3] 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. [4] S. Parrilli, M. Poderico, C. V. Angelino, G. Scarpa, and L. Verdoliva, “A nonlocal approach for SAR image denoising,” in Proc. IGARSS, Jul. 2010, pp. 726–729. [5] H. Zhong, Y. Li, and L. C. Jiao, “SAR image despeckling using bayesian nonlocal means filter with sigma preselection,” IEEE Geosci. Remote Sens. Lett., vol. 8, no. 4, pp. 809–813, Jul. 2011. [6] 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 Sens., vol. 50, no. 2, pp. 606–616, Feb. 2012. [7] F. Argenti, T. Bianchi, A. Lapini, and L. Alparone, “Fast MAP despeckling based on Laplacian-Gaussian modeling of wavelet coefficients,” IEEE Geosci. Remote Sens. Lett., vol. 9, no. 1, pp. 13–17, Jan. 2012. [8] B. Xue, Y. Huang, J. Yang, L. Shi, Y. Zhan, and X. Cao, “Fast nonlocal remote sensing image denoising using cosine integral images,” IEEE Geosci. Remote Sens. Lett., vol. 10, no. 6, pp. 1309–1313, Nov. 2013. [9] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3-D transform-domain collaborative filtering,” IEEE Trans. Image Process., vol. 16, no. 8, pp. 2080–2095, Aug. 2007. [10] C. A. Deledalle, L. Denis, and F. Tupin, “How to compare noisy patches? Patch similarity beyond Gaussian noise,” Int. J. Comput. Vis., vol. 99, no. 1, pp. 86–102, Aug. 2012. [11] M. Beauchemin and K. P. B. Thomson, “The ratio of arithmetic to geometric mean: A first-order statistical test for multilook SAR image homogeneity,” IEEE Trans. Geosci. Remote Sens., vol. 34, no. 2, pp. 604– 606, Mar. 1996. [12] J. Cheng, N. Wang, and C. Tellambura, “Probability density function of logarithmic ratio of arithmetic mean to geometric mean for Nakagami-m fading power,” in Proc. 25th Biennal Symp. Commun., 2010, pp. 348–351. [13] D. Gragnaniello, G. Poggi, and L. Verdoliva, “Classification-based nonlocal SAR despeckling,” in Proc. Tyrrhenian Workshop Adv. Radar Remote Sens., Sep. 2012, pp. 121–125. [14] K. Lengwehasatit and A. Ortega, “Probabilistic partial-distance fast matching algorithms for motion estimation,” IEEE Trans. Circuits Syst. Video Technol., vol. 11, no. 2, pp. 139–152, Feb. 2001. [15] R. Vignesh, B. T. Oh, and C. C. J. Kuo, “Fast non local means (NLM) computation with probabilistic early termination,” IEEE Signal Process. Lett., vol. 17, no. 3, pp. 277–280, Mar. 2010. [16] G. Di Martino, M. Poderico, G. Poggi, D. Riccio, and L. Verdoliva, “SAR image simulation for the assessment of despeckling techniques,” in Proc. IGARSS, Jul. 2012, pp. 1797–1800. [17] G. Di Martino, M. Poderico, G. Poggi, D. Riccio, and L. Verdoliva, “Benchmarking framework for SAR despeckling,” IEEE Trans. Geosci. Remote Sens., 2013, to be published. [Online]. Available: http:// ieeexplore.ieee.org [18] F. Sattar, L. Floreby, G. Salomonsson, and B. Lovstrom, “Image enhancement based on a nonlinear multiscale method,” IEEE Trans. Image Process., vol. 6, no. 6, pp. 888–895, Jun. 1997. [19] E. Moschetti, M. G. Palacio, M. Picco, O. H. Bustos, and A. C. Frery, “On the use of Lee’s protocol for speckle-reducing techniques,” Latin Amer. Appl. Res., vol. 36, pp. 115–121, 2006. [20] R. Touzi, A. Lopes, and P. Bousquet, “A statistical and geometrical edge detector for SAR images,” IEEE Trans. Geosci. Remote Sens., vol. 26, no. 6, pp. 764–773, Nov. 1988.