Hindawi Publishing Corporation International Journal of Biomedical Imaging Volume 2013, Article ID 417491, 23 pages http://dx.doi.org/10.1155/2013/417491

Research Article

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

Hossein Rabbani,1,2 Milan Sonka,2 and Michael D. Abramoff2

1 Biomedical Engineering Department, Medical Image & Signal Processing Research Center, Isfahan University of Medical Sciences, Isfahan 81745, Iran

2 The Iowa Institute for Biomedical Imaging, The University of Iowa, Iowa City, IA 52242, USA

Correspondence should be addressed to Hossein Rabbani; h_rabbani@med.mui.ac.ir Received 7 January 2013; Revised 1 June 2013; Accepted 21 June 2013 Academic Editor: Michael W. Vannier

Copyright © 2013 Hossein Rabbani et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 x 512 x 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 ofGaussian 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], a-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, Mih9ak [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.

(1) Signal transformation of the noisy observation.

(2) Modification of the noisy wavelet coefficients based on some criteria.

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

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

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

(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, 4852]. Figure 1 shows a comparison between subbands of 3D DWT and 3D DCWT.

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 x 512 imagesfrom the Corelimage database [42]. Usually this property can be modeled using proper bivariate pdfs.

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:

Pw(k) (w(k)) = a (k) Pl (w(k)) + (l-a(k)) p2 (w(k))

a (k)e(

2non (k) a12 (k)

(1-a (k)) e(-"2(k)/2o221(k)h(w22(k)/2a222(k))

2no21 (k) o22 (k)

where a(k) e [0,1], <Ju(k), o12(k), o21(k), o22(k) are the mixture model parameters. For each random bivariable, the second component is the parent of the first component; for example, w2(k) represent, the parent of w1(k) at the same spatial position as the kth wavelet coefficient w1 (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.

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

ps( 1) f(w | S = 1)

ps(2) = 1-ps(1) f(w | S = 2)

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

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

E(Wi (k)W2 (k))

= U Wi (k) W2 (k) pw(k) (w (k)) dw (k) = (1 - a (k))

x || w1 (k) w2 (k) p2 (w (k)) dw1 (k) dw2 (k)

+ U(k)UWi (k)W2 (k)Pi (W(k))dWi (k)dW2 (k)

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

Pw(k) (w(k))dw2 (k)

exp (-w2 (k)!2a\i (k))

= a (k)-

+ (\-a(k))

exp (-W2, (к) /2а221 (k)) a21 (к)

a11 (к)

It is easy to see that w1(k), w2(k) are not independent; that is,

Рщк) (w (k)) = pWi(k) (Wi (k)) pWi(k) {w2 (k)). (4) See Appendix B for more explanation.

2.2. Local EM Algorithm. To characterize the parameters in (1), it is necessary to have the parameters a11(k), o21(k), a12(k), o22(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 w(k) for k = 1,...,N, the £-step calculates the responsibility factors for each data as follows:

гг (к)

а(к)р! (w (к))

а(к)Рг (w(k)) + (l-a(k))p2 (w(k)) r2 (к)^-!-Гг (к).

The M-step updates the parameters a(k), агг(к), <J12(k), о2г (к), and o22(k). a(k) is computed by

M £ Гг

jeN(K)

where M is the number of coefficients in the square window N(k) centered at w(k).

The variances <Ju(k), o12(k), o21(k), and o22(k) are computed by [40]:

oL (k)

IjeN(K) ri (j) Wi (k)

^jeN(K)

ri (j)

i, m = 1,2. (7)

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:

* (i) = s (i) g (i),

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:

W (log x (i)) = W (log 5 (i)) + W (log g (i)), (9)

where in this paper W shows 3D DCWT operator. So, we can write

y(k) = w (k) + n(k),

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-homo-morphic 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

W (x (i)) = W(s (i) g (i)) = W(s (i) + s (i) (g (i) - l)) = W(s(0) + W(s(0 Ы0-!)).

Again we can write

y(k) = w(k) +n(k),

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 wJk), y„(k), and np(k) show the parent coefficients of w(k), y(k), and n(k), respectively, we can write

yp (k) = wp (k) + Hp (k).

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:

y(k) = w(k) +n(k),

where w(k) = (w(k),wp(k)), ~y(k) = (y(k), yp(k)),andn(k) = (n(k),np(k)) 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 ^-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:

Pn (n(k)) = exp (-

P„ (n(k)) =

2nol К (k)n2 (k)\

n2 (к) + n22 (к)

п2 (к) + п2 (к) 2а2

where a2 = 2a2 shows the noise variance.

Now our goal is the estimation of w(k) from ~y(k) = w(k) + n(k), where n(k) 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:

(k) = U w (k) pW(k)\y(k) (w (к) | у (к)) dw (к)

,,, Py(k)\w(k) (У(к) 1 w(k)) pw(k)(w (k))

w(k)-j——-dw(k)

Py(k) (y (k))

w(k) pW)Щк) (у (к) I w(k)) pw{k) (w(k))dw(k)

Py(k) (y (k))

w(k) pmmk) (у (к) I w(k)) рЩк) (w(k))dw(k)

Я Py(k)Щк) (У(k) 1 w (k)) рЩк) (w (k)) dw (k) w (k) Pw (У (k) - w (k)) pw{k) (w (k)) dw (k)

Я Pn (y (k) - w (k)) p^{k) (w (k)) dw (k)

3.1. Denoising Based on Modeling Noise-Free Data by Bivari-ate Gaussian PDF with Local Variance. In order to solve (17), we must know the prior distribution of 3D DCWT coefficients, that is, pW(k)(w(k)). Defining Gauss(x, a) := exp(-x2/(2a2))/(aV2n), if we suppose that w(k), wp(k) are independent Gaussian pdf with variances a1 (k) and a2(k), the following bivariate Gaussian pdf with local variances can be proposed for the noise-free wavelet coefficients:

Pw(k) (w(k)) = Pw(k) (w(k)) ■ Pwp(k) (wP (k))

= Gauss (w (k) ,a (k)) ■ Gauss (wp (k), ap (k))

Pw(k) (w(k))

w2 (k) 2a2 (k)

w\ (k)

2aj (k)

x(2na(k)ap (k))

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

w(k) =

JJ w (k) ft (y (k) -w(k)) pMk) (w (k)) dw (k) JJp^ (y(k)-w(k))pmfM) (w(k))dw(k)

a2 (k)

= y (k)

a2 (k) + al '

And so we can write

w(k) =

y (k) a2 (k) yP (k) a2p (k) a2 (k)+a2n ' a2p (k) + a2n

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

(k)=2z(k)^l(2-a^

a2 (k)z2 (k)

x (erfcx (z(k)) - erfcx (-z (k)))

+ (2 + z(k) ^nerfcx(-z(k))

a2 a2 (k)

- z (k) ^n erfcx (z (k)))\

z(k) = l«±J_L

Z(k) ^2/V\ \111„2 1^2 i

a2 (k)\2/a2 +2/a2 (k)

2 (c :(u)=^\ Vn Jo

-t -2tu 1.

And so we can write

Hk)==((2z(k)V2(2-a-2))

,Jn(1 a2 (k)z2 (k)

x (erfcx (z (k)) - erfcx (-z (k)))

y a2 ' a2 (k) x(2 + z(k) Vnerfcx(-z(k))

-z(k) Vn erfcx(z(k))) I

2z(k)V2( 2

, Jn(i at(k)z(k)

x (erfcx (zp (k)) - erfcx (-zp (k)))

1 1 a2 + auk)

x (2 + zp (k) Vn erfcx(-zp (k)) kVn erfc x(zp k)))

Suppose that the input noise variance is known. To implement (20) or (23), we must know the parameter of the prior a(k) (suppose that a(k) = ap(k)). Mih9ak 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 a(k) or in [64] the variable size of the locally adaptive window is obtained using a region-based approach. However, for each data point y(k), a simple estimation of

o(k) can be formed based on a local neighborhood N(k). In simplest case, we can use a square window N(k) centered at ~y(k) and suppose that in this window the variance is approximately constant. Then, an empirical estimate for o(k) can be obtained as follows:

- = ш I (у2 (j) + y2P (J))-°ï>

jeN(k)

where M is the number of coefficients in N(k) and on can be estimated by [4] on = 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

(k) = (jjw(k)pw (y(k)-w(k))

x [a(k)p1 (w(k))

+ (l-a(k))P2 (w(k))]dw(k)

jj^ (y(k)-w(k))

x [a(k)p1 (w(k))

+ (l-a(k))P2 (w(k))]dw(k))

а (к) ^ w (k) p* (у (k) -w(k)) Pl (w (к)) dw (k) a(k)01 (y(k)) + (l-a(k))02 (y(k))

+ ((1-a(k))

X jj ^ (k) ^ (У (k) - W (k)) P2 (W (k)) dW (k)

x(a(k)gi (у(к)) + (1-а(к))д2 (y(k)))-1,

0, (У (k)) = jj Pn (У (k) - w (k)) pt (w (k)) dw (k),

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

g, (y(k)) = exp (-(1/2)

*(У2 (k)/(°2n + 4(k))

+Ур (k) / (°n + 4 (k)))) (27)

i = 1,2.

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

g, (y(k))

exp (-y2 (k)/2afi (к) - y2p (к)/2^2 (к))

8n(l + afi (к) /a2) (l + tû (к) /<*2) (к) аа (к)

x (2 + Zj (к) ^ erfc x(-zt (к))

- Zj (к) ^пerfcx(z{ (к))) x(2 + zip (к) ^п erfc x(-zip (к))

- zip (к) ^ erfc x(zip

(к))),

i = 1,2,

м у(к)

z' (к)-^йк) y

2/a2 + (2/af1(k))'

ОЛ У> (k) Z» (k)=^(k) 1

2/a2 + (2/a2 (k))'

i = 1,2,

i = 1,2.

i = 1,2.

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

w(k) = (o2n (k)l(o2n (k) + ol)

+R (y(k)) (oh (k) | (oh (k) + o2„))) (30) x(l + R(y(k)))-1 y(k),

Figure 4: A shrinkage function produced from BiGaussMixShrink for sample parameters.

R(y(k)) = ( ( (1-a(k))

Ур (к)

+ < (k) °2n + 4 (k)

x( P2n + o\ 1 (k))(o2n + a\2 (к))У / (к)

x(( a(k) exp (^ , M 2\a2 +a2n (k)

Ур (к)

an + <*22 (k)

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

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

x (izi (k)^-^) + ^^ o\i (k)z22 (k)

x (erfc x (z2 (k)) - erfc x (-z2 (k)))

a2 o2i (k)

x(2 + z2 (k) -^nerfcx (-z2 (k)) -z2 (k) ^nerfcx (z2 (k)))

+ R(y(k)) + 1 + R(y(k))

x ( 2z2 (k)V2\2-

n ( л °\2 (k) ¿I (k)

x (erfc x (z2 (k)) - erfc x (-z2 (k)))

a2 o\2 (k)

x (2 + z2 (k) -^n erfc x (-z2 (k))

-z2 (k) ^n erfc x (z2 (k)))\ ,

where R(y(k))

1 - a(k) a(k)

(1+0*2 (к) fa2) (1+^2,2 (к) /<*2) *22 (к) a22 (к) x (1+ oh (к) /а2) (1 + (к) /а2) (к) 22 (к)

exp (-у2 (к) /2о\2 (к) - у2р (к) /2а\2 (к)) x exp (-у2 (к) /2о22 (к) - у2р (к) /2о2и (к))

x((2 + z2 (к) Vл erfcx (-z2 (к)) -z2 (к) Vn erfc x (z2 (к))) x (2 + z2p (k) Vn erfc x (-z2p (k))

-Z2P (k) Vn erfc x(Z2p m

x ((2 + z2 (k) Vn erfcx (-z2 (k)) - z2 (k) Vn erfc x(z2 (k)))

X{2 + Z1P (k) Vn erfc x(-zlp (k))

-Zip (k)^n erfc x(zlp (k)))) \

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 o^ (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),^a'^ + a?.(k) for i, j = 1,2. So, the following

local EM algorithm is used to obtain the parameters.

E-step

a(k)g2 (y (k))

a (k) gi (y(k)) + (1- a (k)) g2 (y (k))' (34) r2(k)^-\-rx (k).

M-step

a(k)^~M ^ ri(j)'

1V1 jeN(K)

IjeN(K) ri (j) / (k) ZjeN(K) ri (J)

IjeN(K) ri (i) y2p(k)

- o„

jiN(K) ri

m = 1,2, (36)

m=1,2, (37)

where M is the number of coefficients in the window N(K) centered at ~y(k). 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.

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

Figure 5: A shrinkage function produced from BiGaussRayMixSh-rink for sample parameters.

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 despeck-ling 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 canbe used for the image-guided retinal therapy. Everyyear, 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 neovascular-ization (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 intravit-real 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

Summarization of the proposed denoising algorithm

Step 1. Complex Wavelet Transformation of Noisy Image to Find ~y(k) and an.

Step 2. Estimation of the Prior (Finding Mixture Model Parameters in Each Subband).

2.1. Initialization for a(k) and all(k),al2(k),a2l(k),a22(k).

2.2. Using (34) to find r1(k),r2(k).

2.3. Using (37) to update a(k) by substituting rl(k),r2(k) from Step 2.2.

2.4. Updating the parameters of prior, all(k),al2(k),a2l (k),a22(k), using (36) and (37).

2.5. Finding gi(y(k)) from (27) for AWGN and (28) for two-sided Rayleigh noise using updated value in Step 2.3.

2.6. Iteration Step 2.2 to 2.4 until parameter convergence.

Step 3. Substituting the Final Parameters in Step 2 in Shrinkage Function (30) for AWGN

(after obtaining R(y(k)) using (31)) and Shrinkage Function (32) for two-sided Rayleigh noise (after obtaining R(y(k)) using (33)).

Step 4. Inverse Complex Wavelet Transformation.

Algorithm 1: Outline of the proposed denoising algorithm.

Figure 6: Macular OCTs and detected SEADs by an expert.

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.

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

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 p1 = (46,20), p2 = (47,21), and p3 = (46,22) we would have h+ (p1,45°) = 2, h+ (p1,225°) = 3 (green dash), h+(p2,45°) = 1, h+(p2,225°) = 3 (orange dash), and h+(p3,45°) = 3, h+(p3,225°) = 3 (red dash).

we can see in this figure, the most important information of OCT data (about retina layers) is located in the center of OCT images.

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 gh e that are obtained using local polynomial approximation (LPA) are employed. The d indicates the direction of filter that is a member of countable set [d1,d2,... ,dL}, 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 d, the length of proposed window is selected from the countable and increasing set {h1,h2,... ,hj}. So, for the noisy observation y(k), we would have the following estimate:

xte (k) = 9h,e (k) *y(k).

Figure 9: Comparison between a circular sector for direction 0 in 2D case (b) with a conical body produced for direction (0, <p) in 3D case

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

For each d 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 hx < h2 < ••• < hj provided that the estimated data using h+ does not have noticeable difference with the estimated data with smaller h's. For this reason, the following confidence intervals are defined:

xht (k) - Rax

xh\ (k) + Ra.

x?„ №

where R is the smoothing parameter and a^n (k) shows the

hs,e( '

variance of x^1 (k) and is obtained as follows:

alt (fc) = | P^ (fc) | (f) df=\py (f) Ghe (f) df, (40)

where P(-) shows the power spectral density function and Ghg(f) is the Fourier transform of ghg, and for a white random process, (40) is simplified to

^ (k) = \ (J) df = °2n 19h,e (k). (41)

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

Ds = ff,

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

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 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,9) = hs+ 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°).

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

Local (L) Nonlocal (NL) Methods Homomorphic (H) Nonhomomorphic (NH) Gaussian noise (G) Two-sided Rayleigh noise (R) MSNRROU MSNRroi2 CNR

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 1311 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

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

Noise parameters

Barbara

Soft :hresh [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

34.61 35.60 27.51 32.44 32.59 33.26 27.94 31.99 32.20 33.56

32.49 33.31 23.22 29.95 30.23 30.75 23.69 29.05 25.81 30.91

29.65 30.52 18.20 26.77 27.60 28.13 18.71 25.81 25.94 27.91

fc0 = 0.05, kl =4 k0 = 0.1, kl = 4 k0 = 0.2, fc, =4

34.13 31.62 27.50

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

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 2j-1h+(M5°) pixel in direction of 45° and 2J-1h+(fc, 225°) pixel in direction of 225° are used to extract the local parameters of kth pixel in this subband (Figure 8).

A similar manner can be proposed in 3D case [67]. However, instead of using 2D direction 0;, we use 3D direction (di,fi). 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 (di,fi), and the sphere is covered (partly) using these cones. Similar to 2D case Q is defined and for each (0;, <p;) the best h called h+(k, dt, <p;) 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 k using the explained LPA-ICI method in this section, and only available data in this window are used for estimating a(k) and a11(k),a12(k),a21(k), and a22(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:

MSNRroi =

CNR = |MSNRroi1 - MSNRROI2|

where ^ROI shows the mean of ROI and a 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 x 512 x 128 Topcon 3D OCT-1000 imaging system using proposed method in [53]. It is clear that the first layer is detected truly after despeckling.

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

(c) (d)

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.

Appendices

A. Directional Selectivity Property of 2D DCWT

Since DWT in 2D domain is produced using separable (row-column) implementation, it has a poor directional selectivity.

15 -1-1-1-1-1-1-r

0 -,-,-,-,-,-,-,-

10 20 30 40 50 60 70 Number of selected ROIs

—•— Original data

Homomorphic nonlocal Gaussian method

- Homomorphic local Gaussian method

Homomorphic nonlocal Rayleigh method

- Homomorphic local Rayleigh method

Nonhomomorphic nonlocal gaussian method

- Nonhomomorphic local gaussian method

Nonhomomorphic nonlocal rayleigh method - Nonhomomorphic local rayleigh method

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

(a) (b) (c)

Figure 14: A comparison between the segmented layers of a 650 x 512 x 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.

B. A Sample Bivariate Gaussian Mixture Model

Figure 16 showsthe pdfofabivariateGaussianmixture 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.

Figure 15: A comparison between subbands of DWT and DCWT. (a) The wavelets in the space domain (LH, HL, and HH). (b) The idealized support of the Fourier spectrum of each wavelet in the 2D frequency domain. We can see the checkerboard artifact of the third wavelet. (c) The complex wavelets in the space domain. (d) The idealized support of the Fourier spectrum of each wavelet in the 2D frequency plane. The absence of the checkerboard phenomenon is observed in both the space and frequency domains [48].

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 x 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 nonstation-ary 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 x 512 grayscale images, namely, Lena, Boat, and Barbara using signal-dependent Gaussian noise with variance og(i) that is defined as a linear function of the pixel intensities s(i) as [54]:

ag (i) = k0s(i) + k1. (C.1)

Since the variance of each noise component is spatially varying with the corresponding content of signal, the nonsta-tionary 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, high-frequency features will be removed and we lose the important information. (LPA-ICI) method for applying shape adaptive window selection around each

Single Gaussian pdf with ff = 1 Univariate gaussian mixture model with ffi = 2, o^ = 8, a = 0.5

- Gaussian mixture pdf

- Single Gaussian with high variance

Single Gaussian with low variance

Single bivariate Gaussian pdf with oi = 1,02 = 2

Bivariate Gaussian mixture model with an =2,012 = 3, 02 = 4, ff22 = 5, a = 5

Figure 16: The pdf of a bivariate Gaussian mixture model for sample parameters and its marginal distribution.

(a) (b)

Figure 17: (a) shows a part of Barbara image denoised using BiShrink [45] for stationary noise with an = 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.

(e) (f) (g) (h)

Figure 18: (a-d) show denoising results for Confocal Microscopy Phantom: from left to right: noise-free image, noisy image, and denoised image with our model and denoised image with Fast TI Haar algorithm. (e-h) show from left to right parts of denoised Barbara image with BiGaussMixShrinkL method and parts of denoised Barbara image with Fast TI Haar algorithm.

References

[1] D. Huang, E. A. Swanson, C. P. Lin et al., "Optical coherence tomography," Science, vol. 254, no. 5035, pp. 1178-1181,1991.

[2] J. M. Schmitt, "Optical coherence tomography (OCT): a review," IEEE Journal on Selected Topics in Quantum Electronics, vol. 5, no. 4, pp. 1205-1215, 1999.

[3] M. A. Choma, M. V. Sarunic, C. Yang, and J. A. Izatt, "Sensitivity advantage of swept source and Fourier domain optical coherence tomography," Optics Express, vol. 11, no. 18, pp. 2183-2189, 2003.

[4] M. E. J. van Velthoven, D. J. Faber, F. D. Verbraak, T. G. van Leeuwen, and M. D. de Smet, "Recent developments in optical coherence tomography for imaging the retina," Progress in Retinal and Eye Research, vol. 26, no. 1, pp. 57-77, 2007.

[5] J. M. Schmitt, S. H. Xiang, and K. M. Yung, "Speckle in optical coherence tomography," Journal ofBiomedical Optics, vol. 4, no. 1, pp. 95-105, 1999.

[6] D. L. Marks, T. S. Ralston, and S. A. Boppart, "Speckle reduction by I-divergence regularization in optical coherence tomography," Journal ofthe Optical Society ofAmerica A, vol. 22, no. 11, pp. 2366-2371, 2005.

[7] Z. Jian, Z. Yu, L. Yu, B. Rao, Z. Chen, and B. J. Tromberg, "Speckle attenuation in optical coherence tomography by curvelet shrinkage," Optics Letters, vol. 34, no. 10, pp. 1516-1518, 2009.

[8] 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," Optics Express, vol. 18, no. 2, pp. 1024-1032, 2010.

[9] S. Chitchian, M. A. Fiddy, and N. M. Fried, "Denoising during optical coherence tomography of the prostate nerves via wavelet shrinkage using dual-tree complex wavelet transform," Journal ofBiomedical Optics, vol. 14, no. 1, Article ID 014031, 2009.

[10] V. Zlokolica, Lj. Jovanov, A. Pizurica, and W. Philips, "Wavelet-based denoising for OCT images, Interaction between image processing, optics and photonics," in Symposium on Optical Science and Technology, Proceedings of SPIE, August 2007.

[11] V. Gupta, C. C. Chan, C.-L. Poh, T. H. Chow, T. C. Meng, and N. B. Koon, "Computerized automation of wavelet based denoising method to reduce speckle noise in OCT images," in Proceedings of the 5th International Conference on Information Technology and Applications in Biomedicine, pp. 120-123, May 2008.

[12] P. Puvanathasan and K. Bizheva, "Speckle noise reduction algorithm for optical coherence tomography based on interval type II fuzzy set," Optics Express, vol. 15, no. 24, pp. 15747-15758, 2007.

[13] A. Pizurica, L. Jovanov, B. Huysmans et al., "Multiresolution denoising for optical coherence tomography: a review and evaluation," Current Medical Imaging Reviews, vol. 4, no. 4, pp. 270-284, 2008.

[14] H. M. Salinas and D. C. Fernandez, "Comparison of PDE-based nonlinear diffusion approaches for image enhancement and denoising in optical coherence tomography," IEEE Transactions on Medical Imaging, vol. 26, no. 6, pp. 761-771, 2007.

[15] K. M. Yung, S. L. Lee, and J. M. Schmilt, "Phase-domain processing of optical coherence tomography images," Journal of Biomedical Optics, vol. 4, no. 1, pp. 125-136,1999.

[16] J. Rogowska and M. E. Brezinski, "Evaluation of the adaptive speckle suppression filter for coronary optical coherence tomography imaging," IEEE Transactions on Medical Imaging, vol. 19, no. 12, pp. 1261-1266, 2000.

[17] A. Ozcan, A. Bilenca, A. E. Desjardins, B. E. Bouma, and G. J. Tearney, "Speckle reduction in optical coherence tomography images using digital filtering," Journal of the Optical Society of America A, vol. 24, no. 7, pp. 1901-1910, 2007.

[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," Optics Letters, vol. 29, no. 24, pp. 28782880, 2004.

[19] S. G. Mallat, A Wavelet Tour of Signal Processing, Academic Press, San Diego, Calif, USA, 1998.

[20] B. Vidakovic, Statistical Modeling by Wavelets, John Wiley & Sons, New York, NY, USA, 1999.

[21] E. P. Simoncelli and E. H. Adelson, "Noise removal via Bayesian wavelet coring," in Proceedings ofthe IEEE International Conference on Image Processing, pp. 379-382, Lausanne, Switzerland, September 1996.

[22] E. P. Simoncelli, "Bayesian denoising of visual images in the wavelet domain," in Bayesian Inference in Wavelet Based Models, P. Muller and B. Vidakovic, Eds., pp. 291-308, Springer, New York, NY, USA, 1999.

[23] A. Achim, A. Bezerianos, and P. Tsakalides, "Novel Bayesian multiscale method for speckle removal in medical ultrasound images," IEEE Transactions on Medical Imaging, vol. 20, no. 8, pp. 772-783, 2001.

[24] J. M. Fadili and L. Boubchir, "Analytical form for a Bayesian wavelet estimator of images using the Bessel K form densities," IEEE Transactions on Image Processing, vol. 14, no. 2, pp. 231240, 2005.

[25] P. A. Khazron and I. W. Selesnick, "Bayesian estimation of bessel K form random vectors in AWGN," IEEE Signal Processing Letters, vol. 15, pp. 261-264, 2008.

[26] M. S. Crouse, R. D. Nowak, and R. G. Baraniuk, "Wavelet-based statistical signal processing using hidden Markov models," IEEE Transactions on Signal Processing, vol. 46, no. 4, pp. 886-902, 1998.

[27] J. Portilla, V. Strela, M. J. Wainwright, and E. P. Simoncelli, "Image denoising using scale mixtures of Gaussians in the wavelet domain," IEEE Transactions on Image Processing, vol. 12, no. 11, pp. 1338-1351, 2003.

[28] H. Rabbani and S. Gazor, "Image denoising employing local mixture models in sparse domains," IET Image Processing, vol. 4, no. 5, pp. 413-428, 2010.

[29] A. Elmzoughi, A. Benazza-Benyahia, and J.-C. Pesquet, "An int-erscale multivariate statistical model for map multicomponent image denoising in the wavelet transform domain," in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '05), pp. II45-II48, March 2005.

[30] H. Rabbani and M. Vafadust, "Image/video denoising based on a mixture of Laplace distributions with local parameters in multidimensional complex wavelet domain," Signal Processing, vol. 88, no. 1, pp. 158-173, 2008.

[31] H. A. Chipman, E. D. Kolaczyk, and R. E. McCulloch, "Adaptive Bayesian wavelet shrinkage," Journal of the American Statistical Association, vol. 92, no. 440, pp. 1413-1421,1997.

[32] M. J. Borran and R. D. Nowak, "Wavelet-based denoising using hidden Markov models," in Proceedings ofthe IEEE Interntional Conference on Acoustics, Speech, and Signal Processing,pp. 39253928, May 2001.

[33] J. K. Romberg, H. Choi, and R. G. Baraniuk, "Bayesian tree-structured image modeling using wavelet-domain hidden Markov models," IEEE Transactions on Image Processing, vol. 10, no. 7, pp. 1056-1068, 2001.

[34] A. Pizurica, W. Philips, I. Lemahieu, and M. Acheroy, "A joint inter- and intrascale statistical model for Bayesian wavelet based image denoising," IEEE Transactions on Image Processing, vol. 11, no. 5, pp. 545-557, 2002.

[35] A. A. Bharath and J. Ng, "A steerable complex wavelet construction and its application to image denoising," IEEE Transactions on Image Processing, vol. 14, no. 7, pp. 948-959, 2005.

[36] L. Çendur and I. W. Selesnick, "Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency," IEEE Transactions on Signal Processing, vol. 50, no. 11, pp. 27442756, 2002.

[37] H. Rabbani, M. Vafadust, S. Gazor, and I. Selesnick, "Image den-oising employing a bivariate cauchy distribution with local variance in complex wavelet domain," in Proceedings of the 12th IEEE Digital Signal Processing Workshop and 4th IEEE Signal Processing Education Workshop, pp. 203-208, Grand Teton National Park, Wy, USA, September 2006.

[38] J. Portilla, "Full blind denoising through noise covariance estimation using gaussian scale mixtures in the wavelet domain," in Proceedings of the International Conference on Image Processing (ICIP '04), pp. 1217-1220, October 2004.

[39] S. Lyu and E. P. Simoncelli, "Modeling multiscale subbands of photographic images with fields of Gaussian scale mixtures," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 31, no. 4, pp. 693-706, 2009.

[40] H. Rabbani, R. Nezafat, and S. Gazor, "Wavelet-domain medical image denoising using bivariate laplacian mixture model," IEEE Transactions on Biomedical Engineering, vol. 56, no. 12, pp. 2826-2837, 2009.

[41] I. Selesnick and K. Li, "Video denoising using 2d and 3d dual-tree complex wavelet transforms," in Wavelet Applications in Signal and Image Processing, vol. 5207 of Proceedings of SPIE, August 2003.

[42] Z. Cai, T. H. Cheng, C. Lu, and K. R. Subramanian, "Efficient wavelet-based image denoising algorithm," Electronics Letters, vol. 37, no. 11, pp. 683-685, 2001.

[43] M. K. Mihçak, I. Kozintsev, K. Ramchandran, and P. Moulin, "Low-complexity image denoising based on statistical modeling of wavelet coefficients," IEEE Signal Processing Letters, vol. 6, no. 12, pp. 300-303, 1999.

[44] H. Rabbani, M. Vafadust, P. Abolmaesumi, and S. Gazor, "Speckle noise reduction of medical ultrasound images in complex wavelet domain using mixture priors," IEEE Transactions on Biomedical Engineering, vol. 55, no. 9, pp. 2152-2160, 2008.

[45] L. Çendur and I. W. Selesnick, "Bivariate shrinkage with local variance estimation," IEEE Signal Processing Letters, vol. 9, no. 12, pp. 438-441, 2002.

[46] H. Rabbani, M. Vafadust, and S. Gazor, "Image denoising based on a mixture of bivariate gaussian distributions with local parameters in complex wavelet domain," in Proceedings of the International Conference on Biomedical and Pharmaceutical Engineering, pp. 174-179, December 2006.

[47] V. Katkovnik, K. Egiazarian, and J. Astola, "Adaptive window size image de-noising based on intersection of confidence intervals (ICI) rule," Journal of Mathematical Imaging and Vision, vol. 16, no. 3, pp. 223-235, 2002.

[48] I. W. Selesnick, R. G. Baraniuk, and N. G. Kingsbury, "The dual-tree complex wavelet transform," IEEE Signal Processing Magazine, vol. 22, no. 6, pp. 123-151, 2005.

[49] B. Wang, Y. Wang, I. Selesnick, and A. Vetro, "Video coding using 3D dual-tree wavelet transform," Eurasip Journal on Image and Video Processing, vol. 2007, Article ID 42761,15 pages, 2007.

[50] W. Jingwei, G. Xinbo, and Z. Juanjuan, "A video watermarking based on 3-D complex wavelet," in Proceedings of the 14th IEEE International Conference on Image Processing, pp. V493-V496, September 2007.

[51] J. Yang, Y. Wang, W. Xu, and Q. Dai, "Image and video denoising using adaptive dual-tree discrete wavelet packets," IEEE Transactions on Circuits and Systems for Video Technology, vol. 19, no. 5, pp. 642-655, 2009.

[52] I. Bayram and I. W. Selesnick, "A dual-tree rational-dilation complex wavelet transform," IEEE Transactions on Signal Processing, vol. 59, no. 12, pp. 6251-6256, 2011.

[53] K. Lee, M. Niemeijer, M. K. Garvin, Y. H. Kwon, M. Sonka, and M. D. Abramoff, "Segmentation of the optic disc in 3-D OCT scans of the optic nerve head," IEEE Transactions on Medical Imaging, vol. 29, no. 1, pp. 159-168, 2010.

[54] W. Y. Lo and I. W. Selesnick, "Wavelet-domain soft-thresholding for non-stationary noise," in Proceedings of the IEEE International Conference on Image Processing, pp. 1441-1444, October 2006.

[55] R. M. Willett and R. D. Nowak, "Fast multiresolution photon-limited image reconstruction," in Proceedings of the 2nd IEEE International Symposium on Biomedical Imaging: Macro to Nano (ISBI '04), vol. 2, pp. 1192-1195, April 2004.

[56] A. Pizurica, W. Philips, I. Lemahieu, and M. Acheroy, "A versatile wavelet domain noise filtration technique for medical imaging," IEEE Transactions on Medical Imaging, vol. 22, no. 3, pp. 323-331, 2003.

[57] S. Gupta, R. C. Chauhan, and S. C. Saxena, "Robust non-homo-morphic approach for speckle reduction in medical ultrasound images," Medical and Biological Engineering and Computing, vol. 43, no. 2, pp. 189-195, 2005.

[58] S. Gupta, L. Kaur, R. C. Chauhan, and S. C. Saxena, "A versatile technique for visual enhancement of medical ultrasound images," Digital Signal Processing, vol. 17, no. 3, pp. 542-560, 2007.

[59] S. Yan, J. Yuan, M. Liu, and C. Hou, "Speckle noise reduction of ultrasound images based on an undecimated wavelet packet transform domain nonhomomorphic filtering," in Proceedings of the 2nd International Conference on Biomedical Engineering and Informatics, pp. 1-5, October 2009.

[60] R. F. Wagner, M. F. Insana, and D. G. Brown, "Statistical properties of radio-frequency and envelope-detected signals with applications to medical ultrasound," Journal of the Optical Society of America A, vol. 4, no. 5, pp. 910-922,1987.

[61] P. M. Shankar, "Ultrasonic tissue characterization using a generalized Nakagami model," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 48, no. 6, pp. 17161720, 2001.

[62] P. M. Shankar, J. M. Reid, H. Ortega, C. W. Piccoli, and B. B. Goldberg, "Use of non-Rayleigh statistics for the identification of tumors in ultrasonic B-scans of the breast," IEEE Transactions on Medical Imaging, vol. 12, no. 4, pp. 687-692,1993.

[63] P. M. Shankar, "A model for ultrasonic scattering from tissues based on the K distribution," Physics in Medicine and Biology, vol. 40, no. 10, pp. 1633-1649, 1995.

[64] I. K. Eom and Y. S. Kim, "Wavelet-based denoising with nearly arbitrarily shaped windows," IEEE Signal Processing Letters, vol. 11, no. 12, pp. 937-940, 2004.

[65] W. Drexler and J. G. Fujimoto, "State-of-the-art retinal optical coherence tomography," Progress in Retinal and Eye Research, vol. 27, no. 1, pp. 45-88, 2008.

[66] M. D. Abramoff, M. K. Garvin, and M. Sonka, "Retinal imaging and image analysis," IEEE Reviews in Biomedical Engineering, vol. 3, pp. 169-208, 2010.

[67] C. Ercole, A. Foi, V. Katkovnik, and K. Egiazarian, "Spatiotemporal pointwise adaptive denoising of video: 3D non-parametric approach," in Proceedings of the 1st International Workshop on Video Process and Quality Metrics for Consumer Electronics (VPQM '05), 2005.

[68] Y. M. Lu and M. N. Do, "Multidimensional directional filter banks and surfacelets," IEEE Transactions on Image Processing, vol. 16, no. 4, pp. 918-931, 2007.

[69] L. Ying, L. Demanet, and E. Candes, "3D discrete curvelet transform," in Wavelets XI, vol. 5914 of Proceedings of SPIE, pp. 1-11, August 2005.

Copyright of International Journal of Biomedical Imaging is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.