3-D compressed sensing optical coherence tomography using predictive coding

: 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%.


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][2][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][14][15]. Many CS approaches for OCT have been proposed that aim to reconstruct the raw interferogram or other hardware-specific signals [16][17][18][19][20][21][22]. Studies combining CS-OCT with alternate imaging methods like full-depth OCT have also been explored [23][24][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][28][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.

Methods
ℓ 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 ∈ R N and y ∈ R M as the full-resolution and undersampled images, respectively. The signal x can be recovered from y by solving the objective function where Ψ is the sparse representation basis (e.g. DFT) and A is an M × 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 ∇x = x t − x t−1 . Using the framework proposed in [39], Eq. (1) can be modified to solve for ∇x as follows. argmin 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 ℓ 1 regularization on the difference image ∇x, it was applied to the denoised version through function D(x, λ) where λ is a denoising parameter. In this case the objective function becomes the following function.
∇x t is in the ℓ 2 term but D(∇x t ) in the ℓ 1 term so it will not be possible to solve using the same approach as in Eq. (2). Instead, the method proposed by Wen, et al. for denoising image restoration [40] can be used to solve Eq. (3) by decoupling the objective function into two subproblems which are solved in an alternating fashion. The subproblems of Eq. (3) as defined as follows.
The iteration begins with an initial guess ∇x i−1 and is penalized to agree with the observation y. The second equation controls the sparsity of the solution via the ℓ 1 norm. Noting that the OCT image is sparse when denoised and transformed to the Fourier basis, ΨD(∇x) is penalized rather than ∇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(∇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 ∇x. Taking the derivative gives A is a matrix of "spikes" corresponding to the sampled a-lines of a given b-scan. We observe that A H A + αI is a diagonal matrix. In particular, the matrix A H A is only non-zero at the diagonal elements k ∈ K that correspond to the sampled entries of ∇x. From this assumption, the solution can be written In this formulation, α is a rough measure of the noise in observation y where α = 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 prox λ ∥ · ∥ 1 is evaluated by performing element-wise soft-thresholding of the argument. The soft-thresholding operation soft() of matrix element u i by threshold β is defined for complexvalued entries as sign(u i ) max(|u i | − β, 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 λ. 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, λ has two values (λ 1 , λ 2 ) which represent the vertical and horizontal standard deviation of the 2-D Gaussian Filter D(·, (λ 1 , λ 2 )). The variability of λ is controlled by setting λ max and λ min such that λ decreases logarithmically over J iterations. λ max and λ min may be set differently for λ 1 and λ 2 . The algorithm is structured to update over an inner and outer iteration. The inner iteration solves subproblems for a fixed value of λ until the update reaches max iteration I or the solution update becomes small (set using τ), while the outer loop iterates J times over λ.

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 × 32 square pixel patches of each b-scan. We empirically tested square patches of different sizes and found 32 × 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 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.

Sampling
In this manuscript, "compression rate" η 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 η 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 η which is a function of the full-resolution b-scan interval I b . If a full-resolution b-scan is acquired every ten b-scans in the volume, then I b = 10. The true compression rate η is calculated as follows (10) 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.

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.

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 ℓ 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.

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 α = 0.1, β = 1, λ max = [3,4], λ min = [0.2, 0.4], J = 20, I = 20, and convergence threshold τ = 10 −3 . We used a filter size of 7 × 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 ρ = 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, Ψ is the level-3 Haar wavelet, the anisotropic TV measure is used, and the remaining parameters are set to µ = 10 4 , β = 20, τ = 0.5, and γ = ( √ 5 + 1)/8.

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 µm axial and 15 µ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 µm, 4.4 µm, 7.1429 µm, 6.3 µm, and 6.7 µm, respectively. All datasets had equal spacing between b-scans except for the retina dataset which had 67 µm b-scan spacing. All OCT volumes were cropped to 512 × 800 × 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].

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 Relative Error = ∥x recon − x∥ 2 ∥x∥ 2 (11) where x is the vectorized original OCT volume and x 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 × 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). 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. 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 η a = 50%, 25%, and 10% were tested using a full-resolution interval 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 η 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. 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 η a = 50%, 25%, and two full-resolution intervals 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 × 800 × 800 pixels and we used a patch size of 32 × 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 η includes the full-resolution b-scans so it is higher than the a-line sampling rate η 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 η. For example, in the case of η a = 25% with staggering on, the relative error improves from 0.28 to 0.27 when I b is lowered from 50 to 10, but at the expense of raising η from 26.5% to 32.5%.

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 µ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

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. 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.
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 η a = 50%, 25%, and 10% using staggering and 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.

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 η a = 50%. Staggering and 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.  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.

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 λ max , λ 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][59][60], cilia and mucus movement [61,62], developmental biology [63], and Radio-Frequency Ablation (RFA) [64][65][66].

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.
Disclosures. The authors declare that there are no conflicts of interest related to this article.