Increasing the field-of-view of dynamic cardiac OCT via post-acquisition mosaicing without affecting frame-rate or spatial resolution.

Optical coherence tomography (OCT) allows imaging dynamic structures and fluid flow within scattering tissue, such as the beating heart and blood flow in murine embryos. For any given system, the frame rate, spatial resolution, field-of-view (FOV), and signal-to-noise ratio (SNR) are interconnected: favoring one aspect limits at least one of the others due to optical, instrumentation, and software constraints. Here we describe a spatio-temporal mosaicing technique to reconstruct high-speed, high spatial-resolution, and large-field-of-view OCT sequences. The technique is applicable to imaging any cyclically moving structure and operates on multiple, spatially overlapping tiled image sequences (each sequence acquired sequentially at a given spatial location) and effectively decouples the (rigid) spatial alignment and (non-rigid) temporal registration problems. Using this approach we reconstructed full-frame OCT sequences of the beating embryonic rat heart (11.5 days post coitus) and compared it to direct imaging on the same system, demonstrating a six-fold improvement of the frame rate without compromising spatial resolution, FOV, or SNR.

1. Introduction

Optical coherence tomography (OCT) is increasingly being adopted as an accurate imaging method to study embryonic cardiovascular anatomy and dynamics. It allows for non-destructive imaging through several millimeters of biological tissue with single cell resolution. OCT has been used for a wide array of biomedical applications, including structural imaging of whole mammalian and avian embryos [1,2], flow measurements [3,4], and 4D imaging of mammalian and avian embryonic hearts [5-7]. High frame rate, signal-to-noise ratio (SNR), and spatial resolution as well as the ability to image over a large field of view are all equally important factors for accurate visualization and analysis.

Improvements in light sources, scanner efficiency, data transmission protocols, and storage capacity have collectively produced improvements in the acquisition frame rate for OCT [8, 9]. The time over which the signal is gathered to produce one image limits the spatial resolution when imaging dynamic samples since sample motion introduces blurring if the image capture time is too long [10]. Increasing the frame rate can be achieved through faster scanning but decreasing the dwell time also reduces the SNR, as fewer photons contribute to the signal. For cardiovascular systems, which exhibit velocities of the order of millimeters per second, framerates of 100 frames per second or more are required to achieve a spatial resolution of 10 μm and below. Thus, there is a need for advanced methods that improve the frame rate while preserving the spatial resolution.

For optical and electron microscopy of static or nearly static samples, tiled imaging, followed by mosaic reconstruction, has become a central tool to extend the limited field of view of high-numerical aperture microscope objectives. Several methods have been proposed to automatically or semi-automatically spatially register (align) individual microscopy tiles to form an image mosaic [11-14] extending methods initially developed for panoramic photography [15, 16]. These technique, however, produce static images and are therefore not directly applicable in our case.

Recently, we demonstrated that dynamic 3D images of the embryonic zebrafish heart can be reconstructed by post-acquisition synchronization of nongated slice sequences acquired on a fast confocal microscope [17, 18]. Building upon our ability to accurately register image sequences in time, we here present a combined (rigid) spatial and (non-rigid) temporal registration method that allows us to extend the scope of static mosaicing methods to improve the field of view and resolution of dynamic 2D OCT by combining tiled image sequences of cyclically moving structures in the mammalian embryonic heart.

2. Method

Our method operates as follows. We assume the object to be imaged has a local scattering intensity that is periodic, that is where T is the period of the heart beat (assuming heart imaging). We acquire N image sequences (see Fig. 1(a) for an example with N = 2) that are spatially tiled such that for 0 ≤ x < L, 0 ≤ y < L, and 0 ≤ t < L, where L and L are the width and height of the tiles, respectively, L is the duration of each sequence and where x is the spatial offset in the x-direction of the nth tile, and t is the time at which recording of the nth sequence started. The duration of the sequence, L, must be at least two periods long, so that one full heart beat period, starting at any arbitrary time in the cardiac cycle, is included in the sequence.

Fig. 1

Tiled acquisition and automatic mosaicing procedure for dynamic images. (a) Overlapping cardiac image series, are acquired sequentially. (b) Each sequence is collapsed to a static image (temporal maximum). (c) Static images are spatially aligned. (d) Image sequences are spatially transformed and padded with zeros. (e) Spatially aligned sequences are temporally synchronized. (f) Synchronized image sequences are blended. (Color online).

We assume that the spatial and temporal offsets, x and t, respectively, are unknown. Our task is to estimate them, which is an image registration (and synchronization) problem (we assume, here, that tiles are aligned in the y direction, but this is not a fundamental requirement of our approach, which could be generalized to include registration in both the x and y directions). The periodicity assumption allows us to decouple the spatial and temporal registration problems as follows. We first collapse each sequence along the temporal direction (Fig. 1(b)) to obtain a spatial signature of the tile that is independent of time We next pair-wise align the N static tiles Ī, that is, we estimate the tiling-offsets x as follows with x0 = 0 (Fig. 1(c)). In words, we find the horizontal position for which the mean-squared difference between the previous and current tiles’ grayscale values (calculated over the region of overlap of the two tiles) is minimal. In practice, we use a rigid-body, pyramid-based, registration implementation proposed by Thévenaz et al. [19], which we constrain to register tiles along the horizontal direction.

We next define zero-padded image sequences in which the original sequences have been placed at their estimated spatial offset (see Fig. 1(d)). These sequences are spatially, but not temporally aligned. Proceeding recursively for n = 1,...,N – 1, and choosing the sequence I′0 as a reference, we align pairs of sequences I′ and I′–1 via a non-uniform time-warping procedure in which we minimize the following cost functions (1 ≤ n < N): where w−1 (with w0(t) = t, the identity function) is the warping function applied to the sequence I′−1, and w is a candidate cost function to warp In so that it matches one period of the warped version of I′−1. The first integral compares the two neighboring, temporally warped sequences (sequences whose time-axis has been deformed), while the second integral keeps the extent of this deformation in check by penalizing candidate warping functions w that stretch or compress the time axis. These two contributions to the cost function are balanced by the parameter 0 < λ ≤ 1 to favor either good matching of the warped and reference sequences (λ = 1) or the temporal integrity of the warped sequence (λ → 0). We constrain the search for warping functions w to the set ℳ = {w ∈ C1([0, T))|0 ≤ w(t) < L and w(t1) < w(t2), t1 < t2} of continuous, non-negative, and strictly increasing functions bounded by L and defined over the interval [0, T). To determine the warping function w that minimizes the cost function, we use a previously described dynamic programming algorithm [20]. The optimal warping function w produced by this algorithm is such that w ∈ {w̄ ∈ ℳ |Q{w̄} = min∈ Q{w}}. Applying this warping function to the sequence I′(x, y,t) results in sequence I′(x, y, w(t)), which is synchronized to all previously warped tiles, I′(x, y, wi(t)), i = 0,..., n – 1 (see Fig. 1(e)).

Once all warping functions are determined, we blend the time-warped sequences according to where, B(x) is a function defined for a < b such that and The final output of our algorithm (see Fig. 1(f)) is a blended, synchronized image series that covers one full period

3. Experiments and Results

To verify our method in practice, we acquired image sequences of a live rat embryo (11.5 days post coitus) using a swept-source OCT (SS-OCT) system as described previously in [3]. Figure 2 shows a representative frame of an SS-OCT image series of the embryonic heart and yolk sac acquired at different frame rates and spatial resolutions. Increasing the frame rate limits the spatial resolution over the region of interest (Fig. 2(a); 512 x 64 pixels, FOV=2.8 mm × 1 mm, 150 frames per second, SNR=15.76dB ). Improving the spatial resolution decreases the frame rate over the same region (Fig. 2(b); 512 × 500 pixels, FOV=2.8 mm × 0.97 mm, 25 frames per second, SNR=15.93dB ). Figure 2(c) shows the series of consecutive image sequences with smaller FOV (each 512 × 64 pixels, FOV=2.8 mm × 0.125 mm, 150 frames per second, 1000 frames (6.7 s, about 13 periods), SNR=15.76dB ) stitched together. The overlaps of the series are marked at the bottom. Figure 2(d) and Media 1 show the result of the synchronization of the same series of sequences. Importantly, the latter reconstructed sequence has the same spatial resolution as the series acquired at low frame rate, Fig. 2(b) (512 × 500 pixels, FOV=2.8 mm × 0.97 mm), but with the same temporal resolution as the low spatial resolution sequence in Fig. 2(a) (150 frames per second) and a better signal to noise ratio (SNR=16.30dB ), thereby demonstrating that we can jointly achieve high spatial and temporal resolutions while preserving the field of view. The increase in SNR is due to the fact that in the areas of overlap the resulting signal is a weighted average of two series.

Fig. 2

Dynamic mosaicing overcomes frame rate, resolution, and FOV interdependence for fast SS-OCT cardiac imaging. (a) High frame rate, low spatial resolution, large FOV: 150 fps, 512×64 pixels (b) Low frame rate, high spatial resolution, large FOV, 25fps, 512×500 pixels. (c) High frame rate, high spatial resolution, small FOV (12 tiles, before dynamic mosaicing, each 150 fps, 512×64 pixels). Marks at bottom of panel indicate lateral position of individual tiles. (d) High frame rate, high spatial resolution, high FOV (after dynamic mosaicing), 150 fps, 512×500 pixels. Inserts represent time-course (500 ms total duration) of pixel intensity (arbitrary units) at location marked by arrows. SNRs were computed over rectangular boxes. Scale bars are 0.2 mm. ( Media 1).

To evaluate the accuracy of our method, we generated test datasets by extracting subsequences from each of N = 11 measured tile image sequence (512 × 64 pixels, FOV=2.8 mm × 0.125 mm, 150 frames per second, 1000 frames; see above and Fig. 2(c)). Specifically, we first split each measured sequence into two overlapping, left and right tile sequences, both of size 512 × 42 pixels and same duration as the original sequence (Figs. 3(b) and 3(c)). We further extracted frames covering two nonoverlapping time intervals (one interval at the beginning and one toward the end of the sequence, both spanning at least two heartbeat periods) from both the original sequence and the left and right sequences. Each measured image sequence thereby resulted in six subsequences labeled A through F (Figs. 3(a)–3(c)).

Fig. 3

Accuracy evaluation of tiled acquisition and automatic mosaicing procedure. (a) Original sequence (one of the tiles in Fig. 2(c)) is split into (b) a left tile sequence and (c) a right tile sequence. (d)–(g) Sequence subsets are registered in space and/or time. (d) Temporal alignment of first period in subsequence A to subsequence B. (e) Spatio-temporal alignment of first period in subsequence C to subsequence F. (f) Spatio-temporal alignmnent of first period in subsequence E to subsequence D. (g) Spatio-temporal alignment of first period in subsequence C to subsequence E.

Given the six subsequences generated from any one of the 11 measured sequences, we registered the first cardiac cycle in subsequence A to subsequence B (an operation we denote by A→B, see also Fig. 3(d)), and similarly, the first cardiac cycle in subsequences C, E, and C to subsequences F, D, and E, respectively (i.e. C→F, E→D, and C→E, see Figs. 3(e), 3(f), and 3(g), respectively).

We first verified that the registration algorithm produced similar reconstructions when merging left or right tiles from an early cardiac cycle (subsequence C or F, respectively) to right or left tiles, respectively, extracted from later cardiac cycles (subsequence F or D, respectively) or when synchronizing non-tiled sequences directly (subsequences A and B). The temporal registration results for C→F and E→D were therefore compared against each other and against those obtained when directly synchronizing non-tiled frames of the same early cardiac cycle (subsequence A) to the same later cardiac cycles (subsequence B), i.e. A→B.

Assuming the notation W, X, Y, or Z for arbitrary subsequence labels A through F, we compared temporal warping function , obtained when synchronizing subsequences W and X, with , obtained when synchronizing subsequences Y and Z (with all four subsequences, W, X, Y, and Z, extracted from the same tile n), by computing a discrete equivalent of which captures the average temporal discrepancy between the two synchronization operations.

Also, we verified that the spatial match of left and right tiles was found correctly (assuming that the sample’s spatial motion was identical from one heartbeat to the next). We calculated the spatial registration error incurred when aligning subsequence W to subsequence X (both extracted from tile n), given the overlap x (known by construction) and the estimated offset , as

The average spatial error and temporal discrepancies were computed in each case for N = 11 tiles, i.e., and Numerical results are reported in the first three columns of Table 1 showing discrepancies of less than a frame in time and spatial registration accuracy of the order of one pixel.

Table 1

Registration Evaluation*

Registration Experiment
(W→X)	C→F	E→D	C→F	C→E
Spatial Registration Evaluation
Average spatial discrepancy
Δx̃W→X (pixels)	0.4 ± 0.1	1.7 ± 0.5	0.4 ± 0.1	0.0 ± 0.00
Temporal Registration Evaluation
Reference transform (Y→Z)	A→B	A→B	E→D	Identity
Average temporal discrepancy
ΔtW→X,Y→Z or ΔtC→E (frames)	0.2 ± 0.1	0.6 ± 0.4	0.6 ± 0.4	0.0 ± 0.00

Each experiment was carried out on N = 11 tiles and the average result is provided. The sequence labels A–F correspond to those in Fig. 3. The notation W→X, stands for “First cardiac cycle of subsequence W is registered to subsequence X”, where symbols W, X, Y, and Z stand for labels in the range A–F.

Finally, as a control to ensure that the spatial registration procedure was accurate and did not affect the accuracy of the temporal registration procedure, we also merged tile-pairs generated from the same cardiac period. To this end, we registered subsequence C to subsequence E (within the same tile, see Fig. 3(g)). The spatial error was obtained as above, while the temporal error was obtained as Both spatial and temporal registration errors, reported in the last column of Table 1, show that registration was perfect in this case for all 11 tiles tested.

4. Conclusion

In summary, we have demonstrated the applicability of an image acquisition protocol and reconstruction algorithm for OCT images of cardio-vascular structures whose motion is cyclical, as dictated by the cyclic contraction of the heart. Dynamic mosaicing overcomes the limitations that arise from the frame rate, resolution, and FOV being interdependent for fast OCT cardiac imaging. Although we believe our space-time registration method is well adapted to the task at hand as it provides a simple solution to separately address the spatial and temporal registration problems (which considerably reduces computational complexity) other methods for space-time registration developed in the context of human cardiac imaging in adults using magnetic resonance imaging, single-photon emission computed tomography, or X-ray computer tomography could be used for similar purposes. In particular, use of registration techniques that are not only nonrigid in time but also in space [21, 22] could likely improve our mosaic reconstructions in cases where the spatial motion of the embryonic heart is not conserved from one heartbeat to the next.

Future work will include investigating possible improvements that could result from using such methods as well as combining our procedure with previously developed noise removal strategies [23] or techniques to reconstruct high frame rate 3D+time volumes from 2D+time series [17].