1. Introduction

Phase-sensitive optical coherence tomography (PhS-OCT) is an emerging branch of optical coherence tomography (OCT), in which the depth-resolved phase angles are obtained in addition to the amplitudes from the reconstructed complex object functions. The retrieved phase angles could be used to sense picometer-level axial displacements via spectral domain phase microscopy (SDPM) [1], to provide better contrast in optical coherence tomography angiography (OCTA) [2–4], to detect the axial flow velocities in Doppler OCT (D-OCT) [5, 6], or to characterize the tissue stiffness through optical coherence elastography [7–9].

Towards the clinical translation of these functional OCT variants, the imaging speed of the system is highly valued in order to reject the motion artifacts [10, 11], facilitate lateral oversampling [12], or generate 4D time-lapse/time-resolved datasets [13, 14]. These make the swept-source (SS) implementation of OCT more favorable than its spectral domain (SD) counterpart due to the much higher A-line acquisition rate [15]. In fact, recent advancements in swept laser technologies have further widened the gap; the sweeping rate of the wavelength-swept laser source has been dramatically increased from 2 kHz to 28 MHz in the past decades [16–21].

However, the phase stability, a key performance specification, of SS-OCT is often lower than that of SD-OCT due to the following factors. Firstly, most wavelength-swept light sources used in SS-OCT utilize mechanical motions to rapidly tune the output wavelength, which leads to a poor repeatability of the wavenumber-versus-time scanning curve (k-t curve). Secondly, it is difficult to exactly synchronize the onset of the data acquisition with the beginning of the wavelength sweeping. The mis-synchronization will cause timing jitters. Both factors add multiplicative noises to the measured complex phase angle [22].

To address these issues, Yasuno et al first suggested that a calibration signal obtained in advance from a reference interferometer could be used to remap the OCT signal [23]. This method was later developed into a simultaneously recorded calibration channel by Gora et al [24], and further modified by Shangguan et al [25]. By using this simultaneously recorded calibration signal, both the scanning variability and the timing jitter could be corrected, and it is actively adopted by various research groups in recent years [14, 26–29].

Unfortunately, digitizing an extra channel will impose a hefty levy on the data bandwidth of the system. For example, a typical 200 kHz SS-OCT system could generate ∼800 MB of data every second, supposing that each A-line is digitized by 2,000 points at 16-bit precision. This figure will be quickly scaled up, if a faster laser is used (data rate larger than $>6$ GB/s needs to be handled for MHz light source [15]) or an oversampling scheme is adopted to improve the phase stability [22]. This inevitably forces SS-OCT designers to make a compromise among the imaging speed, data bandwidth, and the phase stability. In fact, the data transfer/record rate is now recognized as one of the major trade-offs in the SS-OCT system design [13].

To address this constraint, one straightforward approach is to compress the data channels utilizing Compressed Sensing (CS). The theory, which has been proposed in [30,31] and developed thereafter, states that a signal that has a sparse structure can be recovered with high precision from a vector of linear projections, even in the case where the size of the vector of projections is significantly smaller than the size of the original signal. A signal is said to have a sparse structure, if it can be written as the sum of very few vectors, in a given basis.

Previously, CS has been studied in the context of SD-OCT [32, 33], where the authors proved that it is possible to reconstruct the OCT images by using only a subset of the charge-coupled device (CCD) camera pixels. Their methods exploited the 1D-Fourier transform of the acquired OCT tissue, but did not take advantage of the 2D redundancy of the layers. Another intriguing approach has been proposed by Lebed et al, where OCT images are acquired by using an under-sampling grid [34]. However, their mathematical formulation is not related to CS, as requirements like matrix incoherence [31] are not met in their work. Other works have combined CS with SD-OCT, such as [35, 36], but their methods remain post-processing techniques instead of being real time. It is worth noting that all the previous CS researches are conducted on the basis of SD-OCT, which is not really useful in practice, since an addressable CCD camera is not commonly available and the data bandwidth of SD-OCT is not significant. The only exception is the work conducted by Mididoddi et al, in which they experimentally demonstrated a CS-based photonic time-stretch optical coherence tomography system [37]. However, their method is highly specialized and not compatible with most existing SS-OCT system.

In this manuscript, we proposed to apply CS on the simultaneously recorded calibration signal in PhS-SS-OCT: the calibration signal is randomly sub-sampled during the acquisition to reduce the data bandwidth, and later recovered based on the CS theory. The required sample numbers for exact recovery features a logarithmic growth when spectral oversampling is used. Moreover, the 2D redundancies of the calibration channel are also exploited to achieve the highest compression ratio. In this way, we could break the trade-off imposed by the digitization of an extra channel without losing the merits of possessing high phase stability. To our best knowledge, this is the first experimental implementation and demonstration of a CS-based acquisition in Fourier domain OCT systems. The proposed scheme is universal and could be applied to most existing SS-OCT system.

2. Compressed Sensing-enabled SS-OCT

2.1. Compressed Sensing theory

The mathematical theory of CS has been introduced independently in 2006 by Candès, Romberg and Tao [30], and Donoho [31]. The goal of this theory is to provide tools and methods for the reconstruction of signals with sparse structures from a small collection of linear measurements. Mathematically, the problem of CS consists of recovering a signal $x\in {ℂ}^M$ that has a sparse structure, from a linear observation $y=\mathrm{\Phi }x\in {ℂ}^P$, where $\mathrm{\Phi }\in {ℂ}^{P\times M}$ is called the measurement operator, even when the number of observations is drastically lower than the size of the signal to recover ($P\ll M$). The problem $y=\mathrm{\Phi }x$ is strongly ill-posed in this context, while the prior information on the sparsity of the signal x is key. More generally, we consider that the signal x has a sparse representation in a known dictionary $\mathrm{\Psi }\in {ℂ}^{M\times L}$, meaning that there exists a vector $s\in {ℂ}^L$ such that $x=\mathrm{\Psi }s$ with very few non-zero coefficients. In this case, Ψ is called the sparsifying transform, and s is said S-sparse if only S of its coefficients are non-zero.

The CS theory shows that, under some constraints [30], it is possible to recover an estimator $\widehat{x}$ of the true signal x from the observation y by solving the following optimization problem:

In addition, note that, in the ideal case of noiseless acquisition, the minimum number M of measurements required for exact reconstruction is given by [38]:

2.2. Direct domain Compressed Sensing

The calibration signal in SS-OCT has two properties that motivated us to develop a CS-based technique, suited for faster SS-OCT acquisition:

First, the single calibration signal, corresponding to a single A-line, is a chirped sinusoidal function [24]. Such signal has a highly localized spectrum. In other words, these signals have sparse Fourier transform (see Fig. 1(b)), which is the best case scenario for direct domain CS acquisition, meaning that such signal can be reconstructed with a high precision from very few randomly acquired samples.

Second, any consecutive calibration signals (corresponding to consecutive A-lines, during the raster scanning of the whole sample) have a very similar appearance. Up to a few exceptions, where notable pixel shifts are observed, the pixel-wise difference of two consecutive calibration signals is equal to a constant. As a consequence, if we consider the whole stack of calibration signals as a calibration B-scan, the resulting image will present a strong redundancy along the B-scan direction (see Fig. 1(c)). This redundancy is exploited in this work, with CS reconstruction results improved when considering calibration signals as two dimensional images.

Most of state-of-the-art CS applications (such as MRI [40], Astrophysics [41], Radar [42], Holography [43]) exploit the natural sparsity of signals in the direct domain (Ψ being equal to the identity matrix, or to a wavelet transform) or of their gradient (Ψ being the total variation operator), while sampling the data in a very incoherent space, such as the Fourier or the Fresnel domain. Here, the direct domain refers to the space in which the data is acquired. We propose in this paper to design an alternative approach to CS, where the acquisition is done in the direct domain, and the sparsifying transform corresponds to the Fourier transform. Mathematically, we define $\mathrm{\Phi }\in {\{0,1\}}^{P,M}$ a random selection matrix, such that, for a calibration signal x, the resultant signal $y=\mathrm{\Phi }x$ is a random sub-sampling of the measures of x. In details, the matrix Φ contains at most 1 non-zero coefficient on each line, for a total amount of P non-zero coefficients, meaning that y is a selection of P samples among the M that define x. In addition, we define $\mathrm{\Psi }={\mathcal{F}}^{-1}$, which represents the inverse discrete Fourier transform (IDFT) operator. Then, the CS theory states that we can recover with high precision an estimator of x from the acquisition y by solving:

Note that, when x represents a group of N_k clock signals instead of only one clock signal, the approach is exactly the same. The operator ${\mathcal{F}}^{-1}$ then denotes the 2D-IDFT instead of the 1D-IDFT, and the selection matrix Φ belongs to ${\{0,1\}}^{N_{k}P\times N_{k}M}$.

Previous works in direct-domain CS include the spectral Compressive Sensing approach from [44] for general frequency-sparse signals, modeled as a combination of a small number of sinusoids. The spectral iterative hard thresholding (SIHT) algorithm using heuristic approximations was introduced in [44, 45]. However, such approach is not directly expandable to two dimensional objects, which is crucial in our approach, as we consider grouping consecutive calibration A-lines for improved results. For the sake of generality, we use the NESTA algorithm [46] to solve Eq. (3) in this work, as it can be easily adapted to both 1D and 2D cases.

2.3. 1D vs 2D approach

Equation (3) gives an estimator of a calibration signal x from only a fraction of its samples, denoted as y. In this work, we propose two different approaches for the acquisition of y. While the mathematical formulation is the same, the algorithmic implementations are different.

The first option, which is standard in the domain, is to partially acquire one calibration A-line, and reconstruct an estimation of it through CS optimisation. This approach is one-dimensional, exploits the sparsity of the IDFT of each of the clock signals (see Fig. 1(b)), and can be solved using a multitude of state-of-the-art techniques, such as [44, 45] for instance. Then, the algorithm is run on the N calibration A-lines, for a complete reconstruction of the calibration B-scan.

The second option requires a technical improvement of the setup. As in this case, the acquired data y represents the sub-sampling of a group of N_k calibration A-lines, where N_k is an integer between 2 and 8 with our set-up. Consequently, this method exploits the sparsity of the 2D-IDFT of the $M\times N_k$ portion of the calibration B-scan. Then, the algorithm is run on the $N/N_k$ groups of calibration A-lines, leading to the reconstruction of the complete calibration B-scan.

2.4. Optimal sampling rate for CS-SS-OCT

Using the notations in the work of [22], we can model the clock signal as a chirped sine function that is sampled regularily over M data points:

 

Fig. 1 Exemplary calibration signals acquired in a typical SS-OCT. (a) 1D visualization of one calibration A-line. (b) 1D Fourier Transform of (a) (logarithmic view). (c) 2D visualization of a calibration B-scan. The inset shows the magnified view, which manifests server jitter. (d) 2D Fourier Transform of (c) (logarithmic view).

Download Full Size | PDF

Ideally, the k-t curve would be linear, leading to the calibration signal I being a perfect sinusoidal function. However, the actual function k of the laser is nonlinear (see Fig. 4(c)). The consequence of the nonlinear k-t curve is that the calibration signal consists in a sine wave with modulated amplitude, to which are added sine waves of lower frequencies. If we take the IDFT of the calibration signal, we observe two peaks corresponding to the highest frequency contained in the signal, of rather thin widths (less than 100 coefficients overall), as displayed on Fig. 1(b). Considering the other Fourier coefficients as noise terms, we deduce that the sparsity S of the clock signal is given by the number of coefficients contained at these two locations. As a consequence, since the mutual coherence $\mu \left(\mathrm{\Phi },\mathrm{\Psi }\right)$ is close to 1 in the case of direct-domain compressed sensing, the sub-sampling rate required for almost exact reconstruction of the clock signal is given by:

3. Methods

3.1. OCT system

The SS-OCT system used in this study was developed on the basis of the one that was reported previously [47]. Two modifications are highlighted: (1) the original light source is replaced by a MEMS-based swept source (HSL-20-100-M-3-S, Santec, Japan), (2) the custom-made Michelson interferometer, used for calibration, is substituted by an integrated Mach Zehnder interferometer (INT-MZI-1300, Thorlabs, USA). The schematic of the system is shown in Fig. 2(a). The system have two output channels: (1) OCT channel and (2) calibration channel as illustrated in Fig. 2(b) and 2(c), respectively.

The center wavelength of the system is 1317.5 nm, the axial resolution is 16.8 μm, and the lateral resolution is 8.77 μm. The 6 dB coherence length is 6.5 mm, and the sensitivity of the system is 105.4 dB.

3.2. Hardware implementation of sub-sampling

The schematic diagram of the proposed hardware-based signal sub-sampling is illustrated in Fig. 2(d). The OCT channel is digitized by one DAQ board (ATS 9373, AlazarTech, Canada) at full speed, while the simultaneously recorded calibration channel is digitized by another DAQ board of the same specifications but at a reduced rate via random sub-sampling. Specifically, a set of 65,536 registers, which is referred as “skip bitmap”, are used to specify which trigger event to be ignored during the acquisition by assigning binary values to the corresponding registers. Whenever a trigger event is skipped, no data will be generated. The skip bitmap could be made any arbitrary pattern and manipulated via the provided software API (ATS-SDK V7.1.4 (C++), AlazarTech, Canada).

 

Fig. 2 (a) The SS-OCT setup. Two channels are digitized: (b) OCT channel and (c) calibration channel. (d) The schematic of the hardware-based sub-sampling. The OCT channel is fully digitized by the first DAQ board, while the calibration channel is sub-sampled by the second DAQ. Both digitized signals are later transferred to the host PC. (e) The detailed reconstruction work flow for calibration B-scans. Depending on the size of the sub-sampling mask, the original B-scan will be segmented and randomly sub-sampled at the acquisition. Each segment is reconstructed by using NESTA algorithm and concatenated to form the reconstructed calibration B-scan.

Download Full Size | PDF

For example, we can load a vector of 4,000 bits, which consists of alternating 1s and 0s, onto the skipping table. The DAQ board thus digitizes the signal when every other triggers are received during the first 4,000 clock periods. This process is repeated for the next 4,000 sampling points and so forth. We thus generate a sub-sampled calibration signal whose down-sampling factor is 2. The process is a segment-by-segment process as illustrated in Fig. 2(e).

The fully sampled OCT signal and the sub-sampled calibration signal are later transferred to the host computer via a Peripheral Component Interconnect Express (PCI-e) Gen. 3×8 bus.

It is worth noting that the two DAQ boards have to be carefully connected and synchronized to avoid extra jitter. Detailed discussions on this topic are provided later in Section 5.1.

3.3. CS signal reconstruction

In order to reconstruct the calibration signal from the sub-sampled data, we solve the CS problem (Eq. (3)) using the NESTA algorithm [46], that has been shared by Becker et al on https://statweb.stanford.edu/ candes/nesta/. We implemented the algorithm using Matlab, on a PC workstation 2.93 GHz quad-core central processing unit (CPU), with 8 GB of RAM.

To be more precise, we used the ${\mathcal{l}}_1$ version of the proposed NESTA algorithm, adapted to the operators Φ and Ψ defined in Section 2.2. In addition, the regularization parameter ε has been set to 0 in our experiments, as we aim at exact reconstruction. Note that the denoising property of the regularization parameter is out of the scope of this manuscript.

Again, this reconstruction process is performed segment-by-segment as presented in Fig. 2(e). With these conditions, the reconstruction time of one calibration A-line is between 10 and 50 ms, depending on the sampling rate, and the sampling strategy.

3.4. Experimental procedures

3.4.1. Reconstruction quality: numerical simulation

We first numerically studied the impact on reconstruction quality of different sub-sampling mask configurations including sub-sampling rate ($P/M$) and mask width (N_k).

Conventional SS-OCT calibration data were obtained at 200 Mega Samples/s (MS/s) and 1.6 GS/s for the numerical simulation. Random sub-sampling bitmaps were computer-generated and were used to digitally mask the calibration signal. The sub-sampled data was then reconstructed using the NESTA algorithm.

To evaluate the reconstruction quality, we use the correlation coefficient $r_{x,\widehat{x}}$ between the original signal x and the reconstructed signal $\widehat{x}$, which is given by

A $r_{x,\widehat{x}}$ close to 1 indicates a good reconstruction.

3.4.2. Phase stability

We then evaluated the scheme within its intended context by using the reconstructed calibration signal to remap and stabilize the corresponding OCT spectrum. A “standardized test” was conducted to measure the phase stability of the proposed system as described below.

We placed a microscope slide (1-mm thick, Microscope Slides, Fisherfinest, USA) under the sample arm, and blocked the reference arm. The interference pattern between the light reflected from the top surface and the bottom surface of the sample was obtained from the DAQ board $\#1$ at 800 MS/s. The sub-sampled calibration signal was recorded by DAQ board $\#2$ according to a predefined mask ($P/M=0.3$, N_k = 1), while the same calibration signal was fully digitized by DAQ board $\#1$ for comparison purposes. 1,000 M-scans were obtained at a fixed sample location.

We later reconstructed the calibration signal from its sub-sampled copy, and used this reconstructed signal to remap the OCT signal. The instantaneous phase angles over time at the peak location of the OCT A-lines were extracted and their standard deviation was calculated.

 

Fig. 3 (Numerical Simulation) Correlation coefficient ${\overline{r}}_{x,\widehat{x}}$ between the reconstructed signal $\widehat{x}$ and the original signal x at (a) 200 MS/s and (b) 1.6GS/s. The solid blue curves in (a) and (b) give the levels that are corresponding to ${\overline{r}}_{x,\widehat{x}}=0.9$. (c) The evolution of ${\overline{r}}_{x,\widehat{x}}$ against sub-sampling rate $P/M$ for different mask width N_k. (d) The required sub-sampling rate for good reconstruction (to achieve ${\overline{r}}_{x,\widehat{x}}=0.9$) fits linearlyagainst $\log M/M$, which agrees with our prediction made in Eq. 5.

Download Full Size | PDF

3.4.3. Proof-of-concept: vibrational frequency test

The system configuration is identical to that reported in Section 3.4.2. We used SDPM to measure the vibrational frequency of the sub-sampling-interval displacements of a sample to verify the system’s capability of conducting sensitive phase measurements. The test arrangement included a piezoelectric actuator (PZS001, Thorlabs, USA), which was driven by a a sinusoidal voltage signal from a function generator (AFG3022C, Tektronix, USA). The frequency of the driving sinusoidal is 10 kHz and the peak-to-peak voltage is 1V.

To showcase the performance of our proposed system, we used two other remapping schemes (Yasuno et al [23] and Shangguan et al [25]) on the same OCT dataset for comparison purpose. In the following sections, the terms pre-measured calibration and full calibration to will be used refer these two schemes, respectively.

3.4.4. Proof-of-concept: Doppler OCT

To further validate the proposed method, an experimental phantom was constructed to mimic blood flow. Intralipid emulsion (Sigma Aldrich, USA) was diluted to a concentration of $0.25\%$ in de-ionized water and stored within an intravenous fluid bag and tubing setup. The tubing was fitted into an irrigation pump (CoolFlow, Biosense Webster, USA) which allowed for precise manipulation of laminar flow rates. Imaging was performed over the short axis of the tubing for flow rates ranging between 1-3 mL/min. Wethen obtained Doppler OCT images in addition to optical micro-angiography (OMAG) [48]. Moreover, we compared resultant OMAG images by using two other remapping schemes (see Section 3.4.3).

3.4.5. Ex vivo blood flow detection in swine heart

Epicardial ablation is an effective tool for the treatment of complex arrhythmias such as post-myocardial infarction ventricular tachycardia (VT). In such procedures, care must be taken to avoid energy delivery proximal to coronary vasculature which can lead to vessel occlusion and thrombus formation [49, 50]. In this ex vivo experiment, we showcase the capability of the proposed system to detect the blood blow, which could potentially be used for real-time monitoring for the ablation.

Fresh swine hearts were obtained from Green Village Packing (Green Village, New Jersey, USA). Ventricular slabs containing epicardial vasculature were resected and cannulated using 900 μm diameter flexible tubing coupled to an irrigation pump (CoolFlow, Biosense Webster, USA). Whole swine blood was then diluted at 1:10 ratio in phosphate buffered saline and perfused through the vessel at 5 mL/min. OCT volumes (2048×256×1024 pixels, 1.2×4.8×2.5 mm³, long×widthesdepth) were acquired. The sub-sampling rate of the calibration channel was set to 30% and the mask width was 1. We then used the OMAG algorithm to visualize the blood flow.

4. Results

4.1. Reconstruction quality evaluation

The sampling rate of the DAQ board was first set as 200 MS/s so that each A-line contained 1,000 sampling points ($M=1,000$). Result is in Fig. 3(a) show pseudo-color maps of ${\overline{r}}_{x,\widehat{x}}$ under different experimental conditions. Taking into account the randomness of the masks, five simulations were conducted for each sampling rate, and the average ${\overline{r}}_{x,\widehat{x}}$ over these five simulations was computed.

The blue solid curve is the contour line that is given by ${\overline{r}}_{x,\widehat{x}}=0.9$, which we take as a baseline in this study. It shows that for a mask that consists of only one A-line, a sub-sampling rate ($P/M$) $>40\%$ is necessary to achieve a result better than the baseline. Generally, ${\overline{r}}_{x,\widehat{x}}$ decreases exponentially with $P/M$. On the other hand, $P/M$ could be reduced to $\sim 15\%$ if mask width N_k is set to 8, thanks to the redundancy within B-scans.

We repeated the experiment for a different sampling rate (1.6 GS/s), and the result is plotted in Fig. 3(b). The same reconstruction quality ${\overline{r}}_{x,\widehat{x}}$ could be achieved at a lower sub-sampling rate/larger mask width: For example, by sub-sampling only 5% of the data (N_k = 8), the averaged correlation coefficient between the original calibration signal and the reconstructed one could be as high as $0.98$. To better convey the idea, we plotted ${\overline{r}}_{x,\widehat{x}}$ against the sub-sampling percentage for four combinations of sampling rates and mask widths in Fig. 3(c).

 

Fig. 4 The proposed CS-enabled calibration method requires less data bandwidth without compromising phase stability. (a) The Full calibration signal. (b) The experimentally obtained sub-sampled calibration signal (red squares) and its reconstruction (black solid line). (c) The k-t curves extracted from the full calibration signal and the reconstructed signal. The error, whose amplitude is less than 0.1% of that of the original signal over the entire spectrum, is also plotted. The histogram of the instantaneous phase angles obtained from the peak location of OCT images by using full calibration signal and reconstructed signal are presented in (d) and (e). The standard deviation of the distribution is 4.53 mrad and 4.49 mrad, respectively.

Download Full Size | PDF

We further conducted experiments for three other sampling rates: 600 MS/s, 800 MS/s and 1 GS/s. We recorded the $P/M$ (N_k = 1) that is needed to bring ${\overline{r}}_{x,\widehat{x}}$ to the baseline (0.9) value at these sampling rates. Based on the prediction made in Eq. (5), sub-sampling rates should be linear with $\log M/M$, which fully agrees with the experimental results as shown in Fig. 3(d).

4.2. Phase stability

An exemplary full calibration A-line is plotted in Fig. 4(a), while the sub-sampled acquisition and reconstruction are given in Fig. 4(b). The reconstructed calibration signal is almost identical to the fully sampled one in spite of a scaling factor, which is caused by the difference between the two DAQ boards.

The extracted k-t curve from the full and the reconstructed signals are illustrated side-by-side in Fig. 4(c). The error, which is four orders of magnitude smaller than the original signal, is also plotted in the same panel in red solid line.

In the last step, we computed the histogram of the phase angle distribution at the peak location. Results, using the full and the reconstructed calibration signal, are shown in Fig. 4(d) and 4(e), respectively. The standard deviation of the measured phase angles is equal to $4.53$ mrad (signal-to-noise ratio (SNR) = 52.65 dB) for full calibration and $4.49$ mrad (SNR = 50.86 dB) for reconstructed calibration. The phase stability of the proposed system is not compromised, while only 30% of the bandwidth for the extra calibration channel is used in this example.

It is worth noting that reconstructed results were numerically corrected before presentation. Details can be found in Section 5.3.

 

Fig. 5 The results of the vibrational test. (a) The instantaneous phase evolution over 2 ms and (b) the vibrational spectra obtained by the proposed system along with the controls. The low frequency region (blue boxed) in (b) is magnified and replotted in the inset. The results of using pre-measured calibration (red dash lines) show excessive phase noise in both time domain and spectral domain, while the others present comparable performances to that of using the full calibration. It is worth noting that we could achieve a lower sub-sampling rate ($P/M$) by using a larger mask width (N_k) without compromising the performance. The measured noise floor for the proposed system is below 10 pm.

Download Full Size | PDF

4.3. Vibrational frequency test

Instantaneous phase angles extracted from the surface of the piezoelectric actuator over a small time frame are plotted in Fig. 5(a). Four different remapping strategies are used: pre-measured calibration, full calibration, CS-1 ($P/M=0.3,N_k=1$) and CS-2 ($P/M=0.15,N_k=8$). Both CS results are very similar to that of the full calibration case, while the pre-calibration result manifests apparent noises. We then conducted a discrete Fourier transform (DFT) on the acquired phase angles, and obtained the frequency response of the actuator which is illustrated in Fig. 5(b). The dominant frequency peak at 10 kHz matches that of the stimulation exactly, and is well preserved across all methods. However, if we take a closer look at the low frequency region of the spectra as shown in the inset, we observe the elevated background noise in the case where no simultaneous calibration is used (red dashed line). On the other hand, both reduced-data-rate schemes perform very well and are comparable to that of the conventional dual-channel scheme. The noise floor of the proposed scheme is below 10 pm, if we exclude the discrete spikes that are caused by the mechanical motion of the imaging platform.

 

Fig. 6 The results of the flow velocity measure. (a) The original OCT image of the flow phantom. Its 1D DFT in fast axis direction is given in (e), where we could spot the Doppler frequency shift due to the presence of the intralipid flow. The resultant OMAG images by using full calibration, pre-measured calibration, and proposed CS-based remapping schemes are given in (b), (c) and (d), respectively. In panel (c), the artifacts due to the timing jitter is pointed out by yellow arrows and yellow box. (f)The computed Doppler image by using proposed CS-based remapping. The averaged (4 adjacent A-lines) depth profile at the red line’s location is plotted in (g).

Download Full Size | PDF

4.4. Doppler OCT

We set the flow velocity of the irrigation pump to be 2 mL/min. The Tygon tubing (Saint-Gobain, France) we used in the experiment has an inner diameter about 0.89 mm and an outer diameter of 1.56 mm. It should be noted that the tube was positioned at an angle of 86.05° to the vertical and the cross-sectional area of the tube is 87.11% smaller than that of the pump. The original OCT image of the phantom is shown in Fig. 6(a), and its DFT spectrum against the fast axis is shown in Fig. 6(e). The DFT spectrum shows a strong Doppler frequency shift in the inside region of the tube. After high-pass filtering the image, we obtain the OMAG image by following the steps described in Wang et al [48].

Similar to what we have done in previous subsection, three different remapping schemes (full calibration, pre-measured calibration, and proposed CS-based calibration) were used and the results are given in Fig. 6(b), 6(c), and 6(d), respectively. The CS-reconstructed results ($P/M=0.3$, N_k = 1) shows good image quality comparable to using the full calibration. On the contrary, pre-measured calibration presents multiple artifacts due to the phase noise.

We lastly computed the Doppler image by using the CS-reconstructed calibration, which is displayed in Fig. 6(f). To further showcase the velocity profile, the average of 4 adjacent Doppler A-lines as indicated by the bold red line in Fig. 6(f) was calculated and plotted in Fig. 6(g). The average measured flow speed (half of the maximum velocity) is 0.81 mm/s, which agrees with the prediction.

 

Fig. 7 The blood flow detection in ex vivo swine left marginal artery. (a) Volumetric OCT image overlaid with corresponding OMAG image. The cross-sectional views in y − z and x − z planes are shown in (b) and (c), respectively. The white scale bars represent for 1 mm in each direction.

Download Full Size | PDF

4.5. Ex vivo blow flow detection in swine samples

The obtained OCT image and the calculated OMAG mapping are overlaid and presented in Fig. 7(a). Representative sectional views in both y (slow scanning axis)$-z$ and x (fast scanning axis)$-z$ planes are given in Fig. 7(b) and 7(c), respectively. The left marginal artery is segmented in OMAG image with high contrast. It should be noted that the artifacts presented in the lower boundary of the vessel is mainly due to the reduced SNR in those regions.

5. Discussion

Both simulations and experiments conclude that the proposed CS-based system is capable of conducting phase-sensitive SS-OCT measurements with reduced acquisition data rate. Hence, phase stability of the system is comparable to that of using full calibration. Here in this section, we provide some side notes to further help the readers to understand the results and their significances.

5.1. Hardware-based framework suitable for future compression on OCT channel

The proposed hardware configuration could be extended to the OCT channel to randomly decimate the signal. A CS-based algorithm, which utilizes 1D sparsity of the A-lines, could be implemented either in real-time [51] or post-processing [33] to reconstruct the OCT images at a reduced data rate. Time-elapse imaging and Doppler OCT could be good candidate applications: considering that the scenes update slowly in either fast or slow axis, we could take advantage of the redundancy in those directions by combining scenes together and conducting a 2D random sub-sampling.

5.2. Nyquist sampling versus Compressed Sensing

The calibration signal shown in Fig. 1 appears to be band-limited, and we theoretically should be able to reconstruct the original signal by just sampling at the theoretical Nyquist frequency. However, the actual calibration signal acquired at the Nyquist frequency suffers from aliasing caused by spectral leakage toward the high frequencies, which will eventually lead to compromised phase calibration of the system [22]. An over-sampling scheme is a plausible approach to alleviate this effect [22, 24, 29, 52]. Intuitively, a higher sampling rate could better localize the true starting time of the wavelength sweeping.

Furthermore, a much lower data rate than the Nyquist rate is reached with the proposed method if the 2D redundancy of the calibration signal is exploited, as presented in Fig. 3. The extra digitization channel implemented within this study was to experimentally show the equivalence of the compressed clock to the fully digitized clock. As shown within our study, by using a 1D mask only 30% of the clock signal will need to be recorded. For a 2D mask (8 A-lines), this is further reduced to 15%.

5.3. Synchronization between the DAQ boards

One of the biggest challenges in this study was to synchronize the two DAQ boards in the experiment, so that both the sub-sampled signal and the fully digitized OCT signal are perfectly aligned. Misalignment between the two channels will cause an extra jitter that deteriorate the phase stability of the system. The reason why we are not using a single DAQ board to simultaneously perform the full digitization of OCT channel and the sub-sampling of the calibration channel is purely technical: the DAQ board we are using currently only supports one mode at a time.

In our setup, we have to first ensure that the two boards are of the exact same specifications and are connected to the same trigger signal through same length of electrical connections. But there is still a random jitter that is smaller than 1 clock period between the two boards due to the random starting phase of their crystal oscillators. In fact, if we remap the OCT signals by using reconstructed calibration signals as described in Section 4.2, we observe a fixed amount of jump in the phase as shown in in Fig. 8(a). The histogram of the uncorrected phase angles, which is given in Fig. 8(b), shows two peaks and two groups of distributions. One group corresponds to the reconstructed calibration that has no jitter, and one corresponds to the missed one clock period.

This jitter would not exist if we have a device that could support both modes at the same time. For example, a dual channel DAQ board with a customized field-programmable gate array (FPGA). Therefore, in this study, we performed the numerical correction by measuring the value between the two peaks and subtract it from the group that has the high value. The corrected results are shown in Fig. 8(c) along with the results by using pre-measured calibration. A side-by-side comparison with the full calibration is given in Fig. 8(d) for a short excerpt. To make the two comparable, we subtracted the random starting phase from the corrected phase.

 

Fig. 8 The mis-synchronization due to the usage of two DAQ boards. (a) The raw phase angles obtained by CS recalibration over time and (b) their histogram. A fixed phase jump could be observed in (a), and the histogram in (b) shows two distinct peak with similar distribution. This phase jump is caused by the one pixel shift during the simultaneous triggering of the two DAQ baords. The numerically corrected phase angles (CS-based calibration) are plotted against that of (c) pre-measured calibration and (d) full calibration. The CS-based calibration shows a significantly superior phase stability to that of the pre-measured calibration and a comparable result to that of the full calibration.

Download Full Size | PDF

5.4. Optimal sampling rate

The discussion in Section 2.4 shows that the minimum sampling rate required for exact reconstruction of a clock signal by using direct domain CS is linear with $\log M/M$, with a constant that is close to the sparsity level of the signal according to Eq. (5). Prediction of the value of this constant S^′ remains an open question, which will be the focus of our next work. Moreover, techniques other than direct domain CS such as dynamic compressed sensing algorithms [53] will be considered for comparison purpose as well.

5.5. Towards real-time processing

As mentioned before in the Methods section, the CS reconstruction of the calibration signal is done in post-processing for the time being. Currently, the CPU-based processing requires tens of milliseconds to perform one A-line, which is orders-of-magnitude slower than the required processing rate for streaming. But thanks to the independent nature of the OCT dataset, it is possible to incorporate parallel computing to accelerate the processing [51]. Moreover, using a user-programmable Field-programmable gate array (FPGA) might further speed up the CS reconstruction [54]. For the next step, we will explore both pathways and strive to implement real-time OCT image reconstruction and visualization.

6. Conclusion

To summarize, we have proposed a hardware-based compressed sensing SS-OCT system, and we demonstrated within a biological sample ex vivo localization of flow within a coronary artery using phase based analysis. Three sets of proof-of-concept experiments were conducted and the proposed system is shown to be capable of conducting highly phase-sensitive tasks at a reduced data sampling rate, which is especially useful for ultrahigh-speed applications.

We have also theoretically studied the SS-OCT calibration signal under the framework of compressed sensing theory. The results could be used as a guideline for future system design and development.