SPECKLE REDUCTION IN OPTICAL COHERENCE TOMOGRAPHY IMAGES

Abstract

An optical coherence tomography (OCT) image composed of a plurality of A-scans of a structure is analysed by defining, for each A-scan, a set of neighbouring A-scans surrounding the A-slices scan. Following an optional de-noising step, the neighbouring A-scans are aligned in the imaging direction, then a matrix X is formed from the aligned A-scans, and matrix completion is performed to obtain a reduced speckle noise image.

Claims

1. A computer-implemented method of enhancing an optical coherence tomography (OCT) image composed of a plurality of A-scans, each A-scan being indicative of a structure imaged in a predefined direction, the A-scans being grouped as one or more B-scans, the method comprising, for each A-scan: (a) defining a respective neighbourhood surrounding the A-scan; (b) aligning the A-scan in the pre-defined direction with other A-scans in the respective neighbourhood; (c) using the aligned A-scans to form a first matrix X, respective columns of the first matrix being the aligned A-scans; and (d) seeking a second matrix which the minimizes a cost function indicative of the difference between the first matrix X and the second matrix, the second matrix being constrained to obey a complexity constraint.

2. A method according to claim 1 in which the complexity constraint is that the second matrix is constrained to have a rank no higher than a predefined value.

3. A method according to claim 1 in which the complexity constraint is that, representing the second matrix as the sum of a low rank matrix L and a sparse matrix S, the complexity constraints being that the low rank matrix is constrained to have a rank less than a first predetermined value, and the sparse matrix S is constrained to have a number of non-zero elements less than a second predetermined value.

4. A method according to any preceding claim comprising a pre-step of noise reduction in the plurality of A-scans.

5. A method according to claim 4 in which the noise-reduction is a speckle reduction anisotropic diffusion algorithm.

6. A method according to any preceding claim wherein step (c) is performed by: defining a first A-scan group consisting of one or more A-scans including the A-scan; and for each neighbouring A-scan: (i) defining a second respective A-scan group, the second respective A-scan group consisting of one or more A-scans including the neighbouring A-scan; (ii) seeking a respective translation value which, if the one or more A-scans of the respective second A-scan group are translated in the imaging direction by the translation value, minimises a difference between the first A-scan group and the second A-scan group; and (iii) translating the neighbouring A-scan by the respective translation value.

7. A method according to claim 6 in which step (c) is performed by a block matching algorithm.

8. A method according to claim 7 wherein the block matching algorithm is the diamond search algorithm.

9. A computer system comprising a processor and a memory device configured to store program instructions operative, when performed by the processor, to cause the processor to perform a method according to any of claims 1 to 8.

10. A computer program product, such as a tangible recording medium, storing program instructions operative, when performed by a processor, to cause the processor to perform a method according to any of claims 1 to 8.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0023] An embodiment of the invention will now be described for the sake of example only with reference to the following drawings, in which:

[0024] FIG. 1 is a flow chart showing a method which is an embodiment of the invention;

[0025] FIG. 2 illustrates OCT images of a three-dimensional volume, the OCT images being N B-scans B.sub.1, B.sub.2, . . . , B.sub.N;

[0026] FIG. 3 is composed of FIGS. 3(a) and 3(b), which respectively illustrate two methods of defining a neighbourhood of an A-scan A.sub.i,j as (a) a rectangular neighborhood and (b) a circular/elliptical neighborhood;

[0027] FIG. 4 illustrates an A-scan alignment step of the method of FIG. 1;

[0028] FIG. 5 shows results for image reconstruction using a first version of the method of FIG. 1;

[0029] FIG. 6 shows results for image reconstruction using a second version of the method of FIG. 1;

[0030] FIG. 7 is a box plot of CNR comparing three methods of speckle reduction; and

[0031] FIG. 8 is composed of FIGS. 8(a)-(f), in which FIGS. 8(a) and (d) show raw OCT data, FIGS. 8(b) and 8(e) show a moving average of the raw data, and FIGS. 8(c) and 8(f) show data obtained by the method of FIG. 1.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0032] Referring to FIG. 1, a method according to the invention is illustrated. The input to the method is an OCT volume, of the kind which can be produced by a conventional OCT method. The OCT volume is made up of N B-scans (also referred to in the art as “frames” or “slices”), B.sub.i for i=1, 2, . . . N, as illustrated in FIG. 2. Each B-scan B.sub.i consists of MA-scans, A.sub.i,j=1, 2, . . . , M, which are each one-dimensional scans formed by imaging the object with a one-dimensional pencil of light extending in the direction which is vertical in FIG. 1. In other words, the OCT volume consists of a set of A-scan lines A.sub.i,j, i=1, 2, . . . , N, j=1, 2, . . . , M, where N is the number of B-scans in the volume and M is the number of A-scans in each B-scan.

[0033] A basic concept underlying the method of FIG. 1 is to formulate A.sub.i,j as a sum of its underlying clean image l.sub.i,j, and noise g.sub.i,j, i.e.,

A.sub.i,j=l.sub.i,j+g.sub.i,j

The objective is to reconstruct l.sub.i,j from A.sub.i,j. It is impossible to solve this problem without some assumptions about l.sub.i,j and n.sub.i,j, but the embodiment is motivated by the assumption that l.sub.i,j is very similar within a small neighborhood. Based on this assumption, we model the neighboring A-scans using similar model.

A.sub.i+p,j+q=l.sub.i+p,j+q+g.sub.i+p,j+q,

where (p, q)εΩ and Ω denotes a “neighborhood” (or “neighbor region”), that is a set of A-scans which are neighbors of, and surround, the A-scan. Thus, the set of A-scans form a two-dimensional array as viewed in the imaging direction.

[0034] There are many ways to define the neighbor region. A first possible definition is as the rectangular region {(p,q)|−n≦p≦n,−m≦q≦m} as illustrated in FIG. 3(a). An alternative possible definition is as the elliptical region

[00001]$\{\left(p,q\right).Math..Math.\frac{p^2}{n^2}+\frac{q^2}{{m.Math.}^2}\le 1\}$

as illustrated in FIG. 3(b). If n=m then the neighborhoods are respectively square and circular. In the experiments reported below, the neighborhoods were defined as the rectangular neighborhoods of FIG. 3(a), using n=1 and m=2.

[0035] One way to implement this concept would be to define X as the matrix formed by stacking all the A-scans from Ω as columns, i.e., X=[A.sub.i−n,j−m, . . . , A.sub.i+n,j+m], L=[l.sub.i−n,j−m, . . . , l.sub.i+n,j+m], G=[g.sub.i−n,j−m, . . . , g.sub.i+n,j+m]. Each column in L may be assumed to be very similar, so the rank of L is expected to be low, and accordingly the method can use low rank matrix completion to solve L from given X.

[0036] Unfortunately, applying the above technique directly to raw OCT volume may lead to some blur of the A-scans due to content changes among neighboring A-scans. In addition, the inevitable eye movement also leads to some misalignment. For this reason, FIG. 1 includes a step 2 of aligning the A-scans, before the A-scan reconstruction by matrix completion is performed in step 3.

[0037] Furthermore, the original OCT A-scans are affected by large speckle noise, which make it challenging job to align these A-scans in step 2. Therefore, a pre-processing step 1 is used to reduce the noise before step 2 is carried out.

[0038] Many algorithms [1]-[7] can be used for the noise-reduction step 1. We use the SRAD algorithm [5]. SRAD is applied on every scan B.sub.i. Specifically, given an intensity image I.sub.0 (x, y), the output image I(x, y; t) is evolved according to the PDE:

[00002]$\{\begin{array}{c}\frac{\partial I}{\partial t}=\mathrm{div}\left[c\left(q\right).Math.\nabla I\right]\\ I\left(t=0\right)=I_0\end{array},$

where c(q) denotes the diffusion coefficient computed with

[00003]$c\left(q\right)=1/\left(1+\left[q^2\left(x,y;t\right)-q_0^2\left(t\right)\right]/\left[q_0^2\left(t\right).Math.\left(1+q_0^2\left(t\right)\right)\right]\right).Math..Math.\mathrm{or}$$c\left(q\right)=\exp .Math.\left\{-\left[q^2.Math.\left(x,y;t\right)-q_0^2\left(t\right)\right]/\left[q_0^2\left(t\right).Math.\left(1+q_0^2\left(t\right)\right)\right]\right\}.Math..Math.\mathrm{and}.Math.$$q\left(x,y;t\right)=\sqrt{\frac{\left(1/2\right).Math.{\left(.Math.\nabla I.Math./I\right)}^2-\left(1/16\right).Math.{\left({\nabla }^2.Math.I/I\right)}^2}{{\left[1+\left(1/4\right).Math.\left({\nabla }^2.Math.I/I\right)\right]}^2}.}$

In the experimental results presented below, the iteration was performed 50 times.

[0039] In the A-scan alignment step 2, neighboring A-scans A.sub.i+p,j+q, (p,q)εΩ are aligned to A.sub.i,j by a vertical translation A.sub.p,q (i.e., moving up or down) such that the difference in between is minimized, where Q denotes the locations of neighboring A-scans. From our experience, computing the difference from an A-scan alone may lead to a result dominated by local speckle noise. Therefore, we compute the difference between A.sub.i,j and A.sub.i+p,j+q using their neighboring A-scans as well.

[0040] In step 2, for each p and q we define a group (“window”) of A-scans centered at A.sub.i+p,j+q and A.sub.i,j, i.e., G.sub.i+p,j+q=[A.sub.i+p,j+q−w,A.sub.i+p,j+q−w+1, . . . , A.sub.i+p,j+q+w] and G.sub.i,j=[A.sub.i,j−w,A.sub.i,j−w+1, . . . ,A.sub.i,j+w] where w determines the window size, is used to compute the alignment. We find vertical translation Δ.sub.p,q on G.sub.i+p,j+q such that the difference between the translated G.sub.i+p,j+q and G.sub.i,j is minimized Many block matching algorithms can be used to find the alignment. In our numerical investigation of the embodiments, we use the diamond search strategy [8] to find the vertical translation because of its efficiency and easy implementation.

[0041] We now present two methods by which step 3 can be carried out, i.e. two versions of step 3. A first sub-step of both methods is to combine the aligned A-scans as the columns of the matrix X. That is, X=[A.sub.i−n,j−m, . . . , A.sub.i+n,j+m], where the content of the bracket is understood to be the aligned A-scans. Method 1: First, the embodiment decomposes X as a sum of low rank part and noise part, i.e.,

X=L+G,  (1)

where L is the low rank part with rank (L)≦r, and G is the noise. In the experiments reported below, we set r=2.

[0042] The above decomposition is solved by minimizing the decomposition error:

[00004]$\begin{array}{cc}\underset{L}{\min }.Math.{.Math.X-L.Math.}_F^2,s.t..Math..Math.\mathrm{rank}.Math..Math.\left(L\right)\le r& \left(2\right)\end{array}$

[0043] The notation ∥C∥.sub.F of any matrix X denotes the Frobenius norm of the matrix (i.e. the Euclidean norm of a vector formed by stacking the columns of the matrix).

[0044] The optimization problem above is solved by the bilateral random projections [9] algorithm. The algorithm is summarized as follows:

[0045] Algorithm 1 for solving L [0046] Generate a random matrix A.sub.1 [0047] A.sub.1εR.sup.D×r; where D is the number of pixels in an A-scan [0048] Compute {tilde over (L)}=(XX.sup.T).sup.qX; (we use q=2) [0049] Y.sub.1={tilde over (L)}A.sub.1,A.sub.2=Y.sub.1 [0050] Y.sub.2={tilde over (L)}.sup.TY.sub.1,Y.sub.1={tilde over (L)}Y.sub.2 [0051] Using standard methods, calculate QR decomposition Y.sub.1=Q.sub.1R.sub.1, Y.sub.2=Q.sub.2R.sub.2, where Q.sub.1 and Q.sub.2 are orthogonal matrices and R.sub.1 and R.sub.2 are upper triangular matrices;

L=Q.sub.1[R.sub.1(A.sub.2.sup.TY.sub.1).sup.−1R.sub.2.sup.T].sup.1/(2q+1)Q.sub.2.sup.T;

[0052] Method 2: In the second method, the embodiment decomposes the X as a sum of low rank part, sparse part and noise part, i.e.,

X=L+S+G,  (3)

where L is the low rank part with rank (L)≦r, S is the sparse part with card(S)≦k and G is the noise. In the experiments reported below, r was equal to 2, and k was set as 30% of the total elements of the matrix S.

[0053] The above decomposition is solved by minimizing the decomposition error:

[00005]$\hspace{1em}\begin{array}{cc}\underset{L,S}{\min }.Math.{.Math.X-L-S.Math.}_F^2,.Math.s.t..Math.\mathrm{rank}\left(L\right)\le r,.Math.\mathrm{card}\left(S\right)\le k.& \left(2\right)\end{array}$

[0054] The optimization problem above is solved by alternatively solving the following two sub-problems until convergence.

[00006]$\hspace{1em}\{\begin{array}{c}{L}_t=\arg .Math..Math.\underset{\mathrm{rank}\left(L\right)\le r}{\min }.Math.{.Math.X-L_{t-1}-S_{t-1}.Math.}_F^2\\ S_t=\arg .Math..Math.\underset{\mathrm{card}\left(S\right)\le k}{\min }.Math.{.Math.X-L_{t-1}-S_{t-1}.Math.}_F^2\end{array}$

[0055] The algorithm is summarized as follows:

Algorithm 2 for Solving L

[0056] Initialize t:=0; L.sub.0=X; S.sub.0=0;
While ∥X−L.sub.t−S.sub.t∥.sub.F.sup.2/∥X∥.sub.F.sup.2>ε, do
t:=t+1;
{tilde over (L)}=[(X−S.sub.t−1)(X−S.sub.t−1).sup.T].sup.q(X−S.sub.t−1);
Generate a random matrix A.sub.1εR.sup.D×r; here D is the number of pixels in an A-scan.
Compute Y.sub.1={tilde over (L)}A.sub.1,A.sub.2=Y.sub.1;
Y.sub.2={tilde over (L)}.sup.TY.sub.1, Y.sub.1={tilde over (L)}Y.sub.2
Calculate QR decomposition Y.sub.1=Q.sub.1R.sub.1, Y.sub.2=Q.sub.2R.sub.2;
If rank(A.sub.2.sup.TY.sub.1)<r then r=rank(A.sub.2.sup.TY.sub.1), go to the first step; end if;
L.sub.t=Q.sub.1[R.sub.1(A.sub.2.sup.TY.sub.1).sup.−1R.sub.2.sup.T].sup.1/(2q+1)Q.sub.2.sup.T;
S.sub.1=δ.sub.Ω(X−L.sub.t), Ω is the nonzero subset of the first k largest entries of |X−L.sub.t|; δ.sub.Ω is the projection of a matrix to Ω.
End while

[0057] Note that the term “column” is used here with a meaning which includes also a row of the matrix; in other words, in principle, X could be formed by stacking the aligned A-scans as the rows of X, but then the transpose of X would be applied before the equations are solved in methods 1 and 2.

[0058] It can be see that the first method is a special case of the second method with S=0. The main difference is that the second model separates the difference due to content change from the noise.

[0059] Although FIG. 1 shows no further processing steps, optionally various filters can be applied on the reconstructed A-scan volumes. Suitable filters include enhanced Lee, adaptive Wiener filters, etc. For example, an adaptive Wiener filter may be used. In addition, a threshold may be applied to assist the de-noising. Mathematically, we apply the following process:

[00007]$L\left(x,y\right)=\{\begin{array}{cc}{T}_{\mathrm{low}}& L\left(x,y\right)<T_{\mathrm{low}}\\ L\left(x,y\right)& T_{\mathrm{low}}\le L\left(x,y\right)\le T_{\mathrm{high}}\\ T_{\mathrm{high}}& L\left(x,y\right)>T_{\mathrm{high}}\end{array},$

[0060] Where T.sub.low and T.sub.high are two adaptively determined thresholds.

Results

[0061] We first present results of the use of the embodiment for the method 1 (i.e. the first variant of step 3). For convenience, we examine the B-scan formed from the reconstructed A-scans. FIG. 5 shows an example of the image recovery including: [0062] (i) a B-scan of the original data (labelled the “raw image”), [0063] (ii) the recovered image (labelled the “clean part”). Each column has been obtained from a respective performance of method 1. That is, for a given column of the raw image, a corresponding matrix X is formed. From each X, corresponding matrices L and G are formed by method 1. For each column of the raw image, the corresponding column of the “clean part” image is the column of the corresponding matrix L which has the same position in the corresponding matrix L which the raw image column has in the corresponding matrix X. [0064] (iii) a “noise part” image. For each column of the raw image, the corresponding column of the “noise part” image is the column of the corresponding matrix G which has the same position in the corresponding matrix G which the raw image column has in the corresponding matrix X.

[0065] FIG. 6 shows the result when step 3 is performed using the second variant for step 2 (i.e. method 2), modelling X as X=L+S+G. [0066] (i) a B-scan of the original data (labelled the “raw image”), [0067] (ii) the recovered image (labelled the “clean part”). Each column has been obtained from a respective performance of method 1. That is, for a given column of the raw image, a corresponding matrix X is formed. From each X, corresponding matrices L, S and G are formed by method 2. For each column of the raw image, the corresponding column of the “clean part” image is the column of the corresponding matrix L which has the same position in the corresponding matrix L which the raw image column has in the corresponding matrix X. [0068] (iii) a “spare part” image. For each column of the raw image, the corresponding column of the “spare part” image is the column of the corresponding matrix S which has the same position in the corresponding matrix S which the raw image column has in the corresponding matrix X. [0069] (iv) a “noise part” image. For each column of the raw image, the corresponding column of the “noise part” image is the column of the corresponding matrix G which has the same position in the corresponding matrix G which the raw image column has in the corresponding matrix X.

[0070] In both cases, the image L is much sharper than the original image X. Visually, the two results from the two variants of step 3 are very close without much visible difference.

[0071] To evaluate the image quality, the contrast to noise ratio (CNR) metric [6] is used. It measures the contrast between a set of regions of interest and a background region.

[00008]$\mathrm{CNR}=\left(1/R\right).Math.\underset{l}{\overset{R}{.Math.}}.Math.\left({\mu }_r-{\mu }_b\right)/\sqrt{{\sigma }_r^2+{\sigma }_b^2}$

where μ.sub.r and σ.sub.r are the mean and variance of all the region of interest, μ.sub.b and σ.sub.b are the mean and variance of the background region. Ten regions of interest from different retinal layers and one background region from each B-scan image were used to compute CNR.

[0072] Two typical OCT volumes with 256 B-scan 992×512 images are used for testing. Table 1 shows the comparison between (i) the raw volume data, (ii) a conventional method using a moving algorithm, and (iii) the embodiment using method 2. As we can observe, the embodiment has a significantly higher CNR. FIG. 7 uses box plot to show the variation of CNRs from different frames.

TABLE-US-00001 TABLE 1 Comparison of CNR between the raw data, the conventional moving average speckle reduction algorithm, and an embodiment of the invention Raw Volume Moving Average Embodiment CNR 1.30 2.73 5.26

[0073] FIG. 8 shows visual difference between the proposed and baseline methods using two examples of raw data (FIGS. 8(a) and 8(d)). The images produced by the moving average method are shown in FIGS. 8(b) and 8(e), while the images produced by an embodiment of the present invention using method 2 are shown in FIGS. 8(c) and 8(d).

[0074] Compared to the raw frames and their moving average, the embodiment provides much better volume quality for structure examination. Many more details of the retinal layers and structures such as lamina cribrosa can be examined more comfortably.

[0075] Furthermore, the case of performing step 3 by low rank recovery using bilateral random projection is about 20 times faster than traditional algorithms such as robust PCA [8].

[0076] The computational process may be implemented using a computer system having a processor and a data storage device which stores program instructions for implementation by a processor, to cause the processor to carry out the method above substantially automatically (that is, without human involvement, except optionally for initiation of the method). The computer system may further comprise a display device having a screen for displaying the result of the method, and/or a printing device for generating a printed form of the result of the method, and/or an interface for transmitting the result of the method electronically. The computer system may be a general purpose computer, which receives the raw images from an interface. Alternatively, it may be part of an OCT imaging device which first captures the images and then processes them using the embodiment.

[0077] The resulting images may be used by medical personnel as one tool to perform a diagnosis that a patient is suffering from a condition, and/or to select a medical treatment, such as a surgical procedure, to be carried out on the subject.

REFERENCES

[0078] The disclosure of the following references is incorporated herein by reference: [0079] [1] J. S. Lee, “Speckle analysis and smoothing of synthetic aperture radar images,” Computer Graphics Image Processing, vol. 17, pp. 24-32, 1981. [0080] [2] D. T. Kuan, A. A. Sawchuk, T. C. Strand, and P. Chavel, “Adaptive noise smoothing filter for images with signal dependent noise,” IEEE Trans. Pat. Anal. Mach. Intell., vol. 7, pp. 165-177, 1985. [0081] [3] A. Lopes, R. Touzi, and E. Nesby, “Adaptive speckle filter and scene heterogeneity,” IEEE Trans. Geosci. Remote Sens., vol. 28, pp. 992-1000, 1990. [0082] [4] P. Perona and J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Patt. Ana. Mach. Intell., vol. 12, pp. 629-639, 1990. [0083] [5] Y. Yu and S. T. Acton, “Speckle reducing anisotropic diffusion,” IEEE Trans. Image Processing, vol. 11, pp. 1260-1270, 2002. [0084] [6] S. Aja-Fernandez and C. Alberola-López, “On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering,” IEEE Trans. Image Processing, vol. 15, no. 9, 2006. [0085] [7] K. Krissian, C. F. Westin, R. Kikinis, and K. G. Vosburgh, “Oriented speckle reducing anisotropic diffusion,” IEEE Trans. Image Processing, vol. 16(5), pp. 1412-1424, 2007. [0086] [8] S. Zhu and K. Ma, “A new diamond search algorithm for fast block-matching motion estimation,” IEEE Trans. Image Processing, vol. 9, pp. 287-290, 2000. [0087] [9] T. Zhou and D. Tao, “GoDec: Randomized Low-rank & Sparse Matrix Decomposition in Noisy Case”, International Conference on Machine Learning, pp. 33-40, 2011