Open Access
Add to BibSonomy
Abstract
We present a compressed sensing (CS) algorithm and sampling strategy for reconstructing 3-D Optical Coherence Tomography (OCT) image volumes from as little as 10% of the original data. Reconstruction using the proposed method, Denoising Predictive Coding (DN-PC), is demonstrated for five clinically relevant tissue types including human heart, retina, uterus, breast, and bovine ligament. DN-PC reconstructs the difference between adjacent b-scans in a volume and iteratively applies Gaussian filtering to improve image sparsity. An a-line sampling strategy was developed that can be easily implemented in existing Spectral-Domain OCT systems and reduce scan time by up to 90%.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
OCT is capable of acquiring three-dimensional images at micron resolution over a large field of view. A typical OCT image volume can contain over 100 million pixels of information. The data requirements of OCT imaging experiments requiring time-lapse imaging (both 2-D and 3-D in time) [1–3], mosaic imaging [4,5], or real-time acquisition [6] can meet or exceed the data through-put capabilities of image acquisition hardware. These restrictions may prohibit data-intensive experiments or necessitate specialized solutions for handling and storing terabytes of data. Long acquisition times can affect image quality through motion artifacts, particularly for in-vivo imaging [7].
Compressed Sensing (CS) is a technique in sparse representation theory for reconstructing highly undersampled images at full-resolution under the assumption that the signal is sparse in some basis [8]. For a known sampling pattern, the problem can be modeled as a linear relationship $y = Ax$ where $y$ is the observed, undersampled signal and $x$ is the sparse, fully sampled signal. The sensing matrix $A$ provides a mapping between $x$ and $y$. Typically, $A$ represents a random sampling pattern, however, this matrix can also represent more structured sampling patterns that can be emulated by imaging hardware [9]. The linear CS equation is under-determined, but the signal $x$ can be exactly recovered using convex optimization [10]. This simple optimization problem can then be modified and applied to a diverse set of applications by taking advantage of structure in $x$ [11,12].
CS has revolutionized imaging fields like MRI by decreasing image acquisition time and data storage needs by up to 90% [13–15]. Many CS approaches for OCT have been proposed that aim to reconstruct the raw interferogram or other hardware-specific signals [16–22]. Studies combining CS-OCT with alternate imaging methods like full-depth OCT have also been explored [23–25]. This manuscript focuses on methods for reconstructing Spectral-Domain OCT (SD-OCT) volumes as a post-processing approach. Lebed, et al. first demonstrated CS reconstruction of OCT volumes by modifying the scanning pattern to randomly omit full b-scans and a-lines [26]. Xu, et al. have published several studies investigating 3-D CS-OCT by undersampling and reconstructing both the raw interferogram and the image volume in a multi-step reconstruction process [27–29]. Learning-based approaches have been proposed for situations where the sample type is known a-priori using an "energy-guided" approach [30] and Dictionary Learning [31]. These are the first CS-OCT studies to leverage image structure to improve reconstruction accuracy.
Predictive Coding (PC) is a general approach for data compression that leverages structure by representing a signal in terms of the change between successive data points [32]. The underlying assumption of PC is that changes between data points are small and infrequent which makes the difference representation more efficient. Predictive models in CS that assume spatio-temporal structure have been applied to natural images [33], hyperspectral imaging [34], and MRI [15,35,36]. To our knowledge, PC approaches have not previously been applied to CS-OCT.
This study proposes a novel approach to 3-D CS-OCT using a Denoising Predictive Coding (DN-PC) approach that takes advantage of the inherent spatial structure in OCT volumes. We show that by reconstructing the difference between adjacent b-scans in a volume, higher reconstruction accuracy is achieved over traditional methods. Using DN-PC, we provide the first demonstration of CS-OCT reconstruction for a diverse collection of biological samples with complex tissue structure including retina, cardiac tissue, uterine tissue, breast tissue, and ligament. We describe our novel imaging method for rapid acquisition and reconstruction of OCT volumes and demonstrate its success on artificially sub-sampled volumes. Additionally, the role of speckle noise in reconstruction performance is explored. Reconstruction performance is compared with popular CS algorithms through qualitative and quantitative assessment.
2. Methods
$\ell ^{1}$ CS signal recovery is first introduced as the traditional approach, then expanded to predictive coding and the novel approach DN-PC. We begin by defining vectorized images $x \in \mathbb {R}^{N}$ and $y \in \mathbb {R}^{M}$ as the full-resolution and undersampled images, respectively. The signal $x$ can be recovered from $y$ by solving the objective function
where $\mathbf {\Psi }$ is the sparse representation basis (e.g. DFT) and $\mathbf {A}$ is an $M \times N$ matrix which encodes the undersampling pattern. Multiple methods exist for solving an objective function of this form such as Iterative Soft Thresholding (IST) [37] or Alternating Directions Method of Multipliers (ADMM) [38].Suppose $x_t$ and $x_{t-1}$ represent adjacent images in a volumetric scan and the difference image is $\nabla x = x_t - x_{t-1}$. Using the framework proposed in [39], Eq. (1) can be modified to solve for $\nabla x$ as follows.
This objective function has the same form as Eq. (1) so it can be solved using an identical solver. Speckle noise is inevitable source of corruption in OCT images, degrading image quality and potentially hindering accurate CS reconstruction. Consequently, incorporating denoising into the objective function may improve reconstruction performance. Suppose rather than using $\ell ^{1}$ regularization on the difference image $\nabla x$, it was applied to the denoised version through function $D(x,\lambda )$ where $\lambda$ is a denoising parameter. In this case the objective function becomes the following function.
The iteration begins with an initial guess $\nabla x_{i-1}$ and is penalized to agree with the observation $y$. The second equation controls the sparsity of the solution via the $\ell ^{1}$ norm. Noting that the OCT image is sparse when denoised and transformed to the Fourier basis, $\Psi D(\nabla x)$ is penalized rather than $\nabla x$ itself. The change of basis is necessary to ensure incoherence between the representation and measurement domains [8]. We chose to use a Gaussian filter for $D(\nabla x)$, although other denoising methods such as BM3D have been demonstrated [41].
The first subproblem is solved by taking the derivative, setting it equal to zero, and solving for $\nabla x$. Taking the derivative gives
In this formulation, $\alpha$ is a rough measure of the noise in observation $y$ where $\alpha = 0$ corresponds to the noiseless case.
The solution of the second equation is found by using the proximity operator and takes the following form
which for $\textrm {prox}_{\lambda \|\cdot \|_1}$ is evaluated by performing element-wise soft-thresholding of the argument. The soft-thresholding operation soft() of matrix element $u_i$ by threshold $\beta$ is defined for complex-valued entries as $\textrm {sign}(u_i)\max (|u_i|-\beta ,0)$.A summary of the DN-PC method is given in Algorithm 1. A key feature of this method is the use of an adaptive denoising parameter $\lambda$. Similar to the approach used in [41], more of the important image features can be recovered by first denoising strongly, and then iteratively decreasing the degree of denoising. In DN-PC, $\lambda$ has two values $(\lambda _1,\lambda _2)$ which represent the vertical and horizontal standard deviation of the 2-D Gaussian Filter $D(\cdot ,(\lambda _1,\lambda _2))$. The variability of $\lambda$ is controlled by setting $\lambda _{max}$ and $\lambda _{min}$ such that $\lambda$ decreases logarithmically over $\mathcal {J}$ iterations. $\lambda _{max}$ and $\lambda _{min}$ may be set differently for $\lambda _1$ and $\lambda _2$. The algorithm is structured to update over an inner and outer iteration. The inner iteration solves subproblems for a fixed value of $\lambda$ until the update reaches max iteration $\mathcal {I}$ or the solution update becomes small (set using $\tau$), while the outer loop iterates $\mathcal {J}$ times over $\lambda$.
 |
2.1 Compressed sensing pipeline
A flow chart describing the DN-PC image reconstruction process is shown in Figure 1. Undersampling of the OCT volume is simulated by omitting a-lines at a regular interval. The sparsely sampled OCT volume is reconstructed by iterating over $32 \times 32$ square pixel patches of each b-scan. We empirically tested square patches of different sizes and found $32 \times 32$ pixel patches to provide the best trade-off between reconstruction time and accuracy. A given patch $(m,n)$ is reconstructed over all $T$ b-scans before advancing to the next patch, where $m$ and $n$ are the row and column indices of a patch, respectively. To reconstruct patch $(m,n)$ at b-scan $t$, the difference image is acquired by first undersampling then subtracting patch $(m,n)$ at b-scan $t-1$. Patch $t-1$ is undersampled by multiplying it by the sampling matrix $\mathbf {A}$. The DN-PC algorithm produces a reconstruction of the difference patch which is then added to the full resolution patch at $(m,n,t-1)$ to get the reconstructed patch $(m,n,t)$.
One challenge in reconstructing the difference image is that the reconstruction accuracy is dependent on the patch from the previous b-scan because errors from each b-scan can propagate through the entire volume. A unique sampling scheme (Fig. 2) is proposed to mitigate this problem by using staggered sampling and periodic full-resolution acquisitions. Staggered sampling means that the sampling pattern is shifted by one a-line between adjacent b-scans to avoid omitting the same a-lines for the entire volume. Full-resolution b-scans are acquired periodically to "reset" any remaining propagated error.
Fig. 2. Diagram of DN-PC sampling using a-line staggering and periodic full-resolution b-scan acquisitions.
2.2 Sampling
In this manuscript, "compression rate" $\eta _{a}$ refers to the number of a-lines sampled in each image patch, i.e. a 25% compression rate means that one in every four a-lines were acquired. The compression rate of a b-scan $\eta _{b}$ in units of pixels can be calculated as follows
where the operator floor() rounds the argument down to the nearest integer value, and $N_{patch}$ is the total number of pixels per image patch. Compression is defined in a third way for volumes which takes into account the periodically acquired full-resolution b-scans. This is the true compression rate $\eta$ which is a function of the full-resolution b-scan interval $\mathcal {I}_{b}$. If a full-resolution b-scan is acquired every ten b-scans in the volume, then $\mathcal {I}_{b} = 10$. The true compression rate $\eta$ is calculated as follows where $N_{b}$ is the number of pixels per a full-resolution b-scan, $N_{vol}$ is the number of pixels per a full-resolution volume, and $T$ is the number of b-scans in the volume.Random a-line sampling is also employed at times in this study. In this case the randomly chosen a-lines are consistent for patches in the same column to replicate random a-line sampling in hardware. A-lines were chosen pseudo-randomly to constrain the maximum distance between sampled a-lines depending on the sampling rate. Staggering was not employed in this case and instead a new pattern was generated for each b-scan in the volume. Equation (10) applies to random sampling in the same way as the other sampling methods.
3. Experimental methods
Multiple experiments were performed to test and validate the proposed DN-PC algorithm for volumetric CS-OCT reconstruction. We first analyze the effects of denoising and PC by examining pixel decay plots of an OCT image compared with its difference image and a noise-only image. Next, we evaluated sampling strategies by testing reconstruction accuracy at multiple a-line sampling rates and comparing staggered with uniform a-line sampling. The proposed DN-PC method was then evaluated on OCT volumetric datasets of five different tissue samples acquired at full-resolution and then synthetically sub-sampled. The results were quantitatively evaluated and representative images were selected for visual comparison.
3.1 CS algorithms for comparison
The performance of DN-PC was compared with two other algorithms. The first method uses patch-based reconstruction of the raw OCT b-scan and iterates over all b-scans in the volume. The employed algorithm, called YALL1, is an optimized technique for $\ell ^{1}$ minimization [42], i.e. is solves Eq. (1). The second method also uses a PC approach, but with our own implementation of TVL1 reconstruction based on the method from Yang, et al. [43] called RecPF. This algorithm utilizes a Total-Variation (TV) regularization term which promotes smoothness while also preserving edges. Our implementation simply allows the method to be used on OCT volumes and in the Predictive Coding framework, so we refer to it as TVL1-PC. All three algorithms were tested and implemented in MATLAB 2020a using a Windows 10 desktop with an Intel Core i9-9900K CPU at 3.6 GHz and 128 GB of RAM.
3.2 Algorithm parameters
All methods tested in this study rely on reconstruction parameters which affect end performance. We chose these parameters empirically and utilized the same ones in all tests unless otherwise specified. For DN-PC, we chose $\alpha = 0.1$, $\beta = 1$, $\lambda _{\textrm {max}} = [3,4]$, $\lambda _{\textrm {min}} = [0.2,0.4]$, $\mathcal {J} = 20$, $\mathcal {I} = 20$, and convergence threshold $\tau = 10^{-3}$. We used a filter size of $7 \times 9$ for the Gaussian denoising filter, which is rectangular to smooth vertically streaking that can appear as a result of a-line subsampling. A Gaussian filter was chosen over other filters because it reduces image noise and because it is a linear filter. The linearity is important because it means that filtering the difference image is equivalent to filtering its constitutive frames and then taking the difference. For YALL1, the Discrete Cosine Transform (DCT) was chosen as the sparsifying basis. The convergence tolerance was $5*10^{-4}$ and we set parameter $\rho = 5*10^{-4}$. While staggered a-line sampling and periodic full-resolution b-scan sampling are not necessary to use with YALL1, they are both employed in all cases for accurate comparison. For TVL1-PC, $\Psi$ is the level-3 Haar wavelet, the anisotropic TV measure is used, and the remaining parameters are set to $\mu = 10^{4}$, $\beta = 20$, $\tau = 0.5$, and $\gamma = (\sqrt {5} + 1)/8$.
3.3 Image acquisition and datasets
Each of the five datasets used in this study contain OCT volumes of different, structurally complex tissue samples: human right atria [4,44], human uterus [5], human retina [45], bovine Anterior Cruciate Ligament (ACL) [46,47], and human breast [48]. The human retina data is publicly available and managed by Farsiu, et al. [45]. The heart, uterus, ACL, and breast datasets were collected internally using a commercial TELESTO SD-OCT system (Thorlabs, GmbH, Germany) with 6.5 $\mu$m axial and 15 $\mu$m lateral resolution. Each dataset was imaged with some degree of lateral oversampling which is important to consider when analyzing CS reconstruction results. The lateral spacings for the uterus, ACL, breast, heart, and retina datasets are 4 $\mu$m, 4.4 $\mu$m, 7.1429 $\mu$m, 6.3 $\mu$m, and 6.7 $\mu$m, respectively. All datasets had equal spacing between b-scans except for the retina dataset which had 67 $\mu$m b-scan spacing. All OCT volumes were cropped to $512 \times 800 \times 800$ pixels for consistent comparison with the exception of the retina volumes which have only 100 b-scans. Prior to reconstruction, all datasets were converted to double precision and the pixel intensity was scaled to a range of $[0,1]$.
3.4 Metrics
Several quantitative metrics were used to assess and compare CS reconstruction performance. The first is Relative Error which measures the intensity differences between the true and reconstructed volumes. It is defined as
where $x$ is the vectorized original OCT volume and $x_{\textrm {recon}}$ is the reconstructed version. When evaluating the relative error for images, the Frobenious norm is used instead. Structural Similarity Index (SSIM) is another popular metric which uses luminance, contrast, and structure to evaluate the similarity between two images [49]. The SSIM of two images is a value between 0 and 1 where an SSIM of 1 indicates that the two images are identical. Where SSIM is reported for a volume, we provide the average SSIM over all b-scans in the volume. Because we are reconstructing image volumes, we also measure the 3-D Multi-Scale SSIM (MULTI-SSIM 3D) [50]. This is a variant of the SSIM metric for image volumes that applies the same algorithm at multiple scales and produces an aggregate score.While measuring exact reconstruction error is important, we were also interested in analyzing our ability to reconstruct important tissue features independently from speckle noise and other noise sources. 2-D median filtering is a popular and light-weight choice for OCT image denoising. In some cases, we applied a $3 \times 3$ pixel median filter to both the ground truth and the reconstructed volumes before measuring relative error and SSIM to obtain a more honest assessment of the algorithm’s ability to reconstruct important tissue structures. Denoised metrics are reported using the identifier (DN) (see Table 2).
4. Results
Figure 3 demonstrates how image sparsity is affected by denoising and difference image operations using a small image patch of a glass cover slip. The pixel decay plots were generated by vectorizing the image patch and sorting the pixels in descending order of intensity. Plots which decay to zero more quickly correspond to a sparser image. The first column of images (Fig. 3(a), 3(c), and 3(e)) are image patches while the second column (Fig. 3(b), 3(d), and 3(f)) are the corresponding difference images. The rows show an image of the cover slip (Fig. 3(a)), the same image when denoised (Fig. 3(c)), and a patch of only noise (Fig. 3(e)). Figures 3(g) and 3(h) show the pixel decay plots for the six image patches in the image domain and Discrete Cosine domain, respectively. The images in Fig. 3(a)–3(f) show that the difference operation preserves noise, but denoising prior to taking the difference (Fig. 3(c), 3(d)) isolates the structural differences of interest between adjacent b-scans. In the image domain, the noise patch is the least sparse while the denoised difference image is the most sparse. In all cases, the difference operation and denoising created sparser image patches than their counterparts.
Fig. 3. Speckle noise degrades image sparsity. (a) Raw OCT image patch. (b) Raw difference image. (c) Denoised image. (d) Denoised difference image. (e) Speckle noise OCT image patch. (f) Noisy difference image. (g) Pixel decay plot for images a-f. (h) Pixel decay plot for the discrete cosine transformed (DCT) images a-f. g and h show that the denoised difference image is the sparsest of images a-f both in the image and DCT domains. Scale bar = 100 $\mu$m.
The effects of sampling parameters on reconstruction performance were tested to determine an optimal reconstruction configuration. Figure 4 shows the relative error from reconstructing a 50 b-scan subset of an OCT heart volume. In Fig. 4(a), three different a-line sampling rates $\eta _a =$ 50%, 25%, and 10% were tested using a full-resolution interval $\mathcal {I}_b = 25$ b-scans and a-line staggering. The staggering suppresses error as a function of distance from the last full-resolution b-scan, though a small linear increase in the error is visible with 10% sampling. Figure 4(b) demonstrates the effect of staggered sampling by comparing the relative error of the same reconstructed volume using $\eta _a = 50\%$ but with and without staggering. In the "no staggering" case, the same a-lines are omitted every b-scan. In both cases, the b-scan at index 1 is fully sampled. Without staggering, the error increases linearly from 0.3 to 0.33 over 50 b-scans. With staggering, the error dips initially and then plateaus to a value around 0.27. Not only did staggering lower the average error, but it also suppressed the rate of error as a function of distance from the last full-resolution b-scan.
Fig. 4. Compression analysis. (a) Relative reconstruction error for 50 b-scans at 3 a-line sampling rates using the DN-PC method with a-line staggering and a full-resolution interval of 25 b-scans (dashed line). Lower sampling rates cause both the average error and the error rate as a function of distance from the full resolution b-scan to increase. (b) Relative reconstruction error with and without a-line staggering at 50% sampling with a 50 b-scan full-resolution interval. Alternating acquired a-lines between adjacent b-scans mitigates residual error.
The observations from Fig. 4 were quantitatively verified for a full OCT volume from the human cardiac dataset and reported in Table 1. Tweleve use-cases were tested using different sampling methods, two a-line sampling rates $\eta _a = 50\%, 25\%$, and two full-resolution intervals $\mathcal {I}_b = 10, 50$. Three sampling patterns were used to compare the effects of staggering and uniform versus random sampling. For the random sampling cases, the max gap was 3 a-lines for 50% sampling and 6 a-lines for 25% sampling. The OCT volume dimensions were $512 \times 800 \times 800$ pixels and we used a patch size of $32 \times 32$ pixels. The "Full-Res B-Scans" column of the table shows the total number of full-resolution b-scans obtained for the two intervals. Similarly, the column "Sampled A-Lines/B-Scan" shows that 25% and 50% sampling results in acquisition of 200 and 400 a-lines per b-scan, respectively. The true compression rate $\eta$ includes the full-resolution b-scans so it is higher than the a-line sampling rate $\eta _a$ (see Eq. (10)), though the margin of increase is larger for smaller sampling rates. In all cases, staggering improves the relative reconstruction error. The full-resolution interval trades off between relative error and $\eta$. For example, in the case of $\eta _a = 25\%$ with staggering on, the relative error improves from 0.28 to 0.27 when $\mathcal {I}_b$ is lowered from 50 to 10, but at the expense of raising $\eta$ from 26.5% to 32.5%.
Table 1. Quantitative summary of the effects of a-line sampling $\eta _a$, staggering ,and the full-resolution interval $\mathcal {I}_b$ on compression and relative error.
4.1 USAF resolution target
A volume of the USAF 1951 Resolution Target was acquired to assess the impact of different sampling rates on the lateral resolution of the reconstruction. The volume was acquired without lateral or axial oversampling (i.e. with 15 $\mu$m lateral spacing). Figure 5 shows an example en-face image from the full-resolution volume and reconstructed volumes using 50%, 25%, and 10% a-line sampling (Fig. 5(a)–5(d)). Insets showing magnified views of group 2 for each reconstruction are shown in Fig. 5(e)–5(h). Element 5 of group 2 is the smallest resolvable element in the full-resolution acquisition and in the 50% sampling reconstruction. With 25% and 10% sampling, elements 4 and 3 are the smallest resolvable line pairs, respectively.
Fig. 5. USAF Target lateral resolution test. En-Face images from OCT volume of the USAF target at (a) full resolution, (b) 50% a-line sampling, (c) 25% a-line sampling, (d) 10% a-line sampling. (e-h) Insets showing a magnified view of elements in line group 2.
4.2 Multiple tissue type test
DN-PC was used to reconstruct OCT volumes from five different tissue samples: human heart, human uterus, human retina, bovine ACL, and human breast tissue. Example b-scans from each of the reconstructed volumes are shown in Fig. 6. The different tissue types are organized by row and the different sampling rates are organized by column. The different samples and images were chosen to showcase a variety of tissue structures, image textures, and noise environments. Qualitatively, the examples with 50% sampling are nearly indistinguishable from the corresponding full-resolution b-scans, while the 10% samples appear noisier and fine features are blurred. Example en-face images from the same volumes are shown in Fig. 7. Similar degradation of image quality is observed for 10% a-line sampling compared with 50%. Unlike in the b-scan images, horizontal streaking is visible in the en-face images along the fast-scan axis which are artifacts from errors in reconstruction. The retina volumetric scans include only 100 b-scans so en-face images from those samples were omitted as they do not provide valuable information even in the full-resolution volume.
Fig. 6. B-Scan examples from DN-PC reconstructed OCT volumes. Rows correspond to the tissue sample and columns correspond to the a-line sampling rate. Scale bar = 500 $\mu$m.
Fig. 7. En-Face image examples from DN-PC reconstructed OCT volumes. Rows correspond to the tissue sample and columns correspond to the a-line sampling rate. Scale bar = 500 $\mu$m.
Reconstructed volumes from the uterus and ACL datasets were rendered in 3-D to compare volumetric features with the full-resolution volumes. Figure 8 shows images from the uterus volume rendering in the first row and the ACL volume in the second row. Sampling rates are organized by column. Collagen fibers were labelled and identified in the full-resolution volumes (first column) which are visible in the reconstructions at both 50% and 25% a-line sampling. The 3-D perspective shows the ability of DN-PC reconstructed volumes to preserve volumetric features visible in both the en-face and axial image planes.
Fig. 8. 3-D Reconstructions. OCT volumes of human uterus (first row) and bovine ACL (second row) tissue reconstructed using DN-PC. The first column is the fully sampled volume, the second column was reconstructed with 50% a-line sampling, and the third column was reconstructed with 10% a-line sampling. The white arrows point to tissue structures and artifacts of interest. Scale bar = 150 $\mu$m.
DN-PC volumetric reconstruction performance for 5 different tissue types was quantitatively measured and reported in Table 2 which shows the relative error and average SSIM of each reconstruction. A representative OCT volume from each of the 5 tissue types was reconstructed at three a-line sampling rates $\eta _a =$ 50%, 25%, and 10% using staggering and $\mathcal {I}_b = 10$. Relative error and SSIM are reported with and without denoising (labeled DN) following reconstruction. The denoised results are improved over the raw data results across all test cases which suggests strong preservation of tissue structures. DN-PC achieved the best performance for the cardiac volume, while the retina and breast volumes proved the most challenging.
Table 2. Quantitative summary of relative error using DN-PC for 5 different tissue sample types at a-line sampling rates of 50%, 25%, and 10%.
4.3 Algorithm test
DN-PC performance was compared with two other CS reconstruction methods, YALL1 and TV-L1 PC using 100 b-scan sub-volumes of all five tissue samples at $\eta _a = 50\%$. Staggering and $\mathcal {I}_b = 10$ were used in all cases for accurate comparison. In each case, relative error, average SSIM, MULTI-SSIM 3D, and computation time were recorded. Quantitative results are reported in Table 3.
Table 3. Reconstruction algorithm comparison test using 50% a-line sampling.
Figure 9 shows an example b-scan from the heart dataset at full-resolution and reconstructed using each algorithm at 50% a-line sampling. Images in Fig. 9(a)–9(d) are the full-resolution b-scan, YALL1 reconstruction, TVL1-PC reconstruction, and DN-PC reconstruction, in that order. Insets of a magnified portion of the full-resolution myocardial tissue surface are shown for each algorithm in Fig. 9(e)–9(h), where the inset is marked by the rectangle (red). Difference images at the same inset location are shown in Fig. 9(i)–9(l). The insets show that YALL1 is susceptible to streaking artifacts from the a-line sampling. The TVL1-PC reconstruction is of similar quality to the original image, however, the difference image is also has streaking. DN-PC did not reconstruct the original image as precisely as TVL1-PC, but the difference image contains more changes from tissue structure rather than an exact noise pattern. Comparing with Table 3, DN-PC has similar relative error and worse SSIM score then TVL1-PC and takes considerably less time to reconstruct. In the case of the heart sample, the 100 b-scan volume was reconstructed in 19.12 minutes with DN-PC and 616.46 minutes (over 10 hours) with TVL1-PC. Average SSIM tended to have a large discrepancy between the TVL1-PC and DN-PC results despite qualitatively appearing very similar. The MULTI-SSIM 3-D metric gave much better scores to all the reconstructions and reflected the qualitative similarity between TVL1-PC and DN-PC reconstructions as observed in Fig. 9.
Fig. 9. Heart tissue CS reconstruction (b-scans) at 50% a-line sampling. (a) Original. (b) YALL1 reconstruction. (c) TVL1-PC reconstruction. (d) DN-PC reconstruction. Insets indicated by the red square for the original (e), YALL1 (f), TV-L1 (g), and DN-PC (h) images. Insets from the difference images for the original (i), YALL1 (j), TV-L1 (k), and DN-PC (l). Inset location is marked by the rectangle (red) in a-d. Scale bar = 500 $\mu$m.
5. Discussion
This study explored the proposed compressed sensing technique Denoised Predictive Coding (DN-PC) and it ability to reconstruct highly undersampled OCT volumes. This study is first to evaluate a CS-OCT method using five different, clinically relevant tissue types. A critical step for CS to be adopted in clinical OCT systems is to demonstrate that the method provides reliable reconstructions of complex, varied tissue types without prior knowledge of the sample. Quantitative results indicated that tissue type did not affect image reconstruction performance to the same degree as other parameters like sampling rate. However, two sample types, retina and breast, were more challenging to reconstruct than the others. It’s likely that the retina dataset had higher reconstruction error because it was acquired using a different OCT system than the other four datasets. The noise variance in particular higher for the retina dataset, suggesting that denoising parameters like $\lambda _{\textrm {max}}$, $\lambda _{\textrm {min}}$ should be adjusted for image volumes collected with different OCT systems. The source of error in the breast sample is less clear, but one explanation is that adipose tissue could be a difficult feature to reconstruct. Adipose visually look like small bubbles in OCT b-scans and because DN-PC excels at preserving the overall tissue structure it will struggle the most when the tissue is composed of mostly small, fine features. This problem could be mitigated by adjusting the denoising parameters to prevent potential blurring of the adipose edges. While this test laid important groundwork for developing a CS-OCT solution that can generalize to any tissue type, further investigation using a larger pool of samples is needed to understand how reconstruction behavior is related to tissue structure.
Reconstruction performance was characterized in this manuscript by relative error, SSIM, MULTI-SSIM 3D, and computation time. Because no gold standard metric exists to assess reconstruction accuracy, these metrics were chosen assuming that they were the most common and intuitive measures available. One area of ambiguity with regard to performance analysis is the reconstruction of noisy images (which applies to all OCT images). In Table 2 for example, the relative error of the raw reconstructions differed significantly from the denoised (DN) reconstructions. Ultimately, we felt it was important to include both measures because they inherently explain different aspects of the algorithm performance. The denoised results measures the ability to reconstruct important tissue structures independently of noise, while the raw reconstruction results measure how closely the reconstruction exactly matches the raw image, which could be equally important in applications like Speckle Variance imaging that rely on noise statistics [2].
The algorithm comparison test in Table 3 revealed interesting results and behavior of the tested methods. While DN-PC indeed has a consistent disadvantage in SSIM, the SSIM from the YALL1 results are misleading as the images contain streaking artifacts which render the images unusable. Comparing the ability of each method to preserve tissue structure, the MULTI-SSIM 3D gives a more accurate representation. TV-L1 consistently outperforms DN-PC by a small margin, but at the expense of infeasibly long computation time. We recommend that DN-PC be used over the other methods in any situation in which the user intends to use a-line subsampling, would like to reconstruct the volume quickly, and if the b-scan spacing is small. If reconstruction time on the order of hours is no issue for the user, then the TV-L1 approach is an adequate substitute.
A unique challenge in any CS framework is the formulation of a sampling strategy which works with the imaging hardware and enables high accuracy reconstruction of the undersampled data. Recent studies have proposed hardware techniques for undersampling using CCD camera masking materials [21] and masking spectral data through a DAQ [20], however, these methods only compress the signal without improving acquisition time. In order to compress and acquire volumes more quickly, modifications have to be made to the scanning method. Wang, et al. demonstrated this for an OCT endoscope by randomly changing the step-size during pull-back acquisition [22], however, this approach is specific to pull-back endoscopes and cannot be applied to bench-top systems. In a-line subsampling, the desired undersampling rate can be controlled by over-driving the lateral scanning mechanism (e.g. galvo) to the desired speed. This modification can be applied to existing OCT systems with virtually no hardware changes. Furthermore, undersampling this way directly reduces scan time. For example, using DN-PC with 25% a-line sampling and a full-resolution interval of 10 b-scans would reduce a one minute scan to 19.5 seconds. This motivated our choice to design a reconstruction algorithm specifically for a-line subsampling rather than spectral subsampling. Reducing scan time could impact areas of OCT research like whole organ imaging [4,51,52], high-speed endoscopic imaging [53,54], and 4-D imaging [6,55]. With an extension to time-lapse imaging, CS-OCT could impact additional application such as particle tracking [56,57], elastography [58–60], cilia and mucus movement [61,62], developmental biology [63], and Radio-Frequency Ablation (RFA) [64–66].
6. Conclusion
We developed and tested Denoising Predictive Coding (DN-PC), a new method for 3-D CS-OCT reconstruction of highly undersampled image volumes. The unique combination of predictive coding and integrated denoising yielded high accuracy reconstructions and superior computation time over similar methods. We demonstrated that DN-PC was robust to tissue type without a-priori knowledge of the sample. CS-OCT has the potential to enable high data volume experiments with long scan times that were previously infeasible and represents an important step towards commercial and clinical adoption of CS-OCT.
Funding
National Institute of Biomedical Imaging and Bioengineering (4DP2HL127776-02).
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. W. Wieser, W. Draxinger, T. Klein, S. Karpf, T. Pfeiffer, and R. Huber, “High definition live 3D-OCT in vivo: design and evaluation of a 4d oct engine with 1 gvoxel/s,” Biomed. Opt. Express 5(9), 2963–2977 (2014). [CrossRef]
2. Y. Ling, X. Yao, U. A. Gamm, E. Arteaga-Solis, C. W. Emala, M. A. Choma, and C. P. Hendon, “Ex vivo visualization of human ciliated epithelium and quantitative analysis of induced flow dynamics by using optical coherence tomography,” Lasers Surg. Med. 49(3), 270–279 (2017). [CrossRef]
3. J. P. McLean, Y. Ling, and C. P. Hendon, “Frequency-constrained robust principal component analysis: a sparse representation approach to segmentation of dynamic features in optical coherence tomography imaging,” Opt. Express 25(21), 25819–25830 (2017). [CrossRef]
4. T. H. Lye, K. P. Vincent, A. D. McCulloch, and C. P. Hendon, “Tissue-specific optical mapping models of swine atria informed by optical coherence tomography,” Biophys. J. 114(6), 1477–1489 (2018). [CrossRef]
5. J. P. McLean, S. Fang, G. Gallos, K. M. Myers, and C. P. Hendon, “Three-dimensional collagen fiber mapping and tractography of human uterine tissue using OCT,” Biomed. Opt. Express 11(10), 5518 (2020). [CrossRef]
6. J. P. Kolb, W. Draxinger, J. Klee, T. Pfeiffer, M. Eibl, T. Klein, W. Wieser, and R. Huber, “Correction: Live video rate volumetric OCT imaging of the retina with multi-MHz A-scan rates,” PLoS One 14(7), e0220829 (2019). [CrossRef]
7. Y. Chen, Y.-J. Hong, S. Makita, and Y. Yasuno, “Three-dimensional eye motion correction by lissajous scan optical coherence tomography,” Biomed. Opt. Express 8(3), 1783–1802 (2017). [CrossRef]
8. D. L. Donoho and M. Elad, “Optimally sparse representation in general (nonorthogonal) dictionaries via 1 minimization,” Proc. Natl. Acad. Sci. U. S. A. 100(5), 2197–2202 (2003). [CrossRef]
9. M. F. Duarte and Y. C. Eldar, “Structured compressed sensing: From theory to applications,” IEEE Trans. Signal Process. 59(9), 4053–4085 (2011). [CrossRef]
10. J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Trans. Inform. Theory 53(12), 4655–4666 (2007). [CrossRef]
11. S. Oshery, Z. Shiz, and W. Zhuy, “Low dimensional manifold model for image processing,” SIAM J. Imaging Sci. 10(4), 1669–1690 (2017). [CrossRef]
12. E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis?” J. Assoc. Comput. Mach. 58(3), 1–37 (2011). [CrossRef]
13. H. Jung, K. Sung, K. S. Nayak, E. Y. Kim, and J. C. Ye, “K-t FOCUSS: A general compressed sensing framework for high resolution dynamic MRI,” Magn. Reson. Med. 61(1), 103–116 (2009). [CrossRef]
14. R. Otazo, E. Candès, and D. K. Sodickson, “Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components,” Magn. Reson. Med. 73(3), 1125–1136 (2015). [CrossRef]
15. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: the application of compressed sensing for rapid MR imaging,” Magn. Reson. Med. 58, 1182–1195 (2007). [CrossRef]
16. X. Liu and J. U. Kang, “Compressive SD-OCT: the application of compressed sensing in spectral domain optical coherence tomography,” Opt. Express 18(21), 22010 (2010). [CrossRef]
17. D. Xu, N. Vaswani, Y. Huang, and J. U. Kang, “Modified compressive sensing optical coherence tomography with noise reduction,” Opt. Lett. 37(20), 4209 (2012). [CrossRef]
18. N. Zhang, T. Huo, C. Wang, T. Chen, J.-G. Zheng, and P. Xue, “Compressed sensing with linear-in-wavenumber sampling in spectral-domain optical coherence tomography,” Opt. Lett. 37(15), 3075 (2012). [CrossRef]
19. D. Xu, Y. Huang, and J. U. Kang, “Compressive sensing with dispersion compensation on non-linear wavenumber sampled spectral domain optical coherence tomography,” Biomed. Opt. Express 4(9), 1519 (2013). [CrossRef]
20. Y. Ling, W. Meiniel, R. Singh-Moon, E. Angelini, J.-C. Olivo-Marin, and C. P. Hendon, “Compressed sensing-enabled phase-sensitive swept-source optical coherence tomography,” Opt. Express 27(2), 855 (2019). [CrossRef]
21. W. Liao, J. Hsieh, C. Wang, W. Zhang, S. Ai, Z. Peng, Z. Chen, B. He, X. Zhang, N. Zhang, B. Tang, and P. Xue, “Compressed sensing spectral domain optical coherence tomography with a hardware sparse-sampled camera,” Opt. Lett. 44(12), 2955 (2019). [CrossRef]
22. J. Wang, Y. Hu, and J. Wu, “Three-dimensional endoscopic OCT using sparse sampling with a miniature magnetic-driven scanning probe,” Appl. Opt. 57(34), 10056 (2018). [CrossRef]
23. L. Yi and L. Sun, “Full-depth compressive sensing spectral-domain optical coherence tomography based on a compressive dispersion encoding method,” Appl. Opt. 57(31), 9316 (2018). [CrossRef]
24. L. Yi, L. Sun, X. Guo, and B. Hou, “Combination of 2D compressive sensing spectral domain optical coherence tomography and interferometric synthetic aperture microscopy,” Appl. Sci. 9(19), 4003 (2019). [CrossRef]
25. C. K. Mididoddi, F. Bai, G. Wang, J. Liu, S. Gibson, and C. Wang, “High-throughput photonic time-stretch optical coherence tomography with data compression,” IEEE Photonics J. 9(4), 1–15 (2017). [CrossRef]
26. E. Lebed, P. J. Mackenzie, M. V. Sarunic, and F. M. Beg, “Rapid volumetric OCT image acquisition using compressive sampling,” Opt. Express 18(20), 21003 (2010). [CrossRef]
27. D. Xu, Y. Huang, and J. U. Kang, “Real-time compressive sensing spectral domain optical coherence tomography,” Opt. Lett. 39(1), 76 (2014). [CrossRef]
28. D. Xu, Y. Huang, and J. U. Kang, “GPU-accelerated non-uniform fast Fourier transform-based compressive sensing spectral domain optical coherence tomography,” Opt. Express 22(12), 14871 (2014). [CrossRef]
29. D. Xu, Y. Huang, and J. U. Kang, “Volumetric (3D) compressive sensing spectral domain optical coherence tomography,” Biomed. Opt. Express 5(11), 3921 (2014). [CrossRef]
30. S. Schwartz, C. Liu, A. Wong, D. A. Clausi, P. Fieguth, and K. Bizheva, “Energy-guided learning approach to compressive FD-OCT,” Opt. Express 21(1), 329 (2013). [CrossRef]
31. L. Fang, S. Li, R. P. McNabb, Q. Nie, A. N. Kuo, C. A. Toth, J. A. Izatt, and S. Farsiu, “Fast acquisition and reconstruction of optical coherence tomography images via sparse representation,” IEEE Trans. Med. Imaging 32(11), 2034–2049 (2013). [CrossRef]
32. Y. Huang and R. P. Rao, “Predictive coding,” WIREs Cogn. Sci. 2(5), 580–593 (2011). [CrossRef]
33. J. Zhang, D. Zhao, and F. Jiang, “Spatially directional predictive coding for block-based compressive sensing of natural images,” 2013 IEEE International Conference on Image Processing, ICIP 2013 - Proceedings pp. 1021–1025 (2013).
34. F. Rizzo, B. Carpentieri, G. Motta, and J. A. Storer, “Low-complexity lossless compression of Hyperspectral imagery via linear prediction,” IEEE Signal Process. Lett. 12(2), 138–141 (2005). [CrossRef]
35. A. Majumdar, “Causal MRI reconstruction via Kalman prediction and compressed sensing correction,” Magn. Reson. Imaging 39, 64–70 (2017). [CrossRef]
36. U. Sümbul, J. M. Santos, and J. M. Pauly, “A practical acceleration algorithm for real-time imaging,” IEEE Trans. Med. Imaging 28(12), 2042–2051 (2009). [CrossRef]
37. I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Comm. Pure Appl. Math. 57, 1413–1457 (2004). [CrossRef]
38. S. Boyd, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” FNT in Machine Learning 3(1), 1–122 (2010). [CrossRef]
39. A. Majumdar, R. K. Ward, and T. Aboulnasr, “Compressed sensing based real-time dynamic MRI reconstruction,” IEEE Trans. Med. Imaging 31(12), 2253–2266 (2012). [CrossRef]
40. Y.-W. Wen, M. K. Ng, and W.-K. Ching, “Iterative algorithms based on decoupling of deblurring and denoising for image restoration,” SIAM J. Sci. Comput. 30(5), 2655–2674 (2008). [CrossRef]
41. E. M. Eksioglu, “Decoupled algorithm for MRI reconstruction using nonlocal block matching model: BM3D-MRI,” J. Math. Imaging Vis. 56(3), 430–440 (2016). [CrossRef]
42. J. Yang and Y. Zhang, “Alternating direction algorithms for ℓ1-problems in compressive sensing,” SIAM J. Sci. Comput. 33(1), 250–278 (2011). [CrossRef]
43. J. Yang, Y. Zhang, and W. Yin, “A fast alternating direction method for TVL1-L2 signal reconstruction from partial Fourier data,” IEEE J. Sel. Top. Signal Process. 4(2), 288–297 (2010). [CrossRef]
44. Y. Gan, D. Tsay, S. B. Amir, C. C. Marboe, and C. P. Hendon, “Automated classification of optical coherence tomography images of human atrial tissue,” J. Biomed. Opt. 21(10), 101407 (2016). [CrossRef]
45. S. Farsiu, S. J. Chiu, R. V. O’Connell, F. A. Folgar, E. Yuan, J. A. Izatt, and C. A. Toth, “Quantitative classification of eyes with and without intermediate age-related macular degeneration using optical coherence tomography,” Ophthalmology 121(1), 162–172 (2014). [CrossRef]
46. J. P. McLean, Y. Gan, T. H. Lye, D. Qu, H. H. Lu, and C. P. Hendon, “High-speed collagen fiber modeling and orientation quantification for optical coherence tomography imaging,” Opt. Express 27(10), 14457–14471 (2019). [CrossRef]
47. D. Qu, P. J. Chuang, S. Prateepchinda, P. S. Balasubramanian, X. Yao, S. B. Doty, C. P. Hendon, and H. H. Lu, “Micro- and ultrastructural characterization of age-Related changes at the anterior cruciate ligament-to-bone insertion,” ACS Biomater. Sci. Eng. 3(11), 2806–2814 (2017). [CrossRef]
48. D. Mojahed, R. S. Ha, P. Chang, Y. Gan, X. Yao, B. Angelini, H. Hibshoosh, B. Taback, and C. P. Hendon, “Fully automated postlumpectomy breast margin assessment utilizing convolutional neural network based optical coherence tomography image classification method,” Academic Radiol. 27(5), e81–e86 (2020). [CrossRef]
49. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]
50. Z. Wang, E. Simoncelli, and A. Bovik, “Multiscale structural similarity for image quality assessment,” in The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003, vol. 2 (IEEE, 2016), pp. 1398–1402.
51. Y. Gan, T. H. Lye, C. C. Marboe, and C. P. Hendon, “Characterization of the human myocardium by optical coherence tomography,” J. Biophotonics 12(12), e201900094 (2019). [CrossRef]
52. C. P. Hendon, T. H. Lye, X. Yao, Y. Gan, and C. C. Marboe, “Optical coherence tomography imaging of cardiac substrates,” Quant. Imaging Med. Surg. 9(5), 882–904 (2019). [CrossRef]
53. J. Mavadia-Shukla, P. Fathi, W. Liang, S. Wu, C. Sears, and X. Li, “High-speed, ultrahigh-resolution distal scanning oct endoscopy at 800 nm for in vivo imaging of colon tumorigenesis on murine models,” Biomed. Opt. Express 9(8), 3731–3739 (2018). [CrossRef]
54. W. Yuan, D. Chen, R. Sarabia-Estrada, H. Guerrero-Cázares, D. Li, A. Quiñones-Hinojosa, and X. Li, “Theranostic OCT microneedle for fast ultrahigh-resolution deep-brain imaging and efficient laser ablation in vivo,” Sci. Adv. 6(15), eaaz9664 (2020). [CrossRef]
55. L. M. Peterson, M. W. Jenkins, S. Gu, L. Barwick, M. Watanabe, and A. M. Rollins, “4D shear stress maps of the developing heart using Doppler optical coherence tomography,” Biomed. Opt. Express 3(11), 3022 (2012). [CrossRef]
56. K. K. Chu, D. Mojahed, C. M. Fernandez, Y. Li, L. Liu, E. J. Wilsterman, B. Diephuis, S. E. Birket, H. Bowers, G. Martin Solomon, B. S. Schuster, J. Hanes, S. M. Rowe, and G. J. Tearney, “Particle-tracking microrheology using micro-optical coherence tomography,” Biophys. J. 111(5), 1053–1063 (2016). [CrossRef]
57. T. Tang, E. Deniz, M. K. Khokha, and H. D. Tagare, “Gaussian process post-processing for particle tracking velocimetry,” Biomed. Opt. Express 10(7), 3196 (2019). [CrossRef]
58. E. V. Gubarkova, A. A. Sovetsky, V. Y. Zaitsev, A. L. Matveyev, D. A. Vorontsov, M. A. Sirotkina, L. A. Matveev, A. A. Plekhanov, N. P. Pavlova, S. S. Kuznetsov, A. Y. Vorontsov, E. V. Zagaynova, and N. D. Gladkova, “Oct-elastography-based optical biopsy for breast cancer delineation and express assessment of morphological/molecular subtypes,” Biomed. Opt. Express 10(5), 2244–2263 (2019). [CrossRef]
59. B. F. Kennedy, X. Liang, S. G. Adie, D. K. Gerstmann, B. C. Quirk, S. A. Boppart, and D. D. Sampson, “In vivo three-dimensional optical coherence elastography,” Opt. Express 19(7), 6623 (2011). [CrossRef]
60. K. V. Larin and D. D. Sampson, “Optical coherence elastography - OCT at work in tissue biomechanics [invited],” Biomed. Opt. Express 8(2), 1172–1202 (2017). [CrossRef]
61. Y. He, Y. Qu, J. C. Jing, and Z. Chen, “Characterization of oviduct ciliary beat frequency using real time phase resolved doppler spectrally encoded interferometric microscopy,” Biomed. Opt. Express 10(11), 5650–5659 (2019). [CrossRef]
62. H. M. Leung, M. L. Wang, H. Osman, E. Abouei, C. MacAulay, M. Follen, J. A. Gardecki, and G. J. Tearney, “Imaging intracellular motion with dynamic micro-optical coherence tomography,” Biomed. Opt. Express 11(5), 2768–2778 (2020). [CrossRef]
63. S. Bhat, I. V. Larina, K. V. Larin, M. E. Dickinson, and M. Liebling, “4D reconstruction of the beating embryonic heart from two orthogonal sets of parallel optical coherence tomography slice-sequences,” IEEE Trans. Med. Imaging 32(3), 578–588 (2013). [CrossRef]
64. X. Zhao, X. Fu, C. Blumenthal, Y. T. Wang, M. W. Jenkins, C. Snyder, M. Arruda, and A. M. Rollins, “Integrated rfa/psoct catheter for real-time guidance of cardiac radio-frequency ablation,” Biomed. Opt. Express 9(12), 6400–6411 (2018). [CrossRef]
65. X. Yu, R. P. Singh-Moon, and C. P. Hendon, “Real-time assessment of catheter contact and orientation using an integrated optical coherence tomography cardiac ablation catheter,” Appl. Opt. 58(14), 3823–3829 (2019). [CrossRef]
66. C. P. Fleming, H. Wang, K. J. Quan, and A. M. Rollins, “Real-time monitoring of cardiac radio-frequency ablation lesion formation using an optical coherence tomography forward-imaging catheter,” J. Biomed. Opt. 15(1), 1–3 (2010). [CrossRef]
