Home | About | Journals | Submit | Contact Us | Français |

**|**Opt Express**|**PMC2898712

Formats

Article sections

Authors

Related links

Optics Express

Opt Express. 2010 January 18; 18(2): 1024–1032.

Published online 2010 January 7. doi: 10.1364/OE.18.001024

PMCID: PMC2898712

NIHMSID: NIHMS201511

Zhongping Jian,^{
1
,}^{*} Lingfeng Yu,^{
1
} Bin Rao,^{
1
} Bruce J. Tromberg,^{
1
} and Zhongping Chen^{
1
,}^{
2
}

Received 2009 October 21; Revised 2009 December 14; Accepted 2009 December 18.

Copyright ©2010 Optical Society of America

This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-No Derivative Works 3.0 Unported License, which permits download and redistribution, provided that the original work is properly cited. This license restricts the article from being modified or used commercially.

This article has been cited by other articles in PMC.

Optical coherence tomography is an emerging non-invasive technology that provides high resolution, cross-sectional tomographic images of internal structures of specimens. OCT images, however, are usually degraded by significant speckle noise. Here we introduce to our knowledge the first 3D approach to attenuating speckle noise in OCT images. Unlike 2D approaches which only consider information in individual images, 3D processing, by analyzing all images in a volume simultaneously, has the advantage of also taking the information between images into account. This, coupled with the curvelet transform’s nearly optimal sparse representation of curved edges that are common in OCT images, provides a simple yet powerful platform for speckle attenuation. We show the approach suppresses a significant amount of speckle noise, while in the mean time preserves and thus reveals many subtle features that could get attenuated in other approaches.

Optical coherence tomography (OCT) has been undergoing rapid development since its introduction in the early 1990s [1]. It provides high resolution, cross-sectional tomographic images of internal structures of specimens, and therefore gains a wide variety of application in the field of biomedical imaging. Compared with other medical imaging modalities, 3D OCT has advantages in that it is non-invasive and it can acquire and display volume information in real time. However, due to its coherent detection nature, OCT images are accompanied with a significant amount of speckle noise, which not only limits the contrast and signal-to-noise ratio of images, but also obscures fine image features.

Various methods have been developed to minimize the effect of speckle noise. Those methods can generally be classified into two categories: the first one performs noise attenuation by acquiring extra data, such as using spatial compounding and frequency compounding [2, 3]. While effective, this method generally requires extra effort to acquire data and cannot process images from standard OCT systems, and is therefore less preferred than the second category, which uses digital signal processing techniques to process images acquired with standard OCT systems. Different digital signal processing algorithms have been proposed, including for example enhanced Lee filter [4], median filter [4], symmetric nearest neighbor filter [4], adaptive Wiener filter [4], I-divergence regularization [5], as well as filtering in a transform domain such as the wavelet [4, 6–9]. Recently we described a speckle suppression algorithm in a transform domain called curvelets [10]. There we showed the curvelet representation of OCT images is very efficient, and with that, we significantly improved qualities of OCT images in the respects of signal to noise ratio, contrast to noise ratio, and so on.

In almost all those algorithms, however, speckle reduction is performed on each image in a volume individually, and then all despeckled images are put together to form a volume. This process treats images as if they are independent from each other and therefore no relationship among different images is utilized, which is a waste of information provided by 3D OCT data. As many biological structures have layered structures not just in 2D, but also in 3D, and speckle noise is still random in 3D, we would expect that a despeckling algorithm based on 3D processing will be more powerful in attenuating noise and preserving features, especially those fine features across different images.

There are a number of ways to do 3D processing, such as extending those two-dimensional filters mentioned above to three dimensional, or performing a 3D transform followed by processing in the transformed domain. The 3D transform can be, for example, the 3D wavelet transform, the 3D curvelet transform, or a hybrid one, such as a 2D curvelet transform of individual images followed by a one-dimensional wavelet transform along the perpendicular direction. Given the many superior properties of the curvelet transform, here we extend our earlier work of 2D curvelets to 3D, by performing the speckle attenuation in the 3D curvelet domain. We will first introduce some background information of the curvelet transform and its properties, then describe our algorithm in detail, and finally present the curvelet despeckling results tested on three-dimensional Fourier domain OCT images.

The curvelet transform is a recently developed multiscale mathematical transform with strong directional characters [11–13]. It is designed to efficiently represent edges and other singularities along curves. The transform decomposes signals using a linear and weighted combination of basis functions called curvelets, in a similar way as the wavelet transform decomposes signals as a summation of wavelets. Briefly, the curvelet transform is a higher-dimensional extension of the wavelet transform. While the wavelet transform provides structured and sparse representations of signals containing singularities that satisfy a variety of local smoothness constraints, including for example piecewise smoothness, they are unable to capitalize in a similar effective fashion on signals of two and more dimensions. The curvelet transform can measure information of an object at specified scales and locations and only along specified orientations. To achieve that, curvelets first partition the frequency plane into dyadic coronae, and (unlike wavelets) then subpartition the coronae into angular wedges [11]. Curvelets have time-frequency localization properties of wavelets, yet (unlike wavelets) also show a high degree of directionality and anisotropy. The curvelet transform is particularly suitable for noise attenuation, as it maps signals and noise into different areas in the curvelet domain, the signal’s energy is concentrated in a limited number of curvelet coefficients, and the reconstruction error decays rapidly as a function of the largest curvelet coefficients.

The two-dimensional curvelets are, roughly speaking, 2D extensions of wavelets. They are localized in two variables and their Fourier duals (e.g., x-y and fx-fy), and are uniquely decided by four parameters: scale, orientation, and two translation parameters (x,y)). There are several software implementations of the curvelet transform, and the one often used is the wrapping method of Fast Discrete Curvelet Transform (FDCT) [11]. The left of Fig. 1 shows a curvelet partitioning of fx-fy plane, where there are 6 scales, represented by the squares, and going from the inner to outer, the scale is j =1,2,3,…6. Each scale is further partitioned into a number of orientations, and the number doubles every other scale starting from the second (coarsest) scale. That is, going from the inner to the outer, the number of orientations is l=1, n, 2n, 2n, 4n, 4n… where n is the number of orientation at the second (coarsest) scale. This way, the directional selectivity increases for finer scales. The right side of Fig. 1 shows two example curvelets at the specific scales and orientations denoted by A and B in the partition diagram, respectively. Each curvelet oscillates in one direction, and varies more smoothly in the others. The oscillations in different curvelets occupy different frequency bands. Each curvelet is spatially localized, as its amplitude decays rapidly to zero outside of certain region. The directional selectivity of curvelets can be observed, for example, (A) is mainly along the horizontal direction while (B) is in another direction. This property can be utilized to selectively attenuate/preserve image features along certain directions.

Left: A schematic of the curvelet partitioning of fx-fy domain. The number of scales is 6, and the number of orientations at the second scale is 8. Right: two example curvelets, shown for the scale and orientation A and B, respectively. The curvelet A **...**

The three-dimensional (3D) transform is very similar to the two-dimensional transform, except that each scale is defined by concentric cubes in the fx-fy-fz domain, and the division into orientations is performed by dividing the square faces of the cubes into sub-squares. Like in 2D transform, the number of orientations is specified for the second (coarsest) scale, which then determines the number of sub-squares in each direction. For example, a value of 8 orientations would lead to 64 sub-squares on each face. And since there are 6 faces to each cube, there would be a total of 384 orientations at that scale. The number of orientations doubles every other scale for finer scales, the same way as in the 2D transform.

The curvelet-based despeckling algorithm consists of the following steps:

- I.A preprocessing step is first applied to the acquired data to compensate for the motion of the target during the data acquisition process. 3D OCT data is acquired image by image, not obtained at one single shot. Although the scanning time can be quite short, the target can still move during that short time. The motion can have significant impact on acquired images. For example, it can seriously distort the shape of the target, making the edge detection and other image analysis especially challenging. It can also make some continuous features across images not continuous any more, which would make the corresponding 3D curvelet transform coefficients smaller than they should. Those smaller coefficients can then be attenuated during the despeckled process, which in turn, can lead to the loss of image features. To minimize the impact of the motion, those features are first aligned in all directions. For example, for our acquired retina images, data are preprocessed based on the idea that Retinal Pigment Epithelium (RPE) in neighboring images should be continuous, and blood vessels in fundus image should have minimal abrupt changes. The aligned data is then further processed in the next steps.
- II.Take a logarithm operation of the aligned data. This is to convert the multiplicative noise into additive noise, as it is well known that speckles can be well modeled as multiplicative noise. That is, log(s) = log(x) + log(z), where s is the measured data, x is the noise free signals to be recovered, and z is the speckle noise.
- III.Take the 3D forward curvelet transform of the data to produce the curvelet coefficients. The curvelet transform is a linear process, so the additive noise is still additive after the transform: S
= X_{j,l,p}+ Z_{j,l,p}, where S_{j,l,p}, X_{j,l,p}, and Z_{j,l,p}are the coefficients for measured data, speckle-free signals, and speckle noise, respectively;_{j,l,p}*j*,*l*and*p*are parameters used for the curvelet transform,*j*is the scale,*l*is the orientation, and*p*is the spatial coordinates. - IV.Selectively attenuate the obtained curvelet coefficients. A hard threshold T
is applied to each curvelet coefficients S_{j,l}so that ${\overline{S}}_{j,l,p}$ = S_{j,l,p},when |S_{j,l,p}|>T_{j,l,p}, and ${\overline{S}}_{j,l,p}$ =0 when |S_{j,l}|≤T_{j,l,p}._{j,l} - V.Take the inverse 3D curvelet transform of the attenuated curvelet coefficients to reconstruct despeckled data. The obtained data is in logarithm scale, so an exponential calculation of base 10 is applied to convert the despeckled data back to the original linear scale when needed.

In the process, one of the most important steps is the selection of the threshold T_{j}_{,}
* _{l}*, which determines to a large extent the performance of the algorithm. Here we use a simple yet powerful strategy called k-sigma method to set the threshold [14], in which T

To quantify the performance of the algorithm, we compute five quality metrics [6]: contrast-to-noise ratio (CNR), which measures the contrast between image features and noise, and defined to be
$CNR=10\mathrm{log}[({\mu}_{s}-{\mu}_{b})/\sqrt{{\sigma}_{s}^{2}+{\sigma}_{b}^{2}}]$; equivalent number of looks (ENL), which measure the smoothness of areas that should be homogeneous but are corrupted by speckle noise, and defined to be
$ENL={\mu}_{s}^{2}/{\sigma}_{s}^{2}$, where μ_{s} and σ_{s} are the mean and standard deviation of a signal area, and μ_{b} and σ_{b} are the mean and standard deviation of a background noise area, respectively; peak signal to noise ratio (SNR), defined as
$SNR=20\mathrm{log}[\mathrm{max}(x)/\sigma ]$, where x is the amplitude data and σ is the noise variance of the background noise area; crosscorrelation (XCOR), which measures the similarity between the images before and after denoising, and is defined as
$XCOR={\displaystyle \sum _{m,n}{s}_{m,n}{y}_{m,n}}/\sqrt{{\displaystyle \sum _{m,n}{s}_{m,n}^{2}}{\displaystyle \sum _{m,n}{y}_{m,n}^{2}}}$, where s is the intensity data before denoising, y is the intensity data after denoising, and m and n are the indexes of the images; and FWHM, the full width at half maximum, which measures the image sharpness. Both CNR and ENL are computed using log scale data, and are averaged over many areas. SNR and XCOR are computed using linear scale data. The value of XCOR is smaller than 1, and the larger XCOR is, the closer the denoised image is to the original image.

The image data is acquired by a Fourier domain OCT system [15]. The low-coherence light source has a center wavelength of 890nm and an FWHM bandwidth of 150nm. A broadband optical isolator was used to prevent optical feedback before light enters a 2 by 2 broadband fiber- coupler-based interferometer. Light at the reference arm was focused onto a reference mirror. The sample arm was modified from the patient module of a Zeiss Stratus OCT instrument. The detection arm was connected to a high performance spectrometer, which makes the system bench-top sensitivity of 100 dB with 650 μw light out of the sample-arm fiber and 50 μs CCD integration time. A 9 dB of SNR roll-off from 0 mm imaging depth to 2 mm depth was observed. The system speed was set to be 16.7 K A-lines/s, with its CCD A-line integration time being 50 μs and the line period being 60 μs. With the system, we acquired a 3D volume of human retina, with a lateral resolution of 7.8 μm and axial resolution of 4 μm.

We applied our algorithm to the acquired data. Figure 2 shows experimentally acquired cross-sectional images of human retina in three perpendicular planes: (a) x-y (B-scan), (b) x-z, and (c) y-z, respectively, where x is in the depth direction, y is perpendicular to x and is in the B-scan plane, z is perpendicular to both x and y directions and is the third dimension. Figure 3 shows the same images after being denoised by the 3D algorithm. For direct comparison, the images in two figures are shown on the same color scale and no pixel thresholding is applied. The background region, where there are no distinct image features, is the upper region in (a) and (b), as well as the middle and right noise region of (c). In obtaining the despeckled results, we have tested a number of combinations of parameters to perform the 3D curvelet transform, and the used values are: the number of scales is 3, and the number of orientations at the second coarsest scale is 16. A common threshold k=0.42 is used at all scales and orientations.

(color online) acquired cross-sectional retina images before denoising at different planes: (a) x-y plane (B-scan plane), (b) x-z plane along the vertical solid white line in (a), and (c) the cross-section image in the y-z plane along the horizontal solid **...**

Much of the noise in the images has been reduced, which is most obvious in the background regions. To have a better comparison, Fig. 4 shows a one-dimensional cross-section of the image at the indicated white dotted line in Fig. 2, from images (a), (b) and (c), respectively. The despeckled signals are much cleaner than the original ones: the strong noise fluctuation in the original signals is attenuated, not only at the places where only noise resides, but also in other parts where noise is superimposed on the signals. And the attenuation of the speckle noise is achieved when the edge sharpness and image features of the original signal are both well preserved, demonstrating the ability of the algorithm in preserving signals while attenuating noise.

The despeckling process makes some features of the object more obvious. For example, it is challenging, from the original signals (blue dotted lines) in Fig. 4 (a) and (b), to judge where the layered structure of the retina is, but it is much easier to do so from the denoised signals (red solid lines): the denoised signals, with the noised fluctuation removed, provide more distinct peaks and therefore the locations of the layered structure. This is especially useful for further automatic image analysis, as the less the ambiguity there is, the more accurate the results will be.

Often times some image features are not distinct in a single image, but they are continuous across many neighboring images. In 2D despeckling, those weak image features tend to be attenuated with the speckle noise, as their amplitude and therefore transformed coefficients are close to those of noise. They, however, can be better preserved in 3D processing, as a three-dimensional curvelet transform would give relatively large coefficients for those continuous features across images than for randomly appeared speckle noise. An example is the two yellow features indicated by two black arrows in Fig. 3(c). They are easily discernible in the despeckled data, but can be barely observed from the image before despeckling. Another example is the photoreceptor inner and outer segment junction (IS/OS) indicated by the black arrow in Fig. 3(b), which is nicely continuous across images (along the direction of z) and distinct from its neighboring features, but the same feature is less distinct in the image before despeckling.

To see this effect more clearly, Fig. 5 shows the same images in Fig. 2 denoised by 2D despeckling algorithm, where the threshold in the 2D algorithm [10] is chosen so that the crosscorrelation between Fig. 5(a) and Fig. 2(a) is the same as the crosscorrelation between Fig. 3(a) and Fig. 2(a). Not only the features indicated by the black arrows are more distinct and continuous in the 3D despeckling results, but also the layers of tissue where the white arrows reside in Fig. 5 are more preserved in the 3D results. The reason for this preservation difference is that these layers of tissue have the signals that are comparable to those of noise, as a result, when only a single image is despeckled in 2D despeckling, their transformed coefficients are close to those of noise and therefore can be attenuated easily. On the other hand, in 3D despeckling, because of the continuous features, the transformed coefficients are larger than those of noise and therefore are preserved better.

The improvement of the image quality is also reflected in quality metric numbers. Table 1 lists the results of the quality metrics for three different thresholds of 3D method and one threshold for 2D method, and Fig. 6 shows the trend of SNR and crosscorrelation XCOR for more 3D thresholds. Comparing the original signal to the despeckled signal at threshold k=0.5, the signal to noise ratio is significantly increased by 32.59dB, the contrast to noise ratio is increased by 3.17dB, the sharpness, calculated based on the FWHM of the photoreceptor inner and outer segment junction (IS/OS) from the Fig. 4(b), is improved by 1.55 times, and the smooth region is more smooth after despeckling, with the equivalent number of looks increased by more than 3 times. All those are achieved when the crosscorrelation is 0.914. Although the number 0.914 might not seem ideal, as we have seen from Fig. 2, ,3,3, and and5,5, the sharpness and features of the original images are still well preserved in the despeckled images.

SNR and Crosscorrelation as a function of different threshold k in the 3D despeckling algorithm. The algorithm improves the most SNR of 32.59 dB at k=0.5, and the crosscorrelation between the original image and the despeckled image is 0.914. The crosscorrelation **...**

With the increase of threshold k, as expected, SNR, CNR and ENL all increase while XCOR decreases. However, the signal to noise ratio does not always increase, instead it reaches the maximum of 65.54dB at k=0.5, then begins to drop to ~60dB at k=1.0, as shown in Fig. 6; the crosscorrelation decreases initially at small k values, and then it does not change significantly for k between 0.6 and 1.0. This is a very interesting phenomenon, as we would think the crosscorrelation should decrease all the time with increasing thresholds. It, however, is explainable and from another perspective, shows the advantage of processing in the curvelet domain; that is, curvelets provide a sparse representation so that most signal energy is concentrated in a limited number of curvelet coefficients, and the curvelet reconstruction error decays rapidly as a function of maximum curvelet coefficients. As a result, although increasing k leads to zeroing of more curvelet coefficients, so long as the threshold is not large enough to attack those limited number of major curvelet coefficients, an almost the same data can still be reconstructed and therefore the crosscorrelation does not vary much. Of course, increasing the threshold further would, eventually, lead to the loss of image features and the decrease of the crosscorrelation. The existence of the relatively large comfortable range of threshold showcases yet another advantage of the algorithm over many other methods such as those based on the wavelet transform, whose performance is usually more sensitive to the selection of the threshold.

In this manuscript we presented a three-dimensional despeckling algorithm based on the curvelet transform. We demonstrated that the algorithm can significantly suppress speckle noise, and at the same time can preserve image features well. We showed the 3D processing uses inter-information among images, and preserves weak features across images. These encouraging results highlight the power of the curvelet processing in biophotonic and biomedical image processing. Of course, those benefits come at a cost. The current 3D curvelet transform [11] requires more memory and is computationally more demanding, when comparing to other 2D methods such as 2D curvelet and wavelet transforms. This limits the size of data volume that can be processed. For example, the data size can be processed by our dual-processor laptop computer (2.0GHz CPU, 2GB memory) is 256×512×64, and it takes about 20 minutes to do 3D processing, while 2D curvelet processing only takes 2 minutes. However, this size limit is mainly imposed by the Matlab, and we believe the data size can be significantly larger and the time can be less when using the C++ implementation of the curvelet transform. Another issue with the 3D processing is that, although the 3D algorithm performs better in preserving weak features across images than 2D curvelet algorithm, the quality of individual despeckled images seem slightly worse than those processed by the 2D algorithm [10]. This may partially be attributed to the motion of the eye during the data acquisition process, whose artifacts are not completely removed in the preprocessing step. The motion degrades the correlation of image features between different frames, and therefore makes 3D curvelet coefficients smaller than they should. Those smaller coefficients might then be attenuated as if they were from random noise. Better compensation scheme in the preprocessing step, less motion during acquisition process, and future higher speed OCT acquisition system may significantly improve this result.

The results presented in the manuscript still have room for improvement. We believe that more advanced thresholding technique, such as those implemented for wavelets [6, 8, 16], can also be applied for curvelets after proper adaptation. And we have not utilized all the advantages provided by the curvets, such as multi directionality, which can be used to preserve or attenuate image features along certain directions. Further, despeckling is just the beginning of imaging analysis. Given many superior properties of curvelets, we anticipate the curvelet domain provides a promising platform for quantitative image analysis, such as pattern recognition and so on that could play a significant role in future analysis of biophotonic and biomedical images.

This work was supported in part by the National Institutes of Health ( LAMMP-P41RR01192, EB00293, CA 91717).

1. Huang D., Swanson E. A., Lin C. P., Schuman J. S., Stinson W. G., Chang W., Hee M. R., Flotte T., Gregory K., Puliafito C. A., Fujimoto J. G., “Optical Coherence Tomography,” Science 254(5035), 1178–1181 (1991).10.1126/science.1957169 [PubMed] [Cross Ref]

2. Schmitt J. M., “Array detection for speckle reduction in optical coherence microscopy,” Phys. Med. Biol. 42(7), 1427–1439 (1997).10.1088/0031-9155/42/7/015 [PubMed] [Cross Ref]

3. Schmitt J. M., Xiang S. H., Yung K. M., “Speckle in Optical Coherence Tomography,” J. Biomed. Opt. 4(1), 95 (1999).10.1117/1.429925 [PubMed] [Cross Ref]

4. Ozcan A., Bilenca A., Desjardins A. E., Bouma B. E., Tearney G. J., “Speckle reduction in optical coherence tomography images using digital filtering,” J. Opt. Soc. Am. A 24(7), 1901 (2007).10.1364/JOSAA.24.001901 [PMC free article] [PubMed] [Cross Ref]

5. Marks D. L., Ralston T. S., Boppart S. A., “Speckle reduction by I-divergence regularization in optical coherence tomography,” J. Opt. Soc. Am. A 22(11), 2366 (2005).10.1364/JOSAA.22.002366 [PubMed] [Cross Ref]

6. Adler D. C., Ko T. H., Fujimoto J. G., “Speckle reduction in optical coherence tomography images by use of a spatially adaptive wavelet filter,” Opt. Lett. 29(24), 2878–2880 (2004).10.1364/OL.29.002878 [PubMed] [Cross Ref]

7. Gargesha M., Jenkins M. W., Rollins A. M., Wilson D. L., “Denoising and 4D visualization of OCT images,” Opt. Express 16(16), 12313–12333 (2008).10.1364/OE.16.012313 [PMC free article] [PubMed] [Cross Ref]

8. Puvanathasan P., Bizheva K., “Speckle noise reduction algorithm for optical coherence tomography based on interval type II fuzzy set,” Opt. Express 15(24), 15747–15758 (2007).10.1364/OE.15.015747 [PubMed] [Cross Ref]

9. Xiang S. H., Zhou L., Schmitt J. M., “Speckle Noise Reduction for Optical Coherence Tomography,” Proc. SPIE 3196, 79 (1997).10.1117/12.297921 [Cross Ref]

10. Jian Z., Yu Z., Yu L., Rao B., Chen Z., Tromberg B. J., “Speckle Attenuation by Curvelet Shrinkage in Optical Coherence Tomography,” Opt. Lett. 34, 1516 (2009).10.1364/OL.34.001516 [PMC free article] [PubMed] [Cross Ref]

11. Candès E. J., Demanet L., Donoho D. L., Ying L., “Fast Discrete Curvelet Transforms,” SIAM Multiscale Model. Simul. 5(3), 861 (2006).10.1137/05064182X [Cross Ref]

12. E. J. Candès, and D. L. Donoho, “Curvelets–a surprisingly effective nonadaptive representation for objects with edges,” in Curves and Surface Fitting, C. Rabut, A. Cohen, and L. L. Schumaker, eds. (Vanderbilt University Press, Nashville, TN., 2000).

13. Candès E. J., Donoho D. L., “New tight frames of curvelets and optimal representations of objects with piecewise C^{2} singularities,” Commun. Pure Appl. Math. 57, 219 (2003).10.1002/cpa.10116 [Cross Ref]

14. Starck J.-L., Candès E. J., Donoho D. L., “The Curvelet Transform for Image Denoising,” IEEE Trans. Image Process. 11(6), 670–684 (2002).10.1109/TIP.2002.1014998 [PubMed] [Cross Ref]

15. Rao B., Yu L., Chiang H. K., Zacharias L. C., Kurtz R. M., Kuppermann B. D., Chen Z., “Imaging pulsatile retinal blood flow in human eye,” J. Biomed. Opt. 13(4), 040505 (2008).10.1117/1.2967986 [PMC free article] [PubMed] [Cross Ref]

Articles from Optics Express are provided here courtesy of **Optical Society of America**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |