Open Access
5 March 2015 Multiframe denoising of high-speed optical coherence tomography data using interframe and intraframe priors
Author Affiliations +
Abstract
Optical coherence tomography (OCT) is an important interferometric diagnostic technique, which provides cross-sectional views of biological tissues’ subsurface microstructures. However, the imaging quality of high-speed OCT is limited by the large speckle noise. To address this problem, we propose a multiframe algorithmic method to denoise OCT volume. Mathematically, we build an optimization model which forces the temporally registered frames to be low-rank and the gradient in each frame to be sparse, under the constraints from logarithmic image formation and nonuniform noise variance. In addition, a convex optimization algorithm based on the augmented Lagrangian method is derived to solve the above model. The results reveal that our approach outperforms the other methods in terms of both speckle noise suppression and crucial detail preservation.

1.

Introduction

Optical coherence tomography (OCT), which dates back to 1991,1 provides cross-sectional views of biological tissues’ subsurface microstructures2,3 and has become a widely used diagnostic technique in the medical field due to its noninvasive nature. For example, ophthalmology has benefited from OCT, as it can image the retina and thus aid in the diagnosis.4 Different from common CCD imaging, OCT is an interferometric technique, and typically uses near-infrared laser light to penetrate into the scattering medium before capturing the backscattered optical waves for the final imaging.5 With the development of ultrahigh resolution OCT6 and Fourier-domain OCT,7,8 it is feasible to visualize biological tissues at a cellular level and up to the depth of nearly 1 mm below the surface with high sensitivity and image quality. Furthermore, the image acquisition speed of OCT systems has greatly improved with the development of high-speed sensors and tunable lasers with MHz scanning rates, which allow real-time imaging of in vivo tissues.8

One of the important challenges limiting high-speed OCT’s development is its unsatisfying image quality caused by speckle noise. Due to the coherence of optical waves, speckle noise arises under limited spatial-frequency bandwidth of the interference signals.5 The generation mechanism of OCT determines that the properties of speckle noise are related not only to the laser source, but also to the tissue’s structural properties,911 and thus results in nonuniform speckle noise over the entire image. Due to the significance of high precision in medical diagnosis, it is vital to remove speckle noise from OCT images for image quality enhancement.

Much effort has been exerted for denoising OCT images, and various approaches are reported. These methods mainly fall into two categories, namely single-frame methods and multiframe methods. Single-frame methods often assume a prior model (either parametric or nonparametric) for the latent signal and noise, and then remove the noise determinatively or probabilistically from the input single image. Filtering is a widely used strategy, and there are many OCT denoising filters. For example, Ozcan et al.12 apply various digital filters for denoising OCT images, and the results indicate that the nonorthogonal wavelet filter together with the enhanced Lee and the adaptive Wiener filters can significantly reduce speckle noise. Based on the wavelet filter, Yue et al.13 utilize the iterative edge enhancement feature of nonlinear diffusion to improve the denoising results. Similarly, Zhang et al.14 use the nonlinear diffusion in the Laplacian pyramid domain to filter ultrasonic images. The Kovesi Nw filtering technique and the Laplacian pyramid nonlinear diffusion technique are unified together in Ref. 15 to remove both shot noise and speckle noise from OCT images. In addition to the filtering techniques, there are also some other approaches such as regularization and Bayesian inference. The work in Ref. 16 uses a regularization method to minimize the Csiszar’s I-divergence measurement, which would extrapolate additional details from the input noisy images to improve the visual effects. Wong et al.17 and Cameron et al.18 use a statistic Bayesian least square model to reduce OCT speckle noise in logarithmic space. Xie et al.19 take image contents into consideration and propose a salient structure extraction algorithm combined with an adaptive speckle suppression term to enhance ultrasound images. Different from the above filtering or statistic methods, based on sparse coding, the work in Ref. 20 learns an overcomplete dictionary from high signal-to-noise ratio (SNR) images and then utilizes this dictionary to reconstruct low-SNR OCT images and achieve significant noise suppression. In all, making use of the intrinsic redundancy within a single OCT frame can help noise removal to a large extent.

Benefiting from the development of high-speed OCT imaging systems, the correlation between adjacent frames increases, and researchers gradually begin to use multiple OCT frames to attenuate speckle noise and propose various multiframe methods. Technically these methods could be further classified into hardware methods and algorithmic methods. The general idea of hardware methods is to change the parameters of OCT imaging systems to decorrelate speckle noise in different frames, and then compound these frames after image registration to get a noise-free OCT image. For example, Refs. 10 and 2122.23 alter the angle of incident light while capturing different frames,24,25 change the detection angle of backreflected light,26 and change the frequency of the laser beam. Using the dynamic focus OCT facility, the work in Ref. 27 compares the performance of several spatial compounding methods and shows their respective pros and cons. The methods include median filtering, random weighted averaging, and random pixel selection. The main disadvantage of these hardware methods is the complex procedure of data acquisition, which would greatly increase the design complexity of OCT imaging systems.18 In addition, the performance of spatial compounding methods is not satisfactory.

There are several recently reported algorithmic methods for denoising multiple OCT images that make use of the redundancy of latent sharp OCT images in the frequency domain. For example, the work in Ref. 28 performs the three-dimensional (3-D) curvelet transform to the volume data, then thresholds the coefficients, and finally does the inverse 3-D curvelet transform to realize the noise removal. Under the same framework, Mayer et al.4 choose the wavelet domain for coefficient thresholding. In spite of the promising performance, these denoising algorithms run the risk of losing crucial details by directly truncating the coefficients.

Making use of both the intraframe and interframe redundancies of OCT volume data, this paper proposes an algorithmic multiframe optimization method to denoise OCT images. Within each single frame, an OCT image is statistically similar to a natural image and its pixel gradient map tends to be sparse. This serves as the intraframe prior. To avoid piecewise constant artifacts by simply using the total variation constraint, and considering the excellent detail-preservation performance of the low-rank prior in matrix completion29,30 and image reconstruction,31,32 we make use of the interframe redundancy by registering the OCT images along the image count dimension to form a low-rank volume. Subjecting the images to both the image formation model and the nonparametric bound constraint of nonuniform speckle noise, we build a preliminary nonconvex optimization model which jointly minimizes the rank of temporally registered OCT volume and forces the sparsity of its spatial gradient. To solve the above model, we first perform some mathematical transformations and approximations for convexification. Then, considering the superior convergence property of the augmented Lagrange multiplier (ALM) method33 for solving constrained optimization problems, as utilized in Ref. 34, we derive a numeric algorithm based on the ALM method to solve the convexified optimization model. Experiments on a pig eye, human retina, and orange OCT data show that our denoising technique could effectively reduce speckle noise while preserving rich details and that it exhibits superior performance to the other popular methods.

The remainder of this paper is organized as follows: Sec. 2 sequentially describes the preprocessing operations—frame registration, noise variance estimation, model construction, and algorithm derivation. Then in Sec. 3, we apply our method to real-captured OCT images including pig eye, human retina, and orange data and compare our approach with several other popular methods in terms of both visual quality and quantitative evaluation. Finally, we conclude this paper with some conclusions and discussions in Sec. 4.

2.

Method

In this section, we explain the whole operation framework of our approach as diagramed in Fig. 1. Our method mainly includes three steps: (1) preprocessing step including pixel logarithm, frame registration, and noise variance estimation, (2) modeling step, and (3) solving step based on the ALM method. The implementation of each step is detailed in the following subsections.

Fig. 1

Framework of the proposed approach. After taking the logarithm, each OCT frame is represented by a single column in log transformed space, and there is slight misalignment of one frame compared with other frames, as logM shows. Then, by frame registration and noise variance estimation, we get not only the optimization constraints, but also the optimization objective terms—the low rank of L and the sparsity of L. Finally, the model is iteratively solved by a convex optimization algorithm based on the augmented Lagrange multiplier (ALM) method, thus N is separated from L.

JBO_20_3_036006_f001.png

2.1.

Preprocessing

In this section, we conduct three preprocessing operations on the captured OCT frames including pixel logarithm, frame registration, and noise variance estimation. Due to the interference between the reference beam and the backscattering beam in OCT facilities, the speckle noise in OCT images is multiplicative5,9,10 and can be described as

Eq. (1)

M(s)=L(s)×N(s),
where M(s) denotes the captured value at location s, while L(s) and N(s), respectively, denote the ground truth and measurement noise at the same location. To convert the correlation between L and N from multiplication to addition, we take the logarithm to both sides of Eq. (1) and get

Eq. (2)

logM(s)=logL(s)+logN(s).
In the following, we assume that all the variables have been logarithmically transformed.

The second preprocessing operation is frame registration. Although the OCT equipment usually has a high capture speed, there always tends to be a slight misalignment among OCT image sequences in vivo imaging due to the object motion and other system factors.4 Here, we adopt the registration method in Ref. 4, where a powell optimizer is utilized for minimizing the sum of squared distances among multiple registered images. Specifically, the approach applies translations and rotations to warp the pixels describing the same tissue position in different frames to the same image coordinate.

The third requisite preprocessing operation is the variance estimation of speckle noise. Considering the satisfying performance of the median absolute deviation (MAD) method for noise estimation,35 we here utilize the MAD method as described in Ref. 18. Due to similar tissue properties and light directions in the same neighborhood, we assume a uniform noise variance for each pixel within a small patch. Therefore, the MAD within a small neighborhood N of pixel s is first computed in logarithmic space as18

Eq. (3)

σ^(s,N)=1.4826MN(s)(|logM(si)MN(s)(logM)|),
where MN(s) denotes the median value over s’s neighborhood N(s), and siN is the i’th neighboring pixel of s. To make the estimated deviation more precise, we choose a larger neighborhood N2(s) and calculate the local standard deviation σ^ of its subneighborhood N1(s). Then, we regard the mode of these σ^ as the preliminary noise deviation at position s:

Eq. (4)

σ¯(s)=modeN1(s)N2(s)[σ^(s,N1)].

Finally, to force the smoothness of noise variance among adjacent pixels, we conduct a cubic spline fitting process to amend σ¯ and get the final standard deviation estimation of the OCT noise (the corresponding noise variance can be calculated as the square of the estimated standard deviation). Empirical studies in Ref. 18 state that the noise estimation performs best when the pixel numbers of N1 and N2 are 9×9 and 15×15, respectively.

2.2.

Modeling

In this section, we build our optimization model incorporating both of the interframe and intraframe priors. Suppose that there are k frames in the OCT volume, and the pixel number of each frame is m×n. We denote the temporally registered noisy OCT images, their latent sharp version, and the measurement noise as M, L, and N, respectively. Mathematically, the image formation equation can be written as

Eq. (5)

M=L+N.

By representing each frame as a column vector, the dimensions of M, L, and N are all (m×n)×k. After frame registration, theoretically the entries in one specific row of L should be exactly the same, as shown in Fig. 1. Therefore, we treat L as a low-rank matrix and use the minimization of its nuclear norm L*, which calculates the sum of L’s singular values34 as the interframe prior constraint.

According to the statistical studies,36,37 the adjacent pixels in natural images have similar intensities. Thus, the image gradient centers around zero and follows a heavy tailed distribution, i.e., the gradient of natural images is sparse. Although captured via a different imaging mechanism from the usual CCD imaging methods, the OCT images still follow similar statistics, and we impose the gradient sparsity of the latent OCT images as the intraframe prior. Specifically, the l0-norm, which counts the number of nonzero entries in a matrix, can model the sparseness quite well, thus we can minimize L0 as the intraframe constraint, where is the gradient calculating operator. Adopting the same representation in Ref. 38, here we use matrix multiplication for the gradient calculation, namely L0=a=12HaL0, where H1 and H2 are, respectively, the horizontal and vertical gradient operators and are defined as the diagonal matrices of corresponding high-pass filters h1=[1,1] and h2=[1;1].

As mentioned before, speckle noise includes both the image information and the zero-mean noise. Since what we are most concerned with is the removal of the latter component, we treat the first component as a part of the latent sharp image and concentrate on attenuating the zero-mean noise. According to the three sigma rule, which indicates that nearly all (99.73%) of the instances of a random variable lie within three times the standard deviation from its mean, we can approximatively formulate the noise constraint as

Eq. (6)

|N|3σ,
where σ is the standard deviation matrix whose dimension is (m×n)×k. By introducing a non-negative matrix variable ε, we can transform the above inequality into an equality as

Eq. (7)

NN9σσ+ε=0,
in which is the entry-wise product, i.e., for two matrices X and Y, (XY)ij=XijYij.

Based on the above notations, the optimization model for denoising is defined as

Eq. (8)

{L*,N*}=argminL*+λa=12HaL0,s.t.M=L+NNN9σσ+ε=0,
with λ being a positive weighting parameter to balance different objective regulation terms.

2.3.

Solving the Model

In this section, we derive our optimization algorithm based on the ALM method to solve the above model in Eq. (8). The model is obviously nonconvex, so we first conduct several convexification transformations to the model. As shown in Refs. 39 and 40, replacing the l0-norm with the l1-norm is one typical convexification transformation. Here, we replace HaL0 with HaL1, where ·1 denotes the sum of the matrix entries’ absolute values. Further, it is hard to directly utilize the ALM method to solve the convexified objective function due to its high nonlinearity. To address this problem, we replace the variables whose nuclear norm or l1-norm needs minimization with two introduced auxiliary variables S1 and S2. In addition, we pack a=12HaL1 as PL1 for computation simplicity, where P=[H1;H2]. With the above substitutions, we can rewrite the model as

Eq. (9)

min.S1*+λS21,s.t.G1=S1L,G2=S2PL,G3=ML+N,G4=NN9σσ+ε.
Here, G14 are supposed to be 0 in theory.

As stated before, we utilize the ALM method to solve the above model. This method adopts an iterative optimization strategy, and successively updates every variable within each iteration. In the following, we derive the updating rules for each variable. First, the augmented Lagrangian function of Eq. (9) is

Eq. (10)

f=S1*+λS21+j=14(Yj,Gj+θ2GjF2),
where ·,· denotes the inner product, Y defines the Lagrangian multipliers (in matrix form), and ·F refers to the Frobenius norm that calculates the root of all the square entries’ sum in a matrix. Here, θ is a penalty parameter balancing the four equation constraints in Eq. (9) and follows the standard ALM updating rule as θ(k+1)=min(ρθ(k),θmax), where ρ and θmax are both user-defined parameters and k indexes the iteration. The updating rules of the other variables including S, L, N, ε, and Y are derived as follows.

Optimize S. By removing all the items irrelevant to S1 in f, we can get

f(S1)=S1*+θ2S1(L(k)θ1Y1(k))F2.

According to the ALM algorithm, we can get the updating rule of S1 as

Eq. (11)

S1(k+1)=Usθ1(Stemp)VT,
where UStempVT is the singular value decomposition (SVD) of L(k)θ1Y1(k), and
sθ1(x)={xθ1,x>θ1x+θ1,x<θ10,others.

Similarly, keeping only the items related to S2 in f yields

f(S2)=λ[S21+θ2λS2(PL(k)θ1Y2(k))F2],
and we can get the updating rule of S2 as

Eq. (12)

S2(k+1)=sλθ(PL(k)θ1Y2(k)).

Optimize L and N. By keeping only the items related to L, f is simplified as

f(L)=θ2S1L+θ1Y1F2+θ2S2PL+θ1Y2F2+θ2MLN+θ1Y3F2,
and the partial derivative of f(L) with respect to L is
f(L)L=θ(LS1(k)θ1Y1(k))+θ[PTPLPT(S2(k)+θ1Y2(k))]+θ(LM+N(k)θ1Y3(k)).

Similarly, the partial derivative of f(N) with respect to N is

f(N)N=θ(NM+L(k)θ1Y3(k))+θ[NN9σσ+ε(k)+θ1Y4(k)]2N.

Obviously, it is hard to get the closed-form solution to either [f(L)/L]=0 or [f(N)/N]=0, so we resort to the gradient descent method to approximatively update these two variables as

Eq. (13)

L(k+1)=L(k)Δ×f(L)L|L=L(k),

Eq. (14)

N(k+1)=N(k)Δ×f(N)N|N=N(k).
Here, Δ stands for the learning rate.

Optimize ε. The derivative of f with respect to ε is

f(ε)ε=θ(ε+N(k)N(k)9σσ+θ1Y4(k)).

Under the non-negative assumption on ε, we can get its updating rule as

Eq. (15)

ε(k+1)=max(9σσN(k)N(k)θ1Y4(k),0).

All the algorithm parameters are set as follows: λ=0.2, θ(0)=0.01, ρ=1.6, θmax=10, ΔL=0.01, and ΔN=0.05. These constant parameters are statistically set by testing the algorithm on a series of OCT data to obtain the best denoising performance and least running time. In addition, all the parameters are fixed across all the experiments in this paper. We take the average of input preregistered frames as the initialization of the denoised frame. For more clarity, the entire iterative algorithm based on the above derivations is summarized in Algorithm 1.

Algorithm 1

The proposed multiframe algorithm for OCT denoising.

Input: Capturing data M and estimated noise standard deviation σ
Output: Denoised frames L and separated noise N
1 L(0)=M¯, N(0)=ML(0), ε(0)=0; Y14(0)=0;
2 Whilenot convergeddo
3  Update S{1,2}(k+1) according to Eqs. (11) and (12);
4  Update L(k+1) according to Eq. (13);
5  Update N(k+1) according to Eq. (14);
6  Update slack variables ε(k+1) according to Eq. (15);
7  Y{14}(k+1)=Y{14}(k)+θ1G{14};
8  θ(k+1)=min(ρθ(k),θmax);
9  k:=k+1.
10 end

2.4.

Evaluation Criteria

To quantitatively evaluate the denoising performance of different approaches, we utilize three widely used image quality criteria including the peak SNR (PSNR), the structure similarity (SSIM),41 and the figure of merit (FOM)42,43 as evaluation metrics.

  • (a) PSNR: In digital image recovery, PSNR has traditionally been widely used to assess the quality of the processed image Lm×n, with respect to its ground truth Im×n. PSNR is calculated as

    Eq. (16)

    PSNR=10×log10{MAX21mni=1mj=1n[L(i,j)I(i,j)]2},
    where MAX=2b1 is the maximum intensity of b bit images. For example, for the widely used 8 bit images, MAX=255. From the equation, we can see that PSNR intuitively describes the intensity difference between two images and would be smaller for low-quality recovered images. Empirically, the typical PSNR for visually promising images is roughly between 25 and 40 dB.

  • (b) SSIM: The structure similarity criterion is proposed in Ref. 41 to measure the SSIMs between two images. This criterion first selects two corresponding patch sets pL={pLk;k=1K} and pI={pIk;k=1K} from L and I, respectively, with K being the patch number, and then calculates the preliminary SSIM between each patch pair pLk and pIk as

    Eq. (17)

    SSIM(pLk,pIk)=(2μLkμIk+c1)(2σL,Ik+c2)[(μLk)2+(μIk)2+c1][(σLk)2+(σIk)2+c2],
    where μLk and μIk are, respectively, the average pixel intensities of patch pLk and pIk. σL and σI are the patchs’ standard variances, and σL,I is the covariance between pLk and pIk. In addition, c1=(k1MAX)2 and c2=(k2MAX)2 are two constants, with k1 and k2 being two user-defined parameters whose default values are, respectively, 0.01 and 0.03. The final SSIM score between two images is the average of all the patches’ preliminary SSIM scores. The SSIM score ranges from 0 to 1 and is higher when two images have more similar structural information. Compared with traditional metrics such as PSNR, which only reveals the intensity differences between two images, SSIM reflects the similarity in structural information of an image pair, and thus is closer to human perception.

  • (c) Edge preservation: Further processing of denoised OCT images would be likely to involve the segmentation of layers or identification of a particular image feature. Thus, the preservation of edges in denoised OCT images is very important. Here, we also adopt the FOM42,43 to evaluate the edge preservation abilities of various denoising methods. FOM is defined as

    Eq. (18)

    FOM=1max(nL,nI)i=1nL11+γdi2,
    where nL and nI are, respectively, the numbers of detected edge pixels in the reconstructed image and the ground-truth image, di is the Euclidean distance between the i’th detected reconstructed edge pixel and its nearest ground-truth edge pixel, and γ is a scaling constant balancing the penalties to smeared edges and isolated edges, which is typically set to be 1/9.44 In this paper, we use the Canny edge detector under default parameter settings in MATLAB. FOM score ranges from 0 to 1 and is higher when the reconstructed image has clearer edges and is more similar to the ground-truth image.

3.

Experimental Results

In this section, we test our denoising approach on three OCT image sets and compare our denosing results with those of several previously reported popular methods. In addition to the visual comparison, we also perform quantitative comparison in terms of PSNR, SSIM, FOM, and running time.

3.1.

Results on the Pig Eye Data and Performance Comparison

In this experiment, we use the public pig eye OCT dataset in Ref. 45, which is also used as the source data in Ref. 4. The dataset is acquired using a Spectralis HRA & OCT (Heidelberg Engineering) to scan a pig eye in the high-speed mode with 768 A-scans. There are a total of 455 images (each frame contains 768×496pixels) in the dataset including 35 sets, and 13 frames sharing the same imaging position are included in each set. All the 35 positions correspond to a complete 0.384-mm shift in the transversal direction. For each frame, the pixel spacing is 3.87μm in the axial direction and 14μm in the transversal direction. To assess the quality of the recovered images, we need a noise-free benchmark image for reference. However, due to the high-imaging speed, the captured images’ SNR is very low. Therefore, we utilize the same technique as in Ref. 4, in which the averaged image of all the 455 preregistered frames is used as the latent noise-free image.

Since the proposed method is multiframe based, we first investigate the effect of the most important parameter—the number of input frames—on our algorithm. Fixing all the other parameters, we run a MATLAB implementation of our proposed method on an Intel E7500 2.93 GHz CPU computer with 4GB RAM and 64 bit Windows 7 system and compare the performance of different numbers of input frames, ranging from 2 to 13. From the result in Fig. 2, we can see that our algorithm still works using only 2 frames, and the reconstruction quality gradually improves using an increasing number of frames from 2 to 8. This trend is much less obvious with more than 8 frames. Based on this observation, we use 8 input frames in the following experiments to compare our approach with other methods. Note that there exists some slight drop when the frame number grows larger than 8. This is because the frame registration gets more difficult for a longer sequence and thus slightly hampers the reconstruction.

Fig. 2

Reconstruction quality versus input frame number. Input: 2 to 14 frames of the pig eye data. The solid blue line corresponds to the axis on the left ranging from 25 to 32, while the two dashed red lines correspond to the axis on the right ranging from 0.3 to 1.

JBO_20_3_036006_f002.png

Next, we run our algorithm and the other four popular denoising methods on the OCT dataset for comparison. The four methods include the complex diffusion method,46 the Bayesian method,17 the non-stationary speckle compensation method (NSC),18 and the multiframe wavelet OCT denoising method.4 To validate the superior effectiveness of our approach, we compare all the algorithms’ performances visually and quantitatively. What should be noticed is that the complex diffusion method, the Bayesian method, and NSC are all single-frame methods, so we take the average image of the registered input frames as their single-frame input. In addition, the first two methods assume a spatially invariant noise parameter (standard deviation). Correspondingly, we use the maximum in the estimated deviation matrix as the input standard deviation of noise. By random selection, the serial numbers of the input 8 sequential frames are from “35_6” to “35_13.”

The recovered images are shown in Fig. 3, where only the first frame is presented for each method. We can see that the recovered images of the anisotropic diffusion method, the Bayesian method, and NSC still contain undesired noise, which largely degenerates the image quality and makes these three methods less competitive than the other two techniques. On the whole, the multiframe wavelet method and our method are both superior to the three single-frame approaches. This superiority is attributed to the fact that the multiframe methods utilize the interframe correlation and redundancy, which would offer more information of the latent noise-free data. Comparing the results of these two multiframe methods, we can see that in the smooth regions, the wavelet method leaves out more noise than our method. In the textured regions, the result of our method maintains a higher color contrast which would improve the visual effects. A closer comparison is presented in the closeups. For example, in the green-rectangle-highlighted region, the wavelet method nearly blurs out the details of the white spot on the left side, while our method still contains gray value changes which would provide important information for diagnosis.

Fig. 3

Comparison with the other four popular methods. Input: 8 frames of the pig eye data. (a) is the original image in log transformed space, while (b) is the averaged image of 455 registered frames. (c) is the averaged image of the input 8 frames, and (d)–(g) are the recovered results of four popular methods. The result of our method is shown in (h). The two clipped patches on the right of each subfigure are closeups of the regions of interest.

JBO_20_3_036006_f003.png

Numerical assessments are shown in Table 1. For the entire image, we can see that our method could raise the noisy images’s PSNR from 17.19 to 31.74 dB and the SSIM from 0.13 to 0.91. In terms of all the three evaluation criteria including PSNR, SSIM, and FOM, our method consistently has an advantage over the other methods. Comparing the two multiframe methods, namely the multiframe wavelet method and our method, we can see that our approach is superior in PSNR and SSIM by around 1 dB and 0.1, respectively. Our superior performance is mainly attributed to two factors: the combined constraints from both the interframe and intraframe priors and the good convergence of the derived algorithm. Comparing the numerical evaluation results of the two selected regions of interests, we can also see the clearer advantage of our method over the other methods. What is more, the running time comparison is also provided in Table 1 (the noise estimation time is also included for all the four algorithms). We run the MATLAB codes of our algorithm and the other three methods except for NSC on our computer, while the running time of NSC is provided by its proposers who run their MATLAB and C++ implementation on a different platform. We can see that our approach needs around 36 s to process one frame and is of similar efficiency to the Bayesian method which is the fastest algorithm among the popular methods except for NSC.

Table 1

Quantitative comparisons among different denoising methods.

MetricInputAverageDiffusionBayesianNSCaWaveletOurs
Entire imagePSNR (dB)17.1924.5629.1428.3829.8230.7531.74
SSIM0.120.450.730.700.810.860.91
Running time (s)793326036
Red clipPSNR (dB)15.0322.0226.6026.0727.4727.8528.92
SSIM0.060.290.650.630.710.730.81
FOM0.430.460.510.570.600.610.63
Green clipPSNR (dB)15.1321.9126.1425.3526.8327.8328.75
SSIM0.060.250.600.570.660.720.80
FOM0.480.490.510.530.580.580.58

aThe performance of NSC is tested by its proposers on the AMD Athlon X3 II CPU with 8 GB of RAM and Windows 7 64-bit system, using MATLAB and C++ programming for high-computation efficiency.

Note: The bold values represent the best performance in terms of each metric (i.e., highest PSNR, highest SSIM, highest FOM and lowest running time) among all the methods.

3.2.

Analyzing Optical Coherence Tomography Images of Human Retina

To test the practical denoising effectiveness of our method, we conduct a denoising experiment on human retinal OCT images. We use the same public dataset as that used in Ref. 47, which is acquired by an spectral domain (SD)-OCT imaging system from Bioptigen Inc. with 4.5-μm axial resolution, 500 A-scans per B-scan, and 5 azimuthally repeated B-scans in each volume. Referring to the processing progress described in Sec. 2.1, we first register the OCT frames and then use different methods to denoise these frames. Considering that the anisotropic diffusion method, the Bayesian method, and NSC leave too much noise on the recovered images, and thus show little competitiveness compared with the other two multiframe methods, here we only present the denoising results of the multiframe wavelet method and our method for clearer comparison.

The results are shown in Fig. 4, which exhibit a similar performance ranking to that of the pig eye data. Comparing the denoising results produced by the multiframe wavelet method and our proposed method carefully, we can see that the result of the wavelet method contains undesired edge burrs, while our result presents clearer layer boundaries (such as the horizontal layer edges in the two selected regions), which would greatly help in follow-up analysis of the denoised images, such as OCT layer segmentation and diagnosis.

Fig. 4

Denoising results of the human retinal OCT images. Input: 5 frames of the human retina data. (a) shows 1 of the 5 captured frames, and (b) is the average of the 5 frames. (c) and (d) are, respectively, the results of the multiframe wavelet method and our method. Closeups of two selected regions of interest are shown in (e) and (f), which offer a clearer comparison.

JBO_20_3_036006_f004.png

3.3.

Phantom Study Using the Orange Data

To further validate the advantages of the proposed method in preserving image details while removing noise, we use the same orange dataset as in Ref. 48 to conduct a phantom study. The dataset is acquired by an SD-OCT system with the axial resolution being 4μm and the transversal resolution being 12μm. This data include 100 aligned frames and contain many thin structures corresponding to the tangerine pith. Similar to the experiment on human retina, considering the multiframe methods largely outperform the single-frame methods in both noise attenuation and detail preservation, here we only show the results of the multiframe wavelet method as well as ours. We still use 8 input frames and the results are shown in Fig. 5.

Fig. 5

Phantom study on the proposed method: (a) shows one of the captured frames. (b) and (c) are, respectively, the average of the input 8 frames and all 100 frames. (d) and (e) show, respectively, the results of the multiframe wavelet method and our approach.

JBO_20_3_036006_f005.png

For performance comparison between the two methods, we consider the average of 100 frames as benchmark, as shown in Fig. 5(c). From the results in Fig. 5(d), we can clearly see that some fine structural details are smoothed out by the wavelet method, especially in the region highlighted by the red rectangle. This is because the wavelet method truncates the wavelet coefficients of the input frames to attenuate noise, and some entries describing fine structures tend to be treated as speckle noise and removed. On the contrary, our method can preserve the thin structures quite well during noise suppression, as shown in Fig. 5(e). Compared against the results in Fig. 5(c), one can see clearly that our approach is more competitive in such cases with rich thin structures.

4.

Conclusions and Discussions

4.1.

Summary and Conclusions

In this study, we propose a multiframe OCT denoising method utilizing the constraints from both interframe and intraframe priors. Specifically, the interframe prior refers to the low rank of registered OCT frames and the intraframe prior is the sparsity of the image gradient. Benefiting from the proper convexification transformations and usage of ALM, the derived algorithm converges well on different data. In addition, by incorporating a nonparametric and nonuniform noise description, our approach is applicable for different noise models.

On the adopted benchmark data, our approach could improve the OCT image’s quality by raising the image’s PSNR from 17.19 to 31.74 dB and SSIM from 0.12 to 0.91, in around 36 s for each frame. The comparisons with the other four popular methods on the three datasets reveal that our method has advantages mainly in two aspects: (1) being able to attenuate speckle noise effectively and preserve crucial image details; (2) with efficiency comparable with the reported fastest approaches. Such a high performance of the proposed method is mainly the result of the combined prior modeling and effective optimization algorithm.

4.2.

Limitations and Future Extensions

The performance of our method depends on the registration accuracy because the low-rank prior in the objective function does not hold for an unaligned frame stack. This is also a challenge for other multiframe denoising methods and needs to be addressed by the progress of noise robust matching techniques. In addition, a larger frame number is favorable to take advantage of the low-rank prior. Therefore, one needs to set the system’s frame rate to balance the noise level and the number of available frames in practical applications.

Besides, the widely used anisotropic total variation is utilized as the intraframe prior, which penalizes the diagonal gradients more significantly than the horizontal and vertical ones. This means that the utilized nonuniform constraint on image intensity changes along different directions, which may introduce undesired artifacts in the denoised images. Thus, exploring an isotropic intraframe prior would be one of our future extensions.

In addition, the running time of the proposed method can be further reduced. There are two very time-consuming calculations in the current algorithm including the requisite relatively large iteration number and the SVD. To decrease the iterations, we can utilize an accelerated gradient descent scheme to speed up the method such as using adaptive learning rates instead of a constant rate. For a faster SVD, we can replace the currently adopted default implementation in MATLAB with an accelerated SVD algorithm such as the block Lanczos algorithm.49 Furthermore, we can also implement the proposed algorithm with C or C++ and use graphics processing unit accelerating techniques48,50 to further speed up the method.

Acknowledgments

This work was supported by the National Natural Science Foundation of China, Nos. 61171119, 61120106003, and 61327902. The authors thank Professor Yong Huang and Professor Guohai Situ for their valuable suggestions and help, and also wish to thank Professor Alexander Wong for providing their running results to help with the experiments.

References

1. 

D. Huang et al., “Optical coherence tomography,” Science, 254 (5035), 1178 –1181 (1991). http://dx.doi.org/10.1126/science.1957169 SCIEAS 0036-8075 Google Scholar

2. 

J. M. Schmitt, “Optical coherence tomography (OCT): a review,” IEEE J. Sel. Top. Quantum Electron., 5 (4), 1205 –1215 (1999). http://dx.doi.org/10.1109/2944.796348 IJSQEN 1077-260X Google Scholar

3. 

J. Welzel, “Optical coherence tomography in dermatology: a review,” Skin Res. Technol., 7 (1), 1 –9 (2001). http://dx.doi.org/10.1034/j.1600-0846.2001.007001001.x 0909-752X Google Scholar

4. 

M. A. Mayer et al, “Wavelet denoising of multiframe optical coherence tomography data,” Biomed. Opt. Express, 3 (3), 572 –589 (2012). http://dx.doi.org/10.1364/BOE.3.000572 BOEICL 2156-7085 Google Scholar

5. 

J. M. Schmitt, S. Xiang and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt., 4 (1), 95 –105 (1999). http://dx.doi.org/10.1117/1.429925 JBOPFO 1083-3668 Google Scholar

6. 

M. Wojtkowski et al., “Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express, 12 (11), 2404 –2422 (2004). http://dx.doi.org/10.1364/OPEX.12.002404 OPEXFF 1094-4087 Google Scholar

7. 

A. Fercher et al., “Measurement of intraocular distances by backscattering spectral interferometry,” Opt. Commun., 117 (1C2), 43 –48 (1995). http://dx.doi.org/10.1016/0030-4018(95)00119-S OPCOB8 0030-4018 Google Scholar

8. 

M. Wojtkowski et al., “In vivo human retinal imaging by Fourier domain optical coherence tomography,” J. Biomed. Opt., 7 (3), 457 –463 (2002). http://dx.doi.org/10.1117/1.1482379 JBOPFO 1083-3668 Google Scholar

9. 

J. Schmitt and A. Knüttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A, 14 (6), 1231 –1242 (1997). http://dx.doi.org/10.1364/JOSAA.14.001231 JOAOD6 0740-3232 Google Scholar

10. 

M. Bashkansky and J. Reintjes, “Statistics and reduction of speckle in optical coherence tomography,” Opt. Lett., 25 (8), 545 –547 (2000). http://dx.doi.org/10.1364/OL.25.000545 OPLEDP 0146-9592 Google Scholar

11. 

B. Karamata et al., “Speckle statistics in optical coherence tomography,” J. Opt. Soc. Am. A, 22 (4), 593 –596 (2005). http://dx.doi.org/10.1364/JOSAA.22.000593 JOAOD6 0740-3232 Google Scholar

12. 

A. Ozcan et al., “Speckle reduction in optical coherence tomography images using digital filtering,” J. Opt. Soc. Am. A, 24 (7), 1901 –1910 (2007). http://dx.doi.org/10.1364/JOSAA.24.001901 JOAOD6 0740-3232 Google Scholar

13. 

Y. Yue et al., “Nonlinear multiscale wavelet diffusion for speckle suppression and edge enhancement in ultrasound images,” IEEE Trans. Med. Imaging, 25 (3), 297 –311 (2006). http://dx.doi.org/10.1109/TMI.2005.862737 ITMID4 0278-0062 Google Scholar

14. 

F. Zhang et al., “Nonlinear diffusion in Laplacian pyramid domain for ultrasonic speckle reduction,” IEEE Trans. Med. Imaging, 26 (2), 200 –211 (2007). http://dx.doi.org/10.1109/TMI.2006.889735 ITMID4 0278-0062 Google Scholar

15. 

M. Gargesha et al., “Denoising and 4D visualization of OCT images,” Opt. Express, 16 (16), 12313 –12333 (2008). http://dx.doi.org/10.1364/OE.16.012313 OPEXFF 1094-4087 Google Scholar

16. 

D. L. Marks, T. S. Ralston and S. A. Boppart, “Speckle reduction by I-divergence regularization in optical coherence tomography,” J. Opt. Soc. Am. A, 22 (11), 2366 –2371 (2005). http://dx.doi.org/10.1364/JOSAA.22.002366 JOAOD6 0740-3232 Google Scholar

17. 

A. Wong et al., “General Bayesian estimation for speckle noise reduction in optical coherence tomography retinal imagery,” Opt. Express, 18 (8), 8338 –8352 (2010). http://dx.doi.org/10.1364/OE.18.008338 OPEXFF 1094-4087 Google Scholar

18. 

A. Cameron et al., “Stochastic speckle noise compensation in optical coherence tomography using non-stationary spline-based speckle noise modelling,” Biomed. Opt. Express, 4 (9), 1769 –1785 (2013). http://dx.doi.org/10.1364/BOE.4.001769 BOEICL 2156-7085 Google Scholar

19. 

J. Xie et al., “Boundary enhancement and speckle reduction for ultrasound images via salient structure extraction,” IEEE Trans. Biomed. Eng., 53 (11), 2300 –2309 (2006). http://dx.doi.org/10.1109/TBME.2006.878088 IEBEAX 0018-9294 Google Scholar

20. 

L. Fang et al., “Sparsity based denoising of spectral domain optical coherence tomography images,” Biomed. Opt. Express, 3 (5), 927 –942 (2012). http://dx.doi.org/10.1364/BOE.3.000927 BOEICL 2156-7085 Google Scholar

21. 

N. Iftimia, B. E. Bouma and G. J. Tearney, “Speckle reduction in optical coherence tomography by path length encoded angular compounding,” J. Biomed. Opt., 8 (2), 260 –263 (2003). http://dx.doi.org/10.1117/1.1559060 JBOPFO 1083-3668 Google Scholar

22. 

L. Ramrath et al., “Towards multi-directional OCT for speckle noise reduction,” Medical Image Computing and Computer-Assisted Intervention, 815 –823 Springer, New York (2008). Google Scholar

23. 

M. Hughes, M. Spring and A. Podoleanu, “Speckle noise reduction in optical coherence tomography of paint layers,” Appl. Opt., 49 (1), 99 –107 (2010). http://dx.doi.org/10.1364/AO.49.000099 APOPAI 0003-6935 Google Scholar

24. 

A. Desjardins et al., “Speckle reduction in OCT using massively-parallel detection and frequency-domain ranging,” Opt. Express, 14 (11), 4736 –4745 (2006). http://dx.doi.org/10.1364/OE.14.004736 OPEXFF 1094-4087 Google Scholar

25. 

A. Desjardins et al., “Angle-resolved optical coherence tomography with sequential angular selectivity for speckle reduction,” Opt. Express, 15 (10), 6200 –6209 (2007). http://dx.doi.org/10.1364/OE.15.006200 OPEXFF 1094-4087 Google Scholar

26. 

M. Pircher et al., “Speckle reduction in optical coherence tomography by frequency compounding,” J. Biomed. Opt., 8 (3), 565 –569 (2003). http://dx.doi.org/10.1117/1.1578087 JBOPFO 1083-3668 Google Scholar

27. 

M. R. Avanaki et al., “Spatial compounding algorithm for speckle reduction of dynamic focus OCT images,” IEEE Photonics Technol. Lett., 25 (15), 1439 –1442 (2013). http://dx.doi.org/10.1109/LPT.2013.2266660 IPTLEL 1041-1135 Google Scholar

28. 

Z. Jian et al., “Three-dimensional speckle suppression in optical coherence tomography based on the curvelet transform,” Opt. Express, 18 (2), 1024 –1032 (2010). http://dx.doi.org/10.1364/OE.18.001024 OPEXFF 1094-4087 Google Scholar

29. 

E. Candès and B. Recht, “Exact matrix completion via convex optimization,” Commun. ACM, 55 (6), 111 –119 (2012). http://dx.doi.org/10.1145/2184319 CACMA2 0001-0782 Google Scholar

30. 

Y. Deng et al., “Noisy depth maps fusion for multiview stereo via matrix completion,” IEEE J. Sel. Top. Signal Process., 6 (5), 566 –582 (2012). http://dx.doi.org/10.1109/JSTSP.2012.2195472 1932-4553 Google Scholar

31. 

H. Ji et al., “Robust video denoising using low rank matrix completion,” in IEEE Conf. Computer Vision and Pattern Recognition (CVPR), 1791 –1798 (2010). Google Scholar

32. 

J. Suo et al., “Joint non-Gaussian denoising and superresolving of raw high frame rate videos,” IEEE Trans. Image Process., 23 (3), 1154 –1168 (2014). http://dx.doi.org/10.1109/TIP.2014.2298976 IIPRE4 1057-7149 Google Scholar

33. 

D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods, Academic Press, Waltham, Massachusetts (1982). Google Scholar

34. 

Z. Lin, M. Chen and Y. Ma, “Linearized alternating direction method with adaptive penalty for low-rank representation,” 612 –620 (2010). Google Scholar

35. 

A. Leigh et al., “Comprehensive analysis on the effects of noise estimation strategies on image noise artifact suppression performance,” in IEEE Int. Symp. Multimedia (ISM), 97 –104 (2011). Google Scholar

36. 

T. S. Cho et al., “A content-aware image prior,” in IEEE Conf. Computer Vision and Pattern Recognition, 169 –176 (2010). Google Scholar

37. 

T. S. Cho et al., “Image restoration by matching gradient distributions,” IEEE Trans. Pattern Anal., 34 (4), 683 –694 (2012). http://dx.doi.org/10.1109/TPAMI.2011.166 ITPIDJ 0162-8828 Google Scholar

38. 

F. Heide et al., “High-quality computational imaging through simple lenses,” ACM Trans. Graphics., 32 (5), 149:1 –149:14 (2013). http://dx.doi.org/10.1145/2516971 ATGRDF 0730-0301 Google Scholar

39. 

Y. Deng, Q. Dai and Z. Zhang, “An overview of computational sparse models and their applications in artificial intelligence,” Artif. Intell. Evol. Comput. and Metaheuristics, 345 –369 Springer, Berlin Heidelberg (2013). Google Scholar

40. 

H. Lee et al., “Efficient sparse coding algorithms,” Advances in Neural Information Processing Systems, 801 –808 NIPS, San Diego, California (2006). Google Scholar

41. 

Z. Wang et al., “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process., 13 (4), 600 –612 (2004). http://dx.doi.org/10.1109/TIP.2003.819861 IIPRE4 1057-7149 Google Scholar

42. 

Y. Yu and S. Acton, “Speckle reducing anisotropic diffusion,” IEEE Trans. Image Process., 11 (11), 1260 –1270 (2002). http://dx.doi.org/10.1109/TIP.2002.804276 IIPRE4 1057-7149 Google Scholar

43. 

Y. Yue et al., “Nonlinear multiscale wavelet diffusion for speckle suppression and edge enhancement in ultrasound images,” IEEE Trans. Med. Imaging, 25 (3), 297 –311 (2006). http://dx.doi.org/10.1109/TMI.2005.862737 ITMID4 0278-0062 Google Scholar

44. 

W. K. Pratt, Digital Image Processing, 2ndWiley, New York (1991). Google Scholar

45. 

M. A. Mayer et al., “Image denoising algorithms archive,” (2013) http://www5.cs.fau.de/en/research/software/idaa/ ( September ). 2013). Google Scholar

46. 

R. Bernardes et al., “Improved adaptive complex diffusion despeckling filter,” Opt. Express, 18 (23), 24048 –24059 (2010). http://dx.doi.org/10.1364/OE.18.024048 OPEXFF 1094-4087 Google Scholar

47. 

L. Fang et al., “Fast acquisition and reconstruction of optical coherence tomography images via sparse representation,” IEEE Trans. Med. Imaging, 32 (11), 2034 –2049 (2013). http://dx.doi.org/10.1109/TMI.2013.2271904 ITMID4 0278-0062 Google Scholar

48. 

D. Xu, Y. Huang and J. U. Kang, “Real-time compressive sensing spectral domain optical coherence tomography,” Opt. Lett., 39 (1), 76 –79 (2014). http://dx.doi.org/10.1364/OL.39.000076 OPLEDP 0146-9592 Google Scholar

49. 

S. Wei and Z. Lin, “Accelerating iterations involving eigenvalue or singular value decomposition by block Lanczos with warm start,” (2010). Google Scholar

50. 

Y. Huang, X. Liu and J. U. Kang, “Real-time 3D and 4D Fourier domain Doppler optical coherence tomography based on dual graphics processing units,” Biomed. Opt. Express, 3 (9), 2162 –2174 (2012). http://dx.doi.org/10.1364/BOE.3.002162 BOEICL 2156-7085 Google Scholar

Biography

Liheng Bian is currently pursuing a PhD degree with the Department of Automation, Tsinghua University, Beijing, China. He received his BE degree in electronic information and engineering from Xi’dian University, Xi’an, China, in 2013. His research interests include computational imaging and optical bioimaging.

Jinli Suo received her BS degree in computer science from Shandong University, Shandong, China, in 2004, and the PhD degree from the Graduate University of the Chinese Academy of Sciences, Beijing, China, in 2010. She is currently a lecturer with the Department of Automation, Tsinghua University, Beijing. Her research interests mainly include computer vision, computational photography, and statistical learning.

Feng Chen received his BS and MS degrees in automation from St. Petersburg State Polytechnical University, St. Petersburg, Russia, in 1994 and 1996, respectively, and the PhD degree in automation from Tsinghua University, Beijing, China, in 2000. He is a professor at Tsinghua University. His research interests are probabilistic graphical models and computer vision.

Qionghai Dai received his PhD degree in automation from Northeastern University, Shenyang, China, in 1996. He has been a faculty member since 1997 and a professor since 2005 with the Department of Automation, Tsinghua University, Beijing, China. He has published more than 120 conference and journal papers, and holds 67 patents. His current research interests include computational photography, high resolution microscopy, compressed sensing and vision.

© 2015 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2015/$25.00 © 2015 SPIE
Liheng Bian, Jinli Suo, Feng Chen, and Qionghai Dai "Multiframe denoising of high-speed optical coherence tomography data using interframe and intraframe priors," Journal of Biomedical Optics 20(3), 036006 (5 March 2015). https://doi.org/10.1117/1.JBO.20.3.036006
Published: 5 March 2015
Lens.org Logo
CITATIONS
Cited by 27 scholarly publications.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Optical coherence tomography

Denoising

Speckle

Wavelets

Image quality

Optimization (mathematics)

Image restoration

Back to Top