Despeckling optical coherence tomograms from the human retina is a fundamental step to a better diagnosis or as a preprocessing stage for retinal layer segmentation. Both of these applications are particularly important in monitoring the progression of retinal disorders. In this study we propose a new formulation for a well-known nonlinear complex diffusion filter. A regularization factor is now made to be dependent on data, and the process itself is now an adaptive one. Experimental results making use of synthetic data show the good performance of the proposed formulation by achieving better quantitative results and increasing computation speed.
© 2010 OSA
Optical coherence tomography (OCT) is a non-invasive imaging modality with an increasing number of applications and it is becoming an essential tool in ophthalmology  allowing in vivo high-resolution cross-sectional imaging of the retinal tissue. However, as any imaging technique that bases its image formation on coherent waves, OCT images suffer from speckle noise, which reduces its quality.
Speckle noise creates a grainy appearance that can mask diagnostically significant image features (small or low reflectivity features) and reduce the accuracy of segmentation and pattern recognition algorithms [2–4].
The statistical mechanism of laser speckle formation was first presented by Goodman . In addition to the theoretical results, this study also supports the idea that speckle noise could be rejected by linear filtering. Later, Wagner et al. , Burckhardt et al.  and Abbott et al.  concluded that linear filtering, the way it was presented in , suppresses the noise at the cost of smoothing out image details.
Speckling is a common problem for different imaging modalities such as radio astronomy, synthetic aperture radar (SAR), ultrasound and laser imaging. Consequently, it promotes research and results in many speckling reduction techniques .
Schmitt  described the first OCT speckle suppression technique, in which a compounded image was formed from the sum of the signals obtained from a quadrant photodiode detection system, and was followed in his approach in the work of Bashkansky et al.  and Iftimia et al. . Similarly to angular compounding, other methods to reduce speckling applied before image formation (physical techniques), such as frequency compounding or spatial compounding, were also developed.
Skankar et al.  and by Pircher et al.  proposed a frequency compounding technique increasing contrast and improving quality image, without loss of resolution, while Kim et al.  presented a space diversity speckle reduction technique achieving a substantial reduction in speckling at the cost of decreasing transverse resolution.
Another complex-domain processing technique was applied by Yung et al.  using a zero-adjustment procedure (ZAP) to OCT following its successful application in medical ultrasound by Healey et al. .
Since physical techniques require the modification of OCT system design, postprocessing methods were also developed being the CLEAN algorithm one of the first techniques used for OCT despeckling [17,18]. Amongst others are median filtering , homomorphic Wiener filtering , enhanced Lee filter (ELEE) , symmetric nearest neighbor (SNN) filter, adaptive noise smoothing , multiresolution wavelet analysis , filtering techniques based on rotating kernel transformations , Kuwahara filter  and anisotropic diffusion filtering [4,26].
The median filter calculates the median value in a local neighborhood of each pixel. Koozekanani et al.  applied it to reduce speckle of OCT images as a preprocess to a posterior measurement of the retinal thickness. On the other hand, Ishikawa et al. , with the same purpose of measuring the retinal thickness, used a modified mean filter to reduce noise. The gaussian filter is yet another special filter applied to retinal segmentation algorithms [29–32].
In adaptive filtering, the algorithm is modified locally based on the pixel’s neighbors, combining an effective noise reduction and an ability to preserve image edges. The Lee filter is one of these adaptive filters that has been successfully applied to OCT images.
Ozcan et al.  compared the performance of different filters (an enhanced Lee filter, two wavelet-transform-based filters, a hybrid median filter, a symmetric nearest-neighbor filter, a Kuwahara filter and an adaptive Wiener filter), applied independently or in association, to reduce the speckling effect on an OCT tomogram of a bovine retina . A fuzzy thresholding algorithm in the wavelet domain was proposed by Puvanathasan and Bizheva  to remove speckle noise in OCT images of a human fingertip.
The main purpose of the work presented herein is to report on a system to reduce the speckle noise for both the visual assessment and the improved structure segmentation on high-definition spectral domain Cirrus OCT (Carl Zeiss Meditec, Dublin, CA, USA). This retinal imaging system allows the acquisition of volumes of 200x200x1024 or 512x128x1024 voxels, respectively, for the lateral, azimuthal and axial directions (Fig. 1 ). These volumetric data are obtained from a 6000x6000x2000 μm3 volume of the human macula. Additionally, high-resolution B-scan images of 1024x1024 pixels can be obtained.
The OCT working principle is similar to that of ultrasound, and uses some terminology from that field. The volumetric OCT information is composed of a set of A-scans (depth-wise information on refractive index changes). The scanning is performed along a series of coplanar lines (B-scans) covering a 20° field-of-view of the eye fundus.
The aim of the work herein is to improve the process of speckle noise reduction and to improve the preservation of edge and image features. We specifically aim to apply this filtering process to OCT data from the human eye fundus as a preprocessing step for layer segmentation, which explains the need to preserve the retinal tissue. Additionally, the process intends to be applied as a filter for visual inspection, as it preserves features within the tissue.
As a proof-of-concept, the proposed method will be compared with traditional NCDF methods from the literature by means of quantitative measures, such as the mean square error (MSE), contrast-to-noise ratio (CNR) and average effective number of looks (ENL).
2.1 Nonlinear complex diffusion approach
This equation was discretized by a forward in time and centered in space (FTCS) finite difference scheme,
Assuming that, for all , it can be shown (see ) that this explicit method is stable if
In the work of Salinas et al. , the coefficient of diffusion adopted is39], for small θ the imaginary part of I is a smooth function of its second derivatives and the ratiois proportional to the laplacian of I (in the limit when θ→0). In this way, diffusion is facilitated for smooth areas and attenuated at image edges.
The diffusion coefficient can be approximated by
As motivation for these kind of functions, this expression can be seen as a Lorentzian function (Eq. (8) modified to have its maximum equal to 1 and w = k/2.Eq. (7) as shown in Fig. 2 , simply by modifying the value of k.
The choice for the κparameter (Eq. (6) is therefore important, as it modulates the spread of the diffusion coefficient in the vicinity of its maximum, that is, at edges and homogeneous areas, where the image laplacian vanishes. From the plot, it becomes clear the difference in D for ΔI constant from low- (higher k) to high-intensity areas (lower k), thus increasing the diffusion for low-intensity areas and decreasing it for higher-intensity ones.
2.2 Improved coefficient of diffusion
While the formulation for the coefficient of diffusion (Eq. (6) seems to be a good choice to preserve the location of edges, as stated in , it does not specifically address our need to preserve the variation of intensity across the edge and features of the image within the tissue, which might be improved.
By definition, edges are located in between two areas of perceived distinct intensity levels within the image. Additionally, the OCT background is characterized by a low average intensity level. Conversely, the tissue due to the differences in the refractive index, presents a higher average level in the image.
We intend to manipulate Eq. (6) to take advantage of these facts. In this way, we aim to facilitate diffusion at lower intensity level areas (e.g., vitreous, cysts, fluid-filled regions) and to be conservative within the retinal tissue areas (preserving important details for analysis). As an advantageous byproduct, edges are better preserved as diffusion is decreased at the higher level side of the edge.
To this end, we need to locally modify κ according to data, increasing κ at low level areas and making it smaller at higher level areas. This is a clear distinction to the formulation using (Eq. (6) where κ is made constant for the entire image and over iterations.
We proposed to adaptκlocally by the use of the function
Additionally, using a gaussian filter for D is beneficial to increase speckle removal. While in edges, filtering D does not change D significantly as long as σ is small since is smooth . However, at isolated spiky D points, filtering D turns diffusion less conservative and therefore increases speckle reduction without compromising edge preservation.
2.3 Adaptive time step
As opposed to the majority of nonlinear complex diffusion processes that adopt a constant time step (Δt) close to the time step limit of the convergence of the iterative update process, we opted for an adaptive time step (Eq. (11). The rationale behind this decision is based on the fact that the coefficient of diffusion depends on the gradient of the image (Eq. (6) and, due to noise, this gradient is much higher in the initial steps of the diffusion process.Eq. (5) and parameters a and b control the time step (with a + b≤1).
The typical evolution ofover the iterative process is shown in Fig. 3 .
As expected, a small step size is used at the initial iterations in which higher values of D can be found due to the speckle noise. This is, at steady conditions in which changes over time are small (fraction-wise), the time step can be made larger, while at fast changes in time the time step can be made smaller.
2.4 Quality metrics
To perform a quantitative comparison between the performances of the different filters we computed some well-known speckle-reduction performance metrics . The first measure is the mean squared error (MSE), defined by (Eq. (12), where denotes the samples of the original image, denotes the samples of the filtered image and M and P are the number of pixels in row and column directions, respectively.Eq. (13).Eq. (14) and is a measure between an area of image feature (R) and an area of background noise (H). Thus µr and σr are the mean and variance of the regions of interest, including the homogeneous regions and µH and σH are the mean and variance of the background noise regions.
Considering the number of different possibilities for the testing conditions when comparing the proposed formulation with the original one, tests were performed using: 1) synthetic images added with two different noise models; 2) real OCT data using the total diffusion time proposed on the original formulation (12 s) and the diffusion time we have considered for the qualitative results to be presented (3 s); and 3) the required computing time to achieve the same ENL parameter as the original formulation.
A series of tests were performed to quantify the performance of the proposed filter. In comparison to the traditional nonlinear complex diffusion filter (NCDF), the proposed filter differs in making κ to locally adjust to data (Eq. (9), to low-pass filter the coefficient of diffusion D and to adopt an iterative step in time approach (Eq. (11).
A synthetic image composed of four different areas (Fig. 4 , left) was generated: one presents a 2D gaussian profile (top left), one is composed of concentric rings (top right), one presents a step and a ramp edge (bottom left) and, finally, one simulates the foveal area (bottom right).
Two different noise sources were independently tested, adding the respective generated noise to the synthetic image to compute the performance of the filters.
In the first case, a multiplicative noise type was generated using Matlab  speckle “imnoise” function. This noise was added to the synthetic image (I) to become Inoise = I + γI, where γ is a uniformly distributed random noise with zero mean and variance 0.10 (Fig. 4, center).
In the second case, two independent gaussian functions were generated using Matlab’s “randn” function. The Hadamard product of these two functions was calculated and scaled by a factor of 40 (representing ~15% of the signal level) and added to the synthetic image, therefore resulting in the noisy image Inoise = I + 40.randn.randn (Fig. 4, right).
The following results were obtained using for the proposed filter a diffusion time of 3 s, aand(Eq. (9) and(Eq. (6). The parameters for the gaussian kernel of the low-pass filter for g (Eq. (10) were N = 3 and σ = 10, while for the coefficient of diffusion D we considered a low pass filter with N = 3 and σ = 0.5. The time step parameters were a = 0.25 and b = 0.75 (Eq. (11).
The traditional NCDF was applied using 50 iterations,,and a time step (Δt) of 0.24 s, as in Salinas et al. .
The qualitative assessment of the proposed filter performance can be made from Fig. 5 .
The row-wise comparison shows the smoother images on the right resulting from the proposed filter against the traditional NCDF filter (left).
An intensity profile from the bottom left quadrant synthetic noisy image (Fig. 5, bottom left) can be seen in Fig. 6 , presenting step and ramp edges. Overlapped in different colors to facilitate the comparison are the original signal data (black solid line), the noisy signal (dotted-blue line) and the result of the traditional NCDF and of the proposed filter, the solid green and solid red lines, respectively.
The smoothness of the solid red line comparatively to the solid green line is obvious. Simultaneously the former is closer to the original (noise-free) signal. Moreover, the proposed filter allows for smoother filtered signals at low-intensity signal regions compared to higher-intensity ones, therefore showing it to be tuned to preserve retinal tissue and to be less conservative within the vitreous regions when applied to OCT of the human eye fundus.
In addition, on the ramp edge it demonstrates that smooth transitions on signal are also and better preserved as compared to traditional NCDF filter. This illustrates the possibility of better preserving the small transitions occurring in between some of the inner retinal layers.
To quantitate the filter performance measurement, the three aforementioned metrics, the ENL, the MSE and the CNR, were computed for each of 50 runs. Both the ENL and the MSE were separately analyzed in low- and high-intensity regions of the image (dashed-line regions in Fig. 4, left).
When added to the better quantitative performance in despeckling, the proposed formulation allows these results to be obtained in a fraction of the computing time.
The computing time for data shown in the table (with a diffusion time of 12 seconds for the traditional NCDF) is (mean ± SD) 4.2 ± 0.05 seconds for the traditional NCDF for both noise models, because it is nonadaptive. The proposed formulation (with a diffusion time of 3 seconds), however, presents a computing time of 2.1 ± 0.03 seconds for the “Speckle Imnoise” noise and 5.6 ± 0.08 seconds for the “Randn.Randn” noise, showing the adaption required by this noise model.
Tests conducted in the real human eye fundus OCT data (32 B-scans using the 200x200x1024 macular cube protocol from 32 eye fundus scans from 13 healthy volunteers, 3 eyes with choroidal neovascularization, 2 with cystoid macular edema, 9 with diabetic retinopathy and 5 with age-related macular degeneration) allows comparative performances of 3.4 ± 0.06 to 2.1 ± 0.09 seconds, respectively, for the traditional and proposed filter formulation (61.8% of computing time).
To achieve the same ENL (Eq. (13) over a vitreous region as the traditional NCDF filter (12 seconds of diffusion time) for the above test set (32 B-scans) the comparative performance is 3.4 ± 0.05 and 1.1 ± 0.12 seconds, respectively, for the traditional and the proposed formulation, i.e., the same results were achieved in 34% of the computing time.
The application of the proposed filter to an OCT high-resolution B-scan from the human eye fundus is shown in Fig. 7 . For easy of assessment, a region is shown enlarged and using pseudocolor. In the bottom row, both the traditional (left) and the proposed filter (right) results are shown.
4. Discussion and conclusions
In this report, we propose a new formulation of the nonlinear complex diffusion filter and demonstrate its advantages with regard to edge preservation, speckle filtering capabilities and potential to recover the original (uncorrupted) signal, hence the use of synthetic data.
Based on our results the proposed formulation achieved a superior performance while retaining all the advantages of nonlinear complex diffusion filters. One particular advantage to the original formulation is that the κ parameter value does not need to be defined beforehand; instead, it adapts itself to the data.
For parameters needing to be defined, such as θ, the value used by Salinas et al. , determined to be the best choice for the application of the NCDF filter to the human eye fundus OCT data, was used. On the other hand, the proposed formulation achieved better quantitative and qualitative results in 61.8% of computing time; for example, for a diffusion time of 3 seconds, typically corresponding to 22-26 iterations, against the 12 seconds of diffusion time as in  in 50 iterations. Alternatively, similar results could be achieved in 34% of the computing time.
The adaptive time step seems to be particularly important. While for 2D the maximum step in time to achieve stability of the iterative update is 0.25 second, it only assures this stability as a whole, i.e., for the entire image, but cannot assure the stability of the method for individual pixels. By using an adaptive step in time process, making smaller steps in time according to iterative image changes, we can solve this problem.
Although slightly more complex to compute, the new formulation allows better results to be achieved in less total computing time, an additional advantage of the proposed filter in comparison to the original formulation.
We intend to make use of this new NCDF filter both for helping in visual assessment of retinal structures and as a preprocessing step for image/volume segmentation. These two distinct application areas may require different parameters sets, e.g., total diffusion time and κ limits, but this is still unknown and needs further analysis.
In addition to this particular area of application in the fundus of the human eye, this filter may be applied as well to different data sources corrupted with speckle noise, such as medical ultrasound. We will investigate its application in this imaging field.
Preliminary results have shown the potential advantage of using this preprocessing step in OCT retinal layer segmentation, our main reason for beginning this research. This tool may prove to be a powerful asset either in its current capacity or as a preprocessing stage of a multistage image processing system.
This study is supported in part by the Fundação para a Ciência e a Tecnologia (FCT) under the research project PTDC/SAU-BEB/103151/2008. The authors would like to thanks Dr. Melissa Horne and Carl Zeiss Meditec (Dublic, CA, USA) for their support on getting access to OCT data, and Eng. Torcato Santos and AIBILI Clinical Trial Center technicians for their support in managing data, working with patients and performing scans. This study is registered at ClinicalTrials.org (ID: NCT00797524).
References and links
1. W. Drexler, U. Morgner, R. K. Ghanta, F. X. Kärtner, J. S. Schuman, and J. G. Fujimoto, “Ultrahigh-resolution ophthalmic optical coherence tomography,” Nat. Med. 7(4), 502–507 (2001). [CrossRef] [PubMed]
2. A. F. Fercher, Optical Coherence Tomography: Technology and Applications (Springer, New York, 2008), chap. 4.
3. P. Puvanathasan and K. Bizheva, “Speckle noise reduction algorithm for optical coherence tomography based on interval type II fuzzy set,” Opt. Express 15(24), 15747–15758 (2007). [CrossRef] [PubMed]
4. H. M. Salinas and D. C. Fernández, “Comparison of PDE-based nonlinear diffusion approaches for image enhancement and denoising in optical coherence tomography,” IEEE Trans. Med. Imaging 26(6), 761–771 (2007). [CrossRef] [PubMed]
5. J. W. Goodman, “Some fundamental properties of speckle,” J. Opt. Soc. Am. 66, 1145–1150 (1976). [CrossRef]
6. R. Wagner, S. Smith, J. Sandrik, and H. Lopez, “Statistics of speckle in ultrasound b-scans,” IEEE Trans. Sonics Ultrason. 30, 156–163 (1983). [CrossRef]
7. C. Burckhardt, “Speckle in ultrasound b-mode scans,” IEEE Trans. Sonics Ultrason. 25, 1–6 (1978). [CrossRef]
10. M. Bashkansky and J. Reintjes, “Statistics and reduction of speckle in optical coherence tomography,” Opt. Lett. 25(8), 545–547 (2000). [CrossRef]
11. 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]
12. P. Shankar and V. Newhouse, “Speckle reduction with improved resolution in ultrasound images,” IEEE Trans. Sonics Ultrason. SU-32, 537–543 (1985).
13. 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]
14. 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]
15. K. M. Yung, S. L. Lee, and J. M. Schmitt, “Phase-domain processing of optical coherence tomography images,” J. Biomed. Opt. 4, 125–136 (1999). [CrossRef]
16. A. J. Healey, S. Leeman, and F. Forsberg, “Turning off speckle,” Acoust. Imaging 19, 433–437 (1992). [CrossRef]
17. D. L. Fried, “Analysis of the clean algorithm and implications for superresolution,” J. Opt. Soc. Am. A 12, 853–860 (1995). [CrossRef]
18. J. M. Schmitt, “Restoration of optical coherence images of living tissue using the clean algorithm,” J. Biomed. Opt. 3, 66–75 (1998). [CrossRef]
19. R. Bernstein, “Adaptive nonlinear filters for simultaneous removal of different kinds of noise in images,” IEEE Trans. Circ. Syst. 34, 1275–1291 (1987). [CrossRef]
20. G. Franceschetti, V. Pascazio, and G. Schirinzi, “Iterative homomorphic technique for speckle reduction in synthetic-aperture radar imaging,” J. Opt. Soc. Am. A 12, 686–694 (1995). [CrossRef]
21. J. Lee, “Speckle analysis and smoothing of synthetic aperture radar images,” Comput. Graph. Image Process. 17, 24–32 (1981). [CrossRef]
22. D. Kuan, A. Sawchuk, T. Strand, and P. Chavel, “Adaptive noise smoothing filter for images with signal-dependent noise,” IEEE Trans. Pattern Anal. Mach. Intell. 7, 165–177 (1985). [CrossRef] [PubMed]
23. S. H. Xiang, L. Zhou, and J. M. Schmitt, “Speckle noise reduction for optical coherence tomography,” Proc. SPIE 3196, 79–88 (1998). [CrossRef]
24. J. Rogowska and M. E. Brezinski, “Evaluation of the adaptive speckle suppression filter for coronary optical coherence tomography imaging,” IEEE Trans. Med. Imaging 19(12), 1261–1266 (2000). [CrossRef]
25. M. Kuwahara, K. Hachimura, S. Ehiu, and M. Kinoshita, “Processing of riangiocardiographic images”, in Digital Processing of Biomedical Images, (Plenum Publishing, New York, NY, 1976), pp. 187–203.
26. P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. Mach. Intell. 12, 629–639 (1990). [CrossRef]
27. D. Koozekanani, K. Boyer, and C. Roberts, “Retinal thickness measurements in optical coherence tomography using a Markov boundary model”, in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (Hilton Head, SC, 2000), pp. 2363–2370.
28. H. Ishikawa, D. M. Stein, G. Wollstein, S. Beaton, J. G. Fujimoto, and J. S. Schuman, “Macular segmentation with optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 46(6), 2012–2017 (2005). [CrossRef] [PubMed]
29. A. M. Bagci, M. Shahidi, R. Ansari, M. Blair, N. P. Blair, and R. Zelkha, “Thickness profiles of retinal layers by optical coherence tomography image segmentation,” Am. J. Ophthalmol. 146(5), 679–687 (2008). [CrossRef] [PubMed]
30. M. Baroni, P. Fortunato, and A. La Torre, “Towards quantitative analysis of retinal features in optical coherence tomography,” Med. Eng. Phys. 29(4), 432–441 (2007). [CrossRef]
32. O. Tan, G. Li, A. T. H. Lu, R. Varma, D. HuangO. TanG. LiA. T. H. LuR. VarmaD. HuangAdvanced Imaging for Glaucoma Study Group, “Mapping of macular substructures with optical coherence tomography for glaucoma diagnosis,” Ophthalmology 115(6), 949–956 (2008). [CrossRef]
33. B. Sander, M. Larsen, L. Thrane, J. L. Hougaard, and T. M. Jørgensen, “Enhanced optical coherence tomography imaging by multiple scan averaging,” Br. J. Ophthalmol. 89(2), 207–212 (2005). [CrossRef] [PubMed]
34. T. M. Jørgensen, J. Thomadsen, U. Christensen, W. Soliman, and B. Sander, “Enhancing the signal-to-noise ratio in ophthalmic optical coherence tomography by image registration--method and clinical examples,” J. Biomed. Opt. 12(4), 041208 (2007). [CrossRef] [PubMed]
35. 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]
36. D. Cabrera Fernández, H. M. Salinas, and C. A. Puliafito, “Automated detection of retinal layer structures on optical coherence tomography images,” Opt. Express 13(25), 10200–10216 (2005). [CrossRef] [PubMed]
37. K. Abd-Elmoniem, “Feedback coherent anisotropic diffusion for high resolution image enhancement”, in Proceedings of IEEE International Symposium on Biomedical Imaging, (Washington, DC, 2002), pp. 693–696.
38. A. Araújo, S. Barbeiro and P. Serranho, “Stability of finite difference schemes for complex diffusion processes”, DMUC report 10–23, (2010).
39. G. Gilboa, N. Sochen, and Y. Y. Zeevi, “Image enhancement and denoising by complex diffusion processes,” IEEE Trans. Pattern Anal. Mach. Intell. 26(8), 1020–1036 (2004). [CrossRef]
40. Matlab (Matlab – The MathWorks Inc., Natick, MA, USA). http://www.mathworks.com.