Open Access
18 July 2016 Enhanced low-rank + sparsity decomposition for speckle reduction in optical coherence tomography
Author Affiliations +
Abstract
Speckle artifacts can strongly hamper quantitative analysis of optical coherence tomography (OCT), which is necessary to provide assessment of ocular disorders associated with vision loss. Here, we introduce a method for speckle reduction, which leverages from low-rank + sparsity decomposition (LRpSD) of the logarithm of intensity OCT images. In particular, we combine nonconvex regularization-based low-rank approximation of an original OCT image with a sparsity term that incorporates the speckle. State-of-the-art methods for LRpSD require a priori knowledge of a rank and approximate it with nuclear norm, which is not an accurate rank indicator. As opposed to that, the proposed method provides more accurate approximation of a rank through the use of nonconvex regularization that induces sparse approximation of singular values. Furthermore, a rank value is not required to be known a priori. This, in turn, yields an automatic and computationally more efficient method for speckle reduction, which yields the OCT image with improved contrast-to-noise ratio, contrast and edge fidelity. The source code will be available at www.mipav.net/English/research/research.html.

1.

Introduction

Optical coherence tomography (OCT) resolves optical reflections from internal structures in biological tissues by means of noninvasive low-coherence light.1 Quantification of optical properties of the tissue enables discrimination of different tissues or different pathological states of tissue.2,3 This, furthermore, enables characterization of pathological states such as cystoid macular edema,4 central retinal artery occlusion,5 and atherosclerosis plaques.6 However, the large contrast and granular appearance of speckle stands for a major obstacle in quantitative OCT image analysis.79 Speckle is inherent random signal modulation caused by spatial and temporal coherence of the optical waves that at the same time is a basis for interferometry, the measurement technique on which OCT is based.7,9,10 Thus, speckle has a dual role as a source of noise and as a carrier of information about tissue microstructure.9 Hence, complete speckle reduction is not desirable. On the other side, with biological specimens, speckle reduces contrast and makes boundaries between constitutive tissues more difficult to resolve.7,9,11 Speckle reduction techniques generally belong to two groups: physical compounding and digital filtering.7,12 The former group reduces speckle by incoherently summing different realizations of the same OCT image.1315 These strategies achieve OCT image quality improvement proportional to the square root of the number of realizations. Digital filtering techniques aim to reduce speckle through postprocessing of an OCT image, while preserving image resolution, contrast, and edge fidelity (measured by sharpness in this paper).1618 However, as it is demonstrated in Sec. 3, state-of-the-art digital filtering methods, such as median filtering, can even decrease sharpness when reducing speckle [see also Fig. 1(g)]. Here, we propose a low-rank + sparsity decomposition (LRpSD) method to reduce speckle in OCT images. It leverages LRpSD of logarithm of intensity OCT images. Since speckle can be considered as multiplicative noise on a signal,7 logarithm of the original OCT image X yields log(X)=log(L)+log(S), where L and S, respectively, represent a “clean” OCT image and speckle. To simplify further exposition, we shall slightly abuse notation through substitutions: log(X)X, log(L)L, and log(S)S. Hence, it is assumed that an original OCT image is represented in the logarithmic domain as well as that the result of the image enhancement procedure is raised to an exponential. That is, L^exp(L^) and S^exp(S^), where the hat denotes estimation of the corresponding variable. Hence, we represent the OCT image as X=L+S. Due to the random nature of the scattering, the speckle associated with the matrix S has sparse spatial distribution. Since the clean OCT image carries information on tissue microstructure, the matrix L has a structure. Thus, L can be considered as a low-rank approximation of X. Low-rank matrix approximation (LRMA) with or without additional sparsity term is a fundamental problem in many signal processing applications.19 It is a crucial step in many machine learning2025 and signal processing2628 applications. Exact decomposition X=L+S has been known under the name robust principal component analysis (RPCA)29 or rank-sparsity decomposition.30 However, as properly noted in Ref. 20, adding the “noise” term G to the RPCA model, i.e., X=L+S+G, yields a model capable of describing empirical data more realistically. The “noise” term G can also be interpreted as a modeling error, i.e., it partially takes into account imperfections of the original RPCA model. The fundamental issue in low-rank approximations is that, due to discontinuous and nonconvex nature of the rank function, rank minimization is a nondeterministic polynomial-time (NP) hard problem. Thus, the discrete NP-hard rank minimization problem is often replaced by convex relaxation29,31,32 known as nuclear- or Schatten-1 norm.22,33 However, nuclear norm approximates rank with the sum of singular values, and that is known to be inaccurate.3436 In addition to that, since they require a priori information on the rank value, low-rank approximation methods proposed in Refs. 20 and 21 exhibit high computational complexity when the true value of the rank is not known a priori. Several recent studies have emphasized the benefit of nonconvex penalty functions compared to the nuclear norm for the estimation of the singular values.19,31,3436 In particular, it has been presented in Ref. 19 how nonconvex regularization, which promotes more sparse approximation of singular values,37 can be combined into a convex optimization problem related to the estimation of the low-rank matrices. Herein, we combine nonconvex regularization19 with sparsity constraint for LRpSD in the presence of additive white Gaussian noise (AWGN), i.e., X=L+S+G, where G stands for the AWGN with an unknown variance. In addition to yielding more accurate low-rank approximation L of X, which in turn yields the OCT image with improved contrast-to-noise ratio (CNR), signal-to-noise ratio (SNR), contrast, and edge fidelity, the proposed method does not assume a priori information on the rank value. These also are the main distinctions between the proposed LRpSD method and RPCA method in OCT image enhancement.38,39 These distinctions contribute to computational efficiency in comparison with LRpSD algorithms such as Refs. 20 and 21. The proposed method is illustrated in Figs. 1(a) to 1(c). For the sake of visual comparison, we present, in respective order, in Figs. 1(d) to 1(g) the results of OCT image enhancement by algorithms derived in Refs. 20 and 21, as well as by two-dimensional (2-D) bilateral and median filtering (see Sec. 3 for more details).

Fig. 1

(a) to (c) Flow chart of the “low-rank + sparsity” decomposition approach to speckle reduction in OCT images. Information on image quality metrics, such as CNR, SNR in dB, contrast, and sharpness, can be found in Sec. 2.3. (a) Original OCT image: CNR = 3.61, SNR = 26.23, contrast = 1.14, sharpness = 56.90. (b) Enhanced low-rank approximation of OCT image by proposed algorithm: CNR = 4.17, SNR = 32.26, contrast = 1.44, and sharpness = 61.46. (c) Sparse term containing speckle. (d) OCT image enhanced by the GoDec algorithm (rank = 35):20 CNR = 4.59, SNR = 32.52, contrast = 1.71, and sharpness = 49.01. (e) OCT image enhanced by the RNSC algorithm (rank = 35):21 CNR = 4.31, SNR = 30.61, contrast = 1.43, and sharpness = 55.72. (f) OCT image enhanced by bilateral filtering: CNR = 4.17, SNR = 35.82, contrast = 1.65, and sharpness = 59.79. (g) OCT image enhanced by median filtering: CNR = 4.5, SNR = 30.78, contrast = 1.59, and sharpness = 36.14. For visual comparison, OCT images (a) to (g) were mapped to [0 1] interval with the MATLAB mat2gray command from the interval corresponding to minimal and maximal values of each specific case. The best value for each figure of merit is in bold.

JBO_21_7_076008_f001.png

The rest of this paper is organized as follows. The details of the proposed method for LRpSD are presented in Sec. 2. That is followed by an experimental comparative performance analysis in Sec. 3 and the discussion in Sec. 4. The conclusions are presented in Sec. 5.

2.

Materials and Methods

2.1.

Related Works

Let XR0+I1×I2 be one scan of the OCT image with the size of I1×I2  pixels. The speckle, which occurs due to the random scattering of the light on tissues, acts effectively as multiplicative noise.7 That is,x(i1,i2)=l(i1,i2)×s(i1,i2), where (i1,i2) stands for pixel coordinates and x(i1,i2) stands for the intensity value at (i1,i2). By taking the log of x(i1,i2) we obtain

Eq. (1)

logx(i1,i2)=logl(i1,i2)+logs(i1,i2).
With the slight abuse of notation, we rewrite Eq. (1) on the matrix level as

Eq. (2)

X=L+S+G,
where in relation to Eq. (1), the AWGN term G with zero mean and unknown variance σ2 have been added. As discussed previously, G can also be considered as a modeling error that partially takes into account imperfections of the model. Due to the random nature of the scattering, the speckle associated with the matrix S has sparse spatial distribution. Thus, the matrix L represents the “clean” OCT image that contains information on tissue microstructure. Hence, it is justified to assume that L is a low-rank approximation of X.38,39 Thus, reduction of the speckle within the OCT image can be seen as decomposition of the empirical data matrix (OCT image) X into low-rank matrix L and sparse matrix S. The LRpSD problem [Eq. (2)] can be seen as a composition of two separate problems: the LRMA problem X=L+G that appears in many signal processing applications,19,36,4042 and sparsity constrained signal reconstruction corrupted with the AWGN, X=S+G.4346 Thus, estimation of the low-rank matrix L and sparse matrix S is expressed as the following optimization problem:

Eq. (3)

minL,Srank(L)+τS0subject to  X=L+S+G,
where ·0 counts the number of nonzero entries of S and τ>0 is a tuning parameter. The rank minimization problem is NP-hard. Minimization of the number of nonzero entries is an NP-hard problem as well. Thus, the optimization problem [Eq. (3)] is often replaced by convex relaxation29,31

Eq. (4)

minL,Siσi(L)+τS1subject to  X=L+S+G.
The first term is the 1-norm of the vector σ(L) of singular values of L, and it is known as the nuclear- or Schatten-1 norm of L.22,33 It represents convex relaxation of the rank minimization problem.32 The second term is the 1-norm of the matrix S and it represents convex relaxation of the S0 minimization problem.46 The optimization problem [Eq. (4)] is converted into the following optimization equation:

Eq. (5)

minL,S[Ψ(L,S)=12XLSF2+λiσi(L)+τS1],
where λ is a regularization constant that determines relative importance of the rank penalty term. The solution of the nuclear norm minimization problem, when S is fixed, is obtained directly using the singular value decomposition (SVD) of the matrix XS=UΣVT. It is given with

Eq. (6)

L^=Usoft(Σ,λ)VT,
where soft(Σ,λ) is the soft-thresholding function47 applied to the singular values of XS. The solution [Eq. (6)] is known as the “singular value thresholding” (SVT) method.48 Solution of the S1 minimization problem, when L is fixed, is obtained directly as

Eq. (7)

S^=soft(XL,τ),
where soft(XL,τ) is the soft-thresholding function applied entry-wise to the matrix XL. As emphasized in Refs. 20 and 49, the SVT method tends to underestimate the nonzero singular values. Thus, nuclear norm-based solutions of the low-rank approximation problem will exhibit decreased accuracy in estimation of the “clean” OCT image L.

2.2.

Nonconvex Regularization for Low-Rank + Sparsity Decomposition

Several recent studies have emphasized the benefit of nonconvex penalty functions compared to the nuclear norm for the estimation of the singular values.19,31,3437 In particular, it has been presented in Ref. 19 how nonconvex regularization, which promotes more sparse approximation of singular values,40 can be combined into a convex optimization problem related to the estimation of the low-rank matrices. The LRMA problem is formulated as19

Eq. (8)

minL{Ψ(L)=12XLF2+λi=1kφ[σi(L);a]},
where k=min(I1,I2) and φ is the sparsity-inducing regularizer, possibly nonconvex. To estimate nonzero singular values more accurately and induce sparsity more effectively than the nuclear norm, nonconvex penalty functions parameterized by the parameter a0 are used.19,50 Assumption 1 in Ref. 19 defines conditions that φ:RR have to satisfy. An example of a penalty function φ satisfying assumption 1 is the partly quadratic penalty function19,51,52

Eq. (9)

φ(x;a){|x|a2x2,|x|1a12a,|x|1a.

According to definition 1 in Ref. 19, see also Ref. 53, the proximal operator of φ, θ:RR, is defined as

Eq. (10)

θ(y;λ,a)argminxR[12(yx)2+λφ(x;a)].
If 0a<1/λ, then θ is continuous nonlinear threshold function with threshold value λ, i.e.,

Eq. (11)

θ(y;λ,a)=0,  |y|<λ.

The proximal operator of the partly quadratic penalty [Eq. (9)] is the firm threshold function defined as54

Eq. (12)

θ(y;λ,a)min{|y|,max[(|y|λ)/(1aλ),0]}sign(y).

In case of matrix X, notation θ(X;λ,a) implies that the proximal operator is applied element-wise to X. In addition to partly quadratic penalty function [Eq. (9)], other functions, such as logarithmic function,50 can be used as nonconvex penalty function in Eq. (8). However, the partly quadratic function [Eq. (9)] yielded best experimental results presented in Sec. 3. Thus, we shall not elaborate further on nonconvex penalty functions. According to theorem 2 in Ref. 19 the LRMA problem [Eq. (8)] has globally optimal solution

Eq. (13)

L^=U·θ(Σ;λ,a)·VT,
where the threshold function θ is associated with the nonconvex penalty function φ. We now use this result to obtain a more accurate solution to the problem [Eq. (3)]. In this regard, we substitute nuclear norm term in Eq. (5) with the nonconvex penalty from Eq. (8) and that yields

Eq. (14)

minL,S{Ψ(L,S)=12XLSF2+λi=1kφ[σi(L);a]+τS1}.
The optimization problem [Eq. (14)] can be seen as a special case of the more general linearly constrained convex program55

Eq. (15)

minL,Sf(L)+g(S)subject to  A(L)+B(S)=L+S=XG,
where f(L)=λi=1kφ[σi(L);a] and g(S)=τS1. When A(L) and B(S) in Eq. (15) are identity operators, i.e., A(L)=L and B(S)=S, the problem [Eq. (15)] can be solved by the alternating direction method of multipliers (ADMM). For this purpose, the augmented Lagrangian function is formulated56

Eq. (16)

L(L,S,Λ)=λi=1kφ[σi(L);a]+τS1+Λ,L+SX+β2L+SXF2,
where Λ is the matrix of Lagrange multipliers and β is the penalty parameter. The ADMM decomposes the minimization of L with respect to (L,S) into two subproblems that minimize with respect to L and S, respectively55

Eq. (17)

Lt=argminLL(L,St1,Λt1)=argminLλi=1kφ[σi(L);a]+β2L+St1X+Λt1βF2,

Eq. (18)

St=argminSL(Lt,S,Λt1)=argminSτS1+β2Lt+SX+Λt1βF2,

Eq. (19)

Λt=Λt1+β[Lt+StX],
where in Eqs. (17) to (19) t stands for an iteration index. The subproblem [Eq. (17)] is actually the LRMA problem [Eq. (8)] that admits optimal closed form solution given by Eq. (13)

Eq. (20)

L^=U·θ(Σ;λ,a)·VT,
for the SVD of

Eq. (21)

XSt1Λt1β=UΣVT.

The subproblem [Eq. (18)] admits the optimal solution in the form of Eq. (7)

Eq. (22)

S^t=soft(XLtΛt1β,τ).

We name the proposed algorithm enhanced LRpSD (ELRpSD) algorithm. It is summarized in Algorithm 1.

Algorithm 1

The ELRpSD algorithm.

Input: logarithm of acquired OCT image XRI1×I2 with the size of I1×I2  pixels, regularization constant τ related to speckle term S in Eqs. (14)/(16); regularization constant λ related to low-rank approximation term L in Eqs. (14)/(16). Suggested values: τ=0.1; λ=5.
Suggested value for the penalty parameter β in Eq. (16): β=1. Suggested value for constant a in Eqs. (9) to (13): a=0.6/λ.
1. L(0)=X; S(0)=0; Λ(0)=0; t=1.
2. while not converge do
3. Execute SVD Eq. (21).
4. Update L using Eq. (20).
5. Update S using Eq. (22).
6. Update Λ using Eq. (19).
7. tt+1
8. end while
Output:LL(t+1), SS(t+1).

2.3.

Performance Measure

To quantify the performance of speckle reduction algorithms, appropriate measures have to be defined. In the case of an OCT image, the most commonly used figure of merit is CNR.7,9 It corresponds to the inverse of the speckle fluctuation and it is defined as CNR=μl(X)/σl(X), where μl(X) and σl(X), respectively, correspond to the mean and standard deviation in some selected homogeneous part of the image X. Experimental results reported in Sec. 3 were estimated in the region that corresponds with the top most layer in the OCT image of a retina, which is indicated in Fig. 1 by an arrow.5,57 Since the goal of postprocessing algorithms is not only to reduce speckle but also to preserve image resolution, contrast, and edge fidelity,7 we also estimate contrast, sharpness as well as SNR measures directly from the image. Sharpness is the attribute related to the preservation of fine details (edges) in an image. Contrast is defined as the ratio of the maximum and the minimum intensity of the entire image.58 It reflects the strength of the noise or modeling error term G. Up to some extent, it can be considered as an image quality measure that coincides with the SNR quality measure. Technical details on estimation of sharpness and contrast can be found in Refs. 28 and 58. We estimated sharpness in the entire retinal region from the first (top most) to the tenth (bottom most) layer. Contrast was estimated from the whole image. By following Ref. 18, global SNR value was estimated as SNR=10log[max(Xlin)2/σlin2], where Xlin is the OCT image on a linear intensity scale and σlin2, such that the noise variance was estimated on a region between top of the image and the topmost layer.

3.

Experiments and Results

3.1.

Algorithms for Comparison and Optical Coherence Tomography Image Acquisition

We compare the proposed ELRpSD algorithm with the “Go Decomposition” (GoDec) algorithm,20 which solves the optimization problem [Eq. (4)] with the S1 term replaced with S0 and τ standing for a fraction of the nonzero coefficients of S relative to the overall number of coefficients, which is I1×I2; the semisoft version of the GoDec algorithm (SSGoDec), which solves the problem [Eq. (4)]; the rank N soft constraint (RNSC) for RPCA algorithm,21 which is using partial sum of singular values for more accurate approximation, compared with nuclear norm, of a target rank value; that are used to suppress sparsely distributed noise such as speckle. The MATLAB code for the GoDec and SSGoDec algorithms has been downloaded from Ref. 59. The MATLAB code for the RNSC algorithm has been downloaded from Ref. 60. The MATLAB code for the 2-D bilateral filtering algorithm has been downloaded from Refs. 61 and 62. For 2-D median filtering, the MATLAB function medfilt2 has been used. After computational experiments, we have selected for the GoDec algorithm the bound on S0 to be 0.1×(I1×I2). For the SSGoDec algorithm, the sparsity regularization constant has been selected to be τ=0.1. For the ELRpSD algorithm in Eq. (14), respectively, Eqs. (16) to (21), the parameter values were the following: λ=5, τ=0.1, and β=1. The speckle reduction algorithms were comparatively tested on 10 three-dimensional (3-D) macular-centered OCT images of normal eyes acquired with the Topcon 3-D OCT-1000 scanner. Each 3-D OCT image was comprised of 64 2-D scans with the size of 480×512  pixels. These images have been used previously for the study for optical intensity analysis in Ref. 57, where they were segmented into 10 retina layers. We estimated CNR, contrast, SNR, and sharpness values from the original image as well as from the images with reduced speckle. The images were analyzed with software written in the MATLAB (the MathWorks Inc., Natick, Massachusetts) script language on PC with Intel i7 CPU with the clock speed of 2.2 GHz and 16 GB of RAM.

3.2.

Comparative Results

Here, we present the results of the comparative performance analysis among the ELRpSD, GoDec, SSGoDec, RNSC, 2-D bilateral filtering, and 2-D median filtering algorithms. Parameters of bilateral filters have been tuned to yield approximately the same CNR value (the same level of speckle reduction) as the proposed ELRpSD algorithm. The median filtering has been used with the window of the size 3×3  pixels and that yields slightly higher CNR value than the one achieved by the proposed ELRpSD method. The size of the window can be increased to improve the edge fidelity but that would decrease the CNR, contrast, and SNR values. The algorithms were applied to each 2-D OCT scan separately. CNR, SNR, contrast, and sharpness were estimated from each enhanced 2-D scan and the reported values were averaged over 64 scans for each 3-D OCT image. Afterward, they were averaged further over 10 3-D OCT images. Average computation time of the ELRpSD, GoDec, SSGoDec, RNSC, 2-D bilateral filtering, and 2-D median filtering algorithms per one 2-D OCT scan is, respectively, given as: 4.51, 11.56, 9.63, 3.40, 11.28, and 0.22 s. Note, however, that unlike the ELRpSD algorithm, the GoDec, SSGoDec, and RNSC algorithms required a priori information of a targeted rank value. Since the true value of a rank was not known a priori, the GoDec, SSGoDec, and RNSC algorithms had to be run multiple times for the rank value within selected range. Clearly, huge computational complexity makes them not competitive in comparison with the ELRpSD algorithm. We show in Figs. 2, 3, 4, and 5 in respective order error bars of averaged values of CNR, relative SNR, relative contrast, and relative sharpness estimated from 10 3-D OCT images. Thereby, means and standard deviations of relative SNR values were obtained as follows:

Eq. (23)

Relative mean SNR(%)=100×mean(SNR of enhanced image)mean(SNR of original image)mean(SNR of original image),

Eq. (24)

Relative standard deviation SNR(%)=100×std(SNR of enhanced imageSNR of original image)mean(SNR of original image).

Fig. 2

Average CNR values (mean ± standard deviation) estimated from 10 3-D OCT images. The ELRpSD, bilateral filtering, and median filtering do not require a priori information on targeted rank value. Thus, their CNR estimates are shown as straight lines.

JBO_21_7_076008_f002.png

Fig. 3

Values of SNR in percentage (mean ± standard deviation) estimated from enhanced 3-D OCT images relatively to the SNR of original images. The ELRpSD, bilateral filtering, and median filtering do not require a priori information on targeted rank value. Thus, their SNR estimates are shown as straight lines.

JBO_21_7_076008_f003.png

Fig. 4

Values of contrast in percentage (mean ± standard deviation) estimated from enhanced 3-D OCT images relatively to the contrast of original images. The ELRpSD, bilateral filtering, and median filtering do not require a priori information on targeted rank value. Thus, their estimates of relative contrast value are shown as straight lines.

JBO_21_7_076008_f004.png

Relative values of contrast and sharpness are defined analogously. It can be seen that values of CNR, SNR, and contrast decrease when rank is increased. That is because with the increase of rank, influence of noise, which corresponds to small singular values, is increased as well. However, as seen from Fig. 5, the sharpness, which measures the edge fidelity, is increased with the increase of rank. That is because when low-rank approximation L is based on too few singular values, details important for the preservation of edges are lost. That is why the conflicting requirement of having the high CNR, SNR, and contrast values on one side and good edge fidelity on another side makes it difficult to select targeted value of the rank a priori. In this regard, bilateral filtering and median filtering suffer from the same problem. As can be seen from Fig. 2, bilateral filtering achieved the same value of CNR and higher relative SNR value in comparison with the proposed ELRpSD method, but it yielded reduced relative sharpness in comparison with the relative sharpness achieved by the ELRpSD method. The median filtering yielded high values of CNR, relative SNR, and contrast but destroyed edge fidelity compared with the original image. In both cases, bilateral filtering and median filtering reduction of the edge fidelity is caused by the blurring effect when the spatial bandwidth of the filter becomes too narrow and that is necessary to achieve higher values of CNR, SNR, and contrast. Thus, capability of the proposed ELRpSD method to estimate the rank value directly from the image is very valuable. As can be seen, it achieves the highest value of relative sharpness (the best edge fidelity) compared with other algorithms and, in comparison with the original images, also yields increased values of the CNR, relative SNR, and relative contrast. The GoDec, SSGoDec, and RNSC algorithms achieve comparable value of sharpness with the value of rank equal to 35. At this value of rank, GoDec and SSGoDec have slightly better value of CNR and contrast than ELRpSD, while the RNSC is still worse. Thus, presumably the GoDec and SSGoDec could be used for the speckle reduction on the existing OCT scanner with a rank set to a predefined value of 35. However, if the OCT images are to be acquired on a different scanner the GoDec, SSGoDec, and RNSC algorithms would have to be “calibrated” again. To validate stability of the proposed ELRpSD method, we show in Fig. 6 relative values of CNR, SNR, contrast, and sharpness estimated for each of 10 3-D OCT image separately. As can be seen variations of estimated values are within few percentages.

Fig. 5

Values of sharpness in percentage (mean ± standard deviation) estimated from enhanced 3-D OCT images relatively to the sharpness of original images. The ELRpSD, bilateral filtering, and median filtering do not require a priori information on targeted rank value. Thus, their estimates of relative sharpness value are shown as straight lines.

JBO_21_7_076008_f005.png

Fig. 6

Values of CNR, SNR, contrast, and sharpness in percentage (mean ± standard deviation) estimated from ELRpSD enhanced 3-D OCT images relatively to the corresponding values in original 3-D OCT images.

JBO_21_7_076008_f006.png

4.

Discussion

Large contrast and granular appearance of speckle with an OCT image of biological specimens reduce contrast and make boundaries between constitutive tissues more difficult to resolve. Thus, speckle stands for a major obstacle in quantitative OCT image analysis. Since speckle has a dual role as a source of noise and as a carrier of information about tissue microstructure its complete reduction is not desirable. Hence, speckle reduction is a peculiar problem. In particular, it is a challenge to increase the CNR value, which is used as a figure of merit in speckle reduction, and preserve image resolution, contrast and fidelity of edges. In this regard, we have proposed an approach to speckle reduction, which is based on decomposition of 2-D OCT scans into low-rank approximation of the “clean” image and sparse term that takes into account speckle. In particular, we proposed a method capable of estimating rank on data-driven or automatic ways directly from the experimental OCT image. Moreover, the method is using class of nonconvex regularization, which induces sparse approximation of singular values in the related LRMA problem. That, in turn, yields more accurate approximation of a rank than what is achieved by the more often used approximations based on nuclear norm. As a final result, the proposed method yields the low-rank approximation of the original OCT images with simultaneously increased values of CNR, SNR, sharpness, and contrast. That makes the proposed method suitable for speckle reduction in OCT images acquired at different scanners.

5.

Conclusion

We have developed a method for the speckle reduction in OCT images and named it the ELRpSD algorithm. The method, which is applied on individual 2-D OCT scans, was tested on 10 3-D OCT images comprised of 64 scans each. It was able to simultaneously increase, relative to the original OCT images, values of CNR, SNR, contrast, and sharpness (improved fidelity of edges). In particular, the relative improvement, averaged over 10 3-D OCT images, of the CNR, SNR, contrast and sharpness was in respective order 14.71%, 23.08%, 24.54%, and 14.61%. Therefore, we conclude that the ELRpSD method can be used as a preprocessing method for speckle reduction to enable more accurate quantitative analysis of OCT images.

Acknowledgments

This work was supported in part through a bilateral Chinese-Croatian Grant “Decomposition and contrast enhancement of PET/CT and OCT images by means of nonlinear sparse component analysis.” It was also supported in part by the National Basic Research Program of China (973 Program) under Grant No. 2014CB748600, in part by the National Natural Science Foundation of China under Grant Nos. 81371629, 61401293, 61401294, and 81401472, and in part by the Natural Science Foundation of the Jiangsu Province under Grant No. BK20140052.

References

1. 

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

2. 

W. F. Cheong, S. A. Prahl and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron., 26 (12), 2166 –2185 (1990). http://dx.doi.org/10.1109/3.64354 Google Scholar

3. 

B. C. Wilson and S. L. Jacques, “Optical reflectance and transmittance of tissues-principles and applications,” IEEE J. Quantum Electron., 26 (12), 2186 –2199 (1990). http://dx.doi.org/10.1109/3.64355 Google Scholar

4. 

G. R. Wilkins, O. M. Houghton and A. L. Oldenburg, “Automated segmentation of intraretinal cystoid fluid in optical coherence tomography,” IEEE Trans. Biomed. Eng., 59 (4), 1109 –1114 (2012). http://dx.doi.org/10.1109/TBME.2012.2184759 IEBEAX 0018-9294 Google Scholar

5. 

H. Chen et al., “Quantitative analysis of retinal layers’ optical intensities on 3D optical coherence tomography for central retinal artery occlusion,” Sci. Rep., 5 9269 (2015). http://dx.doi.org/10.1038/srep09269 Google Scholar

6. 

C. Xu et al., “Characterization of atherosclerosis plaques by measuring both backscattering and attenuation coefficients in optical coherence tomography,” J. Biomed. Opt., 13 (3), 034003 (2008). http://dx.doi.org/10.1117/1.2927464 JBOPFO 1083-3668 Google Scholar

7. 

G. van Soest et al., “Frequency domain multiplexing for speckle reduction in optical coherence tomography,” J. Biomed. Opt., 17 (7), 076018 (2012). http://dx.doi.org/10.1117/1.JBO.17.7.076019 JBOPFO 1083-3668 Google Scholar

8. 

X. Zhang et al., “Spiking cortical model-based nonlocal means method for speckle reduction in optical coherence tomography images,” J. Biomed. Opt., 19 (6), 066005 (2014). http://dx.doi.org/10.1117/1.JBO.19.6.066005 JBOPFO 1083-3668 Google Scholar

9. 

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

10. 

J. W. Goodman, “Statistical properties of laser speckle patterns,” Laser Speckle and Related Phenomena, Springer Verlag, Berlin (1984). 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. Baghaie, Z. Yu and R. M. D’souza, “State-of-the-art in retinal optical coherence tomography image analysis,” Quantum Imaging Med. Surg., 5 (4), 603 –617 (2015). http://dx.doi.org/10.3978/j.issn.2223-4292.2015.07.02 Google Scholar

13. 

E. Desjardins et al., “Estimation of the scattering coefficients of turbid media using angle-resolved optical frequency-domain imaging,” Opt. Lett., 32 (11), 1560 –1562 (2007). http://dx.doi.org/10.1364/OL.32.001560 OPLEDP 0146-9592 Google Scholar

14. 

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

15. 

J. Kim et al., “Optical coherence tomography speckle reduction by a partially spatially coherent sources,” J. Biomed. Opt., 10 (6), 064034 (2005). http://dx.doi.org/10.1117/1.2138031 JBOPFO 1083-3668 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. 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.001902 JOAOD6 0740-3232 Google Scholar

18. 

D. C. Adler, T. H. Ko and J. G. Fujimoto, “Speckle reduction in optical coherence tomography images by use of a spatially adaptive wavelet filter,” Opt. Lett., 29 (24), 2878 –2880 (2004). http://dx.doi.org/10.1364/OL.29.002878 OPLEDP 0146-9592 Google Scholar

19. 

A. Parekh and I. W. Selesnick, “Enhanced low-rank matrix approximation,” IEEE Signal Proc. Lett., 23 (4), 493 –497 (2016). http://dx.doi.org/10.1109/LSP.2016.2535227 Google Scholar

20. 

T. Zhou and D. Tao, “GoDec: randomized low-rank & sparse matrix decomposition in noisy case,” in Proc. 28th Int. Conf. Machine Learning (ICML), 33 –40 (2011). Google Scholar

21. 

T. H. Oh et al., “Partial sum minimization of singular values in robust PCA: algorithms and applications,” IEEE Trans. Pattern Anal. Mach. Intell., 38 (4), 744 –758 (2016). http://dx.doi.org/10.1109/TPAMI.2015.2465956 ITPIDJ 0162-8828 Google Scholar

22. 

X. Zhang et al., “Schatten-q regularizer for low rank subspace clustering model,” Neurocomputing, 182 36 –47 (2016). http://dx.doi.org/10.1016/j.neucom.2015.12.009 NRCGEO 0925-2312 Google Scholar

23. 

R. Vidal and P. Favaro, “Low rank subspace clustering (LRSC),” Pattern Recogn. Lett., 43 47 –61 (2014). http://dx.doi.org/10.1016/j.patrec.2013.08.006 Google Scholar

24. 

V. M. Patel, H. V. Nguyen and R. Vidal, “Latent space sparse and low-rank subspace clustering,” IEEE J. Sel. Top. Sig. Proc., 9 (4), 691 –701 (2015). http://dx.doi.org/10.1109/JSTSP.2015.2402643 Google Scholar

25. 

G. Liu et al., “Robust recovery of subspace structures by low-rank representation,” IEEE Trans. Pattern Anal. Mach. Intell., 35 (1), 171 –184 (2013). http://dx.doi.org/10.1109/TPAMI.2012.88 ITPIDJ 0162-8828 Google Scholar

26. 

Y. Hu et al., “Anomaly detection in hyperspectral images based on low-rank and sparse representations,” IEEE Trans. Geosci. Remote Sens., 54 (4), 1990 –2000 (2016). http://dx.doi.org/10.1109/TGRS.2015.2493201 Google Scholar

27. 

C. Li et al., “Hyperspectral image denoising using the low-rank tensor recovery,” J. Opt. Soc. Am. A., 32 (9), 1604 –1612 (2015). http://dx.doi.org/10.1364/JOSAA.32.001604 JOAOD6 0740-3232 Google Scholar

28. 

I. Kopriva et al., “Offset-sparsity decomposition for automated enhancement of color microscopic image of stained specimen in histopathology,” J. Biomed. Opt., 20 (7), 076012 (2015). http://dx.doi.org/10.1117/1.JBO.20.7.076012 JBOPFO 1083-3668 Google Scholar

29. 

E. J. Candès et al., “Robust principal component analysis?,” J. ACM, 58 11 (2011). http://dx.doi.org/10.1145/1970392.1970395 Google Scholar

30. 

V. Chandrasekaran et al., “Rank-sparsity incoherence for matrix decomposition,” SIAM J. Opt., 21 572 –596 (2011). http://dx.doi.org/10.1137/090761793 Google Scholar

31. 

R. Chartrand, “Nonconvex splitting for regularized low-rank + sparse decomposition,” IEEE Trans. Signal Proc., 60 (11), 5810 –5819 (2012). http://dx.doi.org/10.1109/TSP.2012.2208955 Google Scholar

32. 

B. Recht, M. Fazel and P. A. Parillo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” SIAM Rev., 52 (3), 471 –501 (2010). http://dx.doi.org/10.1137/070697835 SIREAD 0036-1445 Google Scholar

33. 

K. Mohan and M. Fazel, “Iterative reweighted algorithms for matrix rank minimization,” J. Mach. Learn. Res., 13 (1), 3441 –3473 (2012). Google Scholar

34. 

C. Lu et al., “Generalized nonconvex nonsmooth low-rank minimization,” in Proc. IEEE Conf. Computer Vision and Pattern Recognition, 4130 –4137 (2014). http://dx.doi.org/10.1109/CVPR.2014.256 Google Scholar

35. 

C. Lu et al., “Generalized singular value thresholding,” in Proc. AAAI Conf. Artificial Intelligence, 1805 –1811 (2015). Google Scholar

36. 

C. Liu et al., “Nonconvex nonsmooth low rank minimization via iteratively reweighted nuclear norm,” IEEE Trans. Image Proc., 25 (2), 829 –839 (2016). http://dx.doi.org/10.1109/TIP.2015.2511584 IIPRE4 1057-7149 Google Scholar

37. 

P.-Y. Chen and I. Selesnick, “Group-sparse signal denoising: non-convex regularization, convex optimization,” IEEE Trans. Signal Proc., 62 (13), 3464 –3476 (2014). http://dx.doi.org/10.1109/TSP.2014.2329274 Google Scholar

38. 

A. Baghaie, R. M. S’souza and Z. Yu, “Sparse and low rank decomposition based batch image alignment for speckle reduction of retinal OCT image,” in 2015 IEEE Int. Symp. on Biomedical Imaging, 226 –230 (2015). http://dx.doi.org/10.1109/ISBI.2015.7163855 Google Scholar

39. 

F. Luan and Y. Wu, “Application of RPCA in optical coherence tomography for speckle noise reduction,” Laser Phys. Lett., 10 035603 (2013). http://dx.doi.org/10.1088/1612-2011/10/3/035603 1612-2011 Google Scholar

40. 

M. Fazel et al., “Hankel matrix rank minimization with applications to system identification and realization,” SIAM J. Matrix Anal. Appl., 34 (3), 946 –977 (2013). http://dx.doi.org/10.1137/110853996 SJMAEL 0895-4798 Google Scholar

41. 

H. M. Nguyen et al., “Denoising MR spectroscopic imaging data with low-rank approximation,” IEEE Trans. Biomed. Eng., 60 (1), 78 –89 (2013). http://dx.doi.org/10.1109/TBME.2012.2223466 IEBEAX 0018-9294 Google Scholar

42. 

R. R. Nadakuditi, “Optshrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage,” IEEE Trans. Inf. Theory, 60 (5), 3002 –3018 (2014). http://dx.doi.org/10.1109/TIT.2014.2311661 IETTAW 0018-9448 Google Scholar

43. 

I. Markovsky, “Structured low-rank approximation and its applications,” Automatica, 44 (4), 891 –909 (2008). http://dx.doi.org/10.1016/j.automatica.2007.09.011 ATCAA9 0005-1098 Google Scholar

44. 

I. Daubechies, M. Defrise and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparisty constraint,” Commun. Pure Appl. Math., 57 (11), 1413 –1457 (2004). http://dx.doi.org/10.1002/cpa.20042 Google Scholar

45. 

F. Luisier, T. Blu and M. Unser, “A new sure approach to image denoising: interscale orthonormal wavelet thresholding,” IEEE Trans. Image Process., 16 (3), 593 –606 (2007). http://dx.doi.org/10.1109/TIP.2007.891064 IIPRE4 1057-7149 Google Scholar

46. 

D. L. Donoho and M. Elad, “Optimally sparse representation in general (non-orthogonal) dictionaries via l1 minimization,” Proc. Natl. Acad. Sci. U.S.A., 100 2197 –2202 (2003). http://dx.doi.org/10.1073/pnas.0437847100 Google Scholar

47. 

D. Donoho, “De-noising by soft-thresholding,” IEEE Trans. Inf. Theory, 41 (3), 613 –627 (1995). http://dx.doi.org/10.1109/18.382009 IETTAW 0018-9448 Google Scholar

48. 

J.-F. Cai et al., “A singular value thresholding algorithm for matrix completion,” SIAM J. Optim., 20 (4), 1956 –1982 (2010). http://dx.doi.org/10.1137/080738970 SJOPE8 1095-7189 Google Scholar

49. 

Q. Liu et al., “A truncated nuclear norm regularization method based on weighted residual error for matrix completion,” IEEE Trans. Image Proc., 25 (1), 316 –330 (2016). http://dx.doi.org/10.1109/TIP.2015.2503238 IIPRE4 1057-7149 Google Scholar

50. 

W. Selesnick and I. Bayram, “Sparse signal estimation by maximally sparse convex optimization,” IEEE Trans. Signal Proc., 62 (5), 1078 –1092 (2014). http://dx.doi.org/10.1109/TSP.2014.2298839 Google Scholar

51. 

I. Bayram, “On the convergence of the iterative shrinkage/thresholding algorithm with a weakly convex penalty,” IEEE Trans. Signal Proc., 64 (6), 1597 –1608 (2016). http://dx.doi.org/10.1109/TSP.2015.2502551 Google Scholar

52. 

C.-H. Zhang, “Nearly unbiased variable selection under minimax concave penalty,” Ann. Stat., 38 (2), 894 –942 (2010). http://dx.doi.org/10.1214/09-AOS729 ASTSC7 0090-5364 Google Scholar

53. 

P. L. Combettes and J.-C. Pesquet, “Proximal thresholding algorithm for minimization over orthonormal bases,” SIAM J. Optim., 18 (4), 1351 –1376 (2007). http://dx.doi.org/10.1137/060669498 SJOPE8 1095-7189 Google Scholar

54. 

Y. Gao and A. G. Bruce, “Waveshrink with firm shrinkage,” Stat. Sin., 7 (4), 855 –874 (1997). STSNEO Google Scholar

55. 

Z. Lin, R. Liu and Z. Su, “Linearized alternating direction method with adaptive penalty for low-rank representation,” in Proc. Advances in Neural Information Processing Systems (NIPS) Conf., 612 –620 (2011). Google Scholar

56. 

S. Boyd et al., “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn., 3 (1), 1 –122 (2010). http://dx.doi.org/10.1561/2200000016 Google Scholar

57. 

X. Chen et al., “Quantitative analysis of retinal layers’ optical intensities on 3D optical coherence tomography,” Invest. Ophthalmol. Visual Sci., 54 (10), 6846 –6851 (2013). http://dx.doi.org/10.1167/iovs.13-12062 Google Scholar

58. 

K. Panetta, C. Gao and S. Agaian, “No reference color image contrast and quality measure,” IEEE Trans. Consum. Electron., 59 (3), 643 –651 (2013). http://dx.doi.org/10.1109/TCE.2013.6626251 ITCEDA 0098-3063 Google Scholar

59. 

T. Zhou and D. Taot, “MATLAB code for the GoDec and SSGoDec algorithms,” (2016) https://sites.google.com/site/godecomposition/home April ). 2016). Google Scholar

60. 

T. H. Oh, “MATLAB code for the RNSC algorithm,” (2016) http://rcv.kaist.ac.kr/v2/bbs/member_detail.php?mb_id=thoo April ). 2016). Google Scholar

61. 

D. Lahman, “MATLAB code for 2D bilateral filtering algorithm,” (2016) http://www.mathworks.com/matlabcentral/fileexchange/12191-bilateral-filtering April ). 2016). Google Scholar

62. 

“MATLAB code for proposed LRpSD algorithm,” www.mipav.net/English/research/research.html Google Scholar

Biography

Ivica Kopriva received his PhD in electrical engineering from the University of Zagreb, Croatia, in 1998, with the topic on blind source separation. He is a senior scientist at the Ruđer Bošković Institute, Zagreb, Croatia. He has coauthored over 40 papers in internationally recognized journals, one book, and holds three US patents. His current research is focused on structured decompositions of empirical data with applications in imaging, spectroscopy, and variable selection.

Fei Shi received her PhD in electrical engineering from Polytechnic University, Brooklyn, New York, in 2006. She is an assistant professor at the Soochow University, Suzhou, China. She has coauthored over 20 papers in internationally recognized journals and conferences. Her current research is focused on medical image processing and analysis.

Xinjian Chen received his PhD from the Chinese Academy of Sciences in 2006. He is a distinguished professor at the Soochow University, Suzhou, China. He conducted postdoctoral research at the University of Pennsylvania, National Institutes of Health, and the University of Iowa, USA, from 2008 to 2012. He has published over 80 top international journal and conference papers. He has also been granted four patents.

© 2016 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2016/$25.00 © 2016 SPIE
Ivica Kopriva, Fei Shi, and Xinjian Chen "Enhanced low-rank + sparsity decomposition for speckle reduction in optical coherence tomography," Journal of Biomedical Optics 21(7), 076008 (18 July 2016). https://doi.org/10.1117/1.JBO.21.7.076008
Published: 18 July 2016
Lens.org Logo
CITATIONS
Cited by 29 scholarly publications.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Optical coherence tomography

Speckle

Signal to noise ratio

Digital filtering

3D image processing

Detection and tracking algorithms

Image enhancement

Back to Top