Stochastic speckle noise compensation in optical coherence tomography using non-stationary spline-based speckle noise modelling

Andrew Cameron; Dorothy Lui; Ameneh Boroomand; Jeffrey Glaister; Alexander Wong; Kostadinka Bizheva

doi:10.1364/BOE.4.001769

1. Introduction

Optical Coherence Tomography (OCT) is an interference-based imaging technology that is capable of imaging biological tissue with micrometer scale resolution and in a non-invasive manner. With the advent of ultrahigh resolution OCT (UHROCT) imaging, it is now possible to visualize biological tissue at the cellular level up to depths of approximately 1 mm below the tissue surface, thus allowing morphological details such as individual tissue layers, clusters of specialized cells, small blood vessels, and lipid deposits to be studied and characterized morphometrically in vivo. Furthermore, the data acquisition speed of OCT systems has greatly improved over the past decade with the development of high-speed CCD and CMOS [1] cameras and tunable lasers with MHz sweep rates [2], thus allowing real-time imaging and display of large volumes of biological tissue. All these developments have opened the door for a variety of medical imaging applications such as in vivo imaging of the human eye, gastrointestinal tract, as well as the epidermis and dermis layers of the skin [3].

As an interferometric method based on coherent optical beams, one of the fundamental challenges with OCT imaging is the presence of speckle noise in the tomograms. The characteristics of the speckle pattern in OCT imagery depend not only on the wavelength and coherent properties of the imaging beam, but also on the underlying structural details of the imaged object [4]. All of these factors lead to a speckle pattern with non-stationary characteristics that vary spatially within the OCT imagery and carries both information regarding the structure of the imaged object, as well as a noise component. The presence of speckle noise can often obscure small but important morphological details in the reconstructed OCT images, resulting in an overall grainy appearance that makes it more difficult for clinicians and clinical scientists to distinguish between different types of tissues. This challenge becomes more important in the case of high-speed OCT imaging systems, where there is a fundamental trade-off between acquisition speed and signal-to-noise ratio. As such, the development of speckle noise compensation strategies is of great importance in obtaining high quality, high resolution OCT imagery of biological tissue.

A number of methods exist for speckle noise compensation in OCT imaging, and can be grouped into two main categories: i) hardware compensation, and ii) algorithmic compensation. In general, hardware speckle noise compensation methods are based on the acquisition of multiple frames with uncorrelated noise that can be averaged together to improve the signal-to-noise ratio (SNR). Schmitt demonstrated a method using array detectors to measure the same signal across different spatial coordinates such that each detector element measures light backscattered from the sample at a different angle [5]. This “angular compounding” approach resulted in several other methods, including the sample beam angle adjustment over time [6], splitting of the incident beam and sampling simultaneously at different incident beam angles [7], or use of mirrors to change the incident beam angle [8]. Other methods use a “spatial diversity” approach, for example varying the position of the imaging lens relative to the sample [9], or varying the position of the sample itself for small samples [10]. Fang et al achieved competitive speckle compensation by using custom scan patterns to image a small fraction of a volume at a high SNR, then learning sparse representation dictionaries from the high-SNR scans and denoising the low-SNR remainder of the volume using these sparse dictionaries [11]. Major disadvantages inherent to any hardware speckle noise compensation approaches include longer data acquisition times, a more complex data acquisition procedure, as well as increased complexity of the opto-mechanical design of the OCT imaging system.

Given the technological challenges with hardware speckle noise compensation methods, over the past two decades significant effort has been invested in the development of algorithmic speckle noise compensation methods which require no hardware modifications nor increased acquisition time. Classical algorithmic approaches for speckle noise compensation such as linear least-squares estimation [12–14] and adaptive median filtering [15] were found to not only provide poor speckle noise compensation, but importantly lead to the loss of fine morphological details that are important for clinical diagnosis and scientific analysis [16]. More advanced data-adaptive speckle noise compensation methods [4, 17] have shown improved contrast and preservation of detail compared to the classical approaches, but still suffer from significant loss in detail under high speckle noise contamination situations.

More recently, much of the focus in algorithmic speckle noise compensation has revolved around wavelet-based strategies [18–22] and diffusion-based methods [20, 23–25]. In wavelet-based strategies, the underlying OCT imagery is decomposed across multiple scales using wavelet transforms, and thresholding is applied to the wavelet coefficients which contribute to the appearance of speckle noise in the imagery. Wavelet-based methods have superior speckle noise compensation capabilities when compared to earlier speckle noise compensation strategies. However, such methods may introduce wavelet-related artifacts that can obscure fine morphological details that are important for clinical diagnosis and scientific study [16]. In diffusion-based methods, the speckle-free OCT imagery is estimated by optimizing a diffusion-based partial differential equation (PDE) formulation that adapts non-linearly based on the underlying image detail. Diffusion-based methods have also been shown to provide strong speckle noise compensation capabilities compared to earlier strategies and do not introduce the type of artifacts that may be created by the wavelet-based methods. However, such methods have shown noticeable loss in detail in situations characterized by high speckle noise contamination.

One important trait of speckle noise patterns in OCT imagery of biological tissue that has not been adequately investigated is the non-stationary nature of the speckle noise characteristics that varies spatially within the OCT imagery. Several speckle-specific image restoration methods [26, 27] use sample variance within the local neighborhood to estimate spatially-varying noise variance. More general approaches, such as the Gaussian Scale Mixture model [28, 29] or Non-Local Means [30], use a single global estimate for noise variance. As such, the aim of the proposed method is to further explore the effectiveness of non-stationary speckle noise modeling for the purpose of speckle noise compensation in OCT imaging.

In this study, we investigate and introduce a novel stochastic approach for speckle noise compensation in OCT via non-stationary spline-based speckle noise modeling. Based on the non-stationary speckle noise model learned from the underlying OCT imaging data, a stochastic speckle noise compensation strategy is developed based on a spatially-adaptive Monte Carlo sampling strategy to estimate the noise-free OCT imagery. In this paper, we will denote the proposed method as the Non-stationary Speckle Compensation (NSC) method. The NSC method bears some similarities to a previous method for speckle noise suppression in OCT [16]; notably, both are Bayesian approaches to denoising speckle-contaminated images using Monte Carlo sampling. However, the NSC method makes several contributions:

Non-stationary speckle noise modeling: while the previous method assumes stationary speckle characteristics, the proposed method assumes that speckle variance is spatially-varying and estimates it using local Median Absolute Deviation, found to have superior noise estimation performance [31].
Stochastic acceptance-rejection sampling: while the previous method accepts randomly drawn samples based on a fixed criterion, the NSC method introduces a novel stochastic acceptance-rejection sampling strategy which accepts samples probabilistically based on the local speckle characteristics, thus better adapting to the underlying image information.
Local structure importance-weighted estimation: while the previous method accepts or rejects samples based on first-order statistics, the NSC approach incorporates more local structure information in its importance-weighting function, thus improving estimation of the posterior.

Fessler proposed a method for tomographic reconstruction using spline smoothing, however the proposed method differs by estimating the non-stationary noise variance (rather than relying on calibration frames) [32]. Furthermore, the method proposed by Fessler uses splines to enforce smoothness on image intensity values, while the proposed method uses splines to enforce smoothness on local variance estimates.

The rest of this paper is organized as follows: in Section 2, we explain the problem formulation and methodology of our approach. In Section 3, we explain the protocol used to acquire testing data and the experimental procedure for each of the three experiments conducted. In Section 4, experimental results using the acquired UHROCT imaging data of the human eye retina and limbus are presented and discussed. Finally, conclusions are drawn and potential future avenues of improvements are discussed in the Discussion.

2. Methodology

2.1. Problem formulation

The speckle noise found in OCT images, which is due to scattering of the laser source, is multiplicative in nature, and can be described by

G (\underline{x}) = F (\underline{x}) N (\underline{x}),

where G represents the measurements taken over the domain of x̱, F represents the true reflectance of the sample, and N denotes the speckle noise. Because different tissue types can have very different backscattering properties, some regions naturally produce more speckle. As such, N should be treated as a non-stationary noise process in x̱. By taking the logarithm of the measured data G, the multiplicative relationship between F and N becomes additive, and Eq. (1) becomes

log G (\underline{x}) = log F (\underline{x}) + log N (\underline{x}) .

Looking at Eq. (2), one can reduce the complexity of estimating the true reflectance of the sample by estimating logF(x̱) rather than F(x̱) itself, and then taking the exponential of the logarithmic estimate.

Based on Eq. (2), one can treat the problem of estimating the true reflectance of the sample as an inverse problem, and tackle the problem from a statistical perspective by formulating it as a Bayesian least-squares estimation problem, [33],

\begin{array}{l} log \hat{F} (\underline{x}) & = \underset{log F (\underline{x})}{argmin E} ({(log F (\underline{x}) - log \hat{F} (\underline{x}))}^{2} | log G (\underline{x})) \\ = \underset{log F (\underline{x})}{argmin} \int {(log F (\underline{x}) - log \hat{F} (\underline{x}))}^{2} p (log F (\underline{x}) | log G (\underline{x})) d (log F (\underline{x})) . \end{array}

One issue with the formulation in Eq. (3) is that it requires knowledge of the posterior p(logF(x̱)|logG(x̱)), which is intractable to obtain analytically in this case. To tackle this issue, we estimate the posterior p(logF(x̱)|logG(x̱)) using a non-parametric, spatially-adaptive Monte Carlo sampling method. By adapting the posterior estimation process using the novel spline-based speckle model introduced herein, we can better capture the spatially-varying speckle noise characteristics in OCT imagery to aid in the speckle noise compensation process.

2.2. Non-stationary spline-based speckle noise model

Since the speckle noise process N is a non-stationary noise process due to the varying tissue backscattering properties, one strategy to deal with this is to construct a non-stationary speckle noise model over the domain of x̱. While the speckle noise process N is non-stationary as a whole across the entire set of measurements, it can be reasonably assumed to be stationary within a small region of interest. As such, one can obtain an initial estimate of the speckle noise variance at a pixel x̱ based on the local noise statistics. Furthermore, it can be observed that pixels within very close proximity of each other have similar speckle noise characteristics, given that neighboring pixels within an image typically capture connected tissues with similar tissue backscattering properties. As such, one may want to ensure that the speckle noise variances estimated for neighboring pixels are consistent. Motivated by these insights and observations, we introduce a spline-based non-stationary speckle model for characterizing speckle noise variances over the domain of x̱.

The proposed spline-based non-stationary speckle model can be constructed as follows. First, to estimate local variances, we compute the Median Absolute Deviation (MAD) within the neighborhood ℵ_N of each pixel by taking the difference of the pixel values with the median pixel value in the neighborhood, and then we find the median of the absolute values of these differences. The MAD can be used to estimate the standard deviation in the log space by multiplying by a constant:

σ (\underline{x}) = 1.4826 {median}_{ℵ_{N}} (| log G ({\underline{x}}_{i}) - {median}_{ℵ_{N}} (log G) |) .

The use of MAD for estimating local variance is robust to the presence of outliers, and was shown in a comprehensive analysis of noise estimation strategies to be the most effective at characterizing local noise variance [31].

Following this, we would like to arrive at an initial estimate of the speckle noise variance at pixel x̱ given the set of computed local variances within a neighborhood ℵ_S. To help determine the speckle noise variance based on the set of computed local variances, we observe that the local variances that best approximates the noise variance are those computed in areas that are largely homogeneous. Furthermore, we observe that, in general, the most frequently occurring areas within an image are such homogeneous areas. Motivated by these observations, for each pixel x̱, we compute the initial noise variance estimate M(x̱) as the mode of the set of computed local variances in the neighborhood ℵ_S:

M (\underline{x}) = {mode}_{ℵ_{S}} (σ)

Based on empirical testing, the neighborhoods ℵ_N and ℵ_S were set to 9×9 and 15×15, respectively, for strong speckle noise variance estimation performance.

To enforce smoothness in the speckle noise variance such that pixels within very close proximity of each other have similar speckle noise characteristics, a cubic spline model is fitted using the initial speckle noise estimates M(x̱) to arrive at the final speckle noise variance estimate S(x̱). Therefore, the spline-based speckle noise variance estimate S(x̱) can be defined as

S (\underline{x}) = \underset{S (\underline{x})}{argmin} p \sum_{j = 1}^{n} {(M ({\underline{x}}_{j}) - S ({\underline{x}}_{j}))}^{2} + (1 - p) \int {‖ D^{2} S (\underline{x}) ‖}^{2} d \underline{x},

where p is a smoothing parameter and D² is the second derivative operator. By varying p, the resulting spline will vary between a least-squares straight line fit (p = 0) and a natural cubic spline interpolation (p = 1). Empirical testing led to a choice of 0.004 for p, which provided a good balance between data fidelity and estimation smoothness enforcement.

2.3. Spatially-adaptive Monte Carlo posterior estimation

Based on the spline-based non-stationary speckle noise model S, a spatially-adaptive Monte Carlo approach is then introduced to estimate the posterior p(logF(x̱)|logG(x̱)). We extend upon the concept of importance-weighted Monte Carlo sampling [34] with a spatially-adaptive approach so that the number of samples taken increases in regions with high speckle noise variance (to improve noise suppression) and decreases in areas where there is less speckle noise (to improve detail preservation).

Given the location x̱₀ of a pixel of interest, we use a neighborhood around x̱₀ as the search space, sampling new pixels uniformly at various values of x̱_i. The neighborhoods around these samples are examined to determine their statistical similarity to the neighborhood around x̱₀ before the samples are accepted. All the accepted samples and their importance weights are denoted by Ω. Intuitively, we want to accept a sample with high statistical likelihood of being a realization of the posterior p(logF(x̱)|logG(x̱)). The probability of accepting a sample into Ω is α(x̱_i|x̱₀), and is calculated by examining the neighboring pixels around each of x̱₀ and x̱_i according to

α ({\underline{x}}_{i} | {\underline{x}}_{0}) = \frac{Π_{j} {(2 π S ({\underline{x}}_{0}))}^{- 1 / 2} exp [- \frac{{(h_{i} [j] - h_{0} [j])}^{2}}{2 S ({\underline{x}}_{0})}]}{Π_{j} exp λ_{j}} .

The normalization terms λ_j are defined such that α(x̱_i|x̱₀) = 1 if the neighbors around x̱_i are identical to those around x̱₀. The measured pixel intensities in log-space at the j^th location in the neighborhoods around x̱₀ and x̱_i are h₀[j] and h_i[j], respectively. We assume that neighboring pixels are independent of one another in order to express the total probability as the product of probabilities for each neighbor. It can be observed that Eq. (7) adapts dynamically the likelihood criterion based on S(x̱₀) according to our spline-based speckle noise model.

The sampling process is repeated until a desired number of samples has been accepted into Ω. These samples are used to estimate the posterior distribution by weighting each sample according to its α(x̱_i|x̱₀), computing the weighted histogram, and then normalizing so that the area under the histogram is unity. The result is our estimate, p̂(logF(x̱)|logG(x̱)), which can be used to reconstruct logF̂(x̱), as described in Eq. (3).

3. Experimental setup

To evaluate the performance of the proposed NSC method for speckle reduction, the proposed novel algorithm was applied to OCT images acquired from different parts of the human eye: 1) healthy human retina; 2) healthy human corneo-scleral limbus and 3) human corneo-scleral lim-bus with pinguecula. The images were acquired with a research-grade, high-speed, ultrahigh-resolution Spectral Domain OCT system (SD-OCT) operating at 1060 nm. The SD-OCT system was powered with a super luminescent diode (λc =1060 nm, Δλ =110 nm, Pout=10 mW) to provide 3 μm and 6 μm axial resolution in the corneal and retinal tissue respectively. At the detection end, the system was interfaced with a high performance spectrometer (P&P Optica) and a fast InGaAs linear array CCD camera (Sensors Unlimited Inc.) with 1024 pixels and readout rate of 47 kHz. Volumetric images were acquired from the human retina and limbus with 1.3 mW power of the incident imaging beam, which resulted in 95 dB system SNR near the zero delay line. The imaging procedure was carried out in the biomedical optical imaging group at the University of Waterloo in accordance with the University of Waterloo ethics regulations for research involving human subjects. The images selected for testing were chosen because they contain a variety of tissue types with different morphology and optical properties that result in spatially varying speckle noise characteristics; they are well-suited for evaluating the speckle compensation performance of the proposed NSC method. For simplification purposes, the NSC method does not account for speckle correlation.

3.1. Experimental OCT imaging data

Three types of UHROCT imagery were used to assess the performance of the NSC method: imagery of healthy human retina, imagery of healthy human corneo-scleral limbus, and imagery of human corneo-scleral limbus with pinguecula. Pinguecula is a common type of conjunctival degeneration in the human eye resulting from prolonged exposure to high levels of ultraviolet light and is typically observed in subjects populating the equatorial to tropical regions. Details about these acquisitions are provided below; raw imaging data is shown in Fig. 1 overlaid with bounding boxes used to calculate some of the metrics described in Section 3.2.

Fig. 1 UHROCT imagery of (A) a healthy human retina, (B) a healthy human corneo-scleral limbus and (C) a human limbus with pinguecula, acquired in-vivo. The blue boxes in these images mark the homogeneous regions of interest (ROIs) used for calculation of the SNR value, while the red boxes are pairs of ROIs that were used to obtain the CNR values.

Download Full Size | PDF

3.1.1. Experiment 1: healthy human retina

The retinal imagery used in this experiment measured 512 A-scans × 512 pixels. The SD-OCT system provided imagery with an axial resolution of 5.7 μm and lateral resolution of 15 μm.

3.1.2. Experiment 2: healthy human corneo-scleral limbus

The imagery of the healthy human corneo-scleral limbus measured 512 A-scans × 512 pixels. The SD-OCT system provided imagery with an axial resolution of 3 μm and lateral resolution of 15 μm.

3.1.3. Experiment 3: human limbus with pinguecula

The imagery of the human corneo-scleral limbus with pinguecula measured 1000 A-scans × 512 pixels. The SD-OCT system provided imagery with an axial resolution of 3 μm and lateral resolution of 15 μm.

3.2. Quantitative performance metrics

A variety of quantitative measures were calculated for each processed imaging dataset in order to quantify the performance of the compared algorithms.

3.2.1. SNR analysis

In order to quantitatively evaluate the performance of the proposed method against the other tested speckle noise compensation methods, the Signal-to-Noise Ratio (SNR) was used. SNR acts as an indicator of the ability to suppress speckle noise. The performance is evaluated by determining if the applied method increases the SNR of the original image and can be compared with other methods based on the largest SNR value. SNR is calculated for the original speckle contaminated OCT images (see Figs. 2 and 4(A)) and the noise-compensated OCT images (see Figs. 2 and 4(B–J)). SNR is calculated as follows in decibels [20]:

S N R = \frac{1}{R} [\sum_{r = 1}^{R} 10 {log}_{10} (\frac{μ_{r}^{2}}{σ_{r}^{2}})],

where μ_r and

σ_{r}^{2}

represent the mean and the variance of the r^th local homogeneous region of interest (ROI), marked by blue boxes in Fig. 1 and R denotes the total number of considered ROIs on the OCT image.

Fig. 2 Results of applying different speckle compensation methods on the human retina imagery.

Download Full Size | PDF

Fig. 3 Results of applying different speckle compensation methods on the healthy human limbus imagery.

Download Full Size | PDF

Fig. 4 Results of applying different speckle compensation methods on the human limbus with pinguecula imagery.

Download Full Size | PDF

3.2.2. CNR analysis

Another quantitative measurement used to assess NSC’s execution is Contrast-to-Noise Ratio (CNR). Similarly to SNR, it is computed for the original speckle contaminated OCT images and the noise-compensated OCT images. The quality metric of CNR acts as an indicator of the ability to improve contrast and preserve structure. A higher CNR value indicates there is a greater separation of image features from a region with background noise. CNR is calculated as follows in decibels [20]:

C N R = \frac{1}{R} [\sum_{r = 1}^{R} 10 {log}_{10} \frac{(μ_{r 1} - μ_{r 2})}{σ_{r 1}^{2} + σ_{r 2}^{2}}],

where μ_r₁ and μ_r₂ represent the mean of two different ROIs that are used to obtain the CNR value and marked by numbered red boxes in Fig. 1 and the σ_r₁ and σ_r₂ are the corresponding variances for these arbitrary ROIs.

3.2.3. ENL analysis

A third quantitative metric is the equivalent number of looks (ENL) which is defined as follows [16]:

E N L = \frac{1}{H} [\sum_{h = 1}^{H} \frac{μ_{h}^{2}}{σ_{h}^{2}}]

where

μ_{h}^{2}

represents the mean and

σ_{h}^{2}

represent the variance of the h^th homogeneous region of interest. It is a measure of smoothness in this homogeneous region, where a greater value depicts a more smooth region.

3.2.4. Edge preservation analysis

The fourth quantitative analytic used is the edge preservation measurement. This analytic is a correlation measure that assesses edge deterioration in an image. The desired measurement is a value of 1, which indicates that the speckle-compensated image has very similar edges to those found in the speckle-contaminated image. The metric is calculated as follows [35]:

η = \frac{Σ (\nabla^{2} V - \bar{\nabla^{2} V}) \cdot (\nabla^{2} \hat{G} - \bar{\nabla^{2} \hat{G}})}{\sqrt{Σ {(\nabla^{2} V - \bar{\nabla^{2} V})}^{2} \cdot Σ {(\nabla^{2} \hat{G} - \bar{\nabla^{2} \hat{G}})}^{2}}}

where ∇²V and ∇²Ĝ are the Laplacian operator applied to the original image and the noise-free reconstruction respectively; and

\bar{\nabla^{2} V}

and

\bar{\nabla^{2} \hat{G}}

are the mean values of a 3 × 3 neighborhood around ∇²V and ∇²Ĝ.

4. Results

The results of the proposed NSC method for the speckle noise compensation were compared with six other speckle noise compensation methods.

Lee filter [12]
Frost filter [13]
Improved Adaptive Complex Diffusion (IACD)* [25]
Kuan filter [26]
Detail Preserving Anisotropic Diffusion (DPAD)* [36]
Speckle Reducing Anisotropic Diffusion (SRAD)* [23]
Bayesian Least Squares Gaussian Scale Mixture (BLS-GSM)* [28]
Stationary Bayesian Speckle Compensation (SBSC)

where * indicates use of author’s provided implementation of the algorithm.

All of these speckle noise compensation methods were applied on the selected OCT images shown in Fig. 1. The impact of these methods on speckle suppression in the images was compared quantitatively using well-known metrics such as the SNR, CNR, ENL and edge preservation [20,35] as well as assessed visually, as explained in the following sections. Regions used to calculate SNR and CNR values are shown overlaid on the speckle-contaminated images shown in Fig. 1. The NSC algorithm was implemented in MATLAB and C++, and run on a computer with an AMD Athlon II X3 445 CPU with 12 GB of RAM. Algorithm run times are presented in Section 4.3.

4.1. Processed images

Images were corrected using a number of speckle noise compensation methods including the proposed method. Results for the healthy human retinal image are shown in Fig. 2; results for the healthy human limbus image are shown in Fig. 3; and results for the human limbus with pinguecula image are shown in Fig. 4.

4.2. Quantitative analysis

The quantitative measures described in Section 3 were calculated on the processed images. The results are presented in the following sections.

4.2.1. SNR analysis

SNR values calculated on the ROIs described in Fig. 1 are shown in Table 1. DPAD resulted in a decrease in SNR compared with the original imagery. NSC, BLSGSM, and SRAD had the greatest improvements in SNR, while the other methods showed more moderate improvements. For the limbus imagery, SRAD caused significant loss of detail.

Table 1. Average SNR values taken over ROIs in the retina (R), healthy human limbus (HL), and human limbus with pinguecula (HLwP) imagery. The highest value for each type of imagery is bolded. IACD and SRAD resulted in some high SNR values but at the cost of overcompensation (see Figs. 2 and 4). BLSGSM and NSC had very close improvements to SNR while other methods had less drastic changes.

View Table | View all tables in this article

4.2.2. CNR analysis

The averaged CNR values obtained for the ROIs described in Fig. 1 are shown in Table 2. In each case, the proposed method achieves a much greater improvement in contrast when compared to the other methods ranging from 6–10 dB greater improvement in CNR. This indicates that the method preserves structure and increases contrast better than the other algorithms.

Table 2. Average CNR values taken over ROIs in the retina (R), healthy human limbus (HL), and human limbus with pinguecula (HLwP) imagery. NSC achieved the best contrast improvement in each type of imagery while maintaining structure.

View Table | View all tables in this article

4.2.3. ENL analysis

The ENL values calculated on the processed imagery are shown in Table 3. IACD has very high ENL values for all three types of imagery, but inspection of Figs. 2 and 4 shows that this is due to overcompensation of speckle which leads to oversmoothing of structures and details. NSC achieved competitive ENL values without oversmoothing.

Table 3. Average ENL values taken over ROIs in the retina (R), healthy human limbus (HL), and human limbus with pinguecula (HLwP) imagery. NSC achieved the best improvements to ENL without smoothing structures and detail (see Figs. 2 and 4).

View Table | View all tables in this article

4.2.4. Edge preservation analysis

The edge preservation values calculated from the processed images are shown in Table 4. BLSGSM and SBSC achieved the highest edge preservation; however, analysis of Figs. 2 and 4 shows that this is due to the preservation of edges present due to speckle. The proposed method achieves the best balance of edge preservation with speckle compensation.

Table 4. Edge preservation values for the retina (R), healthy human limbus (HL), and human limbus with pinguecula (HLwP) imagery. Although BLSGSM achieved higher values, NSC did a better job of balancing speckle suppression with edge preservation (see Figs. 2 and 4).

View Table | View all tables in this article

4.3. Run time

The run-times for each of the compared methods are shown for each of the experiments in Table 5. Note that the compared methods were implemented in MATLAB and not optimized for low run-time. The experiments were run using an AMD Athlon X3 445 processor with 12 GB RAM.

Table 5. Run-times for the various compared methods on each type of imagery.

View Table | View all tables in this article

4.4. Visual assessment

In order to visually assess the impact of the aforementioned speckle compensation methods on improving the OCT images, the speckle-contaminated OCT image and the noise-compensated OCT images are shown in Figs. 2 and 4 for the human retina, the normal corneo-scleral limbus and the limbus with pinguecula. As we observe in part (A) of these figures, the speckle-contaminated OCT images have a grainy pattern due to the existence of speckle noise. This makes it difficult to distinguish the finer morphological details within the images, thus leading to potential difficulties in interpreting the OCT images for diagnostic and scientific study purposes. The results of applying Lee, Frost, Kuan, BLSGSM and DPAD noise compensation methods on the tested OCT images showed that they could reduce the speckle noise to some extent, yet in homogeneous regions the speckle noise was still clearly observable. The SBSC method also exhibited speckle in homogeneous regions. In the case of the Lee filter, it does not perform as well as NSC despite its assumption that the image variance changes locally. This is because it still assumes that it has a set noise variance whereas NSC takes into account the non-stationarity of speckle. The results from applying the IACD and SRAD method demonstrated that it could provide improved noise suppression in some cases when compared to the other methods; however, it overcompensates for the speckle and removes important structural details. We note that the proposed NSC method provided noticeably improved speckle noise suppression performance when compared to the other methods while better preserving the fine morphological and original details.

Discussion

In this study, we introduced a spatially-adaptive Monte Carlo speckle noise compensation method based on a spline-based speckle noise modeling approach for the purpose of reducing non-stationary speckle noise in OCT images. This method spatially-adapts to the spatially-varying noise characteristics of the speckle noise found in OCT imagery. Through quantitative analysis, it was shown that the proposed method can provide improved speckle noise suppression performance in terms of SNR, CNR, ENL and edge preservation. For CNR, the proposed NSC algorithm improved the speckle-contaminated images by over 10 dB. This improvement enhanced the visibility of fine morphological details in the UHROCT tomograms, such as clear delineation of the individual retinal layers, as well as the blood vessel walls of the choroidal vasculature in human retina; the layered structure of the human cornea at the limbus, as well as the vasculature pattern of the pinguecula, that are important for the clinical diagnostics of a variety of pathologies. The novel method also improved the SNR by at least 10 dB over the speckle-contaminated OCT images. In some cases, the SRAD method had better SNR results than the proposed NSC method; however in these cases the SRAD method overcompensated for the presence of speckle noise and instead removed important fine details. ENL analysis showed that SRAD and IACD improved the smoothness of homogeneous regions of interest, however visual assessment proved that structures and fine details were not preserved at the cost of smoothing these regions. NSC had the next highest ENL metrics that showed a successful compromise of smoothing homogeneous regions and preservation of structure. Edge preservation analysis supported these results by proving that NSC had high edge preservation but not to the point of speckle preservation as found by the higher edge preservation metrics shown by BLSGSM and SBSC methods.

Visual assessment showed that the Lee, Frost, Kuan, DPAD, BLSGSM, and SBSC methods had limited speckle noise suppression performance compared to the SRAD, IACD and NSC methods. The SRAD and IACD approaches were more successful than these other methods, but removed image details. The speckle compensation of SBSC was limited in regions with greater noise. The non-stationary noise model used by NSC clearly offer better speckle compensation than similar methods using a stationary noise model such as SBSC. The proposed NSC method exhibited strong speckle noise suppression performance while preserving important structures in the image.

Although NSC was applied to two-dimensional B-scans in these experiments, it can be easily applied to three-dimensional OCT data by repeated application to the constituent B-scans. However, applying it in this method would not fully exploit the additional local structure; as such, we propose investigation into the speckle estimation performance of a higher-dimensional spline fit to incorporate speckle characteristics from neighboring B-scans. Higher-dimensional spline fits can further be used to incorporate additional imaging modalities such as functional imaging techniques. Furthermore, we propose additional research into the performance of different regularization methods, for example the use of a piecewise-constant model instead of a spline fit.

As a result of the positive experiments in this study, we propose that further investigation of the speckle noise compensation method for use in OCT reconstruction be explored. Notably, we propose investigation of other speckle models, such as a piecewise-constant model, to see if improved speckle noise compensation performance can be obtained.

Acknowledgments

The study was funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Canadian Institutes of Health Research (CIHR), and the Ontario Ministry of Research and Innovation. AC, DL, AB, JG and AW were involved in designing the study and performing statistical analysis. KB provided the UHROCT images. AC, DL, AB, JG, AW, and KB were involved in the writing and editing. The authors have declared that no competing interests exist.

References and links

1. B. Potsaid, I. Gorczynska, V. J. Srinivasan, Y. Chen, J. Liu, J. Jiang, A. Cable, J. S. Duker, and J. G. Fujimoto, “Ultrahigh speed spectral / Fourier domain opthalmic OCT imaging,” Proc. SPIE 7163, 716307 (2009). [CrossRef]

2. T. Klein, W. Wieser, C. M. Eigenwillig, B. R. Biedermann, and R. Huber, “Megahertz OCT for ultrawide-field retinal imaging with a 1050 nm Fourier domain mode-locked laser,” Opt. Express 19, 3044–3062 (2011). [CrossRef] [PubMed]

3. W. Drexler and J. Fujimoto, “Biological and medical physics, biomedical engineering,” in Optical Coherence Tomography: Technology and Applications (Springer, 2008). [CrossRef]

4. J. Rogowska and M. E. Brezinski, “Evaluation of the adaptive speckle suppression filter for coronary optical coherence tomography imaging,” 19, 1261–1266 (2000).

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

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

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

8. A. E. Desjardins, B. J. Vakoc, W. Y. Oh, S. M. 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]

9. T. M. Jørgensen, L. Thrane, M. Mogensen, F. Pedersen, and P. E. Andersen, “Speckle reduction in optical coherence tomography images of human skin by a spatial diversity method,” in Optical Coherence Tomography and Coherence Techniques III, vol. 6627 of Proceedings of SPIE-OSA Biomedical Optics, P. Andersen and Z. Chen, eds. (Optical Society of America, 2007), pp. 22.

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

11. L. Fang, S. Li, Q. Nie, J. A. Izatt, C. A. Toth, and S. Farsiu, “Sparsity based denoising of spectral domain optical coherence tomography images,” Biomed. Opt. Express 3, 927–942 (2012). [CrossRef] [PubMed]

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

13. 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 Analysis Mach. Intell. 4, 157–166 (1982). [CrossRef]

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

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

16. A. Wong, A. Mishra, K. Bizheva, and D. A. Clausi, “General Bayesian estimation for speckle noise reduction in optical coherence tomography retinal imagery,” Opt. Express 18, 8338–8352 (2010). [CrossRef] [PubMed]

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

18. 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, 2878–2880 (2004). [CrossRef]

19. A. Ozcan, A. Bilenca, A. E. Desjardins, B. E. Bouma, and G. J. Tearney, “Speckle reduction in optical coherence tomography images using digital filtering,” Opt. Lett. 24, 1901–1910 (2007).

20. P. Puvanathasan and K. Bizheva, “Speckle noise reduction algorithm for optical coherence tomography based on interval type II fuzzy set,” Opt. Express 15, 15747–15758 (2007). [CrossRef] [PubMed]

21. M. Gargesha, M. W. Jenkins, A. M. Rollins, and D. L. Wilson, “Denoising and 4D visualization of OCT images,” Opt. Express 16, 12313–12333 (2008). [CrossRef] [PubMed]

22. Z. Jian, L. Yu, B. Rao, B. J. Tromberg, and Z. Chen, “Three-dimensional speckle suppression in optical coherence tomography based on the curvelet transform,” Opt. Express 18, 1024–1032 (2010). [CrossRef] [PubMed]

23. Y. Yu and S. Acton, “Speckle reducing anisotropic diffusion,” Opt. Express 11, 1260–1270 (2002).

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

25. R. Bernardes, C. Maduro, P. Serranho, A. Araujo, S. Barberio, and J. Cunha-Vas, “Improved adaptive complex diffusion despeckling filter,” Opt. Express 18, 24048–24059 (2010). [CrossRef] [PubMed]

26. D. T. Kuan, A. A. Sawchuk, T. C. Strand, and P. Chavel, “Adaptive noise smoothing filter for images with signal-dependent noise,” IEEE Trans. Pattern Analysis Mach. Intell. 7, 165–177 (1985). [CrossRef]

27. R. Touzi, “A review of speckle filtering in the context of estimation theory,” IEEE Trans Geosci. Remote Sens. 40, 2392–2404 (2002). [CrossRef]

28. J. Portilla, V. Strela, M. J. Wainwright, and E. P. Simoncelli, “Image denoising using scale mixtures of Gaussians in the wavelet domain,” IEEE Trans. Image Process. 23, 1338–1351 (2003). [CrossRef]

29. J.-A. Guerrero-Colón, L. Mancera, and J. Portilla, “Image restoration using space-variant Gaussian scale mixtures in overcomplete pyramids,” IEEE Trans. Image Process. 17, 27–41 (2008). [CrossRef] [PubMed]

30. A. Buades, B. Coll, and J. M. Morel, “A non-local algorithm for image denoising,” in Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, Vol. 2 (IEEE, 2005), pp. 60–65.

31. A. Leigh, A. Wong, D. A. Clausi, and P. Fieguth, “Comprehensive analysis on the effects of noise estimation strategies on image noise artifact suppression performance,” in Proceedings of IEEE International Symposium on Multimedia(IEEE, 2011), pp. 97–104.

32. J. A. Fessler, “Tomographic reconstruction using information-weighted spline smoothing,” in Information Processing in Medical Imaging,H. H. Barrett and A. F. Gmitro, eds. (Springer BerlinHeidelberg, 1993), pp. 372–386. [CrossRef]

33. P. Fieguth, Statistical Image Processing and Multidimensional Modeling (Springer Science+Business Media, New York, 2011), chap. 3, p. 65.

34. M.-H. Chen, “Importance-weighted marginal Bayesian posterior density estimation,” J. Am. Stat. Assoc. 89, 818–824 (1994). [CrossRef]

35. F. Sattar, L. Floreby, G. Salomonsson, and B. Lovstrom, “Image enhancement based on a nonlinear multiscale method,” IEEE Trans. Image Process. 6, 888–895 (1997). [CrossRef] [PubMed]

36. S. Aja-Fernández and C. Alberola-López, “On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering,” IEEE Trans. Image Process. 15, 2694–2701 (2006). [CrossRef] [PubMed]

Method	SNR (dB)
Method	R	HL	HLwP
Speckle Contaminated	12.0	15.2	13.5
Frost	17.4	22.1	21.2
Lee	14.9	19.4	18.4
IACD	12.4	28.3	25.8
Kuan	14.9	19.4	18.1
DPAD	9.4	14.0	12.5
SRAD	27.6	34.8	16.9
BLSGSM	26.6	24.9	25.6
SBSC	16.6	21.4	24.6
NSC	22.7	27.3	26.0

Method	CNR (db)
Method	R	HL	HLwP
Speckle Contaminated	0.84	2.35	0.91
Frost	2.70	5.02	2.47
Lee	2.13	3.75	1.78
IACD	3.94	9.14	5.08
Kuan	2.57	6.98	1.471
DPAD	1.63	5.54	1.464
SRAD	1.58	11.53	1.49
BLSGSM	2.83	5.64	5.84
SBSC	4.09	12.68	4.29
NSC	11.53	21.81	12.15

Method	ENL
Method	R	HL	HLwP
Speckle Contaminated	5.9	38.0	8.3
Frost	220.5	278.7	62.8
Lee	75.5	123.5	29.0
IACD	813.4	1715.0	399.7
Kuan	75.5	123.5	10.2
DPAD	25.2	68.6	9.4
SRAD	39.2	7470.0	20.9
BLSGSM	227.6	81.3	109.3
SBSC	71.7	347.8	101.3
NSC	437.4	751.6	284.9

Method	Edge preservation
Method	R	HL	HLwP
Frost	0.431	0.431	0.244
Lee	0.405	0.402	0.209
IACD	0.420	0.411	0.230
Kuan	0.038	0.009	0.018
DPAD	0.239	0.223	0.081
SRAD	0.416	0.249	0.118
BLSGSM	0.819	0.483	0.472
SBSC	0.562	0.518	0.245
NSC	0.435	0.437	0.254

Method	Time (s)
	R	HL	HLwP
Frost	13.1	14.9	33.9
Lee	0.1	0.1	0.2
IACD	1.9	2.4	4.9
Kuan	15.3	17.7	42.0
DPAD	7.0	8.1	19.8
SRAD	8.5	19.8	45.1
BLSGSM	19.0	30.4	48.1
SBSC	48.4	56.9	123.7
NSC	54.9	64.7	140.4

Stochastic speckle noise compensation in optical coherence tomography using non-stationary spline-based speckle noise modelling

Abstract

1. Introduction

2. Methodology

2.1. Problem formulation

2.2. Non-stationary spline-based speckle noise model

2.3. Spatially-adaptive Monte Carlo posterior estimation

3. Experimental setup

3.1. Experimental OCT imaging data

3.1.1. Experiment 1: healthy human retina

3.1.2. Experiment 2: healthy human corneo-scleral limbus

3.1.3. Experiment 3: human limbus with pinguecula

3.2. Quantitative performance metrics

3.2.1. SNR analysis

3.2.2. CNR analysis

3.2.3. ENL analysis

3.2.4. Edge preservation analysis

4. Results

4.1. Processed images

4.2. Quantitative analysis

4.2.1. SNR analysis

4.2.2. CNR analysis

4.2.3. ENL analysis

4.2.4. Edge preservation analysis

4.3. Run time

4.4. Visual assessment

Discussion

Acknowledgments

References and links

Cited By

Figures (4)

Tables (5)

Equations (11)

Biomedical Optics Express