Open Access
Add to BibSonomyAbstract
Speckle resides in OCT signals and inevitably effects OCT image quality. In this work, we present a novel method for speckle noise reduction in Fourier-domain OCT images, which utilizes the phase information of complex OCT data. In this method, speckle area is pre-delineated pixelwise based on a phase-domain processing method and then adjusted by the results of wavelet shrinkage of the original image. Coefficient shrinkage method such as wavelet or contourlet is applied afterwards for further suppressing the speckle noise. Compared with conventional methods without speckle adjustment, the proposed method demonstrates significant improvement of image quality.
© 2017 Optical Society of America
1. Introduction
Optical coherence tomography (OCT) emerges as a powerful biomedical imaging modality providing depth-resolved information based on temporal coherence gating [1]. However due to the nature of temporal coherence, speckle is ubiquitous in OCT imposing a “grainy” appearance onto the images [2]. Such “grainy” features might provide extra information of the imaged subject as the formation of speckle is related to the optical property of the sample. However from the view of image quality, this speckle appearance (also known as speckle noise) degrades the clarity of OCT images and can greatly impede the quantitative image analysis for clinical diagnosis such as edge detection, segmentation, pattern recognition, etc. Hereby we refer the term “speckle” as speckle noise in this work.
Many methods have been reported for reducing speckle noise, most of which can be roughly categorized into system-based approaches and algorithm-based approaches. The system-based approaches include angular compounding [4–7], frequency compounding [8], partially spatial coherent source illumination [9, 10], adjacent frame averaging [11], etc. Some approaches within this category such as angular compounding depends on the different speckle patterns among repeated B-scans that utilize different sample illumination schemes. Hence after frame-averaging, the speckle noise is minimized without the loss of image fidelity. However one drawback is that the lateral resolution is compromised as only part of the objective NA is used for each illumination. Further, for this category, system modification is required which complicates the OCT system. On the other hand, the algorithm-based approaches do not need system modification. This category includes coefficients shrinkage in the transform domain (such as wavelet [12, 13], curvelet [14] and contourlet [15]), anisotropic diffusion filter [16], Bayesian least-square estimation [17], compressive sensing [18], and neural network algorithm [19], etc. Due to the fact that speckle noise can hardly be differentiated individually from the “true” OCT image-formation signals but be studied statistically [2], some fine details of OCT images might be lost or altered during speckle reduction. For example in the coefficients shrinkage methods, higher threshold results in better noise reduction but also causes the loss of image features [13–15].
In this paper we propose a new speckle reduction method based on speckle adjustment for Fourier-domain OCT (FD-OCT) images. Different from conventional speckle reduction methods that operate upon the whole data set, this method first identifies the pixels that impose the “grainy” appearance on the image thus degrade the image quality using the phase information of complex OCT signals [20]. Hereby we denote these pixels as noise pixels while the rest of the pixels as image pixels in this work. Then the intensity values of noise pixels are directly adjusted by the averaged intensity values of surrounding image pixels. After the speckle adjustment, conventional speckle reduction method such as wavelet or contourlet shrinkage is then applied. To the best of our knowledge such a processing method utilizing identified noise pixels for speckle noise reduction has not been reported. Due to the fact that the direct speckle adjustment greatly reduces speckle noise without altering the OCT image signals, this proposed method leads to, as demonstrated below, better results than conventional methods that are not based on speckle adjustment under otherwise the same condition. In this paper, Section 2 details this speckle adjustment method for speckle noise reduction. Section 3 discusses the results of applying this method to a fingertip OCT image and an endoscopic OCT image of guinea pig esophagus. The results are compared with those using the same speckle reduction methods but without speckle adjustment.
2. Method
The speckle-adjustment-based noise reduction method is performed in two steps. In the first step noise pixels are identified and then the intensity value for those pixels adjusted. In the second step, a conventional speckle reduction method is applied to the speckle adjusted image.
2.1 Speckle detection and adjustment
The rationale and approach for pixelwise speckle detection has been detailed in our previous work [20]. In brief, the phase of the OCT signal is linear to depth for a single backscattered event. When multiply backscattered events are detected simultaneously by OCT (which leads to speckle), all the backscattered fields will contribute to the phase of the detected OCT signal. In this case, the linearity of phase along depth is disturbed. By detecting such phase disturbance, the noise pixels can be identified. In practice, the phase of each complex A-scan data is unwrapped over the depth profile after Fourier transform is performed to the raw data. Then the noise pixels are determined by evaluating the first and second derivatives of the phase along depth.
In this work, the images were acquired from two custom-built spectral-domain OCT (SD-OCT) systems. Hereby we termed these two systems as system 1 and system 2, respectively. System 1 used a Superluminescent Diode (SLD) light source centering at 870 nm with an ~180 nm full width at half maximum (FWHM) spectral range. It achieved a measured axial resolution of 3.4 µm (in air, non-optimized) and a lateral resolution of 7.0 µm. The imaging size was 2 mm by 1.23 mm, corresponding to 2048 x 2048 pixels, along the lateral and axial directions, respectively. The system operated at a 70 kHz A-scan rate. System 2 was an endoscopic OCT system using a Ti:Sapphire laser with a central wavelength of 825 nm and an FWHM of ~150 nm. The system operated at a 20 kHz A-scan rate with a resolution of 3.0 µm (axial in air) x 6.2 µm (lateral). Each frame consisted of 2048 × 2048 (lateral × axial) pixels representing one circumferential scan. The axial imaging depth was about 1.23 mm.
Figure 1(a) illustrates a fingertip image acquired by system 1. Figure 1(b) illustrates the detected speckles plotted in red within the area marked by the yellow box in Fig. 1(a), where the speckles are overlaid onto the grayscale structural OCT image for direct comparison.
Fig. 1 The detected speckles from a fingertip OCT image. (a) Original OCT image. (b) Detected speckles overlaid onto the OCT intensity image within the area marked by the yellow box in (a). The speckles are plotted in red.
After the noise pixels are delineated, the value of each noise pixel is then adjusted to the average value of its surrounding image pixels, which were obtained through a 9 × 9 window size. The same area marked by the yellow box in Fig. 1(a) before and after the speckle adjustment was plotted in Figs. 2(a) and (b), respectively.
Fig. 2 The marked area in Fig. 1(a) before speckle adjustment (a) and after speckle adjustment (b). Scale bar: 20 µm.
In this step, the computation is cost effective, and for an image of a 2048 × 2048 pixels, the computation time for speckle adjustment was less than 1 second on an Intel Core i5 desktop computer using Matlab.
2.2 Conventional method applied on the speckle adjusted image
In the second step, after the intensity of the noise pixels on the image is adjusted, we apply the double-density dual-tree complex wavelet transform (DD-DT-CWT) shrinkage as the conventional method on this speckle adjusted fingertip image. DD-DT-CWT has the properties of both double-density discrete wavelet transform and dual-tree discrete wavelet transform [21]. The threshold K is set as 0.2. The result of DD-DT-CWT shrinkage is plotted in Fig. 3, corresponding to the same area as illustrated in Figs. 2(a) and 2(b).
Fig. 3 The resulted image for the area plotted in Fig. 2(b) after being processed with DD-DT-CWT. Scale bar: 20 µm.
3. Results
3.1 Quantitative image quality metrics
We adopt the following commonly used image quality metrics to quantify the performance of this proposed method [12, 22]:
where SNR is the peak signal to noise ratio, and and in Eq. (1) are the linear magnitude image and the variance of the background region in the linear magnitude image; CNR is the contrast-to-noise ratio. and in Eq. (2) represent the mean and variance of the same background region as used for computing SNR. and represent the mean and variance of the region of interest. R is the total number of regions of interest; ENL is the equivalent number of looks, which measures the smoothness of homogeneous region of interest. and in Eq. (3) represent the mean and variance of the homogeneous region of interest, and H is the total number of homogeneous regions of interest.In addition to the above three parameters for assessing image quality, sharpness preservation is another one, and the traditional parameter as denoted in [12] is based on the Laplace operation performed on both the original image and the processed image. One potential issue with this traditional parameter is that it takes all the edges into account when comparing the processed image with the original image in the regions of interest, including those artificial edges caused by speckle. To circumvent this problem and considering the fact that each interface defines two regions with distinct intensities, the change of interface sharpness leads to the intensity change within each region. In biological sample, such an interface can be easily appreciated as each layered tissue structure provides two interfaces. Hereby we evaluate the edge sharpness by measuring the intensity density within a layered tissue structure. In practice, the intensity along a selected direction (either lateral or longitudinal) within a given tissue layer is normalized and then averaged. We denote this result as AREA. An example will be discussed in detail as illustrated in Fig. 5(c) in the following sections. As demonstrated later, if the edges are blurred due to speckle reduction such as wavelet shrinkage, the AREA will increase. Hence the smaller the AREA, the sharper the edges are.
3.2 Fingertip image
To test the performance of this speckle prior based method, the fingertip image as shown in Fig. 1(a) was processed following the above mentioned steps, namely speckle identification and adjustment. The conventional DD-DT-CWT shrinkage was employed in step 2. Figure 4(a) was reproduced from Fig. 1(a) for convenience. The results of applying this method together with the results of conventional DD-DT-CWT shrinkage without speckle adjustment were plotted in Figs. 4(b) and 4(c), respectively. The threshold K for both speckle adjustment based method and DD-DT-CWT shrinkage is set at 0.2.
Fig. 4 The results of speckle reduction algorithms applied to a fingertip image. (a) Original OCT image. (b) Image processed by the proposed method using DD-DT-CWT shrinkage after speckle adjustment. (c) Image processed by conventional DD-DT-CWT shrinkage without speckle adjustment. In (a), the green boxes marked the homogeneous areas. The red boxes marked the heterogeneous areas. The yellow box marked the background area. The dashed blue line marked the pixel values evaluated for calculating AREA. Scale bar: 200 µm.
In Fig. 4(a), the four green boxes marked the representative homogeneous areas used in Eq. (3) for calculating ENL. Hereby the homogeneous areas were visually determined for evaluation ENL. The three red boxes marked the representative heterogeneous areas together with the homogenous areas marked by the four green boxes were used for calculating CNR. The yellow box in Fig. 4(a) marked the background area used in Eqs. (1) and (2). In order to quantify the AREA the pixel values along the dashed blue line as denoted in Fig. 4(a) were examined and illustrated in Figs. 5(a) and 5(b) for K = 0.2 and K = 1.2, representing normal and strong threshold levels respectively.
Fig. 5 The pixel values along the blue dashed line in Fig. 4(a) from the original image, image processed by the proposed method and image processed by DD-DT-CWT shrinkage for K = 0.2 (a) and K = 1.2 (b). (c) The area of interest for evaluating AREA. The maximum pixel value for the proposed method and DD-DT-CWT shrinkage for both K = 0.2 and K = 1.2 is normalized to 1.
In a wavelet transform that is free of noise, most of the coefficients are effectively zero. Such zero coefficients become nonzero and have small values when noise is introduced. In an ideal case, by thresholding the coefficients, only stronger coefficients that represent signals are kept and small nonzero coefficients caused by noise are excluded. Thus larger threshold leads to a higher probability of eliminating the coefficients that represent signals while smaller threshold increases the chance to nosier results.
Figure 5(c) illustrates the area of interest indicated by the yellow box, from which the intensity density is calculated. In Fig. 5(c), the maximum intensity for each case is normalized to 1. From Fig. 5(a) it is clear that both the proposed method and the DD-DT-CWT shrinkage lead to a greatly reduced noise comparing to the original image, while the proposed method shows better noise reduction than the DD-DT-CWT shrinkage using the same threshold level. From Fig. 5(c), it is found that within the same area of interest, the intensity shape for the case of K = 1.2 is flatter than that for the case of K = 0.2, indicating a higher intensity density.
The calculated image quality metrics are summarized in Table 1 including the cases for threshold K = 1.2.
Table 1. The comparison of the image quality metrics calculated from the fingertip image
From this fingertip image, for K = 0.2 the speckle adjustment based method increases CNR by 6.28 dB and SNR by 28.95 dB from the original image. Compared with the DD-DT-CWT shrinkage without speckle adjustment, this speckle adjustment based method leads to the similar AREA indicating similar edge preservation and similar SNR, but 1.69 dB higher in CNR and 226.51 higher in ENL. For K = 1.2, the AREA of both methods has increased comparing with the case of K = 0.2, indicating loss of sharpness due to speckle reduction.
To further explore the effect of the choice of threshold K upon the performance of our proposed method, we simulate the change of CNR and AREA with respect to the change of K in both the proposed method and the DD-DT-CWT shrinkage. The results are illustrated in Fig. 6.
With the increase of K, the CNR from the images processed by both the proposed method and DD-DT-CWT shrinkage increases. The proposed method consistently shows a better performance regarding CNR for all K values. When K = 0, which means no thresholding is applied, the proposed method has 2.2 dB higher in CNR which results from the pixelwise speckle adjustment. Furthermore in Fig. 6 it is noted that in order to reach the same CNR, the threshold applied for the proposed method is smaller than the DD-DT-CWT shrinkage method. By using a smaller threshold, the image details are expected to be better preserved. For AREA, overall it increases with the increase of K, indicating decreased sharpness of the image by applying a higher threshold in DD-DT-CWT shrinkage. It is noted that in this example, for K less than 0.1 in the case of the proposed method and K less than 0.5 in the case of DD-DT-CWT, the speckle noise affects the intensity within the evaluated region and hence the calculation of AREA.
3.3 Endoscopic guinea pig esophagus image
The speckle adjustment based method was also applied to an OCT image of guinea pig esophagus acquired by system 2 as mentioned in Section 2.1. As for the images with layered structures, contourlet [23] or curvelet [24] shrinkages are more often employed than wavelet shrinkage. Hence instead of using DD-DT-CWT shrinkage, hereby the contourlet shrinkage is applied in the second step of the proposed method, and the results of the proposed method are compared with those of contourlet shrinkage without speckle adjustment.
To better demonstrate the performance of this method, we first show the images in a rectangular coordinate before polar conversion as illustrated in Figs. 7(a)-7(c) for the original OCT image, image processed by the proposed method and image processed by the contourlet shrinkage, respectively. The threshold K is set at 0.08, which turns out to provide a good balance between denoising and image blurring. It needs be noted that Figs. 7(a)-7(c) are the B-scan images acquired within one 360° rotation. The scale bar indicates the scale along the depth direction. Each image contains 2048 A-scans. Figures 7(d)-7(f) illustrate the corresponding images after polar conversion. Figures 7(g)-7(i) illustrate the magnified images of Figs. 7(d)-7(f), respectively, corresponding to the area denoted by the black box in Fig. 7(d).
Fig. 7 The results of speckle reduction algorithms applied to an endoscopic OCT image of guinea pig esophagus. (a) Original OCT image before polar convention. (b) Image processed by the proposed method before polar convention. (c) Image processed by contourlet shrinkage before polar convention. (d) Original OCT image after polar convention. (e) Image processed by the proposed method after polar convention. (f) Image processed by contourlet shrinkage after polar convention. (g)-(i) The magnified images of (d)-(f) corresponding to the area denoted by the black box in (d). In (a), the two green boxes marked the representative homogeneous areas. The red boxes marked the representative heterogeneous areas. The yellow box marked the background area. The vertical dashed blue line indicates the direction for calculating AREA.
In Fig. 7(a) the two green boxes denoted the representative homogeneous areas used in Eq. (3) for calculating ENL. The two red boxes denoted the representative heterogeneous areas, together with the homogeneous areas, were used in Eq. (2) for calculating the CNR. The yellow box denoted the background area used in Eqs. (1) and (2). The dashed blue line indicates the direction along which the pixel values of the homogeneous area are evaluated for calculating AREA.
The results of the quantitative image quality metrics are shown in Table 2 for K = 0.08. For this esophagus OCT image, this proposed method increases CNR by 2.93 dB, SNR by 9.91 dB and ENL by 45.11 when compared with contourlet shrinkage; thus the proposed method offers a higher CNR, SNR and ENL with a similar AREA.
Table 2. The comparison of the image quality metrics calculated from the guinea pig esophagus image
4. Discussion and conclusions
In summary, this paper proposes a new speckle reduction method by first identifying speckles and then adjusting the intensity of the noise pixels to the average intensity of neighboring image pixels, followed by traditional de-speckling algorithms such as DD-DT CWT shrinkage or contourlet shrinkage. The results of applying the proposed method to an OCT fingertip image and a guinea pig esophagus image show superior performance of this new method comparing with the results of conventional speckle reduction methods without speckle adjustment. Different from other speckle reduction methods, speckle adjustment in this proposed method reduces the speckle noise without sacrificing the image quality as only the noise pixels are adjusted. This benefit comes from the phase information of the OCT signals which is not usually exploited in traditional de-speckling image processing approaches. Furthermore due to the fact that the speckle noise has been reduced noticeably in speckle adjustment, only a small threshold is needed in the following coefficients shrinkage, hence the side-effects of the transform domain shrinkage technique can be minimized and the fidelity of the processed image can be greatly preserved. It should be noted that since the phase of OCT signal is utilized, this method suffers from large phase noise from low SNR regions of the OCT image. In practice, however, due to the fact that the intensity of each noise pixel is adjusted according to its surrounding intensities, at low SNR regions or background such drawback is not severe as manifested in Figs. 4 and 7.
Compared with contourlet shrinkage and curvelet shrinkage which are more suitable for images with layered structures such as retina [14, 15], the speckle adjustment in our proposed method is not limited by sample structures as demonstrated by the fingertip image and the esophagus image.
Funding
National Institutes of Health (R01 CA153023 and R01 HL121788); The Coulter Foundation.
Acknowledgement
The authors acknowledge the help with the contourlet program from Dr. Kenneth K. Y. Wong and Mr. Jianbing Xu at the Department of Electrical and Electronic Engineering, the University of Hong Kong.
References and links
1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, and C. A. Puliafito, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef] [PubMed]
2. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in Optical Coherence Tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). [CrossRef] [PubMed]
3. A. Mariampillai, B. A. Standish, E. H. Moriyama, M. Khurana, N. R. Munce, M. K. K. Leung, J. Jiang, A. Cable, B. C. Wilson, I. A. Vitkin, and V. X. D. Yang, “Speckle variance detection of microvasculature using swept-source optical coherence tomography,” Opt. Lett. 33(13), 1530–1532 (2008). [CrossRef] [PubMed]
4. J. M. Schmitt, “Array detection for speckle reduction in optical coherence microscopy,” Phys. Med. Biol. 42(7), 1427–1439 (1997). [CrossRef] [PubMed]
5. M. Bashkansky and J. Reintjes, “Statistics and reduction of speckle in optical coherence tomography,” Opt. Lett. 25(8), 545–547 (2000). [CrossRef] [PubMed]
6. N. Iftimia, B. E. Bouma, and G. J. Tearney, “Speckle reduction in optical coherence tomography by “path length encoded” angular compounding,” J. Biomed. Opt. 8(2), 260–263 (2003). [CrossRef] [PubMed]
7. T. Klein, R. André, W. Wieser, T. Pfeiffer, and R. Huber, “Joint aperture detection for speckle reduction and increased collection efficiency in ophthalmic MHz OCT,” Biomed. Opt. Express 4(4), 619–634 (2013). [CrossRef] [PubMed]
8. M. Pircher, E. Gotzinger, R. Leitgeb, A. F. Fercher, and C. K. Hitzenberger, “Speckle reduction in optical coherence tomography by frequency compounding,” J. Biomed. Opt. 8(3), 565–569 (2003). [CrossRef] [PubMed]
9. J. Kim, D. T. Miller, E. Kim, S. Oh, J. Oh, and T. E. Milner, “Optical coherence tomography speckle reduction by a partially spatially coherent source,” J. Biomed. Opt. 10(6), 064034 (2005). [CrossRef] [PubMed]
10. B. Karamata, P. Lambelet, M. Laubscher, P. P. Salathé, and T. Lasser, “Spatially incoherent illumination as a mechanism for cross-talk suppression in wide-field optical coherence tomography,” Opt. Lett. 29(7), 736–738 (2004). [CrossRef] [PubMed]
11. M. Szkulmowski, I. Gorczynska, D. Szlag, M. Sylwestrzak, A. Kowalczyk, and M. Wojtkowski, “Efficient reduction of speckle noise in Optical Coherence Tomography,” Opt. Express 20(2), 1337–1359 (2012). [CrossRef] [PubMed]
12. D. C. Adler, T. H. Ko, and J. G. Fujimoto, “Speckle reduction in optical coherence tomography images by use of a spatially adaptive wavelet filter,” Opt. Lett. 29(24), 2878–2880 (2004). [CrossRef] [PubMed]
13. S. Chitchian, M. A. Mayer, A. R. Boretsky, F. J. van Kuijk, and M. Motamedi, “Retinal optical coherence tomography image enhancement via shrinkage denoising using double-density dual-tree complex wavelet transform,” J. Biomed. Opt. 17(11), 116009 (2012). [CrossRef] [PubMed]
14. Z. Jian, Z. Yu, L. Yu, B. Rao, Z. Chen, and B. J. Tromberg, “Speckle attenuation in optical coherence tomography by curvelet shrinkage,” Opt. Lett. 34(10), 1516–1518 (2009). [CrossRef] [PubMed]
15. J. Xu, H. Ou, E. Y. Lam, P. C. Chui, and K. K. Y. Wong, “Speckle reduction of retinal optical coherence tomography based on contourlet shrinkage,” Opt. Lett. 38(15), 2900–2903 (2013). [CrossRef] [PubMed]
16. R. Bernardes, C. Maduro, P. Serranho, A. Araújo, S. Barbeiro, and J. Cunha-Vaz, “Improved adaptive complex diffusion despeckling filter,” Opt. Express 18(23), 24048–24059 (2010). [CrossRef] [PubMed]
17. A. Cameron, D. Lui, A. Boroomand, J. Glaister, A. Wong, and K. Bizheva, “Stochastic speckle noise compensation in optical coherence tomography using non-stationary spline-based speckle noise modelling,” Biomed. Opt. Express 4(9), 1769–1785 (2013). [CrossRef] [PubMed]
18. D. Xu, N. Vaswani, Y. Huang, and J. U. Kang, “Modified compressive sensing optical coherence tomography with noise reduction,” Opt. Lett. 37(20), 4209–4211 (2012). [CrossRef] [PubMed]
19. M. R. N. Avanaki, P. P. Laissue, T. J. Eom, A. G. Podoleanu, and A. Hojjatoleslami, “Speckle reduction using an artificial neural network algorithm,” Appl. Opt. 52(21), 5050–5057 (2013). [CrossRef] [PubMed]
20. A. Zhang, J. F. Xi, W. X. Liang, T. Gao, and X. D. Li, “Generic pixel-wise speckle detection in Fourier-domain OCT images,” Opt. Lett. 39(15), 4392–4395 (2014). [CrossRef] [PubMed]
21. I. Selesnick, “The double-density dual-tree DWT,” IEEE Trans. Signal Process. 52(5), 1304–1314 (2004). [CrossRef]
22. A. Ozcan, A. Bilenca, A. E. Desjardins, B. E. Bouma, and G. J. Tearney, “Speckle reduction in optical coherence tomography images using digital filtering,” J. Opt. Soc. Am. A 24(7), 1901–1910 (2007). [CrossRef] [PubMed]
23. M. N. Do and M. Vetterli, “The contourlet transform: An efficient directional multiresolution image representation,” IEEE Trans. Image Process. 14(12), 2091–2106 (2005). [CrossRef] [PubMed]
24. E. J. Candès, L. Demanet, D. L. Donoho, and L. Ying, “Fast Discrete Curvelet Transforms,” Multiscale Model. Simul. 5(3), 861–899 (2006). [CrossRef]
25. J. Xi, A. Zhang, Z. Liu, W. Liang, L. Y. Lin, S. Yu, and X. Li, “Diffractive catheter for ultrahigh-resolution spectral-domain volumetric OCT imaging,” Opt. Lett. 39(7), 2016–2019 (2014). [CrossRef] [PubMed]
