|Home | About | Journals | Submit | Contact Us | Français|
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 . 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 , median filter , symmetric nearest neighbor filter , adaptive Wiener filter , I-divergence regularization , 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 . 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 . 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) . 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.
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:
In the process, one of the most important steps is the selection of the threshold Tj , 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 , in which Tj, l=k×σ1×σ2, where k is an adjustable parameter, σ1 is the standard deviation of noise from a background region in the image data, and σ2 is the standard deviation of noise in the curvelet domain at a specific scale j and orientation l. By choosing a background region that does not have image features, one can directly compute the mean value and the standard deviation σ1. σ2, on the other hand, cannot be directly calculated from the forward curvelet transformed data, because the transformed data contain coefficients of not only noises, but also of image features, and it is not easy to separate them in the curvelet domain. One easier way to get σ2 is to simulate the noise data from the mean value and σ1, by assuming the noise has Gaussian distribution. Then the simulated data is transformed into the curvelet domain. The standard deviation σ2 at a specific scale and orientation can then be directly computed . Although the noise in the background region may not be exactly the same as some speckle noises, the adjustable parameter k compensates that and the value of k can vary with scale and/or orientation. The larger k is, the more noise will be removed, and the best of its value can be determined by trial and error.
To quantify the performance of the algorithm, we compute five quality metrics : contrast-to-noise ratio (CNR), which measures the contrast between image features and noise, and defined to be ; 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 , 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 , 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 , 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 . 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.
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  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.
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  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 . 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).