OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editor: Gregory W. Faris
  • Vol. 5, Iss. 8 — Jun. 8, 2010
« Show journal navigation

General Bayesian estimation for speckle noise reduction in optical coherence tomography retinal imagery

Alexander Wong, Akshaya Mishra, Kostadinka Bizheva, and David A. Clausi  »View Author Affiliations


Optics Express, Vol. 18, Issue 8, pp. 8338-8352 (2010)
http://dx.doi.org/10.1364/OE.18.008338


View Full Text Article

Acrobat PDF (1951 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

An important image post-processing step for optical coherence tomography (OCT) images is speckle noise reduction. Noise in OCT images is multiplicative in nature and is difficult to suppress due to the fact that in addition the noise component, OCT speckle also carries structural information about the imaged object. To address this issue, a novel speckle noise reduction algorithm was developed. The algorithm projects the imaging data into the logarithmic space and a general Bayesian least squares estimate of the noise-free data is found using a conditional posterior sampling approach. The proposed algorithm was tested on a number of rodent (rat) retina images acquired in-vivo with an ultrahigh resolution OCT system. The performance of the algorithm was compared to that of the state-of-the-art algorithms currently available for speckle denoising, such as the adaptive median, maximum a posteriori (MAP) estimation, linear least squares estimation, anisotropic diffusion and wavelet-domain filtering methods. Experimental results show that the proposed approach is capable of achieving state-of-the-art performance when compared to the other tested methods in terms of signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), edge preservation, and equivalent number of looks (ENL) measures. Visual comparisons also show that the proposed approach provides effective speckle noise suppression while preserving the sharpness and improving the visibility of morphological details, such as tiny capillaries and thin layers in the rat retina OCT images.

© 2010 Optical Society of America

1. Introduction

Speckle is an inherent characteristic of images acquired with any imaging technique that is based on detection of coherent waves, for example synthetic aperture radar (SAR), ultrasound, coherent optical imaging, etc. Speckle carries information about both the structure of the imaged object as well as a noise component, and the latter is responsible for the grainy appearance of the images. Optical coherence tomography (OCT) is an imaging technique capable of non-contact, high resolution (few micrometers), 3D imaging of the structure of optically semitransparent objects, including biological tissue. Since OCT is based on interferometric detection of coherent optical beams, OCT images contain speckle. Speckle in OCT tomograms is dependent on both the wavelength of the imaging beam and the structural details of the imaged object [1

1. J. Rogowska and M. E. Brezinski, “Evaluation of the adaptive speckle suppression filter for coronary optical coherence tomography imaging,” IEEE Trans. Med. Imaging 19, 1261–1266 (2000). [CrossRef]

].

With the development of ultrahigh resolution OCT (UHROCT), cellular level resolution can be achieved in biological tissue. This enables the visualization of small morphological details in the UHROCT tomograms such as individual tissue layers, small blood and lymph vessels, calcifications, lipid deposits, small clusters of highly specialized cells, etc. Such small structural details can easily be obscured by the presence of speckle noise in unprocessed UHROCT images, or by the blurring and / or image artefacts introduced by speckle denoising algorithms that are currently used in commercial and research grade OCT systems or have been published in the past. Furthermore, the development of high speed UHROCT technology permits fast acquisition of high density, large volume 3D images of biological tissue. The development of fast speckle noise reduction algorithms for UHROCT with very good preservation of boundaries of layered structure or small morphological features is of high importance. Such algorithms can improve the quality of the visual appearance of UHROCT images and allow for zooming on small features in the image without compromising the sharpness of the details or the overall image quality. Furthermore, speckle noise reduction algorithms can potentially improve the precision and overall performance of other image post-processing algorithms such as layer segmentation and pattern analysis that are often applied to UHROCT images for quantitative analysis. However, speckle noise reduction in OCT images is particularly challenging, because of the difficulty in separating the noise and the information components in the speckle pattern.

In addition to algorithmic approaches to speckle noise reduction, design changes to the OCT imaging technique have been proposed to reduce the presence of speckle noise on the imaging side. Kim et al. [18

18. J. Kim, D. Miller, E. Kim, S. Oh, J. Oh, and T. Milner, ”Optical coherence tomography speckle reduction by a partially spatially coherent source,” J. Biomed. Opt. 10, 640349 (2005).

] proposed the use of a partially spatially coherence broadband light source to reduce the presence of speckle noise. Pircher et al. [19

19. M. Pircher, E. Gtzinger, R. Leitgeb, A. F. Fercher, and C. K. Hitzenberger, “Speckle reduction in optical coherence tomography by frequency compounding,” J. Biomed. Opt. 8, 565–569 (2003). [CrossRef] [PubMed]

] proposed a frequency compounding method which makes use of two incoherent interferometric signals to reduce speckle noise while maintaining high spatial resolution. Iftimia et al. [20

20. 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, 260–263 (2003). [CrossRef] [PubMed]

] proposed an angular compounding approach based on path length encoding, which was also found to reduce the presence of speckle noise. More recently, Desjardins et al. [21

21. A. E. Desjardins, B. J. Vakoc, W. Y. Oh, S. M. R. Motaghiannezam, G. J. Tearney, and B. E. Bouma, “Angle-resolved Optical Coherence Tomography with sequential angular selectivity for speckle reduction,” Opt. Express 15, 6200–6209 (2007). [CrossRef] [PubMed]

] was able to achieve speckle noise reduction through the use of angular compounding at multiple backscattering angles. Finally, Jorgensen et al. [22

22. T. Jorgensen, L. Thrane, M. Mogensen, F. Pedersen, and P. Andersen, “Speckle reduction in optical coherence tomography images of human skin by a spatial diversity method,” in Proc. SPIE 6627, Munich, Germany 66270P (2007).

] proposed a spatial diversity approach, where the probe beams focal plane is shifted to achieve speckle reduction. Other spatial diversity based approaches [23–25

23. J. Schmitt, “Array detection for speckle reduction in optical coherence microscopy,” Phys. Med. Biol. 42(7), 1427–1439 (1997). [CrossRef] [PubMed]

], as well as polarization diversity based approaches [26

26. M. Kobayashi, H. Hanafusa, K. Takada, and J. Noda, “Polarization-Independent Interferometric Optical-Time-Domain Reflectometer,” J. Lightwave Technol. 9, 623–628 (1991). [CrossRef]

] have also been investigated. However, such methods require significant changes to the design of the OCT system, longer data acquisition time and in most cases the effect of speckle denoising is limited. Therefore, the development of algorithmic approaches to speckle noise reduction remains an important part of the OCT images post-processing.

Here we present a novel speckle noise reduction algorithm based on nonlinear log-space general Bayesian least squares estimation. The novelty of the algorithm stems from the conditional posterior sampling approach used to estimate the posterior distribution of the noise-free data in a non-parametric fashion. This allows the algorithm to learn the underlying noise distribution of the observed data dynamically to provide significant speckle noise suppression while preserving image details. When applied to UHROCT images, acquired in-vivo from rodent (rat) retinas, the novel algorithm shows significantly improved overall image quality and clear preservation of small structural features in the retinal images.

2. Theory

Let S be a set of sites on a discrete lattice 𝓛 which defines an OCT image and let sS be a site in 𝓛. Let M = {M(s)|sS}, A = {A(s)|sS}, and N = {N(s)|sS} be fields on S. Given the measured data M(s) that we have acquired, let A(s) and N(s) be random variables representing noise-free data and speckle noise of unknown distribution at site s respectively. Let m = {m(s)|sS}, a = {a(s)|sS}, and n = {n(s)|sS} be realizations of M, A, and N respectively. Speckle in OCT imagery arises from the constructive and destructive interference of the backscattered signal from biological issues [16

16. W. Drexler, “Ultrahigh-resolution optical coherence tomography,” J. Biomed. Opt. 9, 47–74 (2004). [CrossRef] [PubMed]

], and can be modeled as multiplicative noise that is dependent on the wavelength of the imaging beam and the structural composition of the imaged object [17

17. J. Schmitt, S. Xiang, and K. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4, 95–105 (1999). [CrossRef]

],

m(s)=a(s)·n(s).
(1)

The data-dependent nature of speckle noise, as well as the fact that the speckled noise distribution can vary depending on the optical design of the imaging system, makes the problem of separating the noise-free data a(s) from the speckle noise n(s) very challenging. We propose to tackle these issues associated with speckle noise reduction by estimating the noise-free data in the logarithm space using a general Bayesian least squares estimation approach based on conditional posterior sampling. The proposed speckle noise reduction algorithm can be described as follows. To handle the data-dependent nature of speckle noise, the noise-free data a(s) and the speckle noise n(s) are decoupled by projecting the measured data m(s) into the logarithm space,

logm(s)=ml(s)=log[a(s).n(s)]=log{a(s)}+log{n(s)}=al(s)+nl(s).
(2)

In the logarithm space, the general Bayesian least squares estimate of al(s) can be defined by the expression,

âl(s)=argminal(s){E[(al(s)âl(s))2ml(s)]}.
(3)

What the general Bayesian least squares estimate does is minimize the average squared error of the noise-free data estimate âl(s) based on the measured data ml(s). Minimizing the expression in Eq. (3) gives,

âl(s)=p(al(s)ml(s))al(s)dal(s).
(4)

What the expression in Eq. (4) shows is that the optimal Bayesian estimate of the noise-free data al(s) is essentially the statistical average based on the measured data ml(s). The posterior distribution p(al(s) |ml(s)), which represents the probability distribution of the noise-free data âl(s) based on the knowledge of the measured data ml(s), can be highly complicated and nonlinear in nature, making it difficult to solve for âl(s) using Eq. (4). Typically, simpler Bayesian linear least squares estimators [2

2. J. Lee, “Speckle suppression and analysis for synthetic aperture radar,” Opt. Eng. 25(5), 636–643 (1986).

] and estimators based on specific parametric posterior distribution models [6

6. A. Lopes, E. Nezry, R. Touzi, and H. Laur, “Structure detection and adaptive speckle filtering in SAR images,” Int. J. Remote Sens. 14(9), 1735–1758 (1993). [CrossRef]

] have been used instead for speckle noise reduction. However, given the complex nature of the OCT images, such estimators provide poor signal-noise separation and hence result in significant loss of structural detail. To address these issues, we propose to employ a conditional posterior sampling approach to estimate p(al(s)|ml(s)).

What conditional posterior sampling does is select information from the measured data ml(s) to estimate the posterior distribution p(al(s)|ml(s)) based on conditions that identify the relevance of that information to accurate estimation. For example, when estimating the posterior distribution of a pixel with strong boundary characteristics, conditional posterior sampling adaptively selects information from the measured data that has similar boundary characteristics to estimate the posterior distribution, as that would allow the boundary characteristics to be preserved while averaging out the speckle noise.

The conditional posterior sampling approach we used to estimate p(al(s)|ml(s)) can be described as follows. In Markov-Chain Monte Carlo density estimation [27

27. W. Hastings, “Monte carlo sampling methods using Markov chains and their applications,” Biometrika 57(1), 97–109 (1970). [CrossRef]

], an unknown target distribution (in our case, p(al(s)|ml(s))) is estimated in an indirect manner by first sampling from a known initial probability distribution Q. Similarly, in the proposed approach, a random site s′ in S is determined based on an initial probability distribution Q(s′|s). Q(s′|s) is defined as a Gaussian distribution centered at s,

Q(s|s)=12πσspatiale(s's22σspatial2),
(5)

where ∥s′-s2 denotes the squared Euclidean distance of a site s′ from s, and σspatial represents the spatial variance of the initial probability distribution Q(s′|s), which is set to 7 pixels as it was shown to provide accurate estimates during testing. Testing under different imaging conditions as well as different resolutions have shown that the use of σspatial = 7 allows for consistently accurate results, thus making that setting suitable for a wide range of imaging scenarios. The reason for the strong performance using this setting is that the algorithm is adapts to the underlying statistics of the image and as such performs well as long as the selected area to select samples from is large enough to obtain good statistics from, irrespective of the resolution of the image. This initial probability distribution Q(s′|s) generates sites that tend to be in closer proximity to site s. Given the drawn site s′, the inclusion of s′ as a realization of p(al(s)|ml (s)) is determined based on the condition,

μ(s)μ(s)<2σ,
(6)

where μ(s) is the local mean of the neighborhood centered at s and σ is the estimated image noise variance. The local mean μ(s) is computed in a 7 × 7 region centered at s in the current implementation, as that provides sufficient information to obtain good statistics from, irrespective of the resolution of the image. Furthermore, in the current implementation of the proposed method, the noise variance σ was estimated by taking the local variance within a 7 × 7 region for the same reason of obtaining good statistics. Equation (6) enforces the inclusion of s′ as a realization of p(al(s)|ml(s)) within two standard deviations, which accounts for the effect of noise variations. The inclusion of s′ as a realization of p(al(s)|ml(s)) is based on the assumption that the local mean is a reasonable initial estimate of al(s). This conditional sampling is repeated until the maximum number of sites used to estimate the original signal, denoted as γ, is drawn.

Given the set of γ sites drawn from Q(s′|s), denoted as Ω = {s1,s2,…,sγ}, the weight associated with each site si in estimating al(s), denoted as w(si|s), is computed using the following Gibbs-based likelihood function based on the local means of si and s,

w(s'i|s)=exp[|μ(s'i)μ(s)2σ2].
(7)

Equation (7) allows for a more reliable estimate of p(al(s)|ml(s)) by weighing sites with local means similar to the local mean of site s higher since they are more likely to be true realizations of p(al(s)|ml(s)).

Given a set of sites Ω = {s1,s2,…,sγ} and the associated set of weights W = {w(s1,s),w(s2,s),…,w(sγ, s)}, the posterior distribution p(al(s)|ml(s)) is then estimated using a weighted histogram approach. Suppose that the discrete range of possible measured data values (ml) be [Lmin, Lmax], where Lmin and Lmax are the minimum and maximum possible values, respectively. Let the discrete range of possible noise-free data values (al) also be [Lmin, Lmax]. Furthermore, let h(rk) be a weighted histogram, defined over [Lmin, Lmax], where rk is the kth possible noise-free data value. For each site si, the weight w(si|s) is accumulated in the histogram bin of the weighted histogram that corresponds to ml(si) (i.e., h(rk = ml(si))). After constructing the weighted histogram, each histogram bin is then divided by the sum of all weights to construct a normalized histogram representing p(al(s)|ml(s)). Therefore, (al(s)|ml(s)) can be formulated as

p̂(al(s)|ml(s))=k=Ωw(sk|s)δ(alml(sk))Z,
(8)

where δ(․) is the Dirac delta function and Z is a normalization term such that alp̂(al(s)ml(s))=1.

Finally, the estimate of a(s) can be found by back-projecting the Bayesian estimate of al(s) computed using the estimated p(al(s)|ml(s)) [Eq. (4)] from the logarithm space using â(s) = exp[âl(s)]. A flowchart detailing a step-by-step breakdown of the proposed despeckling algorithm is shown in Fig. 1.

3. Results and discussion

To evaluate the effectiveness of the proposed log-space general Bayesian estimation approach for speckle noise suppression in OCT tomograms, the method was applied to a set of UHROCT images of the rat retina acquired in-vivo with a research grade UHROCT system. Two representative images acquired at and away from the optic disc of the rat eye (Fig. 2) are discussed.

Fig. 1. Step-by-step flowchart of the proposed despeckling algorithm.

The images were acquired with a state-of-the art, research grade high speed, UHROCT system operating in the 1060nm wavelength region. A detailed description of the system can be found in [28

28. P. Puvanathasan, P. Forbes, Z. Ren, D. Malchow, S. Boyd, and K. Bizheva, “High-speed, high-resolution Fourier-domain optical coherence tomography system for retinal imaging in the 1060 nm wavelength region,” Opt. Lett. 33(21), 2479–2481 (2008). [PubMed]

]. Briefly, the UHROCT system is based on a spectral domain design powered with a super-luminescent diode (λc = 1020nm, δλ = 110nm, Pout = 10mW) and data is acquired with a 47kHz data rate, InGaAs linear array, 1024 pixel camera (SUI, Goodrich) interfaced with a high performance spectrometer (P&P Optica). The UHROCT system provides 3μm axial and 5μm lateral resolution in retinal tissue and 100dB SNR for 1.3mW optical power incident on the rat cornea.

2D and 3D images were acquired from the retinas of living Long Evans rats and the imaging procedure was carried out in accordance with the University of Waterloo ethics regulations. The two selected images, one acquired at the optic disc of the retina [Fig. 2(a)], and one 2mm away from the optic disc [Fig. 2(b)], containing different morphological features, were selected to test the performance of the speckle noise reduction algorithm. Both images show clear visualization of all intraretinal layers, as well as small morphological details such as the tiny capillaries ( 10μm in diameter) in the inner- and outer plexiform layers (red arrows) and the larger choroidal blood vessels (yellow arrows). For comparison purposes, the same images were also processed with some of the high performance wavelet- or diffusion-based image processing algorithms, such as the adaptive median filter proposed by Loupas et.al. [5

5. T. Loupas, W. Mcdicken, and P. Allen, “An adaptive weighted median filter for speckle suppression in medical ultrasound images,” IEEE Trans. Circuits Syst. 36(1), 129–135 (1989). [CrossRef]

], the linear least squares estimator proposed by Frost et al. [3

3. V. Frost, J. Stiles, K. Shanmugan, and J. Holtzman, “A model for radar images and its application to adaptive digital filtering for multiplicative noise,” IEEE Trans. Pattern Anal. Machine Intell. 4(2), 157–166 (1982). [CrossRef]

], the wavelet-based method proposed by Pizurica et al. [9

9. A. Pizurica, W. Philips, I. Lemahieu, and M. Acheroy, “A versatile wavelet domain noise filtration technique for medical imaging,” IEEE Trans. Med. Imag. 22(3), 323–331 (2003). [CrossRef]

], the MAP method proposed by Lopes at el. [6

6. A. Lopes, E. Nezry, R. Touzi, and H. Laur, “Structure detection and adaptive speckle filtering in SAR images,” Int. J. Remote Sens. 14(9), 1735–1758 (1993). [CrossRef]

], the anisotropic diffusion method proposed by Yu and Acton [12

12. Y. Yu and S. Acton, “Speckle reducing anisotropic diffusion,” IEEE Trans. Image Process. 11(11), 1260–1270 (2002). [CrossRef]

], and the Type II Fuzzy anisotropic diffusion method proposed by Puvanathasan and Bizheva [13

13. P. Puvanathasan and K. Bizheva, “Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images,” Opt. Express 17(2), 733–746 (2009). [CrossRef] [PubMed]

]. All tested methods were implemented in as computationally efficient a manner as possible, using the parameters proposed in the associated research literature. For the proposed method, the number of sites used to estimate the original signal at each site, denoted as γ, was set to 64 and 7 + 7 neighborhoods were used, as they were proven to produce good estimates of p(al(s)|ml(s)). Tests performed under different imaging conditions, resolutions, as well as biological tissues other than retina (e.g., corneal and skin) and the testing have shown that selection of these parameters work well for a wide variety of different imaging scenarios. The reason for the strong performance using such parameters is that the proposed algorithm adapts to the underlying statistics of the image and as such performs well as long as the selected area to select samples from is large enough to obtain good statistics from, irrespective of the resolution of the image. As such, the use of the presented parameters should provide strong speckle noise reduction performance for most practical situations. All algorithms were implemented in MATLAB and tested on an Intel Pentium 4, 3 GHz machine with 1 GB of RAM. For direct comparison of the algorithms performance, 3 regions of the retinal image [Fig. 2(b)] were selected (red-line boxes), focusing on specific morphological features such as the retinal capillaries and surface blood vessels (box #1), the boundaries between retinal layers (box #2), the choroidal blood vessels (box #3). Significant efforts were made to ensure fair comparisons between the proposed algorithm and other tested speckle denoising algorithms. Extensive parametric testing showed that the parameters proposed in the associated research literature for the tested algorithms provide the strongest results that can be obtained using these methods for the tested retinal images, as changes to these parameters yield no improvements in terms of the results.

Fig. 2. UHROCT images of the rat retina acquired near the optic disc (A) and away from the optic disc (B). NFL nerve fiber layer; GCL ganglion cell layer; IPL inner plexiform layer; INL inner nuclear layer; OPL outer plexiform layer; ONL outer nuclear layer; ELM external limiting membrane; IS inner segment; OS outer segment of the photoreceptor layer; RPE retinal pigmented epithelium; C- choroid and S sclera. The red arrows mark tiny capillaries imbedded in the retinal OPL, while the yellow arrows mark large blood vessels in the choroid. The red line boxes mark sections of the retinal image that were enlarged for more direct visual comparison of the performance of the speckle denoising algorithms (see Figs. 3,4 and 5).

Figure 3 shows qualitative, visual comparisons of the original and processed UHROCT retina images acquired away from the optic disc. All retinal layers, as well as small features, such as the tiny capillaries imbedded in the inner plexiform layer of the retina (red arrows) and choroidal blood vessels (yellow arrows) are visible in the image. Blue and red coloured, numbered boxes in the original images mark the regions of interest used in the quantitative comparison of the processed images described below. Although morphological features are visible in the images, presence of speckle noise causes the grainy appearance of the UHROCT unprocessed tomogram. As a result, segmentation algorithms applied to this image may fail to identify properly the boundaries of all retinal layers, specifically the thin retinal layers such as the external limiting membrane (ELM) and the inner plexiform layer (IPL), as well as to define well the inner boundaries of the choroidal blood vessels. Precise determination of choroidal and retinal blood vessels diameter is essential for the accurate calculation of retinal blood flow.

Visual comparison of the processed images reveals that the adaptive median filter (AMF) preserved image detail but provided limited noise suppression. The linear least squares estimation method (LLSQ) provided better noise suppression but at the cost of significant loss of image detail. The wavelet noise reduction methods provided both improved noise suppression and detail preservation, but introduced wavelet related artefacts (Fig. 4). The MAP and anisotropic diffusion (AD) methods provided good noise suppression and detail preservation with minimal appearance of image artefacts. In contrast, application of the proposed method resulted in noticeably improved detail preservation and noise suppression compared to the other tested methods. Areas with fairly homogeneous structure, for example the inside of the inner and outer nuclear layers, as well as part of the sclera [see regions marked with red numbered boxes in Fig. 2(b)] show excellent removal of the grainy appearance characteristic of the original image and most of the other post-processed images and associate with presence of speckle noise. While the speckle noise removal with the new algorithm results in blurring of the homogenous regions of the image, image features corresponding to actual morphological details in the retina remain very well preserved and show improved visibility (contrast).

To emphasize the difference in edge preservation and contrast improvement between the proposed method and the other algorithms, three regions in the original image [Fig. 2(b)] were selected and marked with red line, numbered boxes, and enlarged views of these regions are shown in Figs. 4, 5, and 6. These regions were chosen to contain image details such as very thin retinal layers, tiny capillaries and cross-sections of larger retinal and choroidal blood vessels, in order to evaluate the performance of the individual algorithms under different conditions.

Figure 4 shows a section of the original image [Fig. 2(b)] marked with red box #1, containing a cross-section of a large blood vessel imbedded in the NFL, and the cross-sections of three capilaries positioned at the interface between the IPL and INL. The original inset shows that the presence of speckle noise reduces the image quality to the extent that the cross-section of the second (middle) capillary is barely visible, while the outlines of the rest of the blood vessels appear blurred. The application of the AMF algorithm does not improve image quality significantly, while the LLSQ method blurs the image to the extent that the outlines of blood vessel cross-sections cannot be easily distinguished. Both the AD and wavelet methods improve the overall contrast of the image and the visibility of the blood vessel cross-sections, however these methods introduce significant artefacts in the processed images that blur the contours of the blood vessel cross-sections. Both the MAP and the Type II Fuzzy AD approaches appear to significantly improve the contrast and visibility of the blood vessel cross-sections, but some residual blur is still present in the contours of the cross-sections.

The proposed algorithm removes speckle noise from the homogeneous regions of the IPL and INL layers and greatly improves the contrast and visibility of the image features. Furthermore, the edge preservation of this method is effective, as the outlines of all blood vessels appear sharp and well defined.

Fig. 3. An OCT image of the rat retina acquired away from the optic disc and processed with the following filters: Original image, Adaptive median filter, Linear least square estimation, Anisotropic Diffusion, MAP estimation, general Wavelet, Type II Fuzzy AD, and the new proposed algorithm. The blue and red-line boxes in the original image mark regions in the image used for quantitative comparison of the performance of all image processing algorithms applied to the original image.

Figure 5 shows a region of the original image [Fig. 2(b)] marked with red box #2, containing sections of the RPE layer (thin black line) and the IS / OS portion of the photoreceptor layer (pale and dark gray alternating layers). The original inset shows that the presence of speckle noise reduces the overall image quality and blurs the boundaries of the individual layers. The application of other tested algorithms causes significant blur of the layers boundaries and in the case of the AD and Wavelet methods, introduces image artefacts. Overall, the best performance in terms of improved image contrast and visual appearance, and best edge preservation with practically no image artefacts, for the case of layered structures in the UHROCT image is achieved with the proposed algorithm.

Figure 6 shows the performance of the different algorithms in terms of improving the appearance of choroidal blood vessels to allow for more precise segmentation and measurement of the blood vessel diameter. For this purpose, we have selected a region in the original image [Fig. 2(b)] marked with red box #3, that contains a part of the blood vessel longitudinal cross-section. The presence of speckle noise in the original image in Fig. 6 blurs the contours of the blood vessel to the extent that precise determination of the blood vessel diameter is not possible. The application of the AMF algorithm does not improve image quality significantly. The LLSQ, MAP and Fuzzy type II AD algorithms result in significantly improved contrast between the blood vessels and the surrounding choroidal tissue, at the expense of significant blur of the blood vessel contours, which is most extreme in the case of the LLSQ method. Although the Wavelet approach results in the best contrast improvement, it introduces image artefacts that can compromise the overall visual appearance of the processed image. Both significant blur and image artefacts are present in the image processed with the AD method. The proposed approach generates an optimal combination of improved image contrast, visual appearance and sharp boundaries of the blood vessel cross-section.

Figure 7 shows qualitative, visual comparison of the original and processed OCT retina images acquired at the optic disc of the rat retina [Fig. 2(a)]. Again, all intra-retinal layers, as well as small features, such as tiny capillaries imbedded in the inner plexiform layer of the retina, and choroidal blood vessels are visible in the image. The cross-section of the remains of the hyaloids blood vessel is also visible. Blue and red colored, numbered boxes mark the regions of interest used in the quantitative comparison of the processed images described below. In addition to improving the visual appearance of the layered retinal structure, as in Fig. 2, here the novel speckle noise reduction algorithm results in better visualization of the outline of the optic disc funnel, as well as clearly defined inner and outer boundaries of the hyaloid blood vessel. The interface between the retinal nerve fiber layer (RNF) and the IPL is also significantly sharper with better contrast as compared to the original and the rest of the processed images. The sharper appearance of the optic disc funnel and overall layered retinal structure can aid and significantly improve the performance of automatic segmentation algorithms targeting to evaluate the RNF thickness and optic disc diameter in healthy and diseased retinal images. Quantitative analysis of the RNFL thickness and the dimensions of the optic disk are used in clinical studied for improved and early diagnostics of glaucoma.

To compare the tested methods in a quantitative manner, the average SNR, contrast-to-noise ratio (CNR), equivalent number of looks (ENL), and edge preservation (η) over the in-vivo OCT retinal images were computed for each tested method. The image quality metrics used are the same as the metrics used in [7

7. 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]

, 10

10. 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]

, 13

13. P. Puvanathasan and K. Bizheva, “Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images,” Opt. Express 17(2), 733–746 (2009). [CrossRef] [PubMed]

], and are defined as follows:

SNR=10log10[max(A2)/σ2],
(9)
ENL=1H[h1Hμh2σh2],
(10)
Fig. 4. 3X enlargement of the region in the original OCT image marked with the red box #1 in Fig. 1(b) and processed with the following filters: AMF (adaptive median filter); LLSQ (linear least square estimation); AD (adaptive diffusion); general Wavelet; MAP (maximum a posteriori estimation); Type II fuzzy AD (fuzzy rules controlled adaptive diffusion). Features in the original image such as the large blood vessel in the NFL and the three capillaries positioned at the boundary between the IPL and INL appear very sharp and distinct on the image processed with the proposed algorithm. In contrast, the same image features appear with low contrast and blurred boundaries on the images processed with the rest of the algorithms.
CNR=1R[r=1R(μrμb)σr2+σb2],
(11)
η=(2M2M¯)(2A2A¯)(2M2M¯)2.(2A2A¯)2,
(12)

ENL measures the smoothness of a homogenous region of interest, while CNR measure the difference between an area of image feature and an area of background noise. Edge preservation is a correlation measure that indicates how the edges in the image are degrading. A value close to 1 indicates the filtered image is similar to the reference image. In the expression for SNR, A and σ2 represent the linear magnitude image and the variance of the background noise region in the linear magnitude image respectively. In the expression for ENL, μh and σh2 represent the mean and the variance of the hth homogenous region of interests respectively. In the definition for CNR, μb and σb2 represent the mean and the variance of the same background noise region as in SNR and μr and σr2 represent the mean and the variance of the rth region of interest which includes the homogeneous regions. In the edge preservation measure, ∇2M and ∇2A represent the Laplacian operator performed on the original image and the filtered image respectively. Also, 2M¯ and 2A¯ represent the mean value of a 3 × 3 neighbourhood around the pixel location of ∇2M and ∇2A respectively.

Fig. 5. 3X enlargement of the region in the original OCT image marked with the red box #2 in Fig. 1(b) and processed with the following filters: AMF (adaptive median filter); LLSQ (linear least square estimation); AD (adaptive diffusion); general Wavelet; MAP (maximum a posteriori estimation); Type II fuzzy AD (fuzzy rules controlled adaptive diffusion). The RPE, IS and OS of the photoreceptor layers of the retina appear very sharp and distinct on the image processed with the proposed novel algorithm. In contrast, the same image features appear with low contrast and blurred boundaries on the images processed with the rest of the algorithms.

Fig. 6. Magnified view of a choroidal blood vessel, corresponding to the area marked with red box #3 in the original image [Fig. 1(b)] and processed with the different speckle reduction algorithms. AMF adaptive median filter; LLSQ linear least square estimation; AD anisotropic diffusion; general Wavelet; MAP maximum a posteriori estimation; Type II fuzzy AD (fuzzy rules controlled adaptive diffusion). The proposed method results in very clear delineation of the blood vessel walls as compared to the significant blur and / or image artefacts introduced by the other speckle denoising methods.

Table 1. Image quality metrics evaluated for the rat retina image [Fig. 2(b)]. Values are relative to the original image.

table-icon
View This Table

4. Conclusion

Fig. 7. An OCT image of the rat retina acquired at the optic disc and processed with the following filters: Original image, Adaptive median filter, Linear least square estimation, Anisotropic Diffusion, MAP estimation, general Wavelet, Type II Fuzzy AD, and the new proposed algorithm. The blue and red boxes in the original image mark regions in the image used for quantitative comparison of the performance of all image processing algorithms applied to the original image.

Acknowledgment

The authors would like to thank S. Hariri, A.A. Moayed, C. Hyun and S. Shackheel for assistance with the acquisition of the rat retina OCT images. This project was supported in part by NSERC, OCE and University of Waterloo research grants.

References and links

1.

J. Rogowska and M. E. Brezinski, “Evaluation of the adaptive speckle suppression filter for coronary optical coherence tomography imaging,” IEEE Trans. Med. Imaging 19, 1261–1266 (2000). [CrossRef]

2.

J. Lee, “Speckle suppression and analysis for synthetic aperture radar,” Opt. Eng. 25(5), 636–643 (1986).

3.

V. Frost, J. Stiles, K. Shanmugan, and J. Holtzman, “A model for radar images and its application to adaptive digital filtering for multiplicative noise,” IEEE Trans. Pattern Anal. Machine Intell. 4(2), 157–166 (1982). [CrossRef]

4.

D. Kuan, A. Sawchuk, T. Strand, and P. Chavel, “Adaptive restoration of images with speckle,” IEEE Trans. Acoust. Speech Signal Process. 35(3), 373–383 (1987). [CrossRef]

5.

T. Loupas, W. Mcdicken, and P. Allen, “An adaptive weighted median filter for speckle suppression in medical ultrasound images,” IEEE Trans. Circuits Syst. 36(1), 129–135 (1989). [CrossRef]

6.

A. Lopes, E. Nezry, R. Touzi, and H. Laur, “Structure detection and adaptive speckle filtering in SAR images,” Int. J. Remote Sens. 14(9), 1735–1758 (1993). [CrossRef]

7.

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]

8.

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]

9.

A. Pizurica, W. Philips, I. Lemahieu, and M. Acheroy, “A versatile wavelet domain noise filtration technique for medical imaging,” IEEE Trans. Med. Imag. 22(3), 323–331 (2003). [CrossRef]

10.

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]

11.

S. Aja, C. Alberola, and J. Ruiz, “Fuzzy anisotropic diffusion for speckle filtering,” Proc. IEEE ICASSP 2, 1261–1264 (2001).

12.

Y. Yu and S. Acton, “Speckle reducing anisotropic diffusion,” IEEE Trans. Image Process. 11(11), 1260–1270 (2002). [CrossRef]

13.

P. Puvanathasan and K. Bizheva, “Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images,” Opt. Express 17(2), 733–746 (2009). [CrossRef] [PubMed]

14.

D. Fernandez, H. Salinas, and C. Puliafito, “Automated detection of retinal layer structures on optical coherence tomography images,” Opt. Express 13, 10200–10216 (2005). [CrossRef]

15.

K. Yung, S. Lee, and J. Schmitt, “Phase-Domain Processing of Optical Coherence Tomography Images,” J. Biomed. Opt. 4(1), 125–136 (1999). [CrossRef]

16.

W. Drexler, “Ultrahigh-resolution optical coherence tomography,” J. Biomed. Opt. 9, 47–74 (2004). [CrossRef] [PubMed]

17.

J. Schmitt, S. Xiang, and K. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4, 95–105 (1999). [CrossRef]

18.

J. Kim, D. Miller, E. Kim, S. Oh, J. Oh, and T. Milner, ”Optical coherence tomography speckle reduction by a partially spatially coherent source,” J. Biomed. Opt. 10, 640349 (2005).

19.

M. Pircher, E. Gtzinger, R. Leitgeb, A. F. Fercher, and C. K. Hitzenberger, “Speckle reduction in optical coherence tomography by frequency compounding,” J. Biomed. Opt. 8, 565–569 (2003). [CrossRef] [PubMed]

20.

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, 260–263 (2003). [CrossRef] [PubMed]

21.

A. E. Desjardins, B. J. Vakoc, W. Y. Oh, S. M. R. Motaghiannezam, G. J. Tearney, and B. E. Bouma, “Angle-resolved Optical Coherence Tomography with sequential angular selectivity for speckle reduction,” Opt. Express 15, 6200–6209 (2007). [CrossRef] [PubMed]

22.

T. Jorgensen, L. Thrane, M. Mogensen, F. Pedersen, and P. Andersen, “Speckle reduction in optical coherence tomography images of human skin by a spatial diversity method,” in Proc. SPIE 6627, Munich, Germany 66270P (2007).

23.

J. Schmitt, “Array detection for speckle reduction in optical coherence microscopy,” Phys. Med. Biol. 42(7), 1427–1439 (1997). [CrossRef] [PubMed]

24.

M. Bashkansky and J. Reintjes, “Statistics and reduction of speckle in optical coherence tomography,” Opt. Lett. 25, 545–547 (2000). [CrossRef]

25.

D. Popescu, M. Hewkoa, and M. Sowa, “Speckle noise attenuation in optical coherence tomography by compounding images acquired at different positions of the sample,” Opt. Comm. 1(1), 247–251 (2007). [CrossRef]

26.

M. Kobayashi, H. Hanafusa, K. Takada, and J. Noda, “Polarization-Independent Interferometric Optical-Time-Domain Reflectometer,” J. Lightwave Technol. 9, 623–628 (1991). [CrossRef]

27.

W. Hastings, “Monte carlo sampling methods using Markov chains and their applications,” Biometrika 57(1), 97–109 (1970). [CrossRef]

28.

P. Puvanathasan, P. Forbes, Z. Ren, D. Malchow, S. Boyd, and K. Bizheva, “High-speed, high-resolution Fourier-domain optical coherence tomography system for retinal imaging in the 1060 nm wavelength region,” Opt. Lett. 33(21), 2479–2481 (2008). [PubMed]

OCIS Codes
(030.6140) Coherence and statistical optics : Speckle
(100.0100) Image processing : Image processing
(100.2980) Image processing : Image enhancement
(170.4500) Medical optics and biotechnology : Optical coherence tomography
(100.3008) Image processing : Image recognition, algorithms and filters

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: January 21, 2010
Revised Manuscript: March 24, 2010
Manuscript Accepted: March 30, 2010
Published: April 6, 2010

Virtual Issues
Vol. 5, Iss. 8 Virtual Journal for Biomedical Optics

Citation
Alexander Wong, Akshaya Mishra, Kostadinka Bizheva, and David A. Clausi, "General Bayesian estimation for speckle noise reduction in optical coherence tomography retinal imagery," Opt. Express 18, 8338-8352 (2010)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-18-8-8338


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. J. Rogowska and M. E. Brezinski, "Evaluation of the adaptive speckle suppression filter for coronary optical coherence tomography imaging," IEEE Trans. Med. Imaging 19, 1261-1266 (2000). [CrossRef]
  2. J. Lee, "Speckle suppression and analysis for synthetic aperture radar," Opt. Eng. 25(5), 636-643 (1986).
  3. V. Frost, J. Stiles, K. Shanmugan, and J. Holtzman, "A model for radar images and its application to adaptive digital filtering for multiplicative noise," IEEE Trans. Pattern Anal. Machine Intell. 4(2), 157-166 (1982). [CrossRef]
  4. D. Kuan, A. Sawchuk, T. Strand, and P. Chavel, "Adaptive restoration of images with speckle," IEEE Trans. Acoust. Speech Signal Process. 35(3), 373-383 (1987). [CrossRef]
  5. T. Loupas, W. Mcdicken, and P. Allen, "An adaptive weighted median filter for speckle suppression in medical ultrasound images," IEEE Trans. Circuits Syst. 36(1), 129-135 (1989). [CrossRef]
  6. A. Lopes, E. Nezry, R. Touzi, and H. Laur, "Structure detection and adaptive speckle filtering in SAR images," Int. J. Remote Sens. 14(9), 1735-1758 (1993). [CrossRef]
  7. 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]
  8. 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]
  9. A. Pizurica, W. Philips, I. Lemahieu, and M. Acheroy, "A versatile wavelet domain noise filtration technique for medical imaging," IEEE Trans. Med. Imag. 22(3), 323-331 (2003). [CrossRef]
  10. 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]
  11. S. Aja, C. Alberola, and J. Ruiz, "Fuzzy anisotropic diffusion for speckle filtering," Proc. IEEE ICASSP 2, 1261-1264 (2001).
  12. Y. Yu and S. Acton, "Speckle reducing anisotropic diffusion," IEEE Trans. Image Process. 11(11), 1260-1270 (2002). [CrossRef]
  13. P. Puvanathasan and K. Bizheva, "Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images," Opt. Express 17(2), 733-746 (2009). [CrossRef] [PubMed]
  14. D. Fernandez, H. Salinas, and C. Puliafito, "Automated detection of retinal layer structures on optical coherence tomography images," Opt. Express 13, 10200-10216 (2005). [CrossRef]
  15. K. Yung, S. Lee, and J. Schmitt, "Phase-Domain Processing of Optical Coherence Tomography Images," J. Biomed. Opt. 4(1), 125-136 (1999). [CrossRef]
  16. W. Drexler, "Ultrahigh-resolution optical coherence tomography," J. Biomed. Opt. 9, 47-74 (2004). [CrossRef] [PubMed]
  17. J. Schmitt, S. Xiang, and K. Yung, "Speckle in optical coherence tomography," J. Biomed. Opt. 4, 95-105 (1999). [CrossRef]
  18. J. Kim, D. Miller, E. Kim, S. Oh, J. Oh, and T. Milner, "Optical coherence tomography speckle reduction by a partially spatially coherent source," J. Biomed. Opt. 10, 640349 (2005).
  19. M. Pircher, E. Gtzinger, R. Leitgeb, A.F. Fercher, and C. K. Hitzenberger, "Speckle reduction in optical coherence tomography by frequency compounding," J. Biomed. Opt. 8, 565-569 (2003). [CrossRef] [PubMed]
  20. 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, 260-263 (2003). [CrossRef] [PubMed]
  21. A. E. Desjardins, B. J. Vakoc, W. Y. Oh, S. M. R. Motaghiannezam, G. J. Tearney, and B. E. Bouma, "Angle resolved Optical Coherence Tomography with sequential angular selectivity for speckle reduction," Opt. Express 15, 6200-6209 (2007). [CrossRef] [PubMed]
  22. T. Jorgensen, L. Thrane, M. Mogensen, F. Pedersen, and P. Andersen, "Speckle reduction in optical coherence tomography images of human skin by a spatial diversity method," in Proc. SPIE 6627, Munich, Germany 66270P (2007).
  23. J. Schmitt, "Array detection for speckle reduction in optical coherence microscopy," Phys. Med. Biol. 42(7), 1427-1439 (1997). [CrossRef] [PubMed]
  24. M. Bashkansky and J. Reintjes, "Statistics and reduction of speckle in optical coherence tomography," Opt. Lett. 25, 545-547 (2000). [CrossRef]
  25. D. Popescu, M. Hewkoa, and M. Sowa, "Speckle noise attenuation in optical coherence tomography by compounding images acquired at different positions of the sample," Opt. Commun. 1(1), 247-251 (2007). [CrossRef]
  26. M. Kobayashi, H. Hanafusa, K. Takada, and J. Noda, "Polarization-Independent Interferometric Optical-Time-Domain Reflectometer," J. Lightwave Technol. 9, 623-628 (1991). [CrossRef]
  27. W. Hastings, "Monte Carlo sampling methods using Markov chains and their applications," Biometrika 57(1), 97-109 (1970). [CrossRef]
  28. P. Puvanathasan, P. Forbes, Z. Ren, D. Malchow, S. Boyd, and K. Bizheva, "High-speed, high-resolution Fourier domain optical coherence tomography system for retinal imaging in the 1060 nm wavelength region," Opt. Lett. 33(21), 2479-2481 (2008). [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited