Optical Coherence Tomography Noise Reduction Using Anisotropic Local Bivariate Gaussian Mixture Prior in 3D Complex Wavelet Domain.

In this paper, MMSE estimator is employed for noise-free 3D OCT data recovery in 3D complex wavelet domain. Since the proposed distribution for noise-free data plays a key role in the performance of MMSE estimator, a priori distribution for the pdf of noise-free 3D complex wavelet coefficients is proposed which is able to model the main statistical properties of wavelets. We model the coefficients with a mixture of two bivariate Gaussian pdfs with local parameters which are able to capture the heavy-tailed property and inter- and intrascale dependencies of coefficients. In addition, based on the special structure of OCT images, we use an anisotropic windowing procedure for local parameters estimation that results in visual quality improvement. On this base, several OCT despeckling algorithms are obtained based on using Gaussian/two-sided Rayleigh noise distribution and homomorphic/nonhomomorphic model. In order to evaluate the performance of the proposed algorithm, we use 156 selected ROIs from 650 × 512 × 128 OCT dataset in the presence of wet AMD pathology. Our simulations show that the best MMSE estimator using local bivariate mixture prior is for the nonhomomorphic model in the presence of Gaussian noise which results in an improvement of 7.8 ± 1.7 in CNR.

1. Introduction

Optical coherence tomography (OCT) is an optical signal acquisition and processing method that captures 3D images from within optical scattering media such as biological tissues [1-4]. For example, in ophthalmology, OCT is used to obtain detailed images from within the retina [4]. Similar to other optical tomographic techniques, OCT suffers from speckle noise that reduces the ability of image interpretation [5]. So, noise reduction is an essential part of OCT image processing systems. Until now, several techniques for OCT noise reduction have been reported [6-14]. Initial methods perform in complex domain [15], that is, before producing magnitude of OCT interference signal, while most introduced despeckling methods are applied after an OCT image is formed [6-14]. These methods, which usually suppose multiplicative noise for speckled data, also can be categorized into image domain and transform domain methods. As an example for image domain techniques, in [16] the rotating kernel transform (RKT) filters are applied on an image with a set of oriented kernels and keep the largest filter output for each pixel. Other image domain methods based on enhanced Lee filter [17], median filter [17], symmetric nearest neighbor filter [17], and adaptive Wiener filter [17], and I-divergence regularization [6] and PDE-based nonlinear diffusion methods [14] have been reported in the literature.

Transform domain techniques typically outperform the image domain techniques because incorporating speckle statistics in the despeckling process would be facilitated in sparse domains. Such techniques apply a sparse transform (such as wavelet and curvelet transforms) [7–12, 18] directly on data (viz., nonhomomorphic methods) or on log-transformed data (viz., homomorphic methods), and suppose that in the sparse domain noise is converted to additive white Gaussian noise (AWGN) [13] or other models which can be removed using an appropriate shrinkage function. For example, in [18], a spatially adaptive wavelet thresholding method is used for speckle suppression in log-transformed domain. Since actual signal in OCT images consists of horizontal edges arising from reflections at the layer boundaries, most of the edge information is in “low-pass”-“high-pass” (LH) subbands (and some of it in HH subbands). Therefore, an increased threshold in the vertical subbands using a constant multiplier (K = 4) is chosen to decrease further noise with a minimal effect on the edge sharpness. Other transform domain methods based on hard thresholding in 3D curvelet domain [8], soft thresholding in discrete complex wavelet transform (DCWT) domain [9], and temporal and spatial wavelet-based filtering [10] have been reported in other literatures.

In fact denoising is the problem of obtaining the noise-free data from noisy data observation, which may be solved in a deterministic or probabilistic framework. In the first case, each voxel is considered as an unknown deterministic variable, and non-Bayesian techniques are employed to solve this problem. In the second case, the data is modeled as a random field, and Bayesian methods are used for the estimation of clean data from the noisy environment. Therefore, the proposed prior probability distributions for noise-free data and noise (i.e., proposed as speckle for OCT data) play a key role in the noise reduction problem.

1.1. Statistical Properties of Noise-Free Coefficients

Description of the statistical properties of natural signals can be facilitated in the wavelet domain [19] due to sparseness and decorrelation properties of wavelets [20]. The sparseness property states that the marginal pdf of wavelet coefficients in each subband has a large peak at zero and its tails fall to zero slower than the Gaussian pdf (leptokurtic). On this base, some long-tailed pdfs such as generalized Gaussian distribution (GGD) [21, 22], α-stable distributions [23], Bessel K form densities [24, 25], and mixture pdfs [26-31] have been proposed. Although the decorrelation property of wavelets states that coefficients at the same positions in the adjacent scales are uncorrelated, it does not mean that they are independent. The interscale dependency of wavelet states that large/small values of wavelet coefficients tend to propagate across scales [32]. Some researchers have proposed hidden Markov models (HMMs) [33] and Markov random fields (MRFs) [34] to model the interscale dependency [35]. Recently, it has been shown that some non-Gaussian bivariate joint pdfs for each coefficient and its parent, such as circular symmetric Laplacian pdf [36], bivariate Cauchy distribution [37], (multivariate) Gaussian scale mixture (GSM) model [27, 38, 39], and bivariate Laplacian mixture models [40] are able to capture this property easily and produce better denoising results with lower computational complexity.

The dependencies between wavelet coefficients are not restricted to the interscale dependency. There is another dependency between spatial adjacent coefficients in each subband, namely, intrascale dependency [42]. This dependency states that if a particular wavelet coefficient is large/small, then the spatial adjacent coefficients are likely to be large/small too. Usually this property is captured using local parameters for pdfs [37], and it has been shown that denoising algorithms using this property for statistical modeling of wavelets are able to improve the denoising results [43-45]. For example, Mihçak [43] employs local variance for Gaussian pdf to model intrascale dependency. In [44], a mixture of two Laplace pdf with local parameters is proposed for simultaneously capturing heavy-tailed nature and intrascale dependency. Reference [45], using local variance for proposed model in [36], improves the results for noise reduction application because this local pdf models both interscale and intrascale dependencies. In this paper, we extend the proposed pdf in [46] based on a mixture of bivariate Gaussian pdfs with local parameters for noise-free wavelet coefficients. Since the empirically observed distribution of wavelet coefficient pairs in adjacent scales have elliptical symmetry, we use different variances for marginal pdfs that lead to an elliptical symmetric bivariate pdf instead of circular symmetric pdf. Recently, it has been shown that using anisotropic window instead of square window can improve the denoising results [47]. Based on the special structure of OCT data, we choose an anisotropic windowing procedure for local parameters estimation that results in visual quality improvement.

1.2. Discrete Complex Wavelet Transform (DCWT)

The wavelet based image denoising consists of the following steps.

Signal transformation of the noisy observation.

Modification of the noisy wavelet coefficients based on some criteria.

Inverse signal transformation of modified coefficients.

As explained earlier, the second step depends on the type of estimator and for a minimum mean square error (MMSE) estimator, the proposed model for signal and noise (which we propose as a multiplicative model), the proposed pdf of noise-free wavelet coefficients (modeled, in this paper, as a mixture of bivariate Gaussian pdfs with local parameters), and the proposed pdf for noise (with which we test both Gaussian and two-sided Rayleigh distributions) define the performance of the algorithm. However, for the first and last steps of wavelet-based denoising algorithm, the type of transformation plays a key role. In this paper, we use DCWT [48] instead of ordinary discrete wavelet transform (DWT). Despite DWT being a sparse representation that outperforms many signal processing approaches, it does not lead to an optimum performance in all applications and suffers from several fundamental shortcomings (especially in high-dimensional cases), which DCWT avoids them. These shortcomings are as follows.

In the neighborhood of an edge, the DWT produces both large and small wavelet coefficients. In contrast, the magnitudes of DCWT coefficients are more directly related to their adjacency to the edge. The main reason of this phenomenon is using bandpass filters that produce DWT coefficients which oscillate positively and negatively around the singularities, and this subject complicates wavelet-based processing.

DWT is not shift invariant. It means that a small shift in the input signal of DWT makes the total energy of wavelet coefficients in subband completely differ. This shift greatly perturbs oscillation pattern around singularities of the DWT coefficient which complicates wavelet-domain processing.

Since the DWT coefficients in each subband are produced via critical sampling after using nonideal low-pass and high-pass filters, substantial aliasing would be produced. If the wavelet coefficients are not changed, the inverse DWT cancels this aliasing. Applying any processing method on wavelet coefficients (such as thresholding) disarranges this balance between the forward and inverse transforms which leads to artifacts in the reconstructed signal.

The directional selectivity of 2D DCWT has been explained in Appendix A. Similar to the 2D case, the standard 3D data transforms, which are separable multiplication of 1D tensors, do not provide useful representations with good energy compaction property for 3D data. For example, the multi-dimensional standard separable DWT mixes orientations and motions in its subbands and produces the checkerboard artifacts (Figure 1). In contrast, since the spectrum of the (approximately) analytic 1D wavelet is supported on only one side of the frequency axis, the spectrum of the DCWT in 3D domain is supported in only 1/27 of the 3D frequency plane. So, instead of 3D DWT, usually oriented transforms such as 3D DCWT are proposed for 3D data processing [41, 48–52]. Figure 1 shows a comparison between subbands of 3D DWT and 3D DCWT.

Figure 1

A comparison between the idealized support of the Fourier spectrum of each standard and complex wavelet in the 3D frequency domain. (a) Isosurfaces of the 7 3D wavelets for a standard 3D wavelet transform. The blue and red colors have the same amplitude, but their phases are complement. (b) Isosurfaces of 7 of the 28 3D wavelets for a 3D DCWT. Each subband corresponds to motion in a specific direction [41].

1.3. Organization of the Paper

In Section 2, we explain our proposed pdf for noise-free 3D DCWT coefficients, that is, a mixture of bivariate Gaussian distributions with local parameters. In Section 3, at first we obtain a local thresholding function supposing a priori distribution as a bivariate Gaussian pdf with local variance, and then in a Bayesian framework we produce our new shrinkage functions derived from the proposed pdf and using Gaussian/two-sided Rayleigh noise distribution and homomorphic/non-homomorphic model. In Section 4, we explain the proposed anisotropic window selection procedure for local parameter estimation based on special structure of OCT data. In Section 5, we use our model for wavelet-based denoising of several 3D OCT data. We compare our methods visually and in terms of PSNR. Also in this section, we use the proposed method for nonspeckle noise reduction. Finally, in Section 6, we summarize this paper and suggest some future work.

2. Bivariate Gaussian Mixture Model with Local Parameters

One of the primary properties of the wavelet transform is compression. This property means that the marginal distributions of wavelet coefficients are highly kurtotic, and so long-tailed distributions are suitable models for marginal pdf. A zero-mean mixture model could have a large peak at zero and would be long tailed. For example, in [22, 26, 29, 31] a mixture of Gaussian distributions is proposed to model the heavy-tailed nature of wavelet coefficients. Figure 2 shows this model that consists of two zero-mean Gaussian distributions with two different variances. The Gaussian pdf with low variance can model the large peak at zero and the Gaussian pdf with high variance can model tails of distribution. The secondary properties of the wavelet transform are clustering and persistence. The clustering property, that is called intrascale dependency, states that if a particular wavelet coefficient is large/small, then adjacent coefficients are very likely to also be large/small [36], and usually local pdfs are able to model this property. The persistence property, that is called the interscale dependency, states that large/small values of wavelet coefficients tend to propagate across scales [36]. As an example, Figure 3 illustrates the empirical joint parent-child histogram of wavelet coefficients computed from the 200, 512 × 512 images from the Corel image database [42]. Usually this property can be modeled using proper bivariate pdfs.

Figure 2

Zero-mean Gaussian mixture model (left image) and empirical histogram of wavelet in a subband together with the Gaussian mixture model (right image) [26].

Figure 3

Empirical joint parent-child histogram of wavelet coefficients (computed from the Corel image database) [42].

2.1. Description of the Proposed Model

In this paper, we assume a pdf as a mixture of two bivariate Gaussian pdfs with local parameters in order to model the distribution of wavelet coefficients of images as follows: where a(k) ∈ [0,1],  σ 11(k),  σ 12(k),  σ 21(k),  σ 22(k) are the mixture model parameters. For each random bivariable, the second component is the parent of the first component; for example, w 2(k) represent, the parent of w 1(k) at the same spatial position as the kth wavelet coefficient w 1(k) and at the next coarser scale.

Our proposed model in this paper, that is a mixture of bivariate Gaussian pdfs with local parameters, is mixture, bivariate and local. Therefore, it is able to simultaneously capture the heavy-tailed property and inter- and intrascale dependencies.

After substitution of mixture model in the definition of E(w 1(k)w 2(k)), we can see that this pdf represents two uncorrelated random variables as follows:

Interestingly, the marginal pdf of w (k) for i = 1, 2 is the mixture of two univariate Gaussian pdf with local parameters [44],

It is easy to see that w 1(k), w 2(k) are not independent; that is,

See Appendix B for more explanation.

2.2. Local EM Algorithm

To characterize the parameters in (1), it is necessary to have the parameters σ 11(k),  σ 21(k), σ 12(k), σ 22(k), and a(k). For this mixture model, we use an iterative numerical algorithm to estimate these parameters. The expectation maximization (EM) algorithm is most frequently used to estimate such parameters. Usually, the EM algorithm for mixture models employs all data in each subband to obtain the parameters. Using this global EM algorithm, equal parameters are obtained for all data in each subband. However, to model the intrascale dependency, we must incorporate the local statistics and need to have different parameters for each voxel in each subband. So, we introduce a local version of EM algorithm. This local EM algorithm is able to obtain separate parameters for each voxel by the implementation of EM algorithm in each window N(k) centered at w(k). This iterative algorithm has two steps. Assuming that the observed data for k = 1,…, N, the E-step calculates the responsibility factors for each data as follows:

The M-step updates the parameters a(k),  σ 11(k),  σ 12(k),  σ 21(k),  and  σ 22(k). a(k) is computed by where M is the number of coefficients in the square window N(k) centered at w(k).

The variances σ 11(k),  σ 12(k),  σ 21(k),  and  σ 22(k) are computed by [40]:

3. Denoising Using MMSE Estimator

In this section, the denoising of a 3D OCT data is considered. We assume that dominant noise in OCT data is speckle. In this case as a common model, we propose multiplicative model as follows: where i is the index of voxel and is between 1 and number of voxels.

As explained in Introduction, reported transform-based OCT noise reduction methods in the literatures [7–12, 18] usually at first transform data into log domain, and suppose that noise in log domain is AWGN: where in this paper W shows 3D DCWT operator. So, we can write where w(k), y(k), and n(k) are, respectively, the kth noise-free 3D DCWT coefficients, noisy 3D DCWT coefficients, and noise in the 3D DCWT domain.

Recently, it has been reported [56-59] that non-homomorphic techniques that do not use this nonlinear operation and apply wavelet transform directly on speckled data lead to unbiased estimation of the data and decrease the computational complexity. On this base after applying 3D DCWT (directly) on data, we would have

Again we can write where w(k),  y(k), and n(k) are, respectively, the kth noise-free 3D DCWT coefficients, noisy 3D DCWT coefficients, and noise in the 3D DCWT domain. Since speckle noise g can be modeled as a unit-mean random process independent of the noise-free data, we would have E[W(s(g − 1))] = 0, and also it can be easily shown [58] that E[W(s)W(s(g − 1))] = 0 which means that w(k) and n(k) are zero-mean uncorrelated random variables. If w (k), y (k), and n (k) show the parent coefficients of w(k), y(k), and n(k), respectively, we can write

Based on the persistence property, we need to have a bivariate model based on parent-child pairs. So, we can propose the following bivariate model: where , , and are, respectively, the kth parent-child pairs of noise-free 3D DCWT coefficients, noisy 3D DCWT coefficients, and additive noise in the 3D DCWT domain. In the literature, several models such as K-distribution, Rayleigh, Weibull, log-normal, and Nakagami distributions have been proposed [57, 58, 60–63] for speckle in image domain. In this paper, we test both AWGN and two-sided Rayleigh model for noise in wavelet domain as follows: where σ 2 = 2α 2 shows the noise variance.

Now our goal is the estimation of from , where is a Gaussian or two-sided Rayleigh according to some criteria.

If we employ the MMSE estimator for the estimation problem, we get the posterior mean as an optimal solution:

3.1. Denoising Based on Modeling Noise-Free Data by Bivariate Gaussian PDF with Local Variance

In order to solve (17), we must know the prior distribution of 3D DCWT coefficients, that is, . Defining , if we suppose that w(k), w (k) are independent Gaussian pdf with variances σ 1(k) and σ 2(k), the following bivariate Gaussian pdf with local variances can be proposed for the noise-free wavelet coefficients:

In this case, w(k) and w (k) are uncorrelated and independent, and therefore the MMSE estimator of w(k), w (k) yields the shrinkage function corresponding to univariate Gaussian pdf, that is, Wiener filter [21] as follows:

And so we can write

Similarly, if we choose two-sided Rayleigh pdf for noise distribution, the following estimator is obtained [44]: where

And so we can write

Suppose that the input noise variance is known. To implement (20) or (23), we must know the parameter of the prior σ(k) (suppose that σ(k) = σ (k)). Mihçak et al. [43] showed that using local variance (instead of global variance) for Wiener filter leads to a substantial improvement in denoising results (using local variance allows incorporating the local statistics of image into the proposed prior). It has been shown in the literature that the correctness of estimation of variance is an impact factor for denoising [23, 27, 34, 42–46]. Thus, the proposed criteria for estimation of the variance, such as the involved data for estimation (e.g., in some approaches the coarser scales are used as a source of prior), the type of estimator, and the shape and size of the proposed window for the local estimation of the variance, play key roles in the performance of denoising procedure. For example, in [54] a recurrence equation using a local Gaussian pdf is used for estimation of σ(k) or in [64] the variable size of the locally adaptive window is obtained using a region-based approach. However, for each data point , a simple estimation of σ(k) can be formed based on a local neighborhood N(k). In simplest case, we can use a square window N(k) centered at and suppose that in this window the variance is approximately constant. Then, an empirical estimate for σ(k) can be obtained as follows: where M is the number of coefficients in N(k) and σ can be estimated by [4] σ = median{|noisy wavelet coefficients in finest scale|}/0.6745. In this estimation, we propose the coarser scale as a source of prior, but another estimate can be obtained using only spatial adjacent in the same scale. It has been shown [47] that the local features in the edges of images are not isotropic and so can be better modeled in a shape-adaptive window selection manner. We explain in this regard in Section 4 and try to improve the denoising results by using anisotropic window instead of square window for the estimation of local parameters (such as variance in (24)).

3.2. Denoising Based on Modeling Noise-Free Data by a Mixture of Bivariate Gaussian PDFs with Local Parameters

A nonlinear shrinkage function for wavelet-based denoising, which is derived by assuming that the noise-free wavelet coefficients follow a bivariate Gaussian mixture model with local parameters given by (1), is introduced in this section. Substituting (1) in (17), we can write where

In fact, is the 2D convolution of the pdf of (defined in (15) or (16)) and p (defined in (1)). Using (15) as proposed model for noise, both and p are bivariate Gaussian pdfs. So, we obtain

For two-sided Rayleigh noise (16), more computations are needed. After some simplifications, the final formula would be where

Using (19), we can obtain numerators of (25), and finally (25) for AWGN can be written as where

We call the new obtained bivariate local shrinkage function as BiGaussMixShrinkL. Figure 4 shows this shrinkage function with sample constant parameters.

Figure 4

A shrinkage function produced from BiGaussMixShrink for sample parameters.

Similarly, using (21), we can obtain numerators of (25), and finally (25) for two-sided Rayleigh noise is obtained as where

We call this bivariate local shrinkage function as BiGaussRayMixShrinkL. Figure 5 shows this shrinkage function with sample constant parameters.

Figure 5

A shrinkage function produced from BiGaussRayMixShrink for sample parameters.

For implementation of our denoising algorithm, we must estimate the parameters σ (k) for i, j = 1,2, and a(k) (that are for noise-free data) from noisy observation. For AWGN, the noisy observation would be a Gaussian mixture model with parameters a(k), for i, j = 1,2. So, the following local EM algorithm is used to obtain the parameters.

E-step

M-step where M is the number of coefficients in the window N(K) centered at . As discussed in the literatures [40], for non-Gaussian mixture models, which is a case for two-sided Rayleigh noise, using (34)–(36) finally converge to the final results.

Our denoising algorithm is summarized in Algorithm 1.

Algorithm 1

Outline of the proposed denoising algorithm.

4. Shape Adaptive Windows Selection

It has been shown that using anisotropic and shape adaptive window for local parameter estimation can extremely improve the modeling and processing results. For example, in [47] a new image denoising is introduced that proposes an anisotropic window around each pixel of image and obtains the denoised pixel just by using the located data in the window. Comparing with the denoising methods that are based on proposing isotropic window around each pixel (e.g., [23, 27, 34, 42–46]), the proposed method in [47] is able to segment the image to rather smoothed regions before denoising due to anisotropic window selection that leads to improvement of denoising results. As explained before, the mixture model parameters in each subbands are estimated locally using an isotropic window around each voxel. In this section at first we explain the structure of macular OCT then we introduce 3D “linear polynomial approximation-intersection confidence interval” (LPA-ICI) method for applying shape adaptive window selection around each voxel in 3D DCWT domain. So we will try the despeckling results in 3D DCWT domain by choosing an anisotropic window (instead of isotropic) for estimating the parameters of mixture model in each subband locally.

4.1. OCT Structure

To select the shape-adaptive window, we must take a look at the special structure of OCT data. In ophthalmology, the OCT data shows detailed images from within the retina. The automated analysis of OCT images can be used for the image-guided retinal therapy. Every year, many people become blind as a result of age-related macular degeneration (AMD) due to affecting the central retina where our central vision is perceived. The most sight-threatening form of AMD is called exudative or wet AMD. Choroidal neovascularization (CNV) is a common symptom of the degenerative maculopathy wet AMD. A wealth of powerful new treatments for CNV, especially anti-VEGF agents, have become available very recently to restore normal visual function. The risk of ocular adverse events, including the devastating intraocular infection, endophthalmitis, increases with repeated intravitreal treatment injections, and the effects of chronic treatment with anti-VEGF agents on the retina are unknown. Ideally a more cost-effective, patient-specific dosing strategy with the minimally necessary number of anti-VEGF injections is required. With all the promise, these novel treatments will only reach their full potential when objective and early indices of treatment response are developed. Prior to the introduction of retinal OCT imaging, clinical assessment of whether the preservation or restoration of visual function is successful, which indeed is the ultimate goal of treatment, could only be obtained by measuring visual function. Unfortunately, visual function lags structural response and is cumbersome and noisy, and its reproducibility is limited. Two-dimensional OCT imaging of the retina was introduced several years ago, and was rapidly adopted, among others, to qualitatively measure macular structure as an indicator of AMD treatment response and for guidance of retreatment in CNV recurrence. It is now becoming clear that these simplified structural measures though leading indicators of visual function are inadequate, as they are based on simplified interpretation of single transverse slices of the macula, some patients do not recover visual function even though their total macular thickness has become normally thin after treatment, and others paradoxically gain visual acuity while their macula is still thickened.

True 3D spectral OCT imaging, that became available in 2007 is fast (1.5 s per volume scan), allows full 3D retinal coverage at a much higher resolution and offers improved imaging of subtle differences in retinal structure. In the recent years [65, 66], 3D analysis of 3D OCT as an improved measure of subtle macular structure has been proposed motivated by various hypotheses as follows: A model of retinal response to initial anti-VEGF treatment for CNV, based on quantitative 3D OCT-derived measures, can predict the timing of retreatment.

On this base, developing analysis methods and approaches for 3D spectral OCT image analysis in the presence of wet AMD pathology (Symptomatic Exudate Associated Derangements or SEADs, also known as AMD-related cysts, vessel leakages, etc.) and assessing their performance by comparison to expert analyses are of utmost interest. Another interesting subject is determining how well the quantitative SEAD- and layer-derived measures from 3D OCT predict the patient-specific outcome parameters in response to postinduction anti-VEGF treatment in patients with CNV in order to predict the timing of retreatment.

Figure 6 shows several sample macular OCTs and detected SEADs by an expert as the region of interest (ROI). As we can see in this figure, the most important information of OCT data (about retina layers) is located in the center of OCT images.

Figure 6

Macular OCTs and detected SEADs by an expert.

4.2. 3D LPA-ICI for Data between the First and Last Layers

In [47], a new image denoising based on using an anisotropic window around each pixel of image is introduced. To select the anisotropic window, the linear directional filters g that are obtained using local polynomial approximation (LPA) are employed. The θ indicates the direction of filter that is a member of countable set {θ 1, θ 2,…, θ }, where L is the number of directions. A common choice for L is L = 8 that results in the set {0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°}. For each θ, the length of proposed window is selected from the countable and increasing set {h 1, h 2,…, h }. So, for the noisy observation y(k), we would have the following estimate:

For each θ and k, an appropriate value of h called h + is estimated using the nonlinear intersection of confidence intervals (ICI) rule. h + is the largest h from the h 1 < h 2 < ⋯

According to the ICI rule, D is defined using the following formula:

The largest s that leads to a nonempty value is called s +, and finally h +(k, θ) is obtained using h +(k, θ) = h .

Figure 7 shows an example of mentioned anisotropic window selection for a SEAD.

Figure 7

The red line shows the detected SEAD by an expert. The yellow circles show the isotropic windows with various radii. The green line illustrates the obtained anisotropic based on LPA-ICI rule.

Since applying LPA-ICI in each subband is a time consuming process, a fast version of the mentioned algorithm can be based on only applying LPA-ICI to low-pass subbands using L = 12 with an offset of 15° that results in the set {15°, 45°, 75°, 105°, 135°, 165°, 195°, 225°, 255°, 285°, 315°, 345°} in a 2D case. Since, in this case, each subband is extracting the information concentrated in a specific direction corresponding to {15° (195°), 45° (225°), 75° (255°), 105° (285°), 135° (315°), 165° (345°)}, the extracted h +(k, θ) = h that results from applying LPA-ICI to corresponding low-pass subband is used for obtaining the local parameters of kth pixel. For example, suppose that DCWT is used for 3 scales and we want to calculate the local parameters of coarsest scale for the oriented real subband around 45° (225°). For this reason, h +(k, 45°) and h +(k, 225°) are extracted from the results of applying LPA-ICI on the LL subband of real part (or imaginary part) of DCWT. Then, if we are in the jth scale, only 2 h +(k, 45°) pixel in direction of 45° and 2 h +(k, 225°) pixel in direction of 225° are used to extract the local parameters of kth pixel in this subband (Figure 8).

Figure 8

From left to right: imaginary LL subband of one slice of OCT data, the oriented (imaginary) subband around 45° (225°), h +(·, 45°) for LL subband, and h +(·, 225°) for LL subband extracted by applying LPA-ICI to the LL subband of imaginary part of DCWT. As indicated in the second image for p 1 = (46,20), p 2 = (47,21), and p 3 = (46,22) we would have h +(p 1, 45°) = 2, h +(p 1, 225°) = 3 (green dash), h +(p 2, 45°) = 1, h +(p 2, 225°) = 3 (orange dash), and h +(p 3, 45°) = 3, h +(p 3, 225°) = 3 (red dash).

A similar manner can be proposed in 3D case [67]. However, instead of using 2D direction θ , we use 3D direction (θ , φ ). As shown in Figure 9, in 2D case we use a circular sector for each direction while for 3D case a conical body is produced for direction (θ , φ ), and the sphere is covered (partly) using these cones. Similar to 2D case g , is defined and for each (θ , φ ) the best h called h +(k, θ , φ ) is obtained using ICI rule.

Figure 9

Comparison between a circular sector for direction θ in 2D case (b) with a conical body produced for direction (θ, φ) in 3D case (a).

Note that in order to incorporate the anisotropic window selection for each DCWT coefficient in our OCT denoising algorithm explained in Algorithm 1, instead of using a square window for parameter estimation, an anisotropic window is obtained for each coefficient k using the explained LPA-ICI method in this section, and only available data in this window are used for estimating a(k) and σ 11(k),  σ 12(k),  σ 21(k), and σ 22(k).

5. Experimental Results

In this section, we apply the proposed despeckling algorithm to OCT image noise reduction. For this reason, we use 20 three-dimensional OCT datasets in the presence of wet AMD pathology (SEAD) and use mean signal-to-noise ratio (MSNR) and contrast-to-noise ratio (CNR) as two quality measurements for OCT data. To calculate these measurements, we must define the region of interest (ROI). In this paper, we propose this region within the SEAD as illustrated in Figure 10. The base MSNR and CNR are defined as follows: where μ ROI shows the mean of ROI and σ indicates the standard deviation of a large region outside the ROI (noise ROI in Figure 10).

Figure 10

One slice from a sample OCT image and proposed ROIs for computation of MSNR and CNR reported in Table 1.

Table 1 shows the results of MSNR and CNR for proposed ROIs in OCT data using our algorithm. As discussed in Section 3, various shrinkage functions can be obtained using our algorithm based on applying log transformation before applying 3D DCWT (we use homomorphic prefix for this method and non-homomorphic when we do not use log transformation) and proposing AWGN or two-sided Rayleigh pdf for modeling noise in 3D DCWT domain (we name them BiGaussMixShrinkL and BiGaussRayMixShrinkL, resp.). Figures 11 and 12, respectively, show the results of applying non-homomorphic and homomorphic methods for (a slice of) depicted OCT image in Figure 10. In this figure, also in Table 1, we compare the results of nonlocal version of methods to show the effect of using anisotropic window selection technique. In order to show the SNR improvements, CNR curves for 156 selected ROIs have been depicted in Figure 13. It is clear that non-homomorphic BiGaussMixShrinkL method outperforms the others.

Table 1

The results of MSNR and CNR using several ROIs, shown in Figure 12.

Methods			MSNRROI1	MSNRROI2	CNR
Local (L) Nonlocal (NL)	Homomorphic (H) Nonhomomorphic (NH)	Gaussian noise (G)Two-sided Rayleigh noise (R)
L	H	G	7.00	15.76	8.76
NL	H	G	7.56	17.03	9.47
L	NH	G	12.27	27.76	13.49
NL	NH	G	10.77	22.73	11.95
L	H	R	5.89	13.11	7.22
NL	H	R	8.63	19.59	10.95
L	NH	R	10.75	22.55	11.81
NL	NH	R	10.88	23.05	12.17
	Original image		2.56	5.30	2.74

Figure 11

The results of applying homomorphic methods on proposed image in Figure 10. From top-left clockwise: despeckled data using BiGaussMixShrinkL, nonlocal BiGaussMixShrinkL, nonlocal BiGaussRayMixShrinkL, and BiGaussRayMixShrinkL.

Figure 12

The results of applying non-homomorphic methods on proposed image in Figure 10. From top-left clockwise: despeckled data using BiGaussMixShrinkL, nonlocal BiGaussMixShrinkL, nonlocal BiGaussRayMixShrinkL, and BiGaussRayMixShrinkL.

Figure 13

A comparison between CNR curves for 156 selected ROIs from OCT dataset.

Another way for evaluating the effect of our despeckling algorithm is the investigation of the intralayer segmentation results. Figure 14 shows a comparison between the segmented layers of a 650 × 512 × 128 Topcon 3D OCT-1000 imaging system using proposed method in [53]. It is clear that the first layer is detected truly after despeckling.

Figure 14

A comparison between the segmented layers of a 650 × 512 × 128 Topcon 3D OCT-1000 imaging system using proposed method in [53]. From left to right: original image, denoised image by nonlocal homomorphic BiGaussRayMixShrinkL method, and local homomorphic BiGaussRayMixShrinkL method.

6. Conclusion and Future Work

In this paper, we introduced a new noise reduction algorithm for 3D OCT data. We found new shrinkage functions employing a mixture of bivariate Gaussian for modeling wavelet coefficients in each subband of complex wavelets. The parameters of this mixture model are estimated locally using a shape-adaptive manner based on the special structure of OCT data. We also used this model for denoising of other kinds of noise. Experiments show that our model has better results than other methods visually and in terms of PSNR especially for the crowded images. In this paper, we suppose that the parameters of EM algorithm, in extracted windows are constant. It is possible to improve the EM algorithm, for example, by using recurrence equations. It is possible that we only propose the main section of data (between the first and last layers) containing retina layer information and apply our algorithm on the selected data to improve the speed and performance of denoising process.

Using 3D DCWT instead of other transforms such as 3D DWT is a main reason for improvement of the denoising results [30]. In [27], it has been shown that other kinds of oriented transforms such as steerable pyramid decomposition can produce better denoising results. However, for 3D case, 3D transforms that are applied on whole 3D data (not slice by slice) such as surfacelet [68] and 3D discrete curvelet [69] can be investigated.

Table 2

PSNR (in dB) values of test images for different nonstationary noise levels.

Noise parametersσg(i) = k 0 s(i) + k 1	Lena				Boat				Barbara
	Noisy Image	Soft thresh. [54]	Proposed method in [55]	Our method	Noisy image	Soft thresh. [54]	Proposed method in [55]	Our method	Noisy image	Soft thresh. [54]	Proposed method in [55]	Our method
k 0 = 0.05, k 1 = 4	27.72	34.13	34.61	35.60	27.51	32.44	32.59	33.26	27.94	31.99	32.20	33.56
k 0 = 0.1, k 1 = 4	23.48	31.62	32.49	33.31	23.22	29.95	30.23	30.75	23.69	29.05	25.81	30.91
k 0 = 0.2, k 1 = 4	18.49	27.50	29.65	30.52	18.20	26.77	27.60	28.13	18.71	25.81	25.94	27.91

Table 3

Comparison between PSNRs (in dB) of denoised images with Fast TI Haar algorithm [55] and BiGaussMixShrinkL.

	Lena	Boat	Barbara	Confocal Phantom	Shep Logan Phantom	Bowl
Noisy image	27.22	27.05	27.49	35.74	47.68	28.21
Fast TI Haar	32.11	29.30	26.59	44.49	60.63	46.79
BiGaussMixShrinL	39.88	37.57	37.78	47.36	64.65	47.09