Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Optical coherence tomography using physical domain data compression to achieve MHz A-scan rates

Open Access Open Access

Abstract

The three-dimensional volumetric imaging capability of optical coherence tomography (OCT) leads to the generation of large amounts of data, which necessitates high speed acquisition followed by high dimensional image processing and visualization. This signal acquisition and processing pipeline demands high A-scan rates on the front end, which has driven researchers to push A-scan acquisition rates into the MHz regime. To this end, the optical time-stretch approach uses a mode locked laser (MLL) source, dispersion in optical fiber, and a single analog-to-digital converter (ADC) to achieve multi-MHz A-scan rates. While enabling impressive performance this Nyquist sampling approach is ultimately constrained by the sampling rate and bandwidth of the ADC. Additionally such an approach generates massive amounts of data. Here we present a compressed sensing (CS) OCT system that uses a MLL, electro-optic modulation, and optical dispersion to implement data compression in the physical domain and rapidly acquire real-time compressed measurements of the OCT signals. Compression in the analog domain prior to digitization allows for the use of lower bandwidth ADCs, which reduces cost and decreases the required data capacity of the sampling interface. By leveraging a compressive A-scan optical sampling approach and the joint sparsity of C-scan data we demonstrate 14.4-MHz to 144-MHz A-scan acquisition speeds using a sub-Nyquist 1.44 Gsample/sec ADC sampling rate. Furthermore we evaluate the impact of data compression and resulting imaging speed on image quality.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Throughout the past decade, optical coherence tomography (OCT) has proven to be a versatile tool in medical diagnostics allowing, for example, straightforward assessment of the progress of macular degeneration, multiple sclerosis, and glaucoma [13]. Using the spectrally-dependent interference of two light waves, OCT interrogates depth information of a sample, including many different tissues, and can be scanned to collect a three-dimensional($3$D) data cube of an object’s structure. The massive amount of data collected in a volumetric OCT acquisition poses challenges for data throughput, storage, and manipulation, typically requiring data compression prior to any data processing [4,5]. Additionally, in many applications the signal capture should be performed quickly to avoid motion artifacts that can distort the image [6,7].

Without the need for a scanning reference arm, spectral domain OCT (SD-OCT) techniques have increased A-scan rates from a few kHz to a few MHz [8,9] by using CCD or CMOS cameras to sample and digitize the OCT spectral interference. However, techniques that utilize cameras face limitations from the electronic data readout rate, as well as motion blur due to exposure time [10]. To avoid cameras, swept source OCT techniques using Fourier domain mode locked (FDML) sources and dispersed super-continuum sources have been used to produce a spectrum-to-time mapping of the interference signal that can be sampled serially with a single ADC [11,12]. The use of mode locked lasers (MLL) and optical dispersion have become the state of the art in terms of speed, providing stable pulses with extremely large bandwidths, demonstrating A-scan rates from 7.14 MHz to 90.9 MHz [13,14]. However, CCD/CMOS cameras, swept source methods, and dispersed MLL methods all read out the spectral interference signal serially using different methods of parallel to serial conversion. This introduces a hard limit on the total image voxel rate of the OCT system to the sampling rate of available ADC technology.

In the past decade the information theory community introduced the concept of compressed sensing (CS) [1519], suggesting that the sparsity of natural signals can be leveraged to reduce the number of samples required to capture signals of interest. This has been adopted by the medical field and applied to MRI imaging, photo-acoustic imaging, and OCT [2022]. Particularly, extensive work on under sampling OCT data followed by implementing CS algorithms has demonstrated successful reconstruction with less than 20% of measurements required by the Shannon/Nyquist theory [22,23]. By taking advantage of the compressibility of volumetric OCT data [24], the data cube can be under sampled without the loss of image quality [25]. However many of these methods still require the entire data cube to be recorded at the Nyquist rate followed by digital under sampling after acquisition [2225]. Although compressed reconstruction after data acquisition allows for real time visualization, such digital domain sub-sampling fails to address the physical domain bottleneck at the signal acquisition stage due to serial sampling limits of an ADC.

We have previously developed an architecture to implemented CS at high speed in the optical domain using optical signal processing in a technique we termed Continuous High-Rate Photonically Enabled Compressed Sensing (CHiRP-CS) and have leveraged this architecture for high-speed flow microscopy, ultra wide-band radio frequency (RF) sensing, and preliminary work on OCT [2629]. For OCT we directly record optically computed inner products between the interference signals and know binary patterns such that each ADC sample contains information spanning the entire A-scan depth profile. This allows for real-time optical domain data compression of the OCT signal prior to detection by a high-speed ADC. Consequently, we can remove the limits imposed by Nyquist sampling on A-scan rates and thus decouple imaging rate from ADC sampling rate. Similar recent work, investigated the compressive acquisition and reconstruction of A-scans at 1.51-MHz rates using a 66% compression ratio [30]. Here, we demonstrate A-scan acquisition rates from 14.4-MHz to 144-MHz using compression ratios from 26% to 2.6%, respectively. Such low compression ratios are achieved by implementing compressive sampling in the axial A-scan dimension followed by joint reconstruction of the entire volumetric C-scan data utilizing the multi-dimensional sparsity of the full 3D signals of interest.

2. Experiment

2.1 CHiRP-CS sampling system

The CHiRP-CS optical system is shown in the lower half of Fig. 1. The approach begins with ultrafast laser pulses from a MLL source at a center wavelength of 1550 nm and native 90-MHz repetition rate. Pulses from this laser are first temporally multiplexed twice up to a repetition rate of 360 MHz and then sent through a 853 ps/nm dispersion optical fiber module to spread the 22 nm of optical bandwidth over 20 ns. This resulting spectrum-to-time mapping allows the spectrum of the laser pulses to be amplitude modulated in time with an EOM operating at 11.52 Gb/s using a pre-programmed pseudo-random binary sequence (PRBS). With the 2.77-ns repetition period of the laser pulse train and the over 20-ns chirped pulse duration, neighboring pulses overlap greatly, resulting in up to eight pulses modulated simultaneously. However this modulation is imparted on different portions of their spectra resulting in unique spectral patterns on each pulse. Modulating the 360-MHz pulse source at 11.52 Gb/s results in sequential pulses spectrally shifting the PRBS pattern by 32 bits. In order to reach our final 1.44-GHz repetition rate, the patterned source is temporally multiplexed twice more, with large delays of 166 ns and 306 ns. The dispersed pulses are re-compressed in 50 km of SMF resulting in a sequence of temporally distinct short pulses with unique spectral patterns at a 1.44-GHz repetition rate. This approach results in minimum spectral features of 12.53 GHz, over almost 3 THz of optical bandwidth, defined by the dispersion profile mapping of the EOM modulation rate onto optical bandwidth. Notably, given the 11.52-Gb/s modulation rate the 12.53-GHz minimum spectral features are near the modulation limit from interbit spectral distortion [31]. Subsequently passing the spectrally-coded pulses through the OCT interferometer piece-wise multiplies the binary spectral pattern with the OCT interference signal of interest. Finally, detecting the temporally compressed pulse after this process achieves optical integration. Thus this process optically computes the vector inner product between the binary spectral patterns and the OCT interference signal such that only a single ADC sample of each pulse is required to capture each compressed measurement.

 figure: Fig. 1.

Fig. 1. Experimental setup for conventional time-stretch MHz OCT is shown on top. A 90-MHz MLL is pulse picked down to a 18-MHz repetition rate and dispersed to over 8 nanoseconds using SMF. This is sent into the OCT interferometer and the returned pulses are detected with a 20-GHz balanced photo-detector and digitized at 40 Gsamples/s. Our CHiRP-CS MHz OCT system is shown at the bottom. Pulses from a 90-MHz MLL are dispersed in DCF, spectral encoded with a PRBS using an EOM, then temporally compressed in SMF. The pulses are temporally multiplexed four times, before and after the modulation, for a final 1.44-GHz repetition rate. The pulses are sent into the OCT interferometer and detected with a 1.6-GHz balanced photo-detector and digitized at 1.44 Gsamples/s. MLL - mode-locked laser, SMF - single mode fiber, DCF - dispersion compensating fiber, EOM - electro-optic modulator, PRBS - psudeo-random binary sequence, BPD - balanced photo-detector.

Download Full Size | PDF

2.2 OCT interferometer system

The OCT interferometer is constructed with a mirror as a reference arm and a two-dimensional laterally scanning sample arm with a single 75-mm focal length aspheric lens. The input CHiRP-CS source is split by a 80/20 coupler, where 80 percent enters the sample arm and 20 percent reaches the reference mirror. The returned pulses are then recombined in a 50/50 coupler and detected in a balanced configuration by a 1.6-GHz amplified photo-detector. As illustrated in Fig. 1 a single sample of the compressed pulse amplitude yields the inner product between the spectral interference signals and the unique binary spectral pattern imposed on each pulse by the CHiRP-CS system. Each inner product contains information spanning the entire spectrum and can be described mathematically as

$$y = \langle {\boldsymbol {a}}, {\boldsymbol {s}} \rangle + z,$$
where ${\boldsymbol {a}}$ is one row in the known pseudo-random binary pattern ${\boldsymbol {A}}$, and ${\boldsymbol {s}}$ is the interference signal under test a sample location, and $z$ denotes the noise. We can then reconstruct the depth profile using a sequence of such compressed measurements (vector ${\boldsymbol {y}}$) and compressed sensing algorithms necessitating far fewer measurements than a comparable Nyquist sampling system such as a time-stretch system shown in the top of Fig. 1. Specifically, here we reconstruct the OCT spectra onto 384 spectral pixels using 10 to 100 samples yielding compression ratios from 2.6% to 26%. Practically, this compression allows us to use a much lower speed detector and ADC than the comparable time-stretch system, while achieving the same imaging speed or, alternately higher imaging speed using comparable electronics.

In order to investigate the fidelity of our CS approach, we first acquire Nyquist-sampled time-stretch OCT measurement for comparison [13,14]. As shown on top in Fig. 1, this time-stretch system is implemented as follows. The 90-MHz MLL is sent into an EOM to be pulse picked down to 18 MHz. Pulse picking is necessary to avoid pulse overlap after the pulses then propagate through 853 ps/nm dispersion module to achieve sufficient spectrum-to-time mapping for comparable axial dimension to our CS system. This signal detected by a 20-GHz linear balanced photo-detector, then digitized by a 20-GHz bandwidth 40-GS/s oscilloscope. This time-stretch measurement is used as the ground truth in comparison to our compressed sensing approach to evaluate reconstructed image quality.

2.3 CS reconstruction

To scan a $3$D object ${\boldsymbol {S}}$ of size ${N_1} \times {N_2} \times {N_3}$, the compressed measurements are taken along the third, A-scan dimension. We generate ${N_1} \times {N_2}$ number of measurements for a C-scan, with each measurement of length ${{M}}$. Let ${\boldsymbol {y}}^{({i},{j})} \in {\mathbb {R}}^{{{M}}}, {i} = 1 ,2 ,\ldots , {N_1}, {j} = 1, 2,\ldots , {N_2}$ denotes the measurements vector collected from the signal at a spatial location $({i}, {j})$ : ${\boldsymbol {s}}^{({i},{j})} \in {\mathbb {R}}^{{{M}}}$, and ${\boldsymbol {A}}^{({i},{j})} \in {\mathbb {R}}^{{{M}} \times {N_3}}$ is the sensing matrix associated with that measurement vector. For the reconstructed cube in frequency domain, these dimensions are: ${N_1} = 100 , {N_2} =150$ and ${N_3} = 384$. The recovered image is of size $100 \times 150 \times 192$. The depth dimension is half of the one in the frequency domain since along the depth direction, the final image is calculated on the energy from both positive and negative frequencies. In the recovery algorithm, the number of measurements per line ${{M}}$ ranges from $10$ to $100$, corresponding to compression ratios ${{M}}/ {N_3}$ from $2.6\%$ to $26\%$. The mathematical model of the system in matrix form is simply:

$${\boldsymbol {y}}^{({i},{j})} = {\boldsymbol {A}}^{({i},{j})} {\boldsymbol {s}}^{({i},{j}) } +{\boldsymbol {z}}^{ ({i}, {j}) } ,$$
where ${\boldsymbol {z}}^{ ({i}, {j}) }$ represents the noise vector added to the $({i}, {j})$-th noiseless measurement. The noise vector is generated from multiple sources such as the linearization approximation error of the system, noise from the data collection process, as well as pre-processing error, etc.

With OCT signals, we know the spectral modulation frequencies correspond to different depths of the object. Additionally, in many OCT applications, the object is sparse in the number of reflected layers, resulting a small number of cosine tones along the compressed A-scan dimension. Therefore, if we represent the signal in frequency domain, we would expect the coefficients to be sparse. Accordingly, the inverse discrete Fourier transform matrix is used as the sparsifying basis in the recovery algorithm.

The combination of noise processes erode the compressed measurements, having a larger impact on higher frequency due to the reduced signal strength resulting from electro-optic related bandwidth limitations of the modulated spectral patterns. Consequently, the high frequency parts of the signal become less distinguishable from the noise. While applying sparse recovery with the classic $\ell _1$ min norm, we observe that as the sparsity level $\lambda$ is increasing, the power of the high frequency part vanishes much faster than the low frequency part. As a result, there is a bias toward the low frequency interference depths. To resolve this issue, we employ a weighted $\ell _1$ minimization to this problem. The general weighted $\ell _1$ minimization method was proposed in [32].

Our proposed method is composed of two steps. First, we generate a initialization from the classic $\ell _1$ minimization problem. The weight vector ${\boldsymbol {w}}$ is then computed based on the power spectrum of the initialization solution.

$$\widetilde{{\boldsymbol {X}}_0} = {\mathop{\textrm{arg min}}\limits_{{\boldsymbol {X}}_0 \in {\mathbb {C}}^{{N_3} \times{N_1}{N_2}}}} \frac{1}{2} \sum\limits_{k=1}^{{N_1}{N_2}} \| {\boldsymbol {y}}^{({k})} - {\boldsymbol {A}}^{({k})} {\boldsymbol {\Phi}} {\boldsymbol {x}}^{({k})} \|_2^{2} + \lambda_0 \sum\limits_{{k} = 1}^{{N_1}{N_2}} \| {\boldsymbol {x}}_0^{({k})}\|_{1},$$
where $\| {\boldsymbol {x}} \|_{1}$ calculates the sum of absolute values of all entries of ${\boldsymbol {x}}: \| {\boldsymbol {x}} \|_{1} = \sum _{i= 1 }^{N_3} |x_{i}|$ and $\lambda _0 > 0$ is the sparsity balancing ratio. The non-negative regularization parameter $\lambda _0$ balances the ratio of sparsity of the solution and the fitness of the solution with respect to the measurements. The larger value of sparsity level leads to a more sparse solution.

Second, we pass the $\ell _1$ minimization solution to the weighted $\ell _1$ minimization algorithm. For non-negative weights ${\boldsymbol {w}} = (w_1, w_2,\ldots , w_{N_3})^{T}, w_i \geq 0$, the weighted $\ell _{1}$ norm of vector ${\boldsymbol {x}}$ given the weight vector ${\boldsymbol {w}}$ is defined as: $\| {\boldsymbol {x}}\|_{{\boldsymbol {w}} , 1} = \sum _{i=1}^{N_3} w_i | x_i |$. For variable matrix ${\boldsymbol {X}} = ({\boldsymbol {x}}^{(1)}, {\boldsymbol {x}}^{(2)},\ldots , {\boldsymbol {x}}^{( {N_1} {N_2})} )$, we find:

$$\widetilde{{\boldsymbol {X}}} = {\mathop{\textrm{arg min}}\limits_{{\boldsymbol {X}} \in {\mathbb {C}}^{{N_3} \times{N_1}{N_2}}}} \frac{1}{2} \sum\limits_{k=1}^{{N_1}{N_2}} \| {\boldsymbol {y}}^{({k})} - {\boldsymbol {A}}^{({k})} {\boldsymbol {\Phi}} {\boldsymbol {x}}^{({k})} \|_2^{2} + \lambda \sum\limits_{{k} = 1}^{{N_1}{N_2}} \| {\boldsymbol {x}}^{({k})} \|_{{\boldsymbol {w}},1},$$
where the sparisfying transform ${\boldsymbol {\Phi }}$ is the inverse discrete Fourier basis. The weight ${\boldsymbol {w}}$ in (4) is a function associated with the power spectrum of $\widetilde {{\boldsymbol {X}}_0}$ from solving (3). Smaller weights encourage higher amplitudes by reducing the contribution in the penalty function. We estimated the support based on the amplitudes of the $\ell _1$ minimization solution. To suppress noise on non-support locations and hence enhance the image quality, the weight vector ${\boldsymbol {w}}$ is set to be small for locations in the support and large for non-support locations. Essentially, our choice of ${\boldsymbol {w}}$ enforces joint sparsity along the depth dimension. Within the support, the recovered locations with smaller amplitudes have smaller weights in order to boost the energy of the reconstructed signal in high frequency bands. More details on determining the weight vector are explained in the Section 3.

The above two optimization problems (3)(4) can be solved efficiently by several methods such as proximal gradient descent methods [33,34], gradient projection [35], alternating minimization [36], approximate message passing [37], etc. In our implementation, we use gradient projection for sparse reconstruction (GPSR) [35] to solve both $\ell _1$ and weighted $\ell _1$ minimization programs.

The full image reconstruction process is illustrated in Fig. 2, where the inputs are the compressed measurements ${\boldsymbol {y}}$ and psudeo-random binary pattern ${\boldsymbol {A}}$, and the output is the reconstructed image $\widetilde {{\boldsymbol {X}}}$. The iterative $\ell _1$ minimization algorithm using GPSR is fed the measurements and the known binary patterns to reconstruct an initial image. This is sent into a weighted $\ell _1$ minimization algorithm that produces the final image.

 figure: Fig. 2.

Fig. 2. The full process for reconstructing an image, with an initial GPSR reconstruction is used to seed a weight minimization algorithm that produces the final image, illustrated in Fig. 3(b-c).

Download Full Size | PDF

3. Experimental results

We image a test sample composed of an empty PDMS microfluidic channel mounted on a glass slide and acquire a 1-mm by 1.5-mm transverse scan with a 10-um resolution resulting in a $100 \times 150 \times 384$ OCT data cube. To visualize the reconstruction result, we determine the $3$D power spectrum of the recovered spectral modulation signal by calculating the magnitude for each discrete spectral modulation frequency from $\widetilde {{\boldsymbol {X}}}$ by summing the squared coefficients of positive and negative frequencies at each corresponding location. For our 384 reconstructed spectral pixels, this procedure results in 192 axial pixels in the depth image. The total run time for the algorithm took a few hours to recover the entire 3D cube, which is roughly 25 A-scan lines per minute. If we parallelize the recovery procedure for various groups of A-scan lines, the recovery time can be reduced to a few minutes.

The $3$D C-scan of the reference signal captured using time-stretch OCT is shown in Fig. 3(a), in which the axial dimension spans from bottom to top. The initialization of the compressively acquired and reconstructed C-scan from $\ell _1$ minimization is displayed in Fig. 3(b). The solution from $\ell _1$ minimization detects the locations of the layers fairly accurately however we lose a lot of energy towards the top of the C-scan corresponding to the region of high frequency spectral modulation. Also, the energy for the middle slice with curved channel is aliased to these higher frequency bands.

The $1$D power spectrum statistics along the third dimension are used to determine the weights for weighted $\ell _1$ minimization algorithm to resolve these issues. First, we search for peaks from $\ell _1$ minimization solution to determine the support and non-support part of the signal. We use three parameters to determine the peaks: the threshold ratio $\tau$ between $0$ to $1$ on amplitudes to be in the support, the minimum width $\delta \in {\mathbb {Z}}^{+}$ to be a strong peak, and minimum radius from the peak boundaries to the local maximum (peak centers) $\rho \in {\mathbb {Z}}^{+}$. We first eliminate the support locations based on the first two parameters. Then for the detected peaks, we find the local maximum within each local region, treat it as the peak center and extend its boarder to $2\rho$ if the boarder to center distance is less than $\rho$, encouraging a smoother transaction between support and non-support part in the solution. In our experiment, we pick $\tau = 15 \%, \delta = 3, \rho = 2$. We find peaks within top $15\%$ magnitudes with width at least $3$ pixels and radius at least $2$ pixels from the peak center.

 figure: Fig. 3.

Fig. 3. C-scan recovery from $80$ compressed measurements, or an 18-MHz A-scan rate. a) The 18-MHz time-stretch OCT image which we use as the ground truth reference. b) The output of the GPSR reconstruction algorithm at the mid-stage of our image reconstruction process. c) The final C-scan reconstruction from our improved weighted minimization algorithm showing the distinct third layer of the image.

Download Full Size | PDF

Next we compute the weight vector ${\boldsymbol {w}}$ according to the result from peak detection. Within each peak region, the square root of the averaged power are used as weights. Since the $\ell _1$ minimization recovery solution has a higher energy at a lower frequency peak and a lower energy at a high frequency peak, this method assigns a lower weight for a high frequency peak compared to a lower one. For the non-support part, we set their weights to be a constant higher than all the weights for the support. Here, we set them to be two times the square root of the max value of the $1$D power spectrum vector. All parameters in the peak finding and weight functions are determined empirically, but can be set automatically or tuned after data collection. The higher the weight, the magnitude of corresponding location tends to be zero because the weight contributes more to the penalty. Our designed weight vector significantly decreases the reconstruction noise and boosts the energy in high frequency bands.

The ground truth image in Fig. 3(a) shows three distinct layers, representing the air-glass interface, the glass-PDMS interface, and PDMS-air interface from bottom to top, or low to high spectral modulation frequencies, respectively. The PDMS-glass low contrast interface is interrupted by the air gap channel which is smaller than the axial resolution of the system, creating the clear channel in the center layer of the image. Figure 3(b) illustrated the intermediate result from $\ell _1$ minimization, which has severe aliasing on non-support parts and the top layer corresponding to the highest frequency band is missing. Using our CHiRP-CS pulse source matched to the A-scan rate of the time-stretch system (18 MHz), the final reconstruction Fig. 3(c) shows we can clearly reconstruct three distinct layers of the microfluidic channel with only a fraction of the samples required by Nyquist at a low ADC sampling rate of 1.44 GSample/s. In comparison the 18-MHz A-scan rate achieved with time-stretch MHz OCT requires a detector and ADC over a decade faster. Specifically, a 40 GSample/s ADC is employed here for our reference image, Fig. 3(a).

One intriguing aspect of compressive acquisition is the ability to alter the compression rate and correspondingly imaging speed without modifying the hardware system. Figure 4(a) shows four example reconstruction results with different number of measurements per A-scan, resulting in different A-scan rates. For these reconstructions the samples are chosen at random from the total measurements collected per A-scan, random sample selection is analogous to implementing all stages of temporal multiplexing after the electro-optic modulator with sufficient delay. The resulting image quality for A-scan rates of 144 MHz, 48 MHz, 28.8 MHz and 14.4 MHz corresponding to the reconstructions with $10$, $30$, $50$, and $100$ measurements, respectively. The transverse dimensions are scanned in a 100 by 150 grid, and the depth information is reconstructed onto $192$ axial pixels. For the reconstructions with $10$ and $30$ measurements, there are two and one layers of the image missing respectively.

 figure: Fig. 4.

Fig. 4. a) Example C-scan reconstructions of an 100 x 150 x 192 depth image with 10, 30, 50, and 100 measurements, or 144-MHz, 48-MHz, 28.8-MHz, and 1.44-MHz A-scan rates, respectively. b) The PSNR of the CS reconstruction vs the number of compressed measurements used for reconstruction shows an increase in PSNR around 50 measurements where the third layer becomes clearly visible.

Download Full Size | PDF

The impact of compression rate on image quality is illustrated in Fig. 4. As shown in Fig. 4(b), when the number of measurements increases, the algorithm is able to recover more details of the signal and the reconstruction image quality is improved visually and as quantified by Peak Signal-to-Noise Ratio (PSNR). We calculate PSNR by comparing the reference Nyquist-sampled time-stretch measurement of the section of microfluidic channel to our CS reconstructions. We can see the increase in PSNR from $30$ to $50$ measurements corresponds to properly reconstructing the third layer. The $100$ measurement reconstruction shows some noise reduction, but the increase in PSNR is minimal. Although reconstructions using $10$ and $30$ samples fail to reconstruct all of the layers, there is still a clear object reconstructed. This is due to the sparsity enforced by the compressed reconstruction, meaning the sparser the image, the fewer samples need to be used to form an image. The axial and lateral resolutions of both systems are defined by the same optical bandwidth and free space optics, approximately 42 um and 21 um respectively. Although for the compressed axial dimension, the sparsity can affect the image quality. As we see in Fig. 4, the reduction in image quality using fewer samples correlates with fewer of the layers reconstructed in the CS system but maintains the resolution of the layers.

4. Discussion and conclusion

In this paper we demonstrate a compressed sensing OCT system that both addresses the need for high-speed data acquisition and offers physical-domain data compression. Using this approach we show real-time data compression of the OCT signals allowing A-scan acquisition from deeply sub-Nyquist sets of measurements. Additionally, we develop a tailored reconstruction approach that leverages the joint sparsity in both axial and transverse dimensions to efficiently and accurately reconstruct the full C-scan data cube. Experimentally, we show successful C-scan reconstruction of all layers in a microfluidic channel on a glass slide with only $13\%$ of the measurements required by Nyquist sampling, or a 28.8-MHz A-scan rate. Our results demonstrate that the conventional limit on axial pixel rate imposed by the ADC sampling rate in conventional Nyquist-sampling OCT systems can be surpassed using sub-Nyquist CS sampling resulting in faster collection times, less data, and lower cost detectors and ADCs.

Over the last few decades OCT imaging has moved from the research to practice, making OCT systems a stable medical technique. Compressed sensing can take advantage of the natural image sparsity and breadth of previous knowledge for flexible and fast imaging with less expensive hardware. The number of samples required to reconstruct depends on the sparsity, so more complex objects (i.e. those with more layers) will require more samples and take longer to collect. However, the sparsity of any given sample can be increased by using a basis that better represents the image. We use the DCT basis in this paper as a general transform that is effective for most OCT images, but other general transforms or more image-specific transforms can be used to reconstruct various classes of OCT images. Simply, the more a priori information of the sample, the CS system can collect data faster and reconstruct higher quality images.

Funding

National Science Foundation (ECCS- 1609693, ECCS-1254610).

Acknowledgments

Portions of this work were presented at CLEO in 2016, SM2I-1.

References

1. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography-principles and applications,” Rep. Prog. Phys. 66(2), 239–303 (2003). [CrossRef]  

2. P. H. Tomlins and R. K. Wang, “Theory, developments and applications of optical coherence tomography,” J. Phys. D: Appl. Phys. 38(15), 2519–2535 (2005). [CrossRef]  

3. J. S. Schuman, C. A. Puliafito, J. G. Fujimoto, and J. S. Duker, Optical coherence tomography of ocular diseases (Slack, 2004).

4. K. Zhang and J. U. Kang, “Real-time 4D signal processing and visualization using graphics processing unit on a regular nonlinear-k Fourier-domain OCT system,” Opt. Express 18(11), 11772–11784 (2010). [CrossRef]  

5. I. Grulkowski, M. Gora, M. Szkulmowski, I. Gorczynska, D. Szlag, S. Marcos, A. Kowalczyk, and M. Wojtkowski, “Anterior segment imaging with Spectral OCT system using a high-speed CMOS camera,” Opt. Express 17(6), 4842–4858 (2009). [CrossRef]  

6. S. H. Yun, G. J. Tearney, J. F. De Boer, and B. E. Bouma, “Pulsed-source and swept-source spectral-domain optical coherence tomography with reduced motion artifacts,” Opt. Express 12(23), 5614–5624 (2004). [CrossRef]  

7. K. Zhang, W. Wang, J. Han, and J. U. Kang, “A surface topology and motion compensation system for microsurgery guidance and intervention based on common-path optical coherence tomography,” IEEE Trans. Biomed. Eng. 56(9), 2318–2321 (2009). [CrossRef]  

8. R. Leitgeb, C. Hitzenberger, and A. Fercher, “Performance of Fourier domain vs. time domain optical coherence tomography,” Opt. Express 11(8), 889–894 (2003). [CrossRef]  

9. D. Choi, H. Hiro-Oka, K. Shimizu, and K. Ohbayashi, “Spectral domain optical coherence tomography of multi-MHz A-scan rates at 1310 nm range and real-time 4D-display up to 41 volumes/second,” Opt. Express 3(12), 3067–3086 (2012). [CrossRef]  

10. R. A. Leitgeb, L. Schmetterer, W. Drexler, A. F. Fercher, R. J. Zawadzki, and T. Bajraszewski, “Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography,” Opt. Express 11(23), 3116–3121 (2003). [CrossRef]  

11. R. Huber, D. C. Adler, V. J. Srinivasan, and J. G. Fujimoto, “Fourier domain mode locking at 1050 nm for ultra-high-speed optical coherence tomography of the human retina at 236,000 axial scans per second,” Opt. Lett. 32(14), 2049–2051 (2007). [CrossRef]  

12. S. Moon and D. Y. Kim, “Ultra-high-speed optical coherence tomography with a stretched pulse supercontinuum source,” Opt. Express 14(24), 11575–11584 (2006). [CrossRef]  

13. K. Goda, A. Fard, O. Malik, G. Fu, A. Quach, and B. Jalal, “High-throughput optical coherence tomography at 800 nm,” Opt. Express 20(18), 19612–19617 (2012). [CrossRef]  

14. J. Xu, C. Zhang, K. K. Y. Wong, and K. K. Tsia, “Megahertz all-optical swept-source optical coherence tomography based on broadband amplified optical time-stretch,” Opt. Lett. 39(3), 622–625 (2014). [CrossRef]  

15. E. J. Candes and T. Tao, “Decoding by linear programming,” IEEE Trans. Inf. Theory 51(12), 4203–4215 (2005). [CrossRef]  

16. E. J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]  

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

18. R. G. Baraniuk, “Compressive sensing,” IEEE Signal Process. Mag. 24(4), 118–124 (2007). [CrossRef]  

19. E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]  

20. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med. 58(6), 1182–1195 (2007). [CrossRef]  

21. Z. Guo, C. Li, L. Song, and L. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt. 15(2), 021311 (2010). [CrossRef]  

22. 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–22019 (2010). [CrossRef]  

23. D. Xu, Y. Huang, and J. U. Kang, “Volumetric (3D) compressive sensing spectral domain optical coherence tomography,” Biomed. Opt. Express 5(11), 3921–3934 (2014). [CrossRef]  

24. A. B. Wu, E. Lebed, M. V. Sarunic, and M. F. Beg, “Quantitative evaluation of transform domains for compressive sampling-based recovery of sparsely sampled volumetric OCT images,” IEEE Trans. Biomed. Eng. 60(2), 470–478 (2013). [CrossRef]  

25. M. Young, E. Lebed, Y. Jian, P. J. Mackenzie, M. F. Beg, and M. V. Sarunic, “Real-time high-speed volumetric imaging using compressive sampling optical coherence tomography,” Biomed. Opt. Express 2(9), 2690–2697 (2011). [CrossRef]  

26. B. T. Bosworth and M. A. Foster, “High-speed ultrawideband photonically enabled compressed sensing of sparse radio frequency signals,” Opt. Lett. 38(22), 4892–4895 (2013). [CrossRef]  

27. B. T. Bosworth, J. R. Stroud, D. N. Tran, T. D. Tran, S. Chin, and M. A. Foster, “High-speed flow microscopy using compressed sensing with ultrafast laser pulses,” Opt. Express 23(8), 10521–10532 (2015). [CrossRef]  

28. B. T. Bosworth, J. R. Stroud, D. N. Tran, T. D. Tran, S. Chin, and M. A. Foster, “Ultrawideband compressed sensing of arbitrary multi-tone sparse radio frequencies using spectrally encoded ultrafast laser pulses,” Opt. Lett. 40(13), 3045–3048 (2015). [CrossRef]  

29. J. R. Stroud, B. T. Bosworth, D. N. Tran, T. D. Tran, S. Chin, and M. A. Foster, “72 MHz A-scan optical coherence tomography using continuous high-rate photonically-enabled compressed sensing (CHiRP-CS),” In CLEO: Science and Innovations (Optical Society of America, 2016) paper SM2I-1.

30. 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]  

31. J. R. Stroud, B. T. Bosworth, D. N. Tran, T. P. McKenna, T. R. Clark, T. D. Tran, S. Chin, and M. A. Foster, “Continuous 119.2-GSample/s photonic compressed sensing of sparse microwave signals,” In CLEO: Science and Innovations (Optical Society of America, 2015) paper STh4F-2.

32. M. A. Khajehnejad, W. Xu, A. S. Avestimehr, and B. Hassibi, “Weighted l1 minimization for sparse recovery with prior information,” IEEE International Symposium on Inf. Theory, ISIT 2009, iaw023 (2017). [CrossRef]  

33. P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,” Multiscale Model. Simul. 4(4), 1168–1200 (2005). [CrossRef]  

34. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. on Imag. Sci. 2(1), 183–202 (2009). [CrossRef]  

35. M. A. T. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems,” IEEE J. of Sel. Topics in Signal Process. 1(4), 586–597 (2007). [CrossRef]  

36. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” FNT in Machine Learning 3(1), 1–122 (2010). [CrossRef]  

37. D. L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing,” Proc. Natl. Acad. Sci. U. S. A. 106(45), 18914–18919 (2009). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)

Fig. 1.
Fig. 1. Experimental setup for conventional time-stretch MHz OCT is shown on top. A 90-MHz MLL is pulse picked down to a 18-MHz repetition rate and dispersed to over 8 nanoseconds using SMF. This is sent into the OCT interferometer and the returned pulses are detected with a 20-GHz balanced photo-detector and digitized at 40 Gsamples/s. Our CHiRP-CS MHz OCT system is shown at the bottom. Pulses from a 90-MHz MLL are dispersed in DCF, spectral encoded with a PRBS using an EOM, then temporally compressed in SMF. The pulses are temporally multiplexed four times, before and after the modulation, for a final 1.44-GHz repetition rate. The pulses are sent into the OCT interferometer and detected with a 1.6-GHz balanced photo-detector and digitized at 1.44 Gsamples/s. MLL - mode-locked laser, SMF - single mode fiber, DCF - dispersion compensating fiber, EOM - electro-optic modulator, PRBS - psudeo-random binary sequence, BPD - balanced photo-detector.
Fig. 2.
Fig. 2. The full process for reconstructing an image, with an initial GPSR reconstruction is used to seed a weight minimization algorithm that produces the final image, illustrated in Fig. 3(b-c).
Fig. 3.
Fig. 3. C-scan recovery from $80$ compressed measurements, or an 18-MHz A-scan rate. a) The 18-MHz time-stretch OCT image which we use as the ground truth reference. b) The output of the GPSR reconstruction algorithm at the mid-stage of our image reconstruction process. c) The final C-scan reconstruction from our improved weighted minimization algorithm showing the distinct third layer of the image.
Fig. 4.
Fig. 4. a) Example C-scan reconstructions of an 100 x 150 x 192 depth image with 10, 30, 50, and 100 measurements, or 144-MHz, 48-MHz, 28.8-MHz, and 1.44-MHz A-scan rates, respectively. b) The PSNR of the CS reconstruction vs the number of compressed measurements used for reconstruction shows an increase in PSNR around 50 measurements where the third layer becomes clearly visible.

Equations (4)

Equations on this page are rendered with MathJax. Learn more.

y = a , s + z ,
y ( i , j ) = A ( i , j ) s ( i , j ) + z ( i , j ) ,
X 0 ~ = arg min X 0 C N 3 × N 1 N 2 1 2 k = 1 N 1 N 2 y ( k ) A ( k ) Φ x ( k ) 2 2 + λ 0 k = 1 N 1 N 2 x 0 ( k ) 1 ,
X ~ = arg min X C N 3 × N 1 N 2 1 2 k = 1 N 1 N 2 y ( k ) A ( k ) Φ x ( k ) 2 2 + λ k = 1 N 1 N 2 x ( k ) w , 1 ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.