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

JeaBuem Yoo; Irina V. Larina; Kirill V. Larin; Mary E. Dickinson; Michael Liebling

doi:10.1364/BOE.2.002614

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

I (x, y, t) \approx I (x, y, t + T) \forall t \in ℝ,

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

I_{n} (x, y, t) = I (x_{n} + x, y, t_{n} + t), 0 \leq n < N,

for 0 ≤ x < L_x, 0 ≤ y < L_y, and 0 ≤ t < L_t, where L_x and L_y are the width and height of the tiles, respectively, L_t is the duration of each sequence and where x_n is the spatial offset in the x-direction of the nth tile, and t_n is the time at which recording of the nth sequence started. The duration of the sequence, L_t, 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).

Download Full Size | PDF

We assume that the spatial and temporal offsets, x_n and t_n, 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

{\bar{I}}_{n} (x, y) = max_{0 \leq t < L_{t}} I_{n} (x, y, t) .

We next pair-wise align the N static tiles Ī_n, that is, we estimate the tiling-offsets x_n as follows

{\tilde{x}}_{n} = {\tilde{x}}_{n - 1} + arg min_{0 \leq x < L_{x}} \int_{0}^{L_{y}} \frac{1}{L_{x} - x} \int_{x}^{L_{x}} {| {\bar{I}}_{n} (x' - x, y) - {\bar{I}}_{n - 1} (x', y) |}^{2} d x' d y,

with x ₀ = 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

I_{n}^{'} (x, y, t) = {\begin{array}{l} 0 & 0 \leq x < {\tilde{x}}_{n} \\ I_{n} (x - {\tilde{x}}_{n}, y, t) & {\tilde{x}}_{n} \leq x < {\tilde{x}}_{n} + L_{x} \\ 0 & {\tilde{x}}_{n} + L_{x} \leq x < {\tilde{x}}_{N} + L_{x}, \end{array}

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′ ₀ as a reference, we align pairs of sequences I′_n and I′_n _–1 via a non-uniform time-warping procedure in which we minimize the following cost functions (1 ≤ n < N):

\begin{array}{l} Q_{n} {w} = λ \int_{0}^{T} \int_{0}^{L_{y}} \int_{{\tilde{x}}_{n}}^{{\tilde{x}}_{n - 1} + L_{x}} | I_{n - 1}^{'} (x, y, w_{n - 1} (t)) - I_{n}^{'} (x, y, w (t)) | d x d y d t \\ + (1 - λ) \int_{0}^{T} | 1 - \frac{d}{d t} w (t) | d t, \end{array}

where w_n ₋₁ (with w ₀(t) = t, the identity function) is the warping function applied to the sequence I′_n ₋₁, and w is a candidate cost function to warp In so that it matches one period of the warped version of I′_n ₋₁. 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_n to the set ℳ = {w ∈ C ¹([0, T))|0 ≤ w(t) < L_t and w(t ₁) < w(t ₂), t ₁ < t ₂} of continuous, non-negative, and strictly increasing functions bounded by L_t and defined over the interval [0, T). To determine the warping function w_n that minimizes the cost function, we use a previously described dynamic programming algorithm [20]. The optimal warping function w_n produced by this algorithm is such that w_n ∈ {w̄_n ∈ ℳ |Q_n{w̄_n} = min_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

I_{out}^{(n)} (x, y, t) = B_{{\tilde{x}}_{n}, {\tilde{x}}_{n - 1} + L_{x}} (x) I_{out}^{(n - 1)} (x, y, t) + (1 - B_{{\tilde{x}}_{n}, {\tilde{x}}_{n - 1} + L_{x}} (x)) I_{n}^{'} (x, y, w_{n} (t))

where, B_a,b(x) is a function defined for a < b such that

B_{a, b} (x) = {\begin{array}{l} 1, & x < a \\ 1 - \frac{x - a}{b - a}, & a \leq x \leq b \\ 0, & x > b . \end{array}

and

I_{out}^{(0)} = I_{0}^{'} (x, y, t) .

The final output of our algorithm (see Fig. 1(f)) is a blended, synchronized image series that covers one full period

I_{out} (x, y, t) = I_{out}^{(N - 1)} (x, y, t), 0 \leq t < T .

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

Download Full Size | PDF

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.

Download Full Size | PDF

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

Δ t_{n}^{W \to X, Y \to Z} = \frac{1}{T} \int_{0}^{T} | w_{n}^{W \to X} (t) - w_{n}^{Y \to Z} (t) | d t,

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 $Δ {\tilde{x}}_{n}^{W \to 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}^{W \to X}$ , as

Δ {\tilde{x}}_{n}^{W \to X} = | {\tilde{x}}_{n}^{W \to X} - x_{n} | .

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

Δ {\tilde{x}}^{W \to X} = \frac{1}{N} \sum_{n = 1}^{N} Δ {\tilde{x}}_{n}^{W \to X}

and

Δ t^{W \to X, Y \to Z} = \frac{1}{N} \sum_{n = 1}^{N} Δ t_{n}^{W \to X, Y \to Z} .

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^*

View Table

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

Δ t_{n}^{C \to E} = \frac{1}{T} \int_{0}^{T} | w_{n}^{C \to E} (t) - t | d t .

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

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]

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
Δt ^W→X,Y→Z or Δt ^C→E (frames)	0.2 ± 0.1	0.6 ± 0.4	0.6 ± 0.4	0.0 ± 0.00

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

Abstract

1. Introduction

2. Method

3. Experiments and Results

4. Conclusion

Acknowledgments

References and links

Supplementary Material (1)

Cited By

Figures (3)

Tables (1)

Equations (15)

Biomedical Optics Express