Introduction

Since its initial development more than two decades ago [1], optical coherence tomography (OCT) has matured into an important clinical imaging technique. This is not surprising given the numerous advantages it offers including cellular (i.e., micron-level) axial resolution, noninvasive nature, cost-effectiveness, and the ability to be incorporated into catheters [2] and endoscopes [3]. Unlike many other microscopic methods such as fluorescence microscopy [4], OCT does not require exogenous labeling for visualizing biological structures. Recently there has been a growing interest in applying OCT to pre-clinical studies that involve various small animal models. Chicken embryos have been an important animal model for investigating cardiac morphological changes and fluid dynamic stresses [5–7]. A better understanding of these forces could provide deeper insights into remodeling biology of cardiovascular tissues and the mechanisms towards the normal/abnormal development of the embryo’s cardiovascular system. While various imaging modalities, including micro computed tomography (Micro-CT) and ultrasound [8,9], have been employed to image chicken embryos, none of them could compete with OCT in terms of spatial resolution [10,11]. In our previous study, detailed three-dimensional (3D) morphological information of chicken embryo’s aortic arches (AA) was captured using a commercial OCT system, which provided lateral and axial resolutions of around 10 µm [12]. Based on the 3D model, computational fluid dynamic simulations were conducted to estimate the blood flow pattern and the destruction of induced stress. However, the static 3D AA model was obtained using a fixed sample due to the limited imaging speed of the commercial OCT system (up to 20 B-scans per second). For in vivo multi-dimensional imaging of a beating heart, it is desirable to have an imaging speed at least two orders of magnitude higher.

The imaging speed of OCT has improved significantly over the last decade. Early OCT systems employed time-domain (TD) detection in which an optical beam in the sample arm of an interferometer is scanned transversely across the sample and the interference signal is recorded with a photodetector. The axial scan rate of TD-OCT is limited to several kilohertz due to mechanical movement of the reference mirror. TD-OCT also suffers from insufficient sensitivity for accurate diagnosis [13]. Although increasing the optical power may improve the sensitivity, clinical systems are limited in power due to safety regulations that make high-speed TD-OCT with axial scan rates beyond 10 kHz impractical [14]. Spectral domain (SD)-OCT eliminates the axial mechanical scanning. It offers significant advantages in imaging speed and detection sensitivity [15–18]

In existing OCT systems configured for cross-sectional imaging, A-lines are recorded sequentially when the sample beam is shifted to scan along the transverse direction (i.e., point scanning). While a galvanometer is usually employed as the scanning device, the typical resonant frequency of high-end galvanometers is limited to a few kilohertz. It is therefore very challenging for existing OCT systems to go beyond a few thousand B-scans per second. A video-rate volumetric image acquisition was demonstrated by Wieser et al. with a multi-megahertz swept source based OCT setup [19]. However, the reported 25ms volume acquisition time was made possible by the use of complicated multi-spot scanner optics and a high-end laser source. In addition, the acquisition speed of point-scanning OCT systems is limited by the maximum allowed power for a minimum required detection sensitivity, which in turn is proportional to the applied illumination power [20]. For example, the maximum A-scan rate was estimated to be 1 MHz to achieve the minimum sensitivity requirement of 95 dB for clinical imaging at 1060 nm center wavelength [21].

An alternative approach to point scanning OCT is line-scan OCT in which the sample is illuminated with a focal line [22,23]. Since the A-lines are recorded in parallel rather than sequentially, a cross-sectional image (i.e., B-scan) can be recorded in a single shot without the need for 2D scanning mechanism. This significantly relaxes the system’s technical requirements and allows higher maximum permissible exposure (MPE) levels. Equivalent A-line acquisition speeds of 256 kHz and 1 MHz have been reported using line-scan SD-OCT and line-scan swept-source OCT (SS-OCT), respectively [20,24]. Although line-scan SS-OCT systems exhibit reduced scattering and signal crosstalk due to the temporal separation of spectral sample points, their acquisition rates are limited by the maximum laser sweep rates.

One of the main problems with line-scan OCT system designs reported thus far is the non-uniform field of view (FOV) due to the utilization of a cylindrical lens, which results in an illumination line with a Gaussian power distribution.

In this paper, we report a line-scan SD-OCT system for ultrahigh-speed imaging, utilizing a commercially available supercontinuum light source and a high-speed 2D CMOS camera. Unlike previously reported line-scan SD-OCT systems, the illumination line is generated using a Powell lens to obtain a highly uniform FOV. We utilize phase stepping in the reference arm and differential reconstruction to eliminate the autocorrelation interference. An advanced noise reduction algorithm based on nonlocal means (NLM) is applied in the post-processing stage to significantly enhance the morphological details and suppress motion artifacts. We demonstrate in vivo imaging of chicken heart embryo at acquisition rates up to 5,000 B-scans per second with sensitivities > 90 dB and free-space axial resolution <2.0 µm. To the best of our knowledge, this is the fastest line-scan OCT system reported thus far.

2. Methods

2.1 OCT setup

The experimental setup of the proposed OCT system is illustrated in Fig. 1. A supercontinuum laser (SC-Pro from YSL Photonics) is used as a light source to emit a broad output spectrum covering the entire region from 450 to 2400 nm. The light source is followed by a band-pass filter to select the wavelength range between 660 and 840 nm (centered at 750 nm). This wide spectrum allows the system to achieve high axial resolutions (< 2.0 µm in air). A 6°-angle Powell lens is used to fan out the collimated laser beam (beam diameter ~2.0 mm) as illustrated in the dashed yellow lines in Fig. 1(a). The Powell lens resembles a round prism with a curved roof line. This curved roof is a complex 2D aspheric curve that generates a great amount of spherical aberration which redistributes the light along the line by decreasing the intensity in the center while increasing it at the ends [25]. The Powell lens is followed by an achromatic lens L1 with a focal length (f.l.) of 30 mm to focus the group of diverged beamlets into an array of points forming a uniform and symmetric illumination line. Unlike the conventional approach that employs a cylindrical lens, the use of Powell lens results in a highly uniform wavelength-independent illumination line. A simple telescope, consisting of two achromatic lens L2 (60 mm f.l.) and L3 (100 mm f.l.), relays the illumination and provides a magnification of the illumination line from 0.6 to 1.0 cm. As the laser output is linearly polarized, a half-wave plate (HWP) is used to rotate the polarization to a vertical direction so that most of the source beam is transmitted through the polarization beam splitter (PBS), which passes the vertically polarized light but reflect the horizontally polarized light. The power is then split with a 90/10 beam splitter. The larger fraction of power is sent towards the sample arm where an achromatic lens LS (30 mm f.l.) is used to create a vertical illumination light sheet focusing on the sample with a total incident optical power of 15 ± 1 mW and FOV of 2.0 mm to enable the acquisition of depth-resolved cross-sectional images or B-scans in a single shot. To obtain 3D images, the illumination line is scanned transversely across the sample using a one-dimensional galvanometer scanner yielding a transverse FOV up to 2.6 mm in the scanning direction. A similar setup is used in the reference arm as shown in Fig. 1. However, the reference mirror is mounted on a piezoelectric actuator driven by a square wave at a frequency identical to half the frame rate to realize a π-phase shift between consecutive measured spectra as shown in Fig. 2. This reference beam phase stepping is needed for the elimination of autocorrelation artifacts using differential reconstruction. Compared to an electro-optic phase modulator (EOM), the piezo actuator is compact and can be easily integrated into the system. LabVIEW software (National Instruments) was used to synchronize the movements of the galvanometer scanner and piezo actuator with the data acquisition. That was achieved by sharing the start trigger signal of the analog output from a data acquisition (DAQ) device (USB-6212, National Instruments) controlled by the LabVIEW software.

 

Fig. 1 Schematic diagram of the proposed SD-OCT system: (a) side view and (b) top view. The red dotted line illustrates the backscattered beams. BPF: band-pass filter; BS: beam splitter; CMOS: complementary metal-oxide semiconductor; L: achromatic lens; PA: piezo actuator; PL: Powell lens; PBS: polarized beam splitter; RM: reference mirror; SCL: super-continuum laser; λ/2: half wave plate; and λ/4: quarter wave plate. The table inset summarizes the focal lengths of the utilized achromatic lenses. The photo inset shows the chicken embryo and the orientation of the illumination line on the sample (dashed yellow line).

Download Full Size | PDF

 

Fig. 2 (a) Sensitivity roll-off, and (b) axial resolution measured as a function of the imaging depth at a tomographic frame rate of 3.5 kHz. The error bars are based on sample standard deviation of five repeated measurements.

Download Full Size | PDF

The backscattered light from the sample (i.e., sample beam) interferes with the reflected beam from the reference mirror (i.e., reference beam). A quarter wave plate (QWP) rotates the polarization of reference and sample beams by 90 degrees so that both beams are efficiently reflected by the PBS. The combined beam from the reference and sample arms is then directed towards the slit, which rejects stray light such as surface reflections from optical components and bright specular reflections at the surface of the egg yolk. The transmittance through the slit contains an array of combined beamlets, which are spatially distributed vertically. Each combined beamlet results from the interference between a reference beamlet and a sample beamlet corresponding to one A-line. The depth-dependent intensity of the sample beamlet is encoded into the spectrum of the combined beamlet. A near-infrared (NIR) reflectance grating (GR25-1208, Thorlabs) disperses each combined beamlet to project the wavelength-dependent intensity onto the corresponding row of pixels on the image sensor. GR25-1208 is a 750 nm blaze wavelength reflective diffraction grating. Its maximal efficiency around 750 nm is above 70%.

A 16-bit ultra-fast CMOS camera (PHOTRON FASTCAM SA7) is used to capture the spectrum array. This camera can record images at a speed up to 3,500 Hz with a full image resolution of 1,280 by 1,024 pixels, or up to 30,000 Hz with reduced image resolutions. A raw image captured by the camera contains up to 1024 interferograms, each of which has 1280 spectral points. The acquired spectra are further processed using a program written in MATLAB for image reconstruction and processing.

2.2 Autocorrelation artifact removal

The internal sample reflectivity profile can be estimated from the line-to-line IFT of the wavenumber-dependent detector current described by [26],

As seen in Eq. (1), only the cross-correlation terms carry information about the sample morphology. After IFT, the DC and autocorrelation terms results in artifacts that could severely obscure the real structures. To eliminate the undesired terms, we utilized a differential reconstruction approach [15]. This approach takes advantage of the fact that all cross-correlation terms depend on the reference mirror position while the remaining autocorrelation terms do not. Therefore, the unwanted autocorrelation terms can be eliminated by measuring one additional spectral fringe pattern with a phase shift $\varphi =\pi$ introduced into the reference arm

2.3 NLM algorithm for image sequence denoising

Denoising using simple time averaging generally leads to motion artifacts and blurring for fast moving image sequences, which in turn, makes motion compensation a requirement in nearly all state-of-the-art movie filters. However, motion estimation remains a challenging task due to the ambiguity of trajectories. NLM is a recently developed denoising method that can circumvent the motion estimation problem [27–29]. Since the initial work reported by Buades et. al. in 2005 [28], NLM has become increasingly popular in image sequence denoising.

We modified the NLM algorithm from its original form to reduce the processing time. Traditional NLM algorithms denoise a pixel by involving all pixels in the image sequence. In our case, however, we assume that the structural changes and dynamic deformation are confined within a certain range. Therefore, our NLM algorithm denoises a pixel by taking the mean of all pixels in the corresponding cropped image sequences.

If $NLMu(x,y,k)$ denotes the image intensity at a pixel $(x,y)$ in an image frame $k$, the denoising can be applied using the following filter:

where C is a normalization factor and $\Re$ defines a cropped area centered at the pixel $(x,y)$. ${\sigma }_k$ determines the first exponential weighting function associated with the intensity similarity between pixels $(x,y,k)$ and $(x\text{'},y\text{'},i)$. The second exponential weighting function takes into account the distance between the pixel pair and X is a tuning parameter to control the exponential decay with regard to their spatial proximity. The parameter M defines the range of image frames around k that contribute to the nonlocal means. A total number of 2M neighboring frames are involved in the filtering process.

For OCT image sequences acquired at 10 volumes per second and 500 frames per volume, we chose $M=25$. M was reduced to 10 for the image sequences acquired at 40 volumes per second. $\Re$ and X were both set to 5 pixels in all cases, while the parameter ${\sigma }_k$ was adjusted dynamically according to the average image intensity.

2.4 OCT system characterization

The sensitivity roll-off of the system measured as a function of imaging depth is shown in Fig. 2(a). At a frame rate of 3,500 Hz, the sensitivity was >80 dB from 0 to 0.55 mm and remained >70 dB up to the maximum depth (0.85 mm). The depth-dependent sensitivity drop is attributed to the optical resolution limits of the spectrometer, finite pixel width, reduced fringe visibility at high spatial frequencies, and inter-pixel crosstalk in the spectrometer [30].

Figure 2(b) shows the full-width half-maximum (FWHM) axial resolution measured at various imaging depths. At a minimum depth, the axial resolution was measured to be 1.6~2.0 μm. This value is slightly higher than the theoretical prediction of 1.4 μm [26]. As the imaging depth increases in Fig. 3(b), the axial resolution deteriorates marginally but remains below 3.0 μm at depths ≤0.75 mm. However, a higher degradation of resolution is observed at larger depths, and the largest recorded axial resolution was 5.0 μm at 0.854 mm. This behavior is attributed to the sensitivity of the interpolation process to small errors at high fringe frequencies, resulting in resolutions poorer than the transform-limited resolution at the maximum depth [14]. These errors can be minimized by utilizing a linear array with more pixels. Alternatively, the depth range can be increased by retrieving the missing phase information from the acquired spectrum [31,32], or by employing a multiplexing technique e.g., using optical frequency comb so that the spectrum sampling function can be determined by the comb rather than detector pixel distribution [33].

 

Fig. 3 Measured line profile using Powell lens (red) and cylindrical lens (blue).

Download Full Size | PDF

The lateral resolution was measured using a standard 1951 USAF bar target. The smallest elements visualized using the target were measuring 6.2 µm in width. This value is close to the theoretical lateral resolution estimation of 5.56 μm [26].

To demonstrate the advantage of using Powell lens to generate the illumination line in our system, the illumination intensity distribution was measured along the transverse direction with a glass slide placed in the sample arm. From the recorded 2D intensity matrix, a representative column (along the dimension orthogonal to spectral dispersion) was selected and plotted in Fig. 3 (red curve). For comparison purposes, the Powell lens was replaced by a cylindrical lens and the respective intensity profile was also plotted in Fig. 3 (blue curve). The blue curve shows a Gaussian power distribution with the intensity at the center more than ten times stronger than at the periphery of the FOV, which by definition is transverse scanning range. The FOV is significantly more uniform when Powell lens is used to generate the illumination line as indicated by the red curve in Fig. 3.

2.5 Samples and preparation

The imaging performance of the developed OCT system was validated via a series of in vivo imaging experiments using fertilized chicken embryos. Chicken embryos were chosen for this study because they develop similarly to human embryos and therefore serve as a good animal model system. Fertilized Bovans Goldline chicken eggs were obtained from a local farm and incubated with the blunt-end facing upward in a forced draft incubator at 37.5°C and 80% humidity for a total of 72 hours (HH18 to HH20) of embryonic development according to the defined Hamburger and Hamilton (HH) system [34], which is based on somite counts and morphological features of the developing embryo. At stage HH18, the heart has a tubular structure with no valves, and the primitive ventricle is connected to the arterial system via the cardiac outflow tract (OFT) [34,35]. The OFT wall consists of three layers: myocardium; endocardium, in direct contact with blood flow; and cardiac jelly, sandwiched between the other two. Thereafter at stage HH20, endocardial cushions, formed from localized protrusions of cardiac jelly, act as primitive valves to facilitate lumen closure and prevent backflow into the ventricle during diastolic filling. These endocardial cushions will eventually develop into valves and septa of the heart [36].

Before the imaging session, a small part of the blunt-end shell and chorionic membrane were removed to access the embryo as illustrated in Fig. 1. The illumination line of the OCT system was focused on the OFT (dashed yellow line in Fig. 1), which is a crucial cardiac segment since a significant number of congenital heart defects originate in the OFT [37].

3 Results

We first demonstrate the capability of differential reconstruction to eliminate autocorrelation artifacts. Figure 4(a) shows a cross-sectional image acquired at 3,500 Hz (without scanning in the transverse direction) after wavenumber resampling and reconstruction. All cross-sectional OCT images have been so rescaled that the pixel size is uniform in both directions. Therefore, only a single scalar is used in these images. The two horizontal bands indicated by the yellow arrows are stationary artifacts arising from autocorrelation interference. They are likely to be attributed to the residual surface reflections from optical components in the interferometer. More autocorrelation artifacts may be noticeable when there are strong reflective layers inside the sample. These artifacts prevent accurate delineation of sample structures and image segmentation. In contrast, the autocorrelation artifacts are almost entirely eliminated in Fig. 4(b), thanks to the differential reconstruction approach.

 

Fig. 4 Two-dimensional in vivo B-scans of chicken embryonic heart (outflow tract) obtained with (a) standard and (b) differential reconstruction approaches. The images were acquired at a scanning rate of 3,500 Hz. The yellow arrows indicate autocorrelation artifacts.

Download Full Size | PDF

Next, we validate the improvement in image quality using NLM denoising by quantifying the SNR of the obtained images. SNR is one of the most important parameters to characterize the image quality, and it generally tends to decrease when imaging at an increasing speed. Figure 5(a) shows a raw image of an HH18 embryonic heart acquired at 5,000 B-scans per second and 10 volumes per second, in which heavy speckle noise is evident. By comparing the signal intensity in the yellow ellipsis and the background noise in the yellow rectangle, a signal to noise ratio of 28.7dB is measured. For stationary samples, time averaging can improve the SNR. However, this is not necessarily the case when imaging a dynamic sample. Averaging the image in Fig. 5(a) with neighboring 50 frames results in the image shown in Fig. 5(b), which exhibits a marginal improvement in image quality. The SNR defined with the same signal and background regions is estimated to be improved by 4.1 dB only. Figure 5(c) shows the result of applying NLM denoising using the same 51 raw images. The SNR enhancement is much more significant in this case and it is measured to be 21.9 dB. More importantly, the structural details are enhanced without motion blur. Due to its outstanding performance, NLM denoising will be applied to the images acquired by our system hereafter.

 

Fig. 5 (a) A B-scan acquired at 5,000 Hz. (b) The image in (a) enhanced by averaging over 50 neighboring frames. (c) The image in (a) enhanced by NLM with 50 neighboring frames.

Download Full Size | PDF

ImageJ was then used to create 3D models using the imaging data related to the above-mentioned chicken embryonic heart, which was scanned at 5,000 Hz and 10 volumes per second. Each volume consists of 500 frames (or B-scans) of 640 by 640 pixels [38]. The image size was adjusted so that the size of each pixel was 2.66 μm in all dimensions. The differentially reconstructed images were processed with NLM over neighboring 50 frames ($M=25$). Shown in Fig. 6 are 3D projections of the heart in a diastolic phase. Figure 6(a) shows the projection of the volumetric imaging data along the transverse scanning direction. The same image is reproduced in Fig. 6(b) with morphological structures delineated. The compact myocardium (CM) layers, endocardium (Endo) layers, cardiac jelly (CJ), and the flowing blood cells are clearly visible with high spatial resolution and contrast [39–41]. Figure 6(c) and 6(d) show projections of the 3D model after rotation by 30 and 135 degrees, respectively. The 3D projections over a range of 180 degrees are included in Visualization 1.

 

Fig. 6 (a) Projection of the 3D model of an HH18 chicken embryonic heart along the transverse scanning direction. (b) Endocardium (Endo), cardiac jelly (CJ), compact myocardium (CM), and blood cells identified in the image. (c) Projection of the same 3D model rotated by 30 degrees, (d) Projection of the 3D model rotated by 135 degrees.

Download Full Size | PDF

Another HH18 chicken embryo heart was scanned at 5,000 Hz and 40 volumes per second. Consequently, the number of B-scans included in each volume was reduced to 125. Each B-scan was enhanced with NLM over 20 neighboring frames ($M=10$). Shown in Fig. 7 is the time series of the beating heart projections along the transverse scanning direction. They consist of 23 time points that cover an entire heart beating cycle. Time-dependent deformation of the heart morphological structures and blood flow changes become very smooth, which can be better seen in the corresponding movie (Visualization 2). Another movie (Visualization 3) has been created using ImageJ 3D projection to show the beating heart rotating within an angle range of 0-360° (at a 10° increment).

 

Fig. 7 Montage of a beating HH18 chicken embryonic heart within a complete cardiac cycle. Each panel is a 3D projection of the heart with the corresponding timing information. The 4D imaging data set was acquired at 40 volumes per second.

Download Full Size | PDF

4. Discussion

We have been using around 15-mW optical power for sample illumination when imaging chicken embryos in vivo. This is equivalent to 30 µW per A-line, which is fairly low compared to other line-field OCT systems. For example, Barrick et al recently used a total power of 492 mW at the sample (or 0.48 mW per A-line) with their parallel SD-OCT system using a supercontinuum light source [42]. While a higher sample power is usually desirable for maximizing the sensitivity and maintaining a good signal to noise ratio, we decided not to go beyond 20 mW as the chicken embryo heart might beat irregularly when illuminated with such a high-intensity beam.

The sensitivity of our line-scan system was measured 87 dB at an acquisition speed of 3,500 B-scans per second. Such a high sensitivity, however, seemed valid only for stationary samples such as a coverslip. There was a noticeable sensitivity loss when imaging live samples in vivo. Structures moving towards or away from the objective might result in fringe washout, an adverse effect that line-scan OCT systems suffer more than their point-scanning counterparts. Imaging of the rapid deformation of cardiac tissues and fast flow of blood cells in the OFT is definitely more challenging than dealing with other slowly moving tissues or organs. Besides fringe washout, other issues might also contribute to the signal loss. We used a relatively high NA (~0.1) objective to achieve a lateral resolution of 6.2 μm. The corresponding axial scanning range was therefore limited to a couple of hundred microns. Consequently, some deep structures might be out of the depth of focus and light signal reflected from them was not collected efficiently. Dispersion compensation was not included in our system as we expected marginal sample dispersion when the imaging depth is less than 1 mm. However, this might not be true considering the very large spectral width. In our future study, we will quantify the impact of dispersion and find a solution.

We have been using a B-scan rate of 5 kHz for multi-dimensional imaging of chicken embryo heart. The imaging speed is equivalent to 3.2M A-lines per second, or 2.048G voxels per second. However, the system can reach a maximal A-line rate of 3.584 MHz when the camera is operated at 3,500 frames per second, as there will be more A-lines (1024 vs. 640) in a B-scan. It is more than three times faster than previously reported parallel SD-OCT systems. Image post-processing, however, may take a much longer time. NLM is the most time consuming process. It took roughly 12307 secconds (~3.4 hours) to process a volume of 125 B-scans, each of which is 640 by 600 pixels. For other processes including differential reconstruction, the total processing time is less than a couple of seconds.

While the supercontinuum light source provides a very broad output spectrum, the center wavelength chosen for our OCT setup is 750 nm. Compared with OCT systems at longer center wavelengths, e.g., 1300 nm, a short center wavelength allows better transverse resolution. However, the tissue scattering coefficient becomes higher and the achievable imaging depth in the embryonic heart is therefore reduced.

The most expensive parts in our system are the light source (USD30,000) and the high-speed camera (USD40,000). The total cost (around USD 80,000) is slightly more expensive than standard commercially available OCT systems. For example, the Telesto OCT Base Unit (TEL320, Thorlabs) is priced at USD64,000. It provides an up to 146kHz A-line rate and an axial resolution of 5.5 µm (air).

5. Conclusion

In summary, we have demonstrated an ultra-high speed and high-resolution line-scan SD-OCT system by adopting a supercontinuum white light source (center wavelength: 750 nm, spectral bandwidth: 180 nm), superior line focusing mechanism using Powell lens, ultra-fast image sensor, autocorrelation elimination using phase stepping and differential reconstruction, and NLM denoising algorithm. The system can be used to obtain high-quality images at up to 3.584 million A-lines per second. Free-space axial resolution below 2.0 µm and lateral resolution around 6.2 µm have been demonstrated. To the best of our knowledge, this is the fastest line-scan OCT system reported thus far. This ultra-high speed will allow the tracking of individual blood cells and associated functional changes therein in biological organs.