## Abstract

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.

©2011 Optical Society of America

## 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

*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

*x*<

*L*, 0 ≤

_{x}*y*<

*L*, and 0 ≤

_{y}*t*<

*L*, where

_{t}*L*and

_{x}*L*are the width and height of the tiles, respectively,

_{y}*L*is the duration of each sequence and where

_{t}*x*is the spatial offset in the

_{n}*x*-direction of the

*n*th tile, and

*t*is the time at which recording of the

_{n}*n*th 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.

_{t}We assume that the spatial and temporal offsets, *x _{n}* 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

_{n}*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

*N*static tiles

*Ī*, that is, we estimate the tiling-offsets

_{n}*x*as follows

_{n}*x*

_{0}= 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

*n*= 1,...,

*N –*1, and choosing the sequence

*I′*

_{0}as a reference, we align pairs of sequences

*I′*and

_{n}*I′*

_{n}_{–1}via a non-uniform time-warping procedure in which we minimize the following cost functions (1 ≤

*n*<

*N*):

*w*

_{n}_{−1}(with

*w*

_{0}(

*t*) =

*t*, the identity function) is the warping function applied to the sequence

*I′*

_{n}_{−1}, and

*w*is a candidate cost function to warp In so that it matches one period of the warped version of

*I′*

_{n}_{−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

_{n}*ℳ*= {

*w*∈

*C*

^{1}([0,

*T*))|0 ≤

*w*(

*t*) <

*L*and

_{t}*w*(

*t*

_{1}) <

*w*(

*t*

_{2}),

*t*

_{1}<

*t*

_{2}} of continuous, non-negative, and strictly increasing functions bounded by

*L*and defined over the interval [0,

_{t}*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

_{n}*w*produced by this algorithm is such that

_{n}*w*∈ {

_{n}*w̄*∈

_{n}*ℳ*|

*Q*{

_{n}*w̄*} = min

_{n}

_{w}_{∈}

_{ℳ}*Q*{

_{n}*w*}}. Applying this warping function to the sequence

*I*′

*(*

_{n}*x, y,t*) results in sequence

*I′*(

_{n}*x, y, w*(

_{n}*t*)), which is synchronized to all previously warped tiles,

*I′*(

_{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

*B*(

_{a,b}*x*) is a function defined for

*a*<

*b*such that

## 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.

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)).

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
${w}_{n}^{\text{W}\to \text{X}}(t)$, obtained when synchronizing subsequences W and X, with
${w}_{n}^{\text{Y}\to \text{Z}}(t)$, 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

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
$\mathrm{\Delta}{\tilde{x}}_{n}^{\text{W}\to \text{X}}$ incurred when aligning subsequence W to subsequence X (both extracted from tile *n*), given the overlap *x _{n}* (known by construction) and the estimated offset
${\tilde{x}}_{n}^{\text{W}\to \text{X}}$, as

The average spatial error and temporal discrepancies were computed in each case for *N* = 11 tiles, i.e.,

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

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

## Acknowledgments

This research was supported by DOD ONR YIP (K.V.L.), NIH HL077187, HL095586, and T32HL007676, grants from UCSB/Santa Barbara Cottage Hospital (M.L), AHA10SDG3830006 (I. V. L.), and a Hellman Family Faculty fellowship (M.L.). We thank Sandeep Bhat for help in preparing Fig. 2 and Movie 1.

## References and links

**1. **T. Yelbuz, M. Choma, L. Thrane, M. Kirby, and J. Izatt, “A new high-resolution imaging technology to study cardiac development in chick embryos,” Circulation **106**, 2771–2774 (2002). [CrossRef] [PubMed]

**2. **I. V. Larina, E. F. Carbajal, V. V. Tuchin, M. E. Dickinson, and K. V. Larin, “Enhanced OCT imaging of embryonic tissue with optical clearing,” Laser Phys. Lett. **5**, 476–479 (2008). [CrossRef]

**3. **I. Larina, N. Sudheendran, M. Ghosn, J. Jiang, A. Cable, K. Larin, and M. Dickinson, “Live imaging of blood flow in mammalian embryos using Doppler swept-source optical coherence tomography,” J. Biomed. Opt. **13**, 060506 (2008). [CrossRef]

**4. **I. Larina, S. Ivers, S. Syed, M. Dickinson, and K. Larin, “Hemodynamic measurements from individual blood cells in early mammalian embryos with Doppler swept source OCT,” Opt. Lett. **34**, 986–988 (2009). [CrossRef] [PubMed]

**5. **K. V. Larin, I. V. Larina, M. Liebling, and M. E. Dickinson, “Live imaging of early developmental processes in mammalian embryos with optical coherence tomography,” J. Innov. Opt. Health Sci. **2**, 253–259 (2009). [CrossRef] [PubMed]

**6. **M. Gargesha, M. W. Jenkins, D. L. Wilson, and A. M. Rollins, “High temporal resolution OCT using image-based retrospective gating,” Opt. Express **17**, 10786–10799 (2009). [CrossRef] [PubMed]

**7. **A. Liu, R. Wang, K. Thornburg, and S. Rugonyi, “Efficient postacquisition synchronization of 4-D nongated cardiac images obtained from optical coherence tomography: application to 4-D reconstruction of the chick embryonic heart,” J. Biomed. Opt. **14**, 044020 (2009). [CrossRef] [PubMed]

**8. **T. Schmoll, C. Kolbitsch, and R. A. Leitgeb, “Ultra-high-speed volumetric tomography of human retinal blood flow,” Opt. Express **17**, 4166–4176 (2009). [CrossRef] [PubMed]

**9. **W. Wieser, B. R. Biedermann, T. Klein, C. M. Eigenwillig, and R. Huber, “Multi-megahertz OCT: High quality 3D imaging at 20 million A-scans and 4.5 GVoxels per second,” Opt. Express **18**, 14685–14704 (2010). [CrossRef] [PubMed]

**10. **J. Vermot, S. E. Fraser, and M. Liebling, “Fast fluorescence microscopy for imaging the dynamics of embryonic development,” HFSP J . **2**, 143–155 (2008). [CrossRef]

**11. **A. Can, C. Stewart, B. Roysam, and H. Tanenbaum, “A feature-based technique for joint, linear estimation of high-order image-to-mosaic transformations: Mosaicing the curved human retina,” IEEE Trans. Pattern Anal. **24**, 412–419 (2002). [CrossRef]

**12. **P. Thévenaz and M. Unser, “User-friendly semiautomated assembly of accurate image mosaics in microscopy,” Microsc. Res. Tech. **70**, 135–146 (2007). [CrossRef]

**13. **M. Emmenlauer, O. Ronneberger, A. Ponti, P. Schwarb, A. Griffa, A. Filippi, R. Nitschke, W. Driever, and H. Burkhardt, “XuvTools: free, fast and reliable stitching of large 3D datasets,” J. Microsc. **233**, 42–60 (2009). [CrossRef] [PubMed]

**14. **J. Zupanc, A. Dobnikar, D. Drobne, J. Valant, D. Erdogmus, and E. Bas, “Biological reactivity of nanoparticles: mosaics from optical microscopy videos of giant lipid vesicles,” J. Biomed. Opt. **16**, 026003 (2011). [CrossRef] [PubMed]

**15. **H. Sawhney and R. Kumar, “True multi-image alignment and its application to mosaicing and lens distortion correction,” IEEE Trans. Pattern Anal. **21**, 235–243 (1999). [CrossRef]

**16. **H. Shum and R. Szeliski, “Systems and experiment paper: Construction of panoramic image mosaics with global and local alignment,” Int. J. Comput. Vis. **36**, 101–130 (2000). [CrossRef]

**17. **M. Liebling, A. S. Forouhar, M. Gharib, S. E. Fraser, and M. E. Dickinson, “Four-dimensional cardiac imaging in living embryos via postacquisition synchronization of nongated slice sequences,” J. Biomed. Opt. **10**, 054001 (2005). [CrossRef] [PubMed]

**18. **M. Liebling, A. S. Forouhar, R. Wolleschensky, B. Zimmerman, R. Ankerhold, S. E. Fraser, M. Gharib, and M. E. Dickinson, “Rapid three-dimensional imaging and analysis of the beating embryonic heart reveals functional changes during development,” Dev. Dynam. **235**, 2940–2948 (2006). [CrossRef]

**19. **P. Thévenaz, U. E. Ruttimann, and M. Unser, “A pyramid approach to subpixel registration based on intensity,” IEEE Trans. Image Process. **7**, 27–41 (1998). [CrossRef]

**20. **M. Liebling, J. Vermot, A. Forouhar, M. Gharib, M. Dickinson, and S. Fraser, “Nonuniform temporal alignment of slice sequences for four-dimensional imaging of cyclically deforming embryonic structures,” in “Proc. ISBI 2006,”(2006), pp. 1156–1159.

**21. **A. Frangi, D. Rueckert, J. Schnabel, and W. Niessen, “Automatic construction of multiple-object three-dimensional statistical shape models: Application to cardiac modeling,” IEEE Trans. Med. Imaging **21**, 1151–1166 (2002). [CrossRef]

**22. **D. Perperidis, R. H. Mohiaddin, and D. Rueckert, “Spatio-temporal free-form registration of cardiac MR image sequences,” Med. Image Anal. **9**, 441–456 (2005). [CrossRef] [PubMed]

**23. **S. Bhat, I. V. Larina, K. V. Larin, M. E. Dickinson, and M. Liebling, “Multiple-cardiac-cycle noise reduction in dynamic optical coherence tomography of the embryonic heart and vasculature,” Opt. Lett. **34**, 3704–3706 (2009). [CrossRef] [PubMed]