1. Introduction

Optical coherence tomography (OCT) is an emerging and important imaging technique that provides high-resolution tomographic imaging and cross-sectional microstructure of biological tissues [1]. Micrometer scale resolution and layered structure of eye provided by OCT is handy for the detection of disease. Thus, OCT has a variety of clinical applications and is widely used to assist in the diagnosis of vision related diseases. Due to taking advantage of low-coherence interferometry, OCT images are always corrupted by speckle noise [2], which severely degrades the quality of OCT imaging, consequently makes it particulary challenging to identify the fine features of object for clinical examination. Hence, suppressing speckle noise in OCT images is important to improve image quality and to further increase the accuracy of clinical diagnosis. However, removing such noise while preserving the image fine features is a challenging task.

In the past few decades, great efforts have been taken for reducing speckle noise in OCT images, and a large number of approaches have been reported. Generally speaking, these methods can be mainly categorized as hardware-based approaches and software-based approaches [3]. The hardware-based approaches, such as frequency compounding [4,5]and spatial compounding [6–9], are robust to speckle noise suppression because uncorrelated noise properties vary across different illumination angles and wavelengths [1]. However, it will inevitably have to modify the existing imaging system hardware, which is not easily adapted to standard commercial OCT equipments. An efficient and economical way of speckle noise reduction is software-based approaches, which postprocess OCT images. These methods focus on filtering or numerical algorithms in spatial or a certain transform domain, such as the local statistics based filtering methods [10,11], the diffusion based methods [12–16], Bayesian estimation method [17] and the wavelet based methods [18–21]. The local statistics based filtering methods work fast, but blur detailed structures. The diffusion based methods provide better results with respect to speckle noise reduction, yet over-fitting is still a intrinsically problem because of the nature of diffusion equation that is the loss of information during the diffusion process. Wavelet based filters can obtain good performance on speckle noise reduction, however, they may introduce unwanted visual artifacts that will degrade the OCT images. The aforementioned approaches do not take spatial redundancy of images into account or use a fixed basis to represent OCT images, however, a rich amount of fine structures contained in OCT images cannot be well represented by one fixed wavelet basis. In contrast, self-similarity priori based nonlocal means methods [22–24] and adaptive basis based sparse representation methods [25–27] can obtain better speckle noise reduction. Several recent works have also applied the low-rank representation to speckle noise suppression [28–30]. However, these methods vectorize a local image patch, which breaks apart the local correlation between pixels and destroys the fine structures of objects [31].

The development of OCT imaging system increases the correlation between adjacent frames. Therefore, multi-frame OCT images are gradually used to reduce speckle noise and many multi-frame based methods have been reported. For example, Jian et al. [32] propose a 3D curvelet transform method. The coefficients of 3D curvelet transform are thresholded and inverse transform is used to obtain the final despeckled images. In [19], Mayer et al. first perform wavelet transform to each frame, then threshold the coefficients, and finally reconstruct the despeckled image from the weighted wavelet coefficients. Despite the impressive performance, these approaches lose many fine details of OCT images due to directly truncating the transform coefficients. Recently, low-rank matrix decomposition [33] and weighted nuclear norm minimization [34] are used to denoise multi-frame OCT images. However, each 3D block must be vectorized, consequently destroys the structure of the image cubic blocks.

To address the above-mentioned issues, a multi-linear principal component analysis (MPCA) based method is proposed to despeckle multi-frame OCT images. The proposed method takes advantage of inter-frame and intra-frame redundancy of OCT volume data and directly performs on nonlocal similar cubic blocks without vectorization. To well preserve image local features, nonlocal similar 3D blocks are first extracted by an overlapping way and grouped using k-means++ clustering method, which guarantees that only blocks with similar content are used to learn a nonlocal MPCA basis. The grouping procedure is beneficial to extracting the dominant information from the speckled OCT images and removing the noise. When obtaining the nonlocal MPCA basis for a certain group, MPCA transform is performed on the corresponding group and these transform coefficients are shrunk to remove speckle noise. The filtered cubic blocks are then obtained by inverse MPCA transform. Applying this procedure to each group, the depseckled OCT images are finally obtained by aggregating all filtered cubic blocks since each voxel might belong to several groups. Experimental results show that the proposed method can effectively reduce speckle noise while preserving fine details, and outperforms other popular approaches in terms of objective metrics and visual inspection.

The remainder of the paper is structured as follows: In section 2, we introduce some necessary notions and preliminaries. Section 3 briefly reviews the procedure of PCA denoising. Section 4 presents the proposed speckle noise reduction method in details. In section 5, we apply our method to denoise 3D OCT data, and compare it with several state-of-the-art methods from the perspectives of quantitative evaluation and visual quality. Section 6 concludes the paper.

2. Notions and preliminaries

In this section, we briefly introduce some notions and basic multilinear algebra. Throughout the paper, we denote scalars by non-bold letters, e.g., a; vectors by bold lower case letters, e.g., x; matrices by bold upper case letters, e.g., U; and tensors by calligraphic letters, e.g., 𝒜. For a N^th order tensor 𝒜 ∈ ℝ^{I₁×I₂×···×I_N}, we denote its elements as a_{i1···in ···iN}, where 1 ≤ i_n ≤ I_n. The mode-n vectors of a N^th order tensor 𝒜 are the I_n dimensional vectors obtained from 𝒜 by varying index i_n while keeping the others fixed. The mode-n unfolding matrix of tensor 𝒜 is denoted as A_(n) ∈ ℝ^{I_n×(I₁×···×I_n−1×I_n+1×···×I_N)}, which is composed by taking the mode-n vectors of 𝒜 as its columns.

The mode-n product of a tensor 𝒜 ∈ ℝ^{I₁×I₂×···×I_N} by a matrix U ∈ ℝ^J_n×I_n, denoted by 𝒜 ×_n U, is a N^th tensor 𝒞 ∈ ℝ^{I₁×···×J_n×···×I_N}, whose entries are computed by $c_{i_1\cdots i_{n-1}{j}_n{i}_{n+1}\cdots i_N}={\sum }_{i_n}{a}_{i_1\cdots i_{n-1}{i}_n{i}_{n+1}\cdots i_N}{u}_{j_n{i}_n}$. The mode-n product 𝒞 = 𝒜 ×_n U can also be calculated by the matrix multiplication C_(n) = UA_(n), followed by a re-tensorization of undoing the mode-n flattening. The scalar product of two tensors 𝒜, ℬ ∈ ℝ^{I₁×I₂×···×I_N} is defined as $<𝒜,ℬ>={\sum }_{i_1}{\sum }_{i_2}\cdots {\sum }_{i_N}{a}_{i_1{i}_2\cdots i_N}\cdot b_{i_1{i}_2\cdots i_N}$ and the Frobenius norm of a tensor 𝒜 is defined as: ${\Vert 𝒜\Vert }_F=\sqrt{<𝒜,𝒜>}={\left({\sum }_{i_1,\cdots ,i_N}{\left|a_{i_1\cdots i_N}\right|}^2\right)}^{1/2}$.

Following multilinear algebra, tensor 𝒜 can be expressed as the product

3. PCA denoising

Many researches have shown that PCA based denoising methods [35, 36] can learn adaptive transform basis and obtain impressive achievement. In order to better describe the proposed method, in this section, we briefly review the procedure of PCA denoising method [35].

For a vectorized local image patch y₁, we can find the first n − 1 most similar nonlocal image pathes y_i(i = 2, 3, · · · n) in a search window, and group these vectors into a matrix Y denoted as

The PCA transform can fully de-correlate the centralized dataset Ȳ. In the transform domain, the energy of latent clean signal will concentrate on a small subset, while the energy of noise will uniformly spread over the whole dataset. Therefore, the noise and signal can be better distinguished in the PCA domain [35]. Truncating the coefficients can be used to reconstruct the latent clean signal. Unfortunately, such a procedure leads in practice to poor denoising results compared to other alternatives such as the soft thresholding of the coefficients. The main reason is that all projection directions are relevant for modeling image patches [36]. In [35], Zhang et al. shrink the transform dataset to suppress noise by using linear minimum mean square-error estimation (LMMSE). The LMMSE of the kth row of B is defined as B̂^k = w_kB^k, (k = 1, 2, · · · , m), where $w_k=\mathit{max}\left(\left(\left({\sum }_{i=1}^n{\left(b_i^k\right)}^2/n\right)-{\sigma }^2\right)/\left({\sum }_{i=1}^n{\left(b_i^k\right)}^2/n\right),0\right)$, $b_i^k$ is the (k, i)th entry of matrix B. By transforming B̂ back to the time domain, the denoised result of Ȳ is obtained as $\widehat{\overline{\mathbf{Y}}}=\mathbf{U}\widehat{\mathbf{B}}$. Adding the mean value back to $\widehat{\overline{\mathbf{Y}}}$ gives the denoised dataset Ŷ. The final denoised image can be obtained by applying above procedure to each image patch.

PCA has a wide range of applications in the field of image processing, however, PCA based method can only process vectorized multi-dimensional data, which may break apart the structure information of data. To overcome this problem, many methods have been proposed. For example, Zhang et al. [37] proposed a 2D-PCA method, which calculates the bases in the column direction. Generalized 2D-PCA [38] calculates the bases in both row and column directions, therefore, 2D images can be represented not only accurately but efficiently. For multi-dimensional data, multi-linear principal component analysis (MPCA) [39,40] is proposed to calculate the bases in each mode subspace rather than 1D vector space. MPCA can represent the multi-dimensional data more efficiently.

4. Proposed method

Methods based on PCA transform can obtain impressive performance on image denoising. However, image-to-vector manner of PCA destroys the structure of the image patches. To address this problem, a multi-linear principal component analysis (MPCA) based method is proposed to denoise multi-frame images. The proposed algorithm mainly consists of the following steps: cubic blocks grouping, MPCA transform, shrinking transform dataset, inverse MPCA transform and aggregation. However, noise of OCT images is multiplicative problem. Therefore, a logarithmic transformation has to be performed prior to this procedure to convert the multiplicative speckle noise into additive one. On a logarithmic scale, the observed OCT volume data 𝒴 is modeled as latent clean data 𝒳 and additive noise 𝒩: 𝒴 = 𝒳 + 𝒩.

4.1. Cubic blocks grouping

Image denoising is an ill-posed problem, image priors are thus usually used to regularize this problem. It is well known that many repetitive patterns can be found in image domain and many blocks are similar to a local one. A large number of studies have shown that grouping these similar blocks together is beneficial to extracting the main signal and removing noise. Denote H, W and S as the height, width and number of 2D OCT images, therefore, 3D OCT data can be expressed as 𝒴 ∈ ℝ^H×W×S. By sweeping all across the 3D OCT data with overlaps, we can extract a set of cubic blocks $\mathrm{\Omega }={\left\{{𝒴}_i\in {ℝ}^{h\times w\times S}\right\}}_{i=1}^N$, where N = (H − h + 1) × (W − w + 1) is the total number of extracted cubic blocks over the whole 3D data. Then, k-means++ clustering method [41] is used to group Ω into K clusters denoted by ${\left\{P_k\right\}}_{k=1}^K$, where $P_k={\left\{{𝒴}_i^k\right\}}_{i=1}^{n_k}$ and ${\sum }_{k=1}^K{n}_k=N$, ${𝒴}_i^k$ is the i^th 3D block in the k^th cluster. The aim of the grouping procedure is to obtain adaptive tensor bases so that all blocks in one group can be better approximated.

4.2. MPCA transform

Suppose that all 3D blocks in the given group $P_k={\left\{{𝒴}_i^k\in {ℝ}^{h\times w\times S}\right\}}_{i=1}^{n_k}$ are centered by subtracting the average block ${\overline{𝒴}}^k={\sum }_{i=1}^{n_k}{𝒴}_i^k/n_k$, our aim is to compute U⁽¹⁾, U⁽²⁾, U⁽³⁾, and ${𝒮}_i^k$, such that ${𝒮}_i^k{\times }_1{\mathbf{U}}^{(1)}{\times }_2{\mathbf{U}}^{(2)}{\times }_3{\mathbf{U}}^{(3)}$ gives the best approximation of tensor object ${𝒴}_i^k$. That can be done by solving the following minimization problem

Theorem 2 [39]: If the score tensors ${𝒮}_i^k$ are chosen as ${𝒮}_i^k={𝒴}_i^k{\times }_1{\mathbf{U}}^{{(1)}^T}{\times }_2{\mathbf{U}}^{{(2)}^T}{\times }_3{\mathbf{U}}^{{(3)}^T}$, the minimization of the problem (3) is equivalent to the following maximization problem

Theorem 3 [40]: Let {U⁽ⁿ⁾, n = 1, · · · , N} be the solution to Eq. (4). Then, for given all other projection matrices U⁽¹⁾, · · · , U⁽ⁿ⁻¹⁾, U⁽ⁿ⁺¹⁾, · · · , U^(N), the matrix U⁽ⁿ⁾ consists of the eigenvectors of the matrix

The pseudo-code is given in Algorithm 1. Theoretically, the local optimal solutions of the matrices U⁽ⁱ⁾(i = 1, 2, 3) are obtained because of their dependence on the initialization of U⁽²⁾ and U⁽³⁾. The relative reduction of the mean reconstruction error, i.e. $\xi ={\sum }_{i=1}^{n_k}{\Vert {𝒴}_i^k-{𝒮}_i^k{\times }_1{\mathbf{U}}^{(1)}{\times }_2{\mathbf{U}}^{(2)}{\times }_3{\mathbf{U}}^{(3)}\Vert }_F^2/n_k$, is adopted to check the convergence of the algorithm. When |ξ(j) − ξ(j − 1)| ≤ μξ(j − 1) is satisfied, the algorithm is terminated, where ξ(j) is the mean reconstruction error at the jth iteration, μ is a small positive number. In this paper, U⁽²⁾ and U⁽³⁾ are initialized to the identity matrix for the sake of simplicity, and the algorithm usually converges within two iterations in the following experiments.

Algorithm 1. MPCA transform

View Table | View all tables in this article

4.3. Shrinkage, inverse MPCA transform and aggregation

Given ${\left\{{𝒮}_i^k\right\}}_{i=1}^{n_k}$, U⁽¹⁾, U⁽²⁾, U⁽³⁾ by applying Algorithm 1 to the group $P_k={\left\{{𝒴}_i^k\right\}}_{i=1}^{n_k}$, we aim to shrink transform tensor data ${\left\{{𝒮}_i^k\right\}}_{i=1}^{n_k}$ and apply inverse MPCA transform to obtain the filtered OCT volume data. In fact, the decomposition ${𝒴}_i^k={𝒮}_i^k{\times }_1{\mathbf{U}}^{(1)}{\times }_2{\mathbf{U}}^{(2)}{\times }_3{\mathbf{U}}^{(3)}$ can be rewritten as

Inspired by [35], we define the following shrinkage operator

By transforming ${\widehat{𝒮}}_i^k$ back to time domain, we obtain the denoised result of ${𝒴}_i^k$ as

For N = 1, the input image patches are vectors, where only one projection matrix is computed and MPCA is reduced to PCA. This means that MPCA generalizes PCA. Moreover, PCA based methods cannot directly process multi-dimensional data, and the image-to-vector method ignores the local correlation inside the image and destroys the structure information of image. On the contrary, the image-to-tensor method of MPCA can better preserve the structure of image cubic blocks and the features of local OCT cubic blocks can be efficiently exploited by the computed three projection matrices.

4.4. Iterative regularization

In order to obtain better denoising performance, iterative regularization method is used to update the noisy image by adding some method noise in iteration

Algorithm 2. Proposed multi-frame OCT images denoising method

View Table | View all tables in this article

5. Experimental results

In this section, we test the proposed method on two OCT image datasets and compare it with several state-of-the-art methods, such as nonlocal mean method (PNLM) [24], 2D adaptive complex diffusion method [15], Bayesian method [17], sparse representation method (MSBTD) [26] and multi-frame wavelet method [19]. In order to obtain well denoising performance and reduce the computational cost, we keep w and h fixed, i.e. w = h = 9, for all experiments. But other parameters of the proposed method are tentatively chosen as values that obtained better peak signal-to-noise ratio (PSNR). As shown in Fig. 1, the PSNR of human retinal data is insensitive to parameters. However, parameters have a greater impact on the PSNR of pig eye data, which reaches a maximum at η = 0.2 when keeping the scaling factor (γ = 0.65) fixed and reaches a maximum at γ = 0.65 when keeping the relaxation parameter (η = 0.2) fixed. Therefore, in this paper, we set η = 0.2 and γ = 0.65.

 

Fig. 1 The PSNR of the proposed method under different parameters. (a) influence of parameter η. (b) influence of parameter γ.

Download Full Size | PDF

5.1. Image quality metrics

In order to quantitatively evaluate and analysis the despeckled performance of the proposed method, several metrics are used to quantify the despeckled image quality. These metrics include the cross-correlation (XCOR), the PSNR, the average contrast-to-noise ratio (CNR) and the average equivalent number of looks (ENL). The XCOR depicts the similarity between the despeckled image and the reference one (e.g. Figs. 2(a) and 2(b)), the larger value implies the recovered image is more similar to the reference image. The PSNR measures the intensity difference between two images, and would be higher for high quality despeckled images. The CNR measures the contrast between foreground objects (e.g. the areas in the blue rectangles in Fig. 3(a)) and background region (e.g. the areas in the white rectangles in Fig. 3(a)). The ENL accounts for the smoothness in homogeneous areas (e.g. the areas in the green rectangles in Fig. 3(a)). These metrics are defined as follows:

 

Fig. 2 Reference images. (a) pig eye data. (b) human retinal data.

Download Full Size | PDF

 

Fig. 3 The despeckled results of deferent methods. (a) Noisy image. (b) Averaged image (three frames). (c)–(h) Results of MSBTD method, PNLM method, Baysian method, complex diffusion method, multi-frame wavelet method and the proposed method.

Download Full Size | PDF

5.2. Results on the pig eye data

In following experiments, we use the public pig eye OCT dataset [19], which is acquired by scanning a pig eye in high speed mode with 768 A-scans, using a Spectralis HRA & OCT (Heidelberg Engineering). The dataset includes 35 sets at different eye positions, corresponding to a complete 0.384-mm shift in the transversal direction. Each set contains 13 frames (768 × 496 pixels) sharing the same imaging position. The pixel spacing of each frame is 14 μm in the transversal direction and 3.87 μm in the axial direction. To assess the quality of the despeckled OCT images, a reference image (Fig. 2(a)) is obtained by averaging all preregistered frames. For multi-frame based methods, we take 3 frames, from “35_6” to “35_8”, as input. For single-frame based method, we take the average of 3 frames (Fig. 3(b)) as their input.

The despeckled images of the proposed method and the compared methods are shown in Fig. 3. It can be found that the proposed algorithm realizes delicate edge preservation. Stronger comparison is presented in the closeups. The enlarged red rectangle regions are shown in Fig. 4. We can see that sufficient speckle noise is suppressed and more details are preserved. The quantification results of different methods are listed in Table 1, and the best result is marked in bold. As can be seen, our algorithm harvests overall best results. From Fig. 3 and Table 1, we can find that our algorithm is superior in these metrics, while has good speckle reduction and visual effect. However, as shown in Fig. 4(h), some foreground information is to a certain extent oversmoothed. The main reason for this phenomenon is that the noise level of the input image is high. Therefore, the size of cubic block must be set to a relatively larger constant, which will reduce the number of similar blocks in the corresponding group. Consequently, most of the transform coefficients are shrunk to 0, which means only the main structure can be preserved.

 

Fig. 4 Closeups of the red rectangle region in Fig. 3(a). (a) Noisy image. (b) Averaged image (three frames). (c)–(h) Results of MSBTD method, PNLM method, Baysian method, complex diffusion method, multi-frame wavelet method and the proposed method.

Download Full Size | PDF

Table 1. Performance metrics for different methods operating on the pig eye data.

View Table | View all tables in this article

5.3. Results of human retinal data

To further test the despeckling ability of the proposed methods, we have implemented the algorithm on human retinal OCT images. We use the public dataset [25]. Each image is acquired by 840-nm wavelength spectral domain OCT imaging system from Bioptigen, Inc. with an axial resolution of ∼4.5μm. Each OCT volume data contains 5 frames and the averaged image (Fig. 5 (b)) is the average of the 5 fames. More information about the image acquisition process can be found in paper [25].

 

Fig. 5 The despeckled results of different methods. (a) Noisy image. (b) Averaged image (five frames). (c)–(h) Results of MSBTD method, PNLM method, Baysian method, complex diffusion method, multi-frame wavelet method and proposed method.

Download Full Size | PDF

Figure 5 shows the results of different methods. Comparing the despeckling results, we can see that multi-frame wavelet method blur the fine edges, but the proposed method preserve image details. To better observe the denoising effect, further comparison is presented and the enlarged regions are shown in Fig. 6. It can be found that the proposed method not only obtains better performance on speckle noise suppression but also preserves thin structures. The quantification results listed in Table 2 also show the superior performance of the proposed method on human retinal images compared with other speckle removal methods.

 

Fig. 6 Closeups of the red rectangle region in Fig. 5(a). (a) Noisy image. (b) Averaged image (five frames). (c)–(h) Results of MSBTD method, PNLM method, Baysian method, complex diffusion method, multi-frame wavelet method and proposed method.

Download Full Size | PDF

Table 2. Performance metrics for different methods operating on the human retinal data.

View Table | View all tables in this article

6. Conclusion

In this paper, we propose a new method for denoising multi-frame OCT images. The good performance of the proposed method is mainly attributed to three reasons: nonlocal self-similarity prior is fully utilized by the proposed algorithm; MPCA transform works directly on the image blocks in each mode space rather than 1D vector space, therefore, image cubic blocks can be efficiently represented by MPCA bases; the proposed shrinkage operator is a data driven method. From the results of the experiments, we can see the proposed method outperforms other popular methods in terms of speckle noise reduction and fine details preservation.