1.

Introduction

Optical coherence tomography (OCT) is a recently created imaging technique to illustrate different information about the internal structures of an object. As in ultrasound imaging, OCT is based on interferometric techniques, except that OCT uses light beams instead of sound profiles.¹ Use of OCT images is more popular in the field of ophthalmology, where it has become a key diagnostic technology in the areas of retinal diseases.² OCT images usually suffer from speckle noise, which restricts the ability of image interpretation. Hence, noise reduction techniques are an essential process for these images. In a wide division, OCT denoising methods can be classified into two basic groups: methods that reduce the noise by some modifications on OCT device imaging such as Refs. 3 and 4 (based on hardware) and the methods that work on the recorded data (based on software). Because of our limitation in access to hardware of devices, we focus more on the second group of denoising methods. Based on our categorization in Ref. 5, these methods can be applied on the spatial or transform domain. For reducing noise in OCT images in the spatial domain, usually some simple methods like low-pass filters, linear smoothing, mean, median, and Wiener filters are used. However, some developed methods like multiple uncorrelated B-scan averaging,^6,7 nonlinear anisotropic filter,^8,9 complex diffusion,^10,11 directional filtering,^12,13 adaptive vector-valued kernel function,¹⁴ support vector machine approach,¹⁵ and Bayesian estimations¹⁶ are also applied in the spatial domain.

In the transform based group, some data adaptive methods like principle component analysis;¹⁷ dictionary learning;^18,19 some nondata adaptive methods that are based on wavelets,^20,21 dual-tree complex wavelet transform^22,23 or curvelets;^24,25 and a circular symmetric Laplacian mixture model in wavelet diffusion²⁶ have been applied for denoising.

Recently, some groups of investigators²⁷ improved the result of denoising by combining dictionary learning and wavelet thresholding and by defining a proper data-dependent initial dictionary.

In another point of view, the denoising problem may be solved in a deterministic or statistical framework.²³ In the first framework, each pixel is considered an unknown deterministic variable, and non-Bayesian techniques are used for denoising. In the statistical framework, a random field is employed to model the data, and Bayesian methods are applied to estimate noise-free data. Thus, finding a suitable prior probability distribution function (pdf) for noise-free data plays an important role in the denoising problem.²³ Although there is a large body of work on OCT image analysis, few studies have been done on OCT images from a statistical modeling point of view. One of these studies is the work of Grzywacz et al.²⁸ in which a probability model for OCT data was suggested and used for early detection of diabetic retinopathy. The denoising method that is proposed in this paper is the second type (statistical framework), so it first needs to obtain an appropriate statistical model for OCT images.

Another issue is that assumption of the Gaussianity plays a key role in many signal processing methods. For example, a Wiener filter is a linear minimum-mean-square error filter in which the observed and noise-free data are jointly Gaussian. It is proven that, for a Gaussian process, nonlinear filters do not lead to better results compared with a Wiener filter. However, sometimes researchers use this Wiener filtering method without considering this Gaussianity assumption, so they cannot fully benefit from the ability of the algorithm. In this paper, it is illustrated that holding this assumption has a notable effect on the denoising results. Therefore, in this work, a Gaussianization transform (GT) that maps the original distribution of the OCT data to Gaussian components is introduced, and, with this transform, all of the good characteristics of the Gaussian process can be used.

In the rest of the paper, a powerful patch-based denoising method, namely the spatially constrained Gaussian mixture model (SC-GMM), is first explained in Sec. 2.1. Since SC-GMM is based on the Gaussian distribution assumption of noise-free data, which is not satisfied for OCT images, in Sec. 2.2 GT is applied on OCT data. By applying GT, the distribution of transformed OCT data is expressed by a combination of Gaussian pdfs. In this base, the proposed denoising method, which is a combination of Gaussianization and SC-GMM, is demonstrated in Sec. 2.3. The visual and numerical results of the proposed OCT denoising method are presented in Sec. 3, and Sec. 4 presents the discussion of our results in comparison with the state of the art methods. Finally, a summary of this paper is presented in Sec. 5.

2.

Method

Today, patch-based image denoising methods are known as an effective way for denoising natural images.^29–31 The proposed denoising method in this paper is also a patch-based technique; it is a modified version of the SC-GMM algorithm introduced by Niknejad et al.³² In Sec. 2.1, the SC-GMM method that is based on a Gaussian distribution assumption of similar image patches in a neighborhood is described. However, it is clear that this assumption is not satisfied for all kinds of datasets like OCT images. Thus, after investigation of the distribution of the original OCT image, it is necessary to convert it to a Gaussian distribution for which GT is used. These steps are explained in Sec. 2.2 and finally the proposed OCT denoising method, which is a combination of GT and SC-GMM, is introduced in Sec. 2.3.

2.1.

SC-GMM Denoising Method

As mentioned above, SC-GMM denoising is a patch-based method in which the whole of the image is partitioned into overlapping patches and it is assumed that each set of neighbor patches can be modeled using a multivariate Gaussian pdf with a specific mean vector and covariance matrix. For this purpose, some exemplar patches are uniformly chosen from the whole image; then, similar patches in the vicinity of each exemplar patch are grouped together to construct a region using the $K$-nearest neighbors (KNN) algorithm. To determine these clusters of patches and estimate the parameters of Gaussian pdfs (covariance matrices and mean vectors), a two-step approach is applied.³³ This estimation approach includes a clustering step and a denoising step.

Assume $y=x+n$ is a noisy image, where $x$ is noise-free data and $n$ is an independent zero-mean Gaussian noise with known variance ${\sigma }_n^2$. Also, let $x_r$ denote the $r$’th exemplar patch with mean vector ${\mu }_r$ and covariance matrix ${\mathrm{\Sigma }}_r$ for $r=1,\dots ,R$. After initialization, $\hat{x}=y$, the following steps are iteratively implemented:

Clustering step. In this step, for each exemplar patch ${\hat{x}}_r$, KNN patches ${\hat{x}}_i$, $i=1,\dots ,K$, that have minimum dissimilarity (based on $l_2$ norm metric, i.e., $d={\Vert {\hat{x}}_r-{\hat{x}}_i\Vert }_2^2$) with the exemplar patch are chosen. These ${\hat{x}}_r$ and ${\hat{x}}_i\mathrm{s}$ are estimated in the previous iteration or the initialization in the first step.

Denoising step. The goal of this step is to restore the image using the clusters obtained in the clustering step. Thus, first the parameters of Gaussian distribution for each set of patches are estimated using the maximum likelihood method

Eq. (1)

$${\hat{\mu }}_r=\frac{1}{K}\sum _{i\in S_r}{\hat{x}}_i,$$

Eq. (2)

$${\hat{\mathrm{\Sigma }}}_r=\frac{1}{K}\sum _{i\in S_r}({\hat{x}}_i-{\hat{\mu }}_r){({\hat{x}}_i-{\hat{\mu }}_r)}^{\mathrm{T}},$$

where $S_r$ shows the set of KNN patches in the $r$’th region. It should be noted that in our implementation the above sample covariance matrix is not invertible, so as with similar works,³⁴ an eigenvalue regularization is applied using

Eq. (3)

$${\hat{\mathrm{\Sigma }}}_r={\hat{\mathrm{\Sigma }}}_r+\delta I,$$

where $\delta$ is a small constant (here $\delta =0.1$) and $I$ is the identity matrix.

After estimating the parameters of the Gaussian pdf, the denoised patch ${\hat{x}}_i$ in the $r$’th region is found by maximizing a posteriori function $\log p(x|y_i,{\hat{\mu }}_r,{\hat{\mathrm{\Sigma }}}_r)$, i.e.,

Eq. (4)

$${\hat{x}}_i=\underset{x}{\arg \max }\log p(x|y_i,{\hat{\mu }}_r,{\hat{\mathrm{\Sigma }}}_r)=\arg \underset{x}{\min }{\Vert y_i-x\Vert }_2^2+{\sigma }_n^2{(x-{\hat{\mu }}_r)}^{\mathrm{T}}{\hat{\mathrm{\Sigma }}}_r^{-1}(x-{\hat{\mu }}_r),$$

where $y_i$ is the corresponding noisy observed patch and Eq. (4) is derived from the assumption of the Gaussian pdf for noise-free patch and noise, i.e., $x_i\sim \mathcal{N}({\hat{\mu }}_r,{\hat{\mathrm{\Sigma }}}_r)$ and $n\sim \mathcal{N}(0,{\sigma }_n^{2}I)$, respectively.

By setting the derivative of Eq. (4) to zero, this optimization problem can be solved, which leads to the linear Wiener filter as follows:

Eq. (5)

$${\hat{x}}_i={(I+{\sigma }_n^2{\hat{\mathrm{\Sigma }}}_r^{-1})}^{-1}(y_i+{\sigma }_n^2{\hat{\mathrm{\Sigma }}}_r^{-1}{\hat{\mu }}_r).$$

The last step for reconstructing the whole image is averaging the weighted restored patches. In this work, based on Ref. 32, these weights are considered proportional to similarity to Gaussian distributions. To be more explicit, patches that are more likely to be generated from the estimated Gaussian distribution of their allocated clusters benefit from higher weights. Thus, the weight of the patch ${\hat{x}}_i$ in the $r$’th region is obtained by

Eq. (6)

$$w(i,r)=\exp [-\frac{\gamma }{2}{({\hat{x}}_{\mathrm{i}}-{\hat{\mu }}_r)}^T{\hat{\mathrm{\Sigma }}}_r^{-1}({\hat{x}}_{\mathrm{i}}-{\hat{\mu }}_r)],$$

where $\gamma$ is an appropriate scaling constant.

Using these weights and averaging on all patches, the denoised image is found. Figure 1 shows the block diagram of the SC-GMM method for image denoising.

Fig. 1

Block diagram of spatially constrained GMM (SC-GMM) method for image denoising.

2.2.

Gaussianization as a Preprocessing Step

In fact, the assumption of Gaussianity plays an important role in many signal processing methods. For example, in this work, the Gaussian distribution assumption for a dataset causes the derivation of the Wiener filter as an optimal linear filter for denoising. Thus, the best denoising result is reached when the initial assumption of Gaussianity is satisfied. However in Ref. 32 and other related works,³⁴ holding this Gaussianity assumption for the pdf of data is not investigated, and here we demonstrate the importance of this assumption by converting the pdf of OCT intraretinal layers to Gaussian distribution and comparing denoising results with and without Gaussianization. For this purpose, the distribution of the OCT images should first be determined, and then an appropriate GT is designed to convert the OCT original distribution to a GMM and satisfy the initial assumption for our denoising algorithm. The block diagram of the Gaussianization process is displayed in Fig. 2.

Fig. 2

Block diagram of GT.

2.2.1.

OCT distribution estimation

Since the speckle noise in the OCT images is a multiplicative noise and working with the additive noise is simpler and more common, first, a logarithm operator is used to convert this noise to an additive Gaussian noise. Thus, all the processes are implemented in the logarithmic domain. In Ref. 35, we discussed that, because of the layered structure of the retina, a mixture model can be well fitted to a whole OCT B-scan.

Moreover, a monotonically decaying behavior of the OCT intensities in each layer means that leptokurtic distributions such as the Laplace pdf are good candidates for describing the statistical properties of layers.³⁵ On the other hand, this image is corrupted with the additive white Gaussian noise in the logarithmic space, so an appropriate statistical model for an OCT image is a mixture distribution. Each component of this mixture model is composed of a Laplace random variable $x_j$, plus a Gaussian noise $n$. We called this model a normal-Laplace mixture (NLM) model and described it in detail in Ref. 35.

Indeed, for each component of the mixture model, we would have $y_j=x_j+n$, where the pdfs of $x_j$ and $n$ are

Eq. (7)

$$f_{X_j}(x)=\frac{1}{\sqrt{2}{\sigma }_j}\exp (-\frac{\sqrt{2}}{{\sigma }_j}|x-{\mu }_j|),$$

Eq. (8)

$$f_N(n)=\frac{1}{\sqrt{2\pi }{\sigma }_n}\exp (-\frac{n^2}{2{\sigma }_n^2}).$$

${\mu }_j$ shows the mean of the Laplace variable and ${\sigma }_j$ and ${\sigma }_n$ display the scaling parameters of the Laplace and Gaussian distributions, respectively.

Since $y_j=x_j+n$ and the noise signal $n$ is assumed to be independent from the noise-free signal $x_j$, the pdf of the normal-Laplace variable $y_j$ is given by the convolution

Eq. (9)

$$f_j(y)=(f_{X_j}*f_N)(y)=\underset{-\infty }{\overset{+\infty }{\int }}{f}_{X_j}(x)f_N(y-x)\mathrm{d}x.$$

After some mathematical computations, we have³⁵

Eq. (10)

$$f_j(y)=\frac{-1}{2\sqrt{2}{\sigma }_j}\left\{\exp \right[\frac{{\sigma }_n^2}{{\sigma }_j^2}+\frac{\sqrt{2}(y-{\mu }_j)}{{\sigma }_j}\left]\mathrm{erfc}\right[-\frac{(y-{\mu }_j)}{\sqrt{2}{\sigma }_n}-\frac{{\sigma }_n}{{\sigma }_j}]-\exp [\frac{{\sigma }_n^2}{{\sigma }_j^2}-\frac{\sqrt{2}(y-{\mu }_j)}{{\sigma }_j}\left]\mathrm{erfc}\right[-\frac{(y-{\mu }_j)}{\sqrt{2}{\sigma }_n}+\frac{{\sigma }_n}{{\sigma }_j}\left]\right\},$$

where $\mathrm{erfc}(x)=e^{-x^2}\mathrm{erfc}x(x)=\frac{2}{\sqrt{\pi }}\underset{0}{\overset{\infty }{\int }}{\mathrm{e}}^{-x^2-t^2-2tx}\mathrm{d}t$.

Thus, $y_j$ has a normal-Laplace distribution, and we show it by $y_j\sim \mathrm{NL}({\mu }_j,{\sigma }_j,{\sigma }_n)$. Additionally, the pdf of the NLM model is presented by $g_Y(y)$ as follows:

Eq. (11)

$$g_Y(y)=\sum _{j=1}^J{a}_j{f}_j(y),$$

where $f_j(y)$ is the pdf of each component, $a_j\ge 0$, $\sum _{j=1}^J{a}_j=1$, and $J$ is the number of mixture components. In Ref. 35, it has been shown that the choice of $J=5$ causes the best results.

After finding a proper model, its parameters are estimated based on a used dataset. The EM algorithm³⁶ is used to estimate parameters iteratively as follows:

• • E-step: In each iteration, the responsibility factors $r_j$ (as an auxiliary variable that for each observed data $y$ represents how likely that observed data was produced by each component ${\{f_j(y)\}}_{j=1:\mathrm{J}}$) are updated by

  Eq. (12)

  $$r_j\leftarrow \frac{a_j{f}_j(y)}{{\sum }_{j=1}^J{a}_j{f}_j(y)},\text{for}j=1,\dots ,J.$$

• • M-step: Then, the parameters $a_j$ are updated by

  Eq. (13)

  $$a_j\leftarrow \sum _{l=1}^N{r}_j(l)/N,$$
  where $N$ shows the number of pixels in image $y$.

In the M-step, the updated equations for ${\{{\mu }_j,{\sigma }_j\}}_{j=1}^J$ are approximately estimated as³⁷

Eq. (14)

$${\mu }_j\leftarrow \sum _{l=1}^{N}y.r_j(l)/\sum _{l=1}^N{r}_j(l),$$

Eq. (15)

$${\sigma }_j^2\leftarrow \max \{\sum _{l=1}^N{|y-{\mu }_j|}^2{r}_j(l)/\sum _{l=1}^N{r}_j(l),\delta \},$$

where $\delta$ is a small positive constant number used to avoid numerical errors.

2.2.2.

Gaussianization transform

As mentioned before, SC-GMM has been designed based on Gaussianity assumption. However, the distribution of each OCT intraretinal layer can be modeled as a normal-Laplace pdf.³⁵ Thus, if a GT is designed to map the original distribution to a Gaussian distribution, it is expected to obtain better denoising results. The core of this transform is solving the equation that compares the cumulative distribution function (CDF) of original dataset and Gaussian distribution. Assume $F_W(w)=0.5[1+\mathrm{erf}(\frac{w-{\mu }_w}{\sqrt{2}{\sigma }_w})]$ expresses the CDF of a Gaussian variable ($w$) with ${\mu }_w$ and ${\sigma }_w$ as the mean and standard deviation of $w$, respectively, and $\mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}\underset{0}{\overset{x}{\int }}{e}^{-t^2}\mathrm{d}t$. Also, by computing the integral of $f_j(y)$ in Eq. (10), the CDF of a normal-Laplace variable is found as

Eq. (16)

$$F_j(y)=\frac{1}{2}\{1+\mathrm{erf}[\frac{(y-{\mu }_j)}{\sqrt{2}{\sigma }_n}]-\exp [\frac{{\sigma }_n^2}{{\sigma }_j^2}+\frac{-\sqrt{2}(y-{\mu }_j)}{{\sigma }_i}]+\frac{1}{2}{e}^{-{\left[\frac{(y-{\mu }_j)}{\sqrt{2}{\sigma }_n}\right]}^2}\{\mathrm{erfc}x[\frac{(y-{\mu }_j)}{\sqrt{2}{\sigma }_n}+\frac{{\sigma }_n}{{\sigma }_j}]+\mathrm{erfc}x[\frac{(y-{\mu }_j)}{\sqrt{2}{\sigma }_n}-\frac{{\sigma }_n}{{\sigma }_j}]\left\}\right\}.$$

By equating these two CDFs, one can find an equation that expresses each value of variable $y$ in the initial distribution, in terms of the equivalent value of $w$ in Gaussian space

Eq. (17)

$$F_Y(y)=F_W(w)=>w={\mu }_w+\sqrt{2}{\sigma }_w.\mathrm{erfinv}[2F_Y(y)-1],$$

where erfinv is the inverse of the error function.

2.2.3.

Construct a GMM distribution

After Gaussianization, the distribution of each component of the NLM model, $w_j$, $j=1,\dots ,J$, becomes Gaussian; by combining these Gaussian components, an image with the GMM pdf is constructed and can be used as the input of the SC-GMM denoising method. In Ref. 38, it is shown that the averaged maximum a posterior (AMAP) method is a suitable way for defining the weights of summation for mixture models. In this method, the weight of each component is defined as $a_j{f}_j(y)/\sum _{j=1}^J{a}_j{f}_j(y)$, so $\hat{y}$, as the Gaussianized image with GMM distribution, is obtained as follows:

Eq. (18)

$$\hat{y}=\sum _{j=1}^J{a}_j.w_j.f_j(y)/\sum _{j=1}^J{a}_j{f}_j(y).$$

Also, in Ref. 39, it was discussed that instead of using the univariate weighting function $a_j{f}_j(y)/\sum _{j=1}^J{a}_j{f}_j(y)$, employing multivariate weighting function leads to better results. Thus, for an OCT image of size $S\times T$, a local $P\times Q$ window (where $P$ and $Q$ are odd) is placed around each pixel $y_{s,t}$ for $s=1,\dots ,S$, $t=1,\dots ,T$; and the multivariate function ${\overline{f}}_j(y)$ for the central pixel of window ($y_{p,q}$) is calculated by multiplying all univariate $f_j(y)s$ in that window i.e.

Eq. (19)

$${\overline{f}}_j(y_{p,q})=\prod _{k=-\frac{P-1}{2}}^{\frac{P-1}{2}}\prod _{l=-\frac{Q-1}{2}}^{\frac{Q-1}{2}}{f}_j(y_{p+k,q+l}).$$

Figure 3 shows how the multivariate function ${\overline{f}}_j(y)$ is produced.

Fig. 3

Using local window around pixels to make a multivariate weighting function ${\overline{f}}_j(y)$.

Hence, Eq. (18) is replaced with

Eq. (20)

$$\hat{y}=\frac{{\sum }_{j=1}^J{a}_j.w_j.{\overline{f}}_j(y)}{{\sum }_{j=1}^J{a}_j{\overline{f}}_j(y)}.$$

2.3.

Proposed Denoising Method

According to Secs. 2.2.1 and 2.2.2, our proposed OCT denoising method is obtained by combining GT and the SC-GMM method as shown in Fig. 4. In this method, the OCT image in the log domain is first Gaussianized using the block diagram of Fig. 2, and then the Gaussianized image is denoised by the SC-GMM method using the block diagram of Fig. 1.

Fig. 4

Block diagram of the proposed combined denoising method, in this figure, third block “GT” includes all the blocks in Fig. 2 (except the Log block) and fourth block “SC-GMM denoising method” is equivalent to the block diagram of Fig. 1.

The proposed denoising algorithm can be summarized in the following steps:

• 1. Transform image to the logarithmic domain

• 2. Fit an NLM model to the data using the EM algorithm:

  • 2.1 Initialization: choose initial values for $a_j$, ${\mu }_j$, ${\sigma }_j$

  • 2.2 E-step: compute responsibility factors using Eq. (12)

  • 2.3 M-step: update parameters $a_j$, ${\mu }_j$, ${\sigma }_j$ using Eqs. (13)–(15)

  • 2.4 Iteration: substitute the updated parameters in the previous step for calculating $f_j(y)$

• 3. Calculate the CDF of each normal-Laplace component by substituting final values of ${\mu }_j$, ${\sigma }_j$, and ${\sigma }_n$ in Eq. (16)

• 4. Gaussianize each normal-Laplace component using Eq. (17)

• 5. Use AMAP method in Eq. (20) to obtain the image $\hat{y}$ with GMM distribution

• 6. Apply SC-GMM algorithm on image $\hat{y}$:

  • 6.1 Partition image into overlapping patches and choose exemplar patches

  • 6.2 For each exemplar patch:

    • 6.2.1 Select KNN patches (a region) based on $l_2$ norm metric in a finite sized window around exemplar patch

    • 6.2.2. Estimate Gaussian parameters in each region by Eqs. (1)–(3)

    • 6.2.3. Denoise the patches in each region by the Wiener filter Eq. (5)

  • 6.3. Obtain reconstructed image by weighted average of denoised patches using Eq. (6)

• 7. Apply exponential operator to find the denoised image.

3.

Results

In this section, some simulation results performed on our dataset are presented. This dataset consists of thirteen three-dimensional (3-D) macular spectral domain OCT images obtained from eyes without pathologies using a Topcon 3DOCT-1000 imaging system in the Ophthalmology Department of Feiz Hospital, Isfahan, Iran. The $x$, $y$, $z$ size of the obtained volumes was $512\times 650\times 128\text{voxels}$, $7\times 3.125\times 3.125{\mathrm{mm}}^3$, voxel size $13.67\times 4.81\times 24.41\mu {\mathrm{m}}^3$.

From this dataset, 130 B-scans were selected randomly (10 B-scans from each subject), and these randomly selected B-scans were denoised with the methods described previously in the paper: once images were denoised only using the SC-GMM method and they were then despeckled by a combination of GT and the SC-GMM according to our method.

Considering the constant values and parameters, for each region, $K=40$ nearest neighbor patches were accumulated. The exemplar patches were selected every 10 pixels along both the column and row directions of the image. The size of the constraining window around each exemplar patch was set to $30\times 30$, and the size of patches was selected as $10\times 10$. Also, the size of local window around each pixel for obtaining multivariate ${\overline{f}}_j(y)$ was $3\times 3$. All of these values were selected based on previous related work and our implementation experiences.

Both the qualitative and quantitative performance of the proposed combined method and SC-GMM were provided. The results of denoising with a simple Wiener filter was also prepared to compare with the proposed approach and prove the need and benefit of GT in the proposed method. Contrast to noise ratio (CNR), edge preservation (EP), and texture preservation (TP) as quantitative measures were computed for each denoising method. The mentioned values were computed based on methods discussed in Ref. 40, and regions for analysis of each criterion are shown in Fig. 5.

Fig. 5

Selected ROIs for calculation of CNR, EP, and TP. The larger elliptical region outside the retinal layers is used as background ROI and other circles represent foreground ROIs for CNR calculation. Rectangular ROIs are used for TP and EP.

CNR measures the contrast between a feature of interest and background noise. EP computes the correlation between edges in the denoised image and the noisy image over the locally selected region of interest (ROI). TP shows the degree of preserving texture in an ROI. Both the TP and EP measures range between 0 and 1 and, in the best condition, go to 1.

In addition, to illustrate the power of the proposed method, its ability in OCT denoising was compared with the BM3D as a benchmark denoising algorithm and recently reported state-of-the-art algorithm in this field.²⁷ Kafieh et al.²⁷ proposed a strong method in speckle reduction of OCT datasets in two-dimensional and 3-D by the combination of dictionary learning and wavelet thresholding. They improved the performance of simple dual-tree complex wavelets by taking advantage of adaptability in dictionary learning methods and introducing complex wavelet-based K-SVD. They also compared their algorithm by other previous outstanding researches in this filed like Fang et al.,¹⁸ Yang et al.,⁴¹, and Luan and Wu¹⁷ and demonstrated the superiority of their method. Thus, this complex wavelet-based K-SVD method (3DCWT-KSVD)²⁷ is chosen as a proper reference to evaluate our proposed method. Figures 6 and 7 show some samples of datasets after application of each one of the mentioned denoising methods (simple Wiener, SC-GMM, 3DCWT-KSVD, BM3D, and the proposed method). Also, the value of the quantitative measurements for each denoising method on all randomly selected slices is summarized in Table 1. Furthermore, to have a comparison between time complexity of the evaluated denoising methods, the mean of the needed computation time for each slice on a PC with Microsoft Windows X32 edition, Intel core2 Duo CPU at 3.00 GHz, 4 GB RAM is reported in Table 1.

Fig. 6

Denoising results in a sample B-scan. The second and third columns show an enlargement of the specified regions of first column images. (a), (g), (m) Original noisy image; (b), (h), (n) simple Wiener; (c), (i), (o) SC-GMM; (d), (j), (p) BM3D; (e), (k), (q) 3DCWT-KSVD; (f), (l), (r) proposed method.

Fig. 7

Denoising results in a sample B-scan. The second and third columns show an enlargement of the specified regions of first column images. (a), (g), (m) Original noisy image; (b), (h), (n) simple Wiener; (c), (i), (o) SC-GMM; (d), (j), (p) BM3D; (e), (k), (q) 3DCWT-KSVD; and (f), (l), (r) proposed method.

Table 1

Evaluation measure results; averaged over 130 OCT B-scans from the Topcon device.

Before analyzing the results of Table 1, it should be noted that visual results of a simple Wiener filter (like those shown in Figs. 6 and 7) indicate the weakness of it in OCT image denoising, and the high value of TP and EP for the Wiener filter in Table 1 is because of the definition of these two measures in which the similarity of the denoised image to the original one plays an important role. However, the low value of CNR for the Wiener filter emphasizes the weakness of this denoising method. Regarding the other mentioned denoising algorithms, according to Table 1, the performance of the proposed method is considerably better than the other methods in TP, which is obviously visible in Figs. 6 and 7. Performance of all of the studied methods is close to each other in EP, and the SC-GMM has the highest EP, which is obtained because of the least difference appearing in the resulting image visually. However, it cannot be considered a positive point lonely, since the values of CNR and TP for this method are very low and indeed this method does not show the proper features as a powerful method for OCT denoising. Regarding the CNR measure, the proposed method and 3DCWT-KSVD have remarkably better results than others. Although 3DCWT-KSVD has the highest CNR among all studied methods, it cannot be the winner in this denoising competition because of its weak results in TP and EP; visual results in Figs. 6 and 7 also confirm this matter.

4.

Discussion

The main objective of this paper is to prove the effect of Gaussianization in the result of the OCT denoising by a Wiener-based denoising method. As our results in Table 1 and Figs. 6 and 7 show, it is obviously clear that using the proposed combined method has been much more successful than the SC-GMM method, which assumes Gaussianity for the noise-free image without any justification for proposing this assumption on OCT datasets.

Our second goal in this work is to generally propose a significant method for OCT denoising. To investigate the power of the proposed method, a comparison with the 3DCWT-KSVD algorithm as a recent state of the art in this field and BM3D as a powerful denoising method was done. Visual results in Figs. 6 and 7 demonstrate better noise reduction using the proposed method while desired textural features are preserved simultaneously. Also, in Table 1, this matter is emphasized by checking the high value of TP in the proposed method compared with the 3DCWT-KSVD and BM3D algorithms. Regarding the CNR measures, 3DCWT-KSVD could provide higher CNR because of the very good speckle removing in the background area, but it does not mean that this method is better for next OCT processing because the main goal behind the OCT denoising is providing an appropriate dataset for the next OCT processing, such as layer segmentation⁴² or lesions identification, while the denoised image obtained by the 3DCWT-KSVD method does not show this ability very well.

To demonstrate this superiority in our proposed method, A-scans of each denoised image were analyzed because many OCT segmentation methods perform based on gradient information from B-scans or analysis of A-scans in OCT datasets.⁴³ The high level of speckle noise in the original noisy OCT image does not allow features of the image to be detectable in segmentation processes; however, applying a good denoising algorithm that suppresses much of the speckle noise and preserves features of the image can provide more reliable segmentation. Figure 8 shows the ability of the proposed method to reduce the speckle noise and protect identifiable feature peaks and valleys. In fact, because of producing the clean and almost noiseless B-scan images in the proposed method [Fig. 8(f)], its A-scans obviously include a peak or valley for each layer of the retina, and the signal amplitude at the boundary layers can be greatly altered. In the other words, an important benefit of the proposed method is using contrast enhancement and denoising together. As depicted in our previous work,³⁵ applying GT on OCT images leads to contrast enhancement of the image and combination of GT and SC-GMM is a proper algorithm for simultaneously denoising and contrast enhancement.

Fig. 8

The cross-section signal along the white solid line. (a) Original noisy image, (b) simple Wiener, (c) SC-GMM, (d) BM3D, (e) 3DCWT-KSVD, (f) proposed method. Red stars on the A-scans show the place of each of 12 boundaries of retinal layers. The name of these layers is ILM, internal limiting membrane; NFL, nerve fiber layer; GCL, ganglion cell layer; IPL, inner plexiform layer; INL, inner nuclear layer; OPL, outer plexiform layer; ONL, outer nuclear layer; IS, inner segment of rods and cones; CL, connecting cilia layer; OS, outer segment of rods and cones, RPE, retinal pigment epithelium, respectively, from top to down. RPE includes two surfaces and three last boundaries in our figures. Green indicators in (f) display some examples of the superiority of the proposed method in reducing the speckle noise and protecting identifiable feature peaks and valleys, compared with 3DCWT-KSVD [orange indicators in (e)].

Another way to evaluate the performance of the proposed method is by analyzing the intraretinal layer segmentation results. For this purpose, a diffusion map-based method introduced in Ref. 43 for OCT boundaries segmentation was used. This segmentation method was tested once on original OCT B-scans and another time on images denoised by each of the aforementioned denoising algorithms.

For each case, the error values were calculated as the difference of computed layer position by a diffusion map segmentation method and manually labeled values by an expert observer.

Tables 2 and 3 show the mean signed and unsigned border positioning errors for each boundary, which is obtained for 130 randomly selected B-scans from the Topcon dataset. The numbering of retinal surfaces in these tables is based on Fig. 8.

Table 2

Summary of mean signed border positioning errors on 12 boundaries from the Topcon dataset.

Table 3

Summary of mean unsigned border positioning errors on 12 boundaries from the Topcon dataset.

Based on Tables 2 and 3, the proposed method seems superior to the other mentioned algorithms in the segmentation of some boundaries (8th to 12th boundaries). However, to illustrate the statistically significant improvement of the proposed method over the other algorithms, Table 4 shows the obtained $p$-values. Based on Table 4, it can be seen that the 8th to 12th boundaries, which are more difficult to identify because of lack of contrast, are significantly better localized in the proposed method compared with the SC-GMM algorithm. Regarding the BM3D and 3DCWT-KSVD algorithms, this superiority is significant just for the 11th and 12th boundaries ($p$-value $<0.05$).

Table 4

Improvement of the proposed method compared with SC-GMM, BM3D, and 3DCWT-KSVD algorithm, on the 8th to 12th boundaries of retinal OCT.

Furthermore, given that the ultimate goal of all the works in the medical image processing field is to help physicians determine better and more accurate diagnosis and treatment of disease, ranking mentioned denoising methods by an expert was considered another acceptable criterion for evaluating our method. Accordingly, an ophthalmologist was asked to score the denoised images by each method (raw data, denoised images by SC-GMM, 3DCWT-KSVD, and the proposed method) between one and four; the score of one to four is assigned “best,” “very good,” “good,” and “worst,” respectively. This ranking was done on all 130 randomly selected B-scans (some sample images are illustrated in Figs. 6 and 7) and the average score on all B-scans for our method was 1.2846 (for 110 slices of all 130 B-scans, the proposed method achieved score one), while this value for raw images was 1.8538, for SC-GMM was 2.8231, and for 3DCWT-KSVD was 3.9615. Based on this result, the superiority of the proposed method in clinical applications is also obvious.

Another noticeable issue is that in this research we worked on single-frame images provided by Topcon imaging system; hence the proposed method may become a useful step in OCT image processing approaches, where single frame detection is a must (where the data acquisition time or system restrictions do not allow for multiframe acquisition). However, to make a good judgment on the power of the method in general, it should be compared with results from a multiframe approach (as angular compounding for example), and hopefully this research can be used as a starting point for further and deeper discussions in this field in the OCT community.

5.

Conclusion

In this paper, a developed denoising method for speckle reduction in OCT images was introduced. Our main issue was the effect of holding the assumptions in different statistical methods. Particularly, in this work, the employed denoising algorithm (SC-GMM) was based on the Wiener filter, and in this method both noise-free data and noise are assumed to be Gaussian; however, holding this assumption has not been investigated in previous research. Since OCT datasets do not satisfy this assumption and an NLM model is fitted on these datasets, a transform to Gaussianize each component of this mixture model was needed. After constructing an image underlying GMM, SC-GMM was used for denoising. The visual and numerical results show the superiority of this combined method and encourage us to compare the proposed method with BM3D and 3DCWT-KSVD algorithms as the state-of-the-art studies in this field. Indeed, the advantage of our method, in addition to appropriate noise reduction, is in preserving and highlighting the features of the image like textures and edges that provide a suitable dataset for the next OCT processing such as layer segmentation or object detection. This matter was demonstrated by investigating the behavior of A-scans in all mentioned denoising methods. In addition, the method was shown to improve segmentation accuracy of some intraretinal layers. Furthermore, as another way to prove the superiority of the proposed method in clinical applications, the result of ranking different methods by an ophthalmologist shows predominance of the proposed method.