## Abstract

Acquiring three dimensional image volumes with techniques such as Optical Coherence Tomography (OCT) relies on reconstructing the tissue layers based on reflection of light from tissue interfaces. One B-mode scan in an image is acquired by scanning and concatenating several A-mode scans, and several contiguous slices are acquired to assemble a full 3D image volume. In this work, we demonstrate how Compressive Sampling (CS) can be used to accurately reconstruct 3D OCT images with minimal quality degradation from a subset of the original image. The full 3D image is reconstructed from sparsely sampled data by exploiting the sparsity of the image in a carefully chosen transform domain. We use several sub-sampling schemes, recover the full 3D image using CS, and show that there is negligible effect on clinically relevant morphometric measurements of the optic nerve head in the recovered image. The potential outcome of this work is a significant reduction in OCT image acquisition time, with possible extensions to speeding up acquisition in other imaging modalities such as ultrasound and MRI.

©2010 Optical Society of America

## 1. Introduction

Volumetric Optical Coherence Tomography (OCT) is emerging as a dominant medical imaging modality for diagnostic ophthalmology. Morphometric analysis and clinical applications require high resolution OCT images which necessitate dense sampling leading to high scan times. Clinical imaging with commercial OCT systems requires volumetric acquisitions lasting several seconds. During such an OCT scan session, the subject is required to stay still, fixating gaze on a point without blinking. Research applications and morphometric analysis often require higher resolution OCT images which necessitate denser sampling leading to even longer duration for volume acquisition. Long scan times increase chances of unavoidable motion artifacts that can corrupt the image. They also cause discomfort and aggravate soreness particularly for those with eye problems. Reducing the total scan time while preserving image quality is a much desired goal due to benefits of reduced discomfort, reduced motion artifacts and possible higher quality, higher resolution images.

Recent advances have been made in image acquisition speed, and leading edge systems have been demonstrated at line rates which are 4-10× faster (see [1] for example) than commercial systems. However, the American National Standards Institute (ANSI) sets strict limits on the maximum permissible optical exposure in the eye, and image quality decreases as acquisition speed increases [2]. Previous efforts to reduce the scan time have predominantly concentrated on increasing the line rate of the OCT engine [1, 2], and designing multi-channel optical configurations [3], but raster scan acquisition geometry itself places an upper limit on the speed at which the image samples can be acquired. For this reason, it is highly desirable to develop methods that can reduce the number of image samples acquired without degrading the overall image quality. One attractive option is to acquire a subset of the B-scans that make up the 3D volumetric field of view followed by filling in the missing data by some suitable ‘recovery’ method. In this work we demonstrate a principled method in which OCT acquisition of a small subset of the entire image volume followed by exploiting sparsity of image coefficients in a carefully chosen transform domain enables the recovery of missing samples with high fidelity.

Exploiting transform sparsity is motivated by the success of data compression in imaging. It is well known that natural, real life images are highly compressible under various algorithms such as the wavelet transform or the discrete-cosine transform. The popular JPEG-2000 image compression algorithm implements a version of the decimated wavelet transform [4] with great success. During the compression stage, many of the transform coefficients are discarded, and only the ones that capture most of the signal's energy are retained. Compressive Sampling (CS) addresses the issues of reconstructing the signal accurately when only a few of the significant coefficients in a transform domain have been given ([5–7] to name a few). In this report, we employ these results to develop a method that can recover the full 3D OCT image volume given only a sub-sampled volume. The theoretical justification for the choice of subsampling required for high fidelity CS recovery comes from the work of Candès and Romberg [8] which suggests sampling the object at randomly spaced vertical and horizontal lines skipping a large portion of the image during acquisition. This scan pattern, schematically illustrated in Fig. 1, is compatible with typical implementations of OCT scanning hardware, first acquiring a series of B-scans with irregular spacing in the horizontal orientation, followed immediately by an acquisition over the same area with irregularly spaced vertical B-scans. Following this randomly subsampled acquisition we show how CS-based reconstruction using an Iterative Soft Thresholding (IST) technique [9] can ‘recover’ the full 3D image data with high fidelity.

## 2. Problem Formulation

Let **f** denote the 3D tissue image (defined on a 3D Cartesian grid) that is to be imaged by OCT using a typical raster scanning geometry. The image **f** is therefore a 3D matrix with values corresponding to the OCT image intensity at each point in the volume. Let **R** be a restriction mask that denotes a subset of the samples of **f** that are actually acquired. Therefore **R** is also a 3D matrix, consisting of binary values corresponding to the scan pattern. The forward model of image acquisition describes the relationship with the acquired volume, **y**, which is a product of the actual volume data modified by the sampling pattern (restriction mask) [10]. Mathematically, this is **y** = **RMf** + **n**, where **y** is the observed incomplete signal containing a fraction of the full 3D image measurements from **f** and **n** is system noise. Here, the matrix **M** contains the basis of the domain in which **f** is acquired. Since OCT images are typically acquired in 3D space with Cartesian geometry, **M** is the Dirac basis; so **M** := **I**.

Let **S** be a transform domain offering good compressibility i. e. a small number of coefficients **x = Sf** can describe the 3D image **f** in that domain. Given the coefficients **x**, the signal **f** can be reconstructed as **f** = **S**
^{H}**x** where **S**
* ^{H}* is the synthesis operator transforming the coefficients back to the image domain. Given incomplete measurements

**y**of the image

**f**, instead of recovering the unknown image

**f**directly, the key idea from CS is instead to find the unknown coefficients

**x**parameterizing

**f**such that the reconstructed signal

**S**

^{H}**x**followed by known masking

**R**is close to the observed incomplete data

**y**i.e. the

*ℓ*

_{2}mean-squared error ‖

**y**–

**RS**

^{H}**x**‖

_{2}is small. From the many potential solutions to this problem, the particular

**x̃**chosen is the one with smallest

*ℓ*

_{1}norm ‖

**x̃**‖

_{1}, representing a measure of sparsity. Hence, CS-based recovery of sparsely sampled OCT images will be achieved via solving the following optimization problem:

Our experiments show that shift invariant wavelets [4] offer good compressibility of OCT images. Hence, they are a natural candidate for the sparse transform domain (**S**) required to solve the above CS-based recovery. The constrained optimization problem stated above is converted to an unconstrained optimization problem and rewritten as:

where λ is a positive parameter that controls the trade-off between the sparsity imposed on the coefficients and the commensurate data mismatch that would result. The IST solver [9] is an attractive algorithm to use to solve this problem as it only consists of matrix-vector operations and converges to a solution in a timely manner. Using this approach, the iterative update **x** → *T _{λ}*(

**x**+

**A**

*(*

^{T}**y**−

**Ay**)), where

**A**is defined to be

**RS**

*, converges to the solution for a particular value of λ as the number of iterations goes to infinity. In practice, only a finite number of iterations are required to achieve convergence. The final recovered image is given by*

^{H}**f̃**=

**S**

^{H}**x̃**. As a measure of recovery quality, signal to noise ratio (SNR), defined by −20log‖

*f*−

*f*̃ ‖

_{2}/log‖

*f*‖

_{2}can be used with high numbers representing better recovery.

## 3. Numerical Results

A full 3D OCT original image volume (**f**) was acquired by raster-scanning the optic nerve head of a healthy male subject using OCT [13]. The spectrometer based OCT system used operated at a line rate of 20kHz. The axial resolution of the system was nominally 4*μ*m in tissue, and the lateral resolution was 25*μ*m. The cross sectional B-scans consisted of 512 A-scans, and the volume consisted of 400 B-scan elevations over an area of 3.7mm × 1.8mm. An example of a transversal slice, or C-scan, extracted from the OCT volume is shown in Fig. 2(a). To simulate the sparsely sampled OCT image volume (**y**) that would be acquired by our proposed method of skipping horizontal and vertical raster-scan lines, we apply to this image a binary restriction mask **R**. This mask consists of vertical and horizontal lines with irregular (but fixed) placing generated through a random number generator in software. Five levels of restriction mask giving 23, 35, 53, 61, and 75 percent missing data were created, and were applied to the original image volume to generate five different sub-sampled volumes **y**. An example of mask that discards 53% of the original samples is shown in Fig. 2(b) where the red lines denote data that was discarded.

The full 3D OCT image volume for each level of subsampling was recovered using our proposed CS based method employing IST solver (detail: using 6 inner and 40 outer loop iterations). SNR based on recovered images corresponding to 23, 35, 53, 61, and 75 percent missing data was found to be 14.9, 13.0, 10.6, 9.8, and 7.5 db respectively showing decreasing SNR with increasing degradation due to missing data. The SNR for bilinear interpolation results was found to be 11.5, 9.2, 6.1, 4.2 and 3.3 dB, respectively, for the above-described data volumes, showing a consistent 3-4 dB difference between the two interpolation schemes. The CS-recovery result for 53% missing data and the error |**f** – **f̃**| image for this cross-section is shown in Fig. 2(c) and 2(d), respectively. Fig. 3 shows two vertical and horizontal B-scans recovered from a volume with 53% missing data compared with bilinear interpolation of the same B-scans. One can observe a distinct, undesirable “staircase” pattern in the bilinear interpolation results. Figure 4 shows the same typical slice that was recovered from volumes with 23, 35, 61 and 75 percent missing data, from a region where data was not sampled. One can see that the overall image quality begins to deteriorate as more data is discarded, but never the less, the fidelity is still acceptable even for a recovered volume from 75% missing data.

#### 3.1. Validation

We performed experiments to quantify how the degradation in recovered images affected the measurement of anatomically meaningful parameters such as the optic nerve cup shape. Morphometric analysis of changes in the shape of the optic nerve head have been proposed as a possible early diagnostic indicator of glaucomatous change. This was done by manually segmenting the clinically relevant Inner Limiting Membrane (ILM) surface for the reconstructed image volume for each of the five levels of degradation using an internal protocol developed in collaboration with clinical author (PJM). These segmentations were performed on each B-scan slice of every volume. The extracted ILM surface from one of the segmentations of the original 3D image is shown in Fig. 5(a) where the color map represents the vertical height (to better present the 3D nature of the cup). The ILM surface extracted from the recovered volume with 53% missing data, shown in Fig. 5(b), closely resembles that of the original ILM.

The recovered image will be deemed to have passed the fidelity test if the location of the ILM surface segmented from the recovered image falls within the variability of the surface segmentations performed multiple times on the original image volume by the same rater. To quantify the segmentation variability due to the manual rater, the ILM in the original 3D image was segmented thrice, each time by the same rater but without consulting the previous segmentations. Using these, an average ILM surface was created representing baseline or ground truth. The standard deviation of each point's position as observed in the three manual segmentations of the original volume was calculated. The surface shown in Fig. 5(c) is the average ILM where the colormap represents the variance in position in the three manual segmentations of the original image volume.

For each point on the ILM surface segmented from the reconstructed data, the relative position to the average ILM surface (the baseline) was found. This process generates the so-called Topographical Change Analysis (TCA) map, [14], representing the variability in position of the ILM surface from a known baseline/ground truth. Points on recovered ILM surface lying more than one standard deviation away from the average baseline surface are shown as the TCA overlaid on top of the summed voxels projection image of the original volume in Fig. 6. One can see that the quality of the summed voxel projection images begins to deteriorate at the 61% missing data volume, but is still acceptible even at 75% missing data. The red and green colors represent posterior and anterior change in position, respectively. For the volume that was recovered from 23% missing data, nowhere did the segmentations fall outside the one standard deviation regions. For the other cases with more missing data, most points in the image were reconstructed faithfully so that the errors in position exceeded one standard deviation in only a few locations. These were typically in slices with missing data that were located along the steeply sloping part of the ILM surface. Given that only a few locations in recovered surfaces deviated by more than one standard deviation for each of the five cases tested (as an example, <8% points on the recovered ILM surface were farther than one standard deviation in the case of 75% missing data), clinical measurements such as cup volume and area (see Table 1) did not significantly differ from the measurements made on the original volume. The cup volume was measured as the space enclosed by a reference plane near the top of the volume and the manual segmentation of the ILM in the cup. The surface area of the ILM closed by the reference plane and the cup was also calculated. Changes in surface area and volume of the nerve cup from recovered images were assessed to be clinically negligible by the clinical author (PJM).

The reduction in scan time, presented as a function of percent of missing data is shown in Fig. 7. This function is computed by counting the number of B-scans required for a particular level of subsampling and dividing this by the number of B-scans required in a full raster scan. The fly-back time of the galvanometer mirrors is implicitly taken into account into the time it takes to perform one B-scan. The benefits in reduced scan time of the proposed scanning pattern begin to be observable when one discards approximately 25% of the data, or more. Using the subsampling corresponding to 53% or 61% missing data roughly halved the acquisition time.

## 4. Discussion

Compressive Sampling (CS), pioneered by works like [9, 11, 12] to name a few, provides a theoretical basis for the recovery of the full 3D imaging volume from only a subset of the acquired data. Given a subsampled volume where random data is missing and using the idea of transform sparsity, there is usually enough information in the available data to reconstruct the original image volume with an accuracy that is at least as good as the best compressed representation of the object [6].

In this work we have shown that the theory of CS can be gainfully applied to reconstruct a full 3D OCT image with high fidelity from only a subset of acquired OCT samples. This type of recovery has clear advantages to using the standard bilinear interpolation, as shown by results in Fig. 3. An obvious “staircase” pattern is visible in the bilinear interpolation result (fifth row from the top in Fig. 3), which is undesirable as it corrupts the inherently smooth nature of the image, while CS-based recovery preserves the smoothness.

The sub-sampling scheme we have chosen is physically motivated by OCT raster-scanning, and is one where the imaging volume is acquired at specific vertical and horizontal B-scans only, thereby shortening the image acquisition time. The ILM surfaces obtained from manual segmentations of the recovered images do not differ markedly from the surfaces obtained from segmentations of the original image. The differences primarily occur in regions where the ILM cup has the steepest topographical changes. The small number of points and the localized nature of the deviations of the recovered ILM surface from the original surface are found to have negligible effect on downstream morphometric measurements. The change in metrics such as surface area and volume calculated from ‘recovered’ images was found to be insignificant for the clinical assessment of diseases such as glaucoma. Acquisition of fewer image samples allows the reduction in total OCT scan time, which in turn leads to lower chances of motion artifacts and less subject discomfort. For scan patterns corresponding to 53% or 61% missing data, for which the recovery is found to of very high fidelity, the acquisition time is roughly halved providing valuable saving in scan time. Data acquisition sparsity was investigated through the selective acquisition of B-scan images. The number of lines acquired in a B-scan could also be reduced through the acquisition of randomly spaced A-scans. However, due to the limited response time of the galvanometer mirrors, and complexity of coordinating the camera capture with settled mirror position, capturing an entire line is more economical. Other sub-sampling patterns, which may reduce scan time by an even greater ratio without degrading the overall image quality, are a topic of further study.

## 5. Conclusions

We have described a subsampling scan pattern consisting of a sequence of irregularly spaced horizontal and vertical B-scans that works well with CS-based recovery. Fewer total B-scans are required for CS-based recovery of the OCT volume, reducing the acquisition time with minimal degradation in image quality. The reduction in the OCT scan time limits the discomfort of a subject during the scan and reduces the amount of motion artifact. Image reconstruction with CS preserved the smoothness of sharp features in the optic cup, whereas standard bilinear interpolation resulted in staircase effect and loss of detail. The Compressed Sensing scan pattern presented in this report is compatible with existing OCT acquisition hardware, and has the potential to reduce the acquisition time in applications such as retinal imaging where the A-scan rate is limited by the maximum permitted optical exposure [2]. The degradation in the CS-recovered image quality is small, even at 75% missing data, and does not significantly affect clinical measurements like cup surface area and volume.

## Acknowledgments

Funding for this work is generously provided by Michael Smith Foundation for Health Research (MSFHR), Natural Sciences and Engineering Research Council of Canada (NSERC) and Canadian Institutes of Health Research (CIHR).

## References and links

**1. **B. Považay, B. Hofer, C. Torti, B. Hermann, A. R. Tumlinson, M. Esmaeelpour, C. A. Egan, A. C. Bird, and W. Drexler, “Impact of enhanced resolution, speed and penetration on three-dimensional retinal optical coherence tomography,” Opt. Express17, 4134–4150 (2009), http://www.opticsexpress.org/abstract.cfm?URI=oe-17-5-4134. [CrossRef] [PubMed]

**2. **T. Schmoll, C. Kolbitsch, and R. A. Leitgeb, “Ultra-high-speed volumetric tomography of human retinal blood flow,”, Opt. Express 17, 4166–4176 (2009), http://www.opticsexpress.org/abstract.cfm?URI=oe-17-5-4166.

**3. **M. K. K. Leung, A. Mariampillai, B. A. Standish, K. K. C. Lee, N. R. Munce, A. Vitkin, and V. X. D. Yang, “High-power wavelength-swept laser in Littman telescope-less polygon filter and dual-amplifier configuration for multichannel optical coherence tomography,” Opt. Lett.34, 2814–2816 (2009), http://ol.osa.org/abstract.cfm?URI=ol-34-18-2814 [CrossRef] [PubMed]

**4. **S. Mallat, *A Wavelet Tour of Signal Processing, Second Edition* (Academic Press, New York, 1999).

**5. **E. J. Candès and D. L. Donoho, “New tight frames of curvelets and optimal representations of objects with piecewise-C^{2} singularities,” Comm. Pure Appl. Math.57, 219–266 (2004), http://www.acm.caltech.edu/∼emmanuel/papers/CurveEdges.pdf. [CrossRef]

**6. **E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory **52**, 489–509 (2006). [CrossRef]

**7. **D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory **52**, 1289–1306 (2006). [CrossRef]

**8. **E. Candès and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Probl.23, 969–985 (2007), http://stacks.iop.org/0266-5611/23/969. [CrossRef]

**9. **I. Daubechies, M. Defrise, and C. De Mol “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Comm. Pure Appl. Math.11, 1413–1457 (2004), http://dx.doi.org/10.1002/cpa.20042. [CrossRef]

**10. **F. J. Herrmann and G. Hennenfent, “Non-parametric seismic data recovery with curvelet frames,” Geophys. J. Int. **173**, 233–248 (2008). [CrossRef]

**11. **M. Holschneider, R. Kronland-Martinet, J. Morlet, and P. Tchamitchian, *Wavelets, Time-Frequency Methods and Phase Space* (Springer-Verlag, Berlin, 1989).

**12. **M. Elad, J. L. Starck, P. Querre, and D.L. Donoho, “Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA),” Appl. Comput. Harmon. Anal.19, 340–358 (2005), http://www.sciencedirect.com/science/article/B6WB3-4GWC29F-2/2/61d7afc314d50b27968d84ff4a16acce. [CrossRef]

**13. **M. Young, S. Lee, E. Gibson, K. Hsu, M. F. Beg, P. J. Mackenzie, and M. V. Sarunic, “Morphometric analysis of the optic nerve head with optical coherence tomography.” In proceedings of OCT and Coherence Domain Optical Methods in Biomedicine XIV7554, (2004), http://link.aip.org/link/?PSI/7554/75542L/1.

**14. **B. C. Chauhan, J. W. Blanchard, D. C. Hamilton, and R. P. LeBlanc, “Technique for Detecting Serial Topographic Changes in the Optic Disc and Peripapillary Retina Using Scanning Laser Tomography,”, Invest. Ophthalmol. Vis. Sci. **41**, 775–782 (2000). [PubMed]