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.


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][2][3][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][7][8][9][10][11][12][13][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][7][8][9][10][11][12][13][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][8][9][10][11][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 "lowpass"-"high-pass" (LH) subbands (and some of it in HH subbands). Therefore, an increased threshold in the vertical subbands using a constant multiplier ( = 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 noisefree 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.

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][27][28][29][30][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][44][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.

Discrete Complex Wavelet Transform (DCWT).
The wavelet based image denoising consists of the following steps.
(1) Signal transformation of the noisy observation.
(2) Modification of the noisy wavelet coefficients based on some criteria.
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 noisefree 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 waveletbased 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.
(1) In the neighborhood of an edge, the DWT produces both large and small wavelet coefficients. In contrast, the magnitudes of DCWT coefficients are International Journal of Biomedical Imaging 3 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.
(2) 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.
(3) Since the DWT coefficients in each subband are produced via critical sampling after using nonideal lowpass 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.
(4) 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][49][50][51][52]. Figure 1 shows a comparison between subbands of 3D DWT and 3D DCWT.

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.

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.

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: ( ) ( ( )) = ( ) 1 ( ( )) + (1 − ( )) 2 ( ( )) = ( ) ( 2 1 ( )/2 2 11 ( ))−( 2 2 ( )/2 2 12 ( )) 2 11 ( ) 12 ( ) where ( ) ∈ [0, 1], 11 ( ), 12 ( ), 21 ( ), 22 ( ) are the mixture model parameters. For each random bivariable, the second component is the parent of the first component; for example, 2 ( ) represent, the parent of 1 ( ) at the same spatial position as the th wavelet coefficient 1 ( ) 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. International Journal of Biomedical Imaging After substitution of mixture model in the definition of ( 1 ( ) 2 ( )), we can see that this pdf represents two uncorrelated random variables as follows: Interestingly, the marginal pdf of ( ) for = 1, 2 is the mixture of two univariate Gaussian pdf with local parameters [44], It is easy to see that 1 ( ), 2 ( ) are not independent; that is, See Appendix B for more explanation.

Local EM Algorithm.
To characterize the parameters in (1), it is necessary to have the parameters 11 ( ), 21 ( ), 12 ( ), 22 ( ), and ( ). 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 International Journal of Biomedical Imaging algorithm is able to obtain separate parameters for each voxel by the implementation of EM algorithm in each window ( ) centered at ( ). This iterative algorithm has two steps. Assuming that the observed data ( ) for = 1, . . . , , the -step calculates the responsibility factors for each data as follows: The -step updates the parameters ( ), 11 ( ), 12 ( ), where is the number of coefficients in the square window ( ) centered at ( ). The variances 11 ( ), 12 ( ), 21 ( ), and 22 ( ) are computed by [40]:

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 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][8][9][10][11][12]18] usually at first transform data into log domain, and suppose that noise in log domain is AWGN: where in this paper shows 3D DCWT operator. So, we can write where ( ), ( ), and ( ) are, respectively, the th noisefree 3D DCWT coefficients, noisy 3D DCWT coefficients, and noise in the 3D DCWT domain.
Recently, it has been reported [56][57][58][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 ( ), ( ), and ( ) are, respectively, the th noise-free 3D DCWT coefficients, noisy 3D DCWT coefficients, and noise in the 3D DCWT domain. Since speckle noise can be modeled as a unit-mean random process independent of the noise-free data, we would have [ ( ( − 1))] = 0, and also it can be easily shown [58] that [ ( ) ( ( − 1))] = 0 which means that ( ) and ( ) are zero-mean uncorrelated random variables. If ( ), ( ), and ( ) show the parent coefficients of ( ), ( ), and ( ), 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 th 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 -distribution, Rayleigh, Weibull, log-normal, and Nakagami distributions have been proposed [57,58,[60][61][62][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: International Journal of Biomedical Imaging  (17), we must know the prior distribution of 3D DCWT coefficients, that is, ( ) ( ( )). Defining Gauss( , ) := exp(− 2 /(2 2 ))/( √ 2 ), if we suppose that ( ), ( ) are independent Gaussian pdf with variances 1 ( ) and 2 ( ), the following bivariate Gaussian pdf with local variances can be proposed for the noise-free wavelet coefficients: In this case, ( ) and ( ) are uncorrelated and independent, and therefore the MMSE estimator of ( ), ( ) yields the shrinkage function corresponding to univariate Gaussian pdf, that is, Wiener filter [21] as follows: And so we can writê Similarly, if we choose two-sided Rayleigh pdf for noise distribution, the following estimator is obtained [44]: where And so we can writê Suppose that the input noise variance is known. To implement (20) or (23), we must know the parameter of the prior ( ) (suppose that ( ) = ( )). 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][43][44][45][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 ( ) 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 8 International Journal of Biomedical Imaging ( ) can be formed based on a local neighborhood ( ). In simplest case, we can use a square window ( ) centered at ( ) and suppose that in this window the variance is approximately constant. Then, an empirical estimate for ( ) can be obtained as follows: where is the number of coefficients in ( ) 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)).
Similarly, using (21), we can obtain numerators of (25), and finally (25) for two-sided Rayleigh noise is obtained aŝ where We call this bivariate local shrinkage function as BiGauss-RayMixShrinkL. Figure 5 shows this shrinkage function with sample constant parameters.
For implementation of our denoising algorithm, we must estimate the parameters ( ) for , = 1, 2, and ( ) (that are for noise-free data) from noisy observation. For AWGN, the noisy observation would be a Gaussian mixture model with parameters ( ),√ 2 + 2 ( ) for , = 1, 2. So, the following local EM algorithm is used to obtain the parameters.

M-step
where is the number of coefficients in the window ( ) 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.

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][43][44][45][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.

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. Step 3. Substituting the Final Parameters in Step 2 in Shrinkage Function (30) for AWGN (after obtaining ( ( )) using (31)) and Shrinkage Function (32) for two-sided Rayleigh noise (after obtaining ( ( )) using (33)).   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.

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 ℎ, 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 is the number of directions. A common choice for is  Table 1.
For each and , an appropriate value of ℎ called ℎ + is estimated using the nonlinear intersection of confidence intervals (ICI) rule. ℎ + is the largest ℎ from the ℎ 1 < ℎ 2 < ⋅ ⋅ ⋅ < ℎ provided that the estimated data using ℎ + does not have noticeable difference with the estimated data with smaller ℎ's. For this reason, the following confidence intervals are defined: where is the smoothing parameter and est ℎ , ( ) shows the variance of est ℎ , ( ) and is obtained as follows: where (⋅) shows the power spectral density function and ℎ, ( ) is the Fourier transform of ℎ, , and for a white random process, (40) is simplified to According to the ICI rule, is defined using the following formula: The largest that leads to a nonempty value is called + , and finally ℎ + ( , ) is obtained using ℎ + ( , ) = ℎ + . Figure 7 shows an example of mentioned anisotropic window selection for a SEAD.
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 = 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 ℎ + ( , ) = ℎ + that results from applying LPA-ICI to corresponding lowpass subband is used for obtaining the local parameters of th 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, ℎ + ( , 45 ∘ ) and ℎ + ( , 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 th scale, only 2 −1 ℎ + ( , 45 ∘ ) pixel in direction of 45 ∘ and 2 −1 ℎ + ( , 225 ∘ ) pixel in direction of 225 ∘ are used to extract the local parameters of th pixel in this subband (Figure 8). 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 ℎ, , , is defined and for each ( , ) the best ℎ called ℎ + ( , , ) is obtained using ICI rule.
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 using the explained LPA-ICI method in this section, and only available data in this window are used for estimating ( ) and 11 ( ), 12 ( ), 21 ( ), and 22 ( ).

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

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 16 International Journal of Biomedical Imaging 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.
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.

A. Directional Selectivity Property of 2D DCWT
Since DWT in 2D domain is produced using separable (rowcolumn) implementation, it has a poor directional selectivity.  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.
For example, the HH wavelet is the product of the highpass functions along the first and second dimensions. Because DWT uses real filters, the HH wavelet mixes +45 ∘ and −45 ∘ orientations that results in the checkerboard artifact because it fails to isolate these orientations. 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 2D domain is supported in only one quadrant of the 2D frequency plane. Figure 15 illustrates a comparison between subbands of 2D DWT and 2D DCWT. Figure 16 shows the pdf of a bivariate Gaussian mixture model for sample parameters. We can see the marginal distribution produced from the model in this figure.

C. Other Kinds of Noise
In this appendix, we briefly explain the abilities of the proposed denoising algorithm in this paper for other kinds of noise.

18
International Journal of Biomedical Imaging C.1. Stationary Noise. We tested the shrinkage function BiG-aussMixShrinkL for stationary noise and compared it with other methods such as ProbShrink [34], BiShrink [45], and BLS-GSM [27] and found that our algorithm outperforms these techniques visually and in terms of peak signal-to-noise ratio (PSNR) for various levels of noise. For example, our algorithm for crowded images preserves details of images while BiShrink [45] in some cases produces blurred images (e.g., compare the area around the corresponding arrows in Figure 17). In addition, in [40] using other bivariate mixture models such as bivariate Laplacian was examined. We compared the results of using our model with method in [40], and our simulations show that our algorithm is faster and outperforms the reported shrinkage functions in [40] for some images. For example, for a 512 × 512 8-bit grayscale Barbara image, an improvement of 0.6 dB is obtained for noise level of 30, and BiGaussMixShrinkL is two times faster.

C.2. Nonstationary Noise. This section presents nonstationary noise reduction examples in complex wavelet domain.
Although the stationary noise model is able to simplify the implementation of denoising algorithms, the statistical properties of the noise are not always accurately described with this assumption. For example, in some applications, the noise statistics are spatially varying and the noise power varies between pixels or samples. In these cases, the nonstationary noise assumption is more reasonable and can improve the denoising results. For example, we contaminate three 512 × 512 grayscale images, namely, Lena, Boat, and Barbara using signal-dependent Gaussian noise with variance ( ) that is defined as a linear function of the pixel intensities ( ) as [54]: Since the variance of each noise component is spatially varying with the corresponding content of signal, the nonstationary processes are able to model the statistical properties of this noise. A comparison between the denoised image using soft thresholding, proposed method in [54], and the denoised image using a mixture of two bivariate Gaussian pdfs with local parameters (BiGaussMixShrinkL) for different noise levels can be seen in Table 2. In this table, the highest PSNR value is bolded. We can see from the table that our proposed algorithm has the better results compared to others especially for the Barbara image (which contains details) in the high-level noise.
In [54], it has been shown that the proposed algorithm based on the signal-dependent Gaussian noise can also be effective for the reduction of Poisson noise. On this base, we use our algorithm for noise reduction of images corrupted by Poisson noise generated using corresponding voxel intensities. For this reason, we use the software provided on http://willett.ece.wisc.edu/software.html to compare our method with Fast TI Haar algorithm [55]. The PSNR of grayscale images Lena, Boat, Barbara, Confocal Microscopy Phantom, Bowl, and Shepp Logan Phantom can be seen in Table 3. We can see that our algorithm outperforms the Fast TI Haar algorithm. Figure 18 illustrates a comparison between the denoised images produced from two algorithms. It is clear that our algorithm has better results especially for the crowded images. In fact, the Fast TI Haar algorithm has reasonable performance for soft images such as Confocal Microscopy Phantom, but since this algorithm blurs the produced images, for images with details such as Barbara, highfrequency features will be removed and we lose the important information. (LPA-ICI) method for applying shape adaptive window selection around each International Journal of Biomedical Imaging  Figure 17: (a) shows a part of Barbara image denoised using BiShrink [45] for stationary noise with = 40 and (b) shows denoised image using our method. Comparing the area around the corresponding arrows, we understand that our method is able to better preserve the details of images.