PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Opt Lett. Author manuscript; available in PMC 2010 April 28.
Published in final edited form as:
PMCID: PMC2860949
NIHMSID: NIHMS138502

Speckle attenuation in optical coherence tomography by curvelet shrinkage

Abstract

We describe an algorithm based on shrinkage in the curvelet domain to attenuate speckles in optical coherence tomography (OCT) images. The algorithm exploits the curvelet transform’s sparse representation of edge discontinuities that are common in OCT images and its ability to map signals and noise into different areas in the curvelet domain. The speckle attenuation is controlled by a single parameter that determines the threshold in the curvelet domain. Applying the algorithm to OCT images shows significant improvement of image quality.

Recent years have seen increasing interest in the use of optical coherence tomography (OCT) for a variety of imaging applications, especially in the field of biomedical imaging [1]. Many issues remain to be addressed in order for OCT to be widely used in clinical applications. One of these is the degradation of OCT image quality owing to speckles [2]. Speckles make it particularly challenging to perform further quantitative image analysis such as edge detection, segmentation, and pattern recognition. Consequently, much research has been carried out to develop methods to reduce speckles. The methods developed include spatial compounding [2,3], frequency compounding [2], and digital filtering such as an enhanced Lee filter [4], median filter [4], a symmetric nearest neighbor filter [4], an adaptive Wiener filter [4], an I-divergence regularization [5], as well as filtering in a transform domain such as the wavelet [4,6-9]. Among them, filtering in a transform domain is one of the most promising approaches.

The basic principle of transform domain filtering is to shrink transformed coefficients according to some threshold values. This method is based on the assumption that, in the transform domain, signals are usually more concentrated and represented by some large coefficients, while noise is more spread out and coded in small ones. As a result, by zeroing out small coefficients and keeping large ones, it is possible to effectively attenuate noise without significantly damaging signal components. Two elements are therefore essential for the method: a domain that provides sparse representation of signals and an optimal threshold value. The wavelet is one such domain. It provides a sparse representation for signals containing singularities that satisfy a variety of local smoothness constraints including, for example, piecewise smoothness [10]. Different thresholding strategies in the wavelet domain have been applied for OCT images previously [6,8,11].

OCT images, however, are usually characterized by edge discontinuities rather than individual point discontinuities as many biological tissues have layered structures. The wavelet representation of this kind of feature is not optimally sparse and therefore limits the application of the thresholding strategy. Recently, another multiscale method called curvelet transform was developed [12]. The basis functions of the curvelet transform, so called curvelets, can represent curved edges much more efficiently than wavelets. This concept is illustrated in Fig. 1, which shows that, for the same image, the signal energy is concentrated in a fewer number of curvelet coefficients than of wavelet coefficients. Importantly, the error of reconstructing original signals based on curvelet coefficients as a function of the largest coefficients decays rapidly, which makes denoising via coefficient thresholding quite robust. Curvelets are also directionally sensitive, which offers many degrees of freedom to manipulate signals along certain directions. In contrast, the 2D wavelet transform decomposes data in only three directional subbands. Given the many superior properties of curvelets, we present in this work a method based on the curvelet transform to attenuate speckle noise in OCT images. Our results indicate that image quality is significantly improved.

Fig. 1
(Color online) Energies of the curvelet and wavelet coefficients of the OCT image shown in Fig. 2(a). The coefficients are sorted in descending order of energy from left to right. Curvelets provide a sparser representation of the image, as the signal ...

There are several software implementations of curvelet transform. We use the wrapping method of fast discrete curvelet transform to perform both the forward and inverse transformations [13]. The forward transformation involves several steps, and one of the major steps is to partition the frequency plane into dyadic annuli and trapezoidal regions. The dyadic annuli provide scale information, and trapezoidal regions provide the orientation information of image features. The curvelets are localized in both frequency plane and space domain, and hence the transformed coefficient is a function of the scale j, the orientation l, and the spatial coordinates p and q.

The denoising procedure of our algorithm closely mimics those in the wavelet domain and consists of four steps: logarithm transformation, forward curvelet transform, thresholding, and inverse curvelet transform. As speckles are usually well modeled as multiplicative noises, a logarithm transformation is first applied to convert them to additive noises: log(s) = log(x) + log(z), where s is the acquired data, x is the noise-free signal to be recovered, and z is the speckle. Then a forward curvelet transform is performed on the converted data to produce curvelet coefficients Sj,l,p,q = Xj,l,p,q + Zj,l,p,q, where S, X, and Z are the coefficients for measured data, noise-free signals, and speckle noise, respectively. Following that, a hard threshold Tj,l is applied to each curvelet coefficient Sj,l,p,q so that Sj,l,p,q = Sj,l,p,q when |Sj,l,p,q| > Tj,l, and Sj,l,p,q = 0, otherwise. Then an inverse curvelet transform is performed to Sj,l,p,q to reconstruct the denoised image, and a simple exponential calculation of base 10 would convert the denoised image in the logarithm scale back to the original linear scale.

To obtain the threshold, we use a method called k sigma [14]: Ti,j = 1σ2. Here k is an adjustable parameter, σ1 is the standard deviation of speckles in a background region, and σ2 is the standard deviation of speckles in the curvelet domain at a specific scale j and orientation l. By choosing a background region, which contains only speckles, one can compute the mean value and the standard deviation σ1. Based on those values, σ2 is computed using the Monte Carlo simulation method [14]. We note that, in the computation of σ2, speckles after logarithm operation are assumed to have a Gaussian distribution, but extensive data analysis shows that they generally do not fully satisfy this assumption, so the computed value will be theoretically slightly off. However, since k is adjustable, that effect turns out not to be critical. The value of k can be scale and orientation dependent. It is usually obtained by trial and error.

Figure 2 shows a Fourier domain OCT image of the optical nerve head (a) before and (b) after curvelet denoising. The details of the OCT instrument and image acquisition have been previously described [15]. Much of the speckle in the original image has been reduced, making some image features hidden in the original image more obvious, as shown at the places, for example, indicated by the two white arrows in (b). To better appreciate the performance of the algorithm, Fig. 2(c) shows a cross section of the images at the indicated dashed line in Figs. 2(a) and 2(b). The denoised signal is much cleaner than the original one: the noise fluctuation at the beginning of the original signal is attenuated to close to zero and the noise superimposed on the signal components is also attenuated significantly. Most importantly, all of those are achieved when the image features and the edge sharpness of the original signal are both well maintained, demonstrating the ability of the algorithm to preserve signal while attenuating noise.

Fig. 2
(Color online) Fourier domain OCT images of the optical nerve head (a) before and (b) after curvelet despeckling. For direct comparison, the images are shown on the same color scale. The SNR, CNR, and ENL are all increased, with a high similarity between ...

To quantify the performance of the algorithm, four quality metrics are computed: contrast-to-noise ratio (CNR) [6], which measures the contrast between image features and noise; average equivalent number of looks (ENL) [6], which measure the smoothness of homogeneous areas; peak-signal-to-noise ratio (SNR), defined as SNR = 20 log10{[max(X)]/σ}, where X is the amplitude data and σ is the noise variance; and cross correlation (XCOR), which measures the similarity between the images before and after denoising and is defined as XCOR=m,nsm,nym,n/(m,nsm,n2ym,n2), where s is the intensity data before denoising, y is the intensity data after denoising, and m and n are the indices of the images. Both CNR and ENL are computed using logarithmic scale data, while SNR and XCOR are computed using linear scale data. The value of XCOR is smaller than 1, and the larger is the XCOR, the closer the denoised image is to the original image.

Table 1 lists the results of the quality metrics for different thresholds. With a small increase in the threshold, as expected, the cross correlation decreases, since more data are zeroed, while the values of SNR, CNR, and ENL all increase. For comparison, we have also performed wavelet-based thresholding. The wavelet filtering is based on undecimated discrete wavelet transform and the threshold is chosen to be 4.2 times the noise variance, obtained from the robust median estimator of the highest subband of the transform [16], with the number 4.2 obtained by trial and error to give the best results. Comparing curvelet- with wavelet-based methods, for a similar cross-correlation value, our curvelet method (k = 0.7) further improves the SNR by 7.84 dB.

Table 1
Image Quality Metrics

While the curvelet transform provides an optimally sparse representation of OCT images, it is still unavoidable that image features can be attenuated in the process of thresholding. This can be the result of setting too large a threshold and/or some image noise may have curvelet coefficients that are comparable to those of signals. To minimize the impact, one can take advantage of the direction selectivity of the curvelets. For example, when signals are mainly along certain directions, a small threshold should be used for those directions to preserve signals and a large threshold can be used for other directions to attenuate more noise. This property of the curvelet transform makes it superior to wavelets, as it can distinguish image events in more directions than the 2D wavelet transform, which only decomposes data into three directions: horizontal, vertical, and diagonal. Another advantage is that here the whole process is computationally very efficient, taking only a few seconds in a typical dual-processor laptop computer; by contrast, while many simple wavelet filters take about a few seconds, the spatially adaptive wavelet filter takes about 7 min [6].

To summarize, a curvelet based algorithm is demonstrated that significantly suppresses speckles in OCT images. The image quality is significantly improved: the SNR is increased by 27.7 dB, and the contrast ratio as well as smoothness measure are enhanced, with a small 7.4% loss of similarity. These results highlight the power of the curvelet transform in OCT image analysis and processing. We anticipate that this approach can be applied to other problems in biophotonic and biomedical imaging where noise poses a serious challenge.

Acknowledgments

This work was supported in part by the National Institutes of Health (NIH) (grants LAMMPP41RR01192, EB-00255, and CA 91717).

References

1. Huang D, Swanson EA, Lin CP, Schuman JS, Stinson WG, Chang W, Hee MR, Flotte T, Gregory K, Puliafito CA, Fujimoto JG. Science. 1991;254:1178. [PubMed]
2. Schmitt JM, Xiang SH, Yung KM. J Biomed Opt. 1999;4:95. [PubMed]
3. Schmitt JM. Phys Med Biol. 1997;42:1427. [PubMed]
4. Ozcan A, Bilenca A, Desjardins AE, Bouma BE, Tearney GJ. J Opt Soc Am A. 2007;24:1901. [PMC free article] [PubMed]
5. Marks DL, Ralston TS, Boppart SA. J Opt Soc Am A. 2005;22:2366. [PubMed]
6. Adler DC, Ko TH, Fujimoto JG. Opt Lett. 2004;29:2878. [PubMed]
7. Gargesha M, Jenkins MW, Rollins AM, Wilson DL. Opt Express. 2008;16:12313. [PMC free article] [PubMed]
8. Puvanathasan P, Bizheva K. Opt Express. 2007;15:15747. [PubMed]
9. Xiang SH, Zhou L, Schmitt JM. Proc SPIE. 1997;3196:79.
10. Mallat S. A Wavelet Tour of Signal Processing. Academic; 1998.
11. Choi H, Milner TE, Bovik AC. Proceedings of the 25th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. 2003
12. Candès EJ, Donoho DL. In: Curves and Surface Fitting. Rabut C, Cohen A, Schumaker LL, editors. Vanderbilt U Press; 2000.
13. Candès EJ, Demanet L, Donoho DL, Ying L. Multiscale Model Simul. 2006;5:861.
14. Starck J-L, Candes EJ, Donoho DL. IEEE Trans Image Process. 2002;11:670. [PubMed]
15. Rao B, Yu L, Chiang HK, Zacharias LC, Kurtz RM, Kuppermann BD, Chen Z. J Biomed Opt. 2008;13:040505. [PMC free article] [PubMed]
16. Donoho DL, Johnstone IM. J Am Stat Assoc. 1995;90:1200.