High temporal resolution OCT using image-based retrospective gating

High temporal resolution OCT imaging is very advantageous for analyzing cardiac mechanics in the developing embryonic heart of small animals. An image-based retrospective gating technique is presented to increase the effective temporal resolution of an OCT system and to allow visualization of systolic dynamics in 3D. The gating technique employs image similarity measures for rearranging asynchronously acquired input data consisting of a time series of 2D images at each z position along the heart volume, to produce a time sequence of 3D volumes of the beating heart. The study includes a novel robust validation technique, which quantitatively evaluates the accuracy of the gating technique, in addition to visual evaluations by 2D multiplanar reformatting (MPR) and 3D volume rendering. The retrospective gating and validation is demonstrated on a stage 14 embryonic quail heart data set. Using the validation scheme, it is shown that the gating is accurate within a standard deviation of 4.7 ms, which is an order of magnitude shorter than the time interval during which systolic contraction ( ∼ 50 ms) occurs in the developing embryo. This gating method has allowed, for the first time, clear visualization of systolic dynamics of the looping embryonic heart in 3D.


Introduction
Studying abnormalities during heart development in the embryo is crucial for the understanding of congenital heart disease. Our group has been developing specialized OCT-based instruments and software and applying them to the study of the developing heart in avian and mouse models [1][2][3][4][5]. The developing morphology is of interest, but even more interesting is the development of cardiac function, the ability of the heart to contract and pump blood. Several groups have shown that mechanical forces influence the morphogenesis of the developing heart [6][7][8][9]. In this paper, we report results on avian heart, a model which has contributed greatly to the understanding of human congenital heart disease. OCT is uniquely well suited for imaging fast cardiac dynamics in living hearts because of its microscopic resolution and sufficient depth of penetration [1][2][3][4][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. However, a real-time, OCT volume scan (∼117,000 A-scans/s) of an embryonic quail heart using an ultrahigh-speed OCT system [1,3] produces only ∼4 volumes per heartbeat when the heart beats at physiological rates, which is clearly insufficient to analyze rapid cardiac dynamics in 3D. A possible remedy to this problem is retrospective cardiac gating, where one obtains data over multiple heartbeats and resynchronizes the data after the fact to reconstruct 3D volumes.
There are gating techniques in the literature for acquiring 4D (3D volume over time) cardiac data. First, in the case of prospective gating, one uses a physiological signal (e.g., ECG) to trigger data acquisition, typically during multiple phases of the cardiac cycle. In subsequent heart beats, more data at each phase is acquired, so that after several heart beats, one can create a beating volumetric heart. This approach has been used in various medical imaging modalities (MRI, CT, nuclear medicine, etc.) in addition to OCT [2,4,25,26]. With OCT imaging, we have previously demonstrated prospective gating of an excised heart using a pacing signal and in vivo prospective gating using a Laser Doppler Velocimetry (LDV) signal [2,4]. Prospective gating suffers if the physiological signal is noisy giving inconsistent triggers, as can be the case when imaging small embryonic structures. This usually results in a reduced temporal resolution. Depending on the imaging modality, additional processing needs to be performed to correct for this inconsistency. In addition, in our previous prospective gating study using LDV [4], significant set up time was required due to positioning requirements of the LDV probe, greatly limiting high-throughput imaging (i.e., multiple embryos in a single session). With regards to temporal resolution, with LDV gating, 8 volumes were acquired in a cardiac cycle (∼100 ms/volume) from a stage 14 quail embryo whose heart rate was slightly depressed (75 beats/minute) due to the imaging environment. This would be clearly inadequate to observe systolic contraction (∼50 ms) occurring within the developing embryo. Second, one can record a physiological signal while acquiring image data as fast as possible over several heartbeats at different spatial locations. One then uses timing from the signal to collect the data into appropriate phases of the cardiac cycle creating a 4D beating heart. This approach is called signal-based retrospective gating [26][27][28][29][30]. In a recent study, a Doppler optical cardiogram signal has been used for this purpose [16]. However, in early embryonic developmental studies with small hearts, obtaining these signals adds to system complexity and cost. As with the case of prospective gating, setup time can be prohibitive when performing studies on multiple embryos. Therefore, a third category of gating techniques called image-based retrospective gating, which rearranges out-of-order data based on image similarity, offers an alternative approach [31][32][33][34]. In particular, Liebling et al. [33] applied image-based retrospective gating on a series of confocal microscopy images of living zebrafish embryos to obtain 4D reconstructed data.
In this paper, we have developed an image-based retrospective gating algorithm which is more robust than our previous prospective gating methods [2,4] and greatly simplifies experimental set up. We have used it to observe 3D/4D cardiac systolic dynamics of embryonic quail hearts at ∼90 volumes per heartbeat (∼270 volumes/s), giving detailed images of cardiac contraction. Our image processing approach has been influenced by the report of Liebling et al. [33]. However, due to the higher SNR of our imaging setup, we use image feature matching in the spatial domain rather than the wavelet transform domain, as done by them. By considering large numbers of pixels rather than relatively few wavelet coefficients, accuracy should be improved and straightforward quality metrics can be computed. In addition, we created a novel quantitative validation technique whereby we obtained a real-time reference scan of the same heart from which out-of-order data was collected for retrospective gating. Subsequently, crosscorrelation image matching was performed between frames from this real-time, lower temporal resolution reference scan and those from the gated data set so as to evaluate the accuracy of retrospective gating. In addition, this provides a means of establishing limits on algorithmic parameters such as spatial and temporal decimation of data, a novel study not present in current literature [31][32][33][34]. We believe that the high temporal resolution afforded by image-based retrospective gating will help us develop more accurate models of cardiac dynamics than previously possible. In an other study by our group [3] that used real-time volume imaging, we had measured the heart tube diameter and computed the maximum wall velocity for the beating heart. The maximum wall velocity in turn had helped us determine the theoretical imaging speeds that would be required to capture fast dynamics of the beating heart with little displacement error. Our current study represents a significant step towards lowering this displacement error through an effective increase in imaging speed, allowing us to visualize systolic dynamics in 3D for the first time with ultra-high temporal resolution and very small errors in gated reconstructions (Section 4). Our gating and validation algorithms, in general, could be adapted to other imaging systems which are slower or faster than our current OCT setup. Although demonstrated on OCT images of the embryonic heart, our technique is applicable to any study involving periodic data acquisition, e.g., respiratory studies.
The rest of the paper is organized as follows. We describe methods for retrospective gating and volume-based cross-validation in Section 2. Section 3 presents results from retrospective gating of a stage 14 quail embryonic heart in the form of 3D volume renderings and movies, followed by results from our validation scheme. In Section 4, we make some observations and discuss limitations of our proposed image-based retrospective gating algorithm. Finally, Section 5 concludes the paper.

Imaging protocol
Fertilized quail eggs were incubated in a humidified, forced draft incubator at 102 °F. At approximately 48 hours of development the eggs were taken out of the incubator, the eggshell was removed, and the contents were cultured in a sterilized 3.5 cm diameter Petri dish. Once in the Petri dish, surviving embryos were placed in our environment-controlled OCT imaging chamber [35] and the hearts were imaged according to the following scanning protocol. Our OCT system employed a buffered Fourier Domain Mode Locked (FDML) laser [3,36,37] operating at a 117 kHz line rate. The axial and lateral resolution of the system was ∼10 μm in air. Our protocol consisted of capturing 120 B-scans (500 lines/B-scan or 4.27 ms/B-scan) spaced 2.72 μm apart over a 512 ms time period (∼1.5 heartbeats for a stage 14 quail embryo) at one slice position. After cross-sectional images were captured in this first slice position, the data (206 MB) was quickly transferred (2-3 s) and images were recorded at the next slice position. A total of 80 slices spaced 9 μm apart were recorded for the whole volume. In order to accurately reconstruct the images, it was important to make the slice step smaller than the lateral resolution of the system. This process was repeated until an entire 4D dataset was recorded. For our real-time reference volume scans, we captured images of the same heart in real-time with a reduced spatial sampling rate (150 lines/B-scan or 1.28ms/B-scan) resulting in about 4 vols/heartbeat or 10 vols/s.

Image-based retrospective gating
Our image-based retrospective gating algorithm consists of multiple steps ( Fig. 1). They are: (i) data decimation, (ii) cardiac cycle period estimation, (iii) data reassembly according to estimated cardiac period, (iv) linear interpolation (in time) of reassembled data, (v) relative shift estimation, (vi) absolute shift estimation, and (vii) reordering of data to create 4D data suitable for multiple visualization platforms.
We first spatially decimate B-scan images in x and y. We typically take every 4 th sample (see Section 4) in x and y, a decimation factor of 4, giving 1/16 the amount of image data We next determine cardiac cycle time or period. Let I m denote image data and x, y, z, τ denote the image x coordinate, image y coordinate, the z-slice, and the frame time, respectively. Let T / be an estimate of the cardiac time period in milliseconds and N τ be the total number of frames acquired. First, any frame time in the original data which is greater than T / is reassigned a time within the interval [0, T / ] using the Eq. (1) below [33] where ⌊ ⌋ refers to the floor operator. (1) After this reassignment, frames after T / are wrapped around to the beginning of the sequence at the new times, τ / . Let i = i(j) denote the function that maps indices from the original sequence, j, into indices in the "wrapped-around" sequence, i, containing images at unequally sampled time points. We formed an objective function D [Eq. (2)] which consists of two terms (i) a sum of squared differences term involving spatial image data from consecutive frames in the wrapped-around sequence, and (ii) a sum of squared differences term involving the difference of frame times in the wrapped-around sequence. The first term helps in maximizing image similarity while the second term ensures a wrapped-around sequence with frame times that are close to each other. The objective function for cardiac period estimation is provided below [33]: (2) In Eq 2, m = ( x, y ) ∈ Z 2 denotes the set of x and y coordinates of image frames. The optimal period value was obtained by determining a value of T / at each slice position that minimizes the objective function in Eq. (2). To optimize the objective function, we used the golden search algorithm [38,39], which was constrained to find a value of T / that lies in a time interval [T min ,T max ], with the interval chosen based on our experience with quail hearts. The result of this period determination is an array of cardiac time period estimates, one for each z position.
Next, we reassembled the temporal image sequence at each slice based on the optimal period, T / , using Eq. (1) and the wrap-around described above. For example, assuming a temporal image sequence of 120 frames and a 4.27 ms frame interval time, the original frame times for each image in the sequence would be 0, 4.27 ms, 8.54 ms, ..., 508.13 ms. If the estimated cycle time is 455 ms, all frames from the original sequence with frame times > 455 ms were reassigned values as per Eq. (1), and wrapped around as described above, resulting in a sequence with the following frame times: 0, 1.89 ms, 4.27ms, 6.16ms, 8.54ms,..., 448.35 ms, 452.62 ms. This resulted in an unevenly spaced (in time) sequence of images at each slice position in the volume.
A linear interpolation was then applied on the reassembled data to produce a predetermined number (L) of equally spaced in time image frames for each z position. The original data set had a 4.27 ms frame interval and 120 frames at each slice position over approximately a beat and a half. We set L = 90 in our experiments resulting in an even spacing of 4.9 ms per frame after interpolation.
Next, we determined the correct relative time shift between the image sequences at two adjacent slice positions. We start with a slice corresponding to the middle of the tubular heart and move outward in each direction. Using interpolated data from each slice position, and assuming periodicity of data, an image cross-correlation-based similarity measure Q is computed between the two time sequences of images at the current and the immediately adjacent slice in z, the volume scan direction. (3) where x, y, z, and t denote the image x & y coordinates, the z-slice, and the frame time respectively. k, k / denote respectively neighboring slice positions, s is the relative shift in number of time frames, and R is an operator that computes the 2D correlation coefficient between two images. The correct relative time shift between image sequences at the two slice positions is determined by maximizing the similarity measure in Eq. (3). It is also possible to formulate a more generalized similarity measure by including image data from a set of neighboring slices spaced different distances apart from the current slice: (4) where k (i) denotes neighboring slice spaced i slices apart from current slice. This measure is potentially more robust since it employs multiple slices to determine relative shifts. However, following some experiments in which we performed visual evaluations of reconstructions followed by a quantitative verification (see Section 2.3 and Section 4), we have determined that the measure in Eq. (3) suffices.
Next, in order to obtain absolute shift values for time sequences of images at each slice position, we first set the absolute shift [S m in Eq. (5) below] of the slice corresponding to the middle of heart tube [m in Eq. (5) below] to a constant value λ. Next, we progressively apply the relative shifts [estimated by maximizing Eq. (3)] to each consecutive slice position by moving outward in both directions from the middle slice. This gives us an array of correct starting offsets, one for each slice position. The mathematical representation of absolute shift computation is given below. (5) Last, the retrospectively gated 4D data is produced from full-resolution data using the cardiac cycle time estimates and the absolute shifts obtained at each slice position within the volume. The full-resolution data is reassembled and wrapped around using the optimal cardiac cycle estimates [Eq. (2)]. Then, a linear interpolation is applied on full-resolution data, similar to the one used for decimated data. The absolute time shifts computed using Eq. (5) are then applied to the interpolated data at each slice position. Each 3D volume in the 4D gated data set is generated from the set of slices corresponding to a given frame time.

Validation
A quantitative evaluation of our gating algorithm was performed using a reference real-time volume scan of the same embryonic heart. In our real-time acquisition, we repeatedly scanned the volume as fast as possible in z along the volume scan direction. These reference volumes provided a mechanism for validating the gating algorithm based on comparison of acquisition time difference between two slices of a given reference scan with the time difference between matching volumes from the gated data [ Fig. 2(a)]. For the real-time volume scan, we collected 150 lines/B-scan with a ∼9 μm sampling step in B-scan direction, and a 9 μm sampling step in z along the volume scan direction, resulting in ∼4 vols/heartbeat (or 10 vols/s assuming a healthy heart rate of 2.5 beats/s). For quantitative validation of gating, we compared the differences between frame acquisition times at any two slice positions in a given reference volume, with the acquisition time difference between matching frames from the gated data volumes at the same two slice positions. The reference and gated B-scans were matched using image registration with cross-correlation as a quality metric. Typically, motion of cardiac tissue is more significant when compared to change in speckle across time. Therefore, speckle decorrelation (with time) has a minor effect on our image registration and matching process, a fact that has been verified visually using the image matching technique of Fig. 2. Since the reference volume scan was acquired at a constant line rate of 150 lines/ B-scan, we expected to see a linear behavior when best matching volume indices from gated data were plotted against spatial slice positions. Generally, a better numerical agreement of the slope of gated data with the slope of acquired reference volume scan implies better gated reconstruction. The degree of conformity to linearity was additionally verified using the R-squared or "goodness of fit" measure in linear regression analysis.
A brief description of our image comparison scheme for validation follows next. Let the slice positions be denoted by z m , m = 1,2,...., N z , where N z is the total number of slices. Let τ m denote frame acquisition time for slice z m in a predetermined 3D reference volume scan. Let L be the total number of time frames at each slice position in gated data. Then, for each slice position, an image frame from the reference volume denoted by I r (x,y,z m ,τ m ) was compared with the temporal sequence of images for the same slice position in the gated data denoted by I g (x,y,z m ,t k ), k = 1,2,....L to create a 2D correlation coefficient (R). In each case, the R estimate was maximized by performing rigid image registration between the gated and reference frames assuming only minor x and y translations between the two frames. For the entire volume, this resulted in a one-dimensional array of correlation coefficients denoted by Q below. The matching image and the corresponding frame time t / at each slice position was determined by finding the maximum value in Q. The frame time t / was then used to compute an index into the time series of volumes constituting the 4D gated data set. This matching volume index was plotted against the slice number. Mathematical equations for validation are provided below: (6) In Eq. (6), x, y, z, and t denote the image x & y coordinates, the z-slice, and the frame time respectively. The subscripts r and g denote the reference and gated data sets respectively.

3.1Image-based retrospective gating
We applied the proposed image-based retrospective gating algorithm to out-of-order data collected from a stage 14 embryonic quail heart. The input data consisted of 120 frames (frame interval = 4.27 ms, or 234 frames/s) collected at 80 slice positions along the heart volume with our FDML-OCT system. With retrospective gating, 90 3D volumes across a single cardiac cycle were generated making it possible to visualize systolic dynamics in a 2D en face slice over the cardiac cycle taken 313 μm deep in the tissue. The movie shown in Fig. 3 was created by applying noise reduction on the series of 2D images across the heartbeat using the nonorthogonal wavelet and adaptive Laplacian pyramid non-linear diffusion image denoising algorithm [1] previously proposed by our group. The resulting movie clearly shows the dynamics of the myocardium, the endocardium, and the cardiac jelly. Figs. 4(a)-4(c) show three movies of the beating heart with cutaway views in the three orthogonal planes that have been named with respect to the embryo coordinates -sagittal (en face), coronal (OCT imaging plane), and transverse (along the heart tube). The high temporal resolution afforded by gating enabled us to clearly visualize systolic dynamics of the beating avian embryo heart in 3D for the first time. Visually perceptible differences in contraction at inflow and outflow close to systole have been observed for the first time [ Fig. 4(c), movie], which was not previously possible with lower temporal resolution scans. In the movie of Fig. 5, a digitally reconstructed radiograph (DRR) volume rendering technique has been applied to a region-of-interest corresponding to the heart tube. It involves summing of voxel intensities along rays cast through the volume perpendicular to the viewing plane. A movie of the beating and rotating heart volume clearly shows the extent of the heart tube, direction of blood flow and the interaction of the heart wall with surrounding tissue. The 2D and 3D visualizations described in this section have been created using Amira, a popular, commercially available 3D visualization software [40].

Validation
Next, we show results from validation of our gating algorithm. First, Fig. 2 shows a visual example demonstrating the accuracy of image-based gating. After gating, we found that the matching volume index in gated data for slice 45 from the reference scan A is 24 [ Fig. 2(a)]. This fact was visually confirmed by comparing slice 45 (middle of heart tube) from both data sets. A high image similarity was found [Figs. 2(d) and 2(e)]. When we compared a different slice from the real-time reference volume to the same gated volume, the images looked dissimilar since they belonged to different cardiac phases [ Figs. 2(b) and 2(c)].
We next demonstrate our quantitative validation experiment (Section 2.3) using two real-time, reference full volume scans (denoted by A and B) at different points of cardiac cycle acquired as per the protocol described in Section 2.1. These two volume scans were chosen because they corresponded to phases of cardiac cycle around which significant tissue movement occurred. Figure 6 shows the plot of matching volume numbers from gated data against spatial slice position for two reference volumes scans denoted by A and B. A linear fit was obtained through data points corresponding to slices containing the tubular heart, thereby validating our gating algorithm. The slope of the plot for reference scan A was 1.364 ms/slice and that for reference scan B was 1.120 ms/slice. We observe that these values closely match the ideal slope of the reference scan (1.275 ms/slice), thereby proving accuracy of our gated reconstructions.

Discussion
Experiments with retrospective gating of embryonic quail heart data were very encouraging. We found that both visual (Figs. [3][4][5] and quantitative assessment (Fig. 6) were important for evaluating our gating algorithm. A few remarks about our gating algorithm follow.
We investigated the effect of applying spatial decimation (along the x and y directions of the B-scan frame) to reduce computation time. For example, to spatially decimate by 4, we pick up every 4 th pixel in both x and y. We often found that cropping the B-scan images to include only the region corresponding to the heart also helps in reducing computation time. We compared cardiac period estimates at all slice positions obtained from original and spatially decimated data (with various decimation factors, e.g., 2, 4, 8, and 16). We observed that up to a factor of about 4, the cardiac period estimates agreed numerically with those obtained from native resolution data. This suggested a lower limit on the sampling rate equal to ¼ of the native sampling rate used for out-of-order data acquisition. We then produced gated reconstructions from original resolution and each of the decimated data sets and applied the validation algorithm described in Section 2.3 in each case. The plots obtained from validation of data with different decimation factors (Fig. 7) helped us establish a lower limit on the spatial decimation of input data to ensure accurate gated reconstructions. A spatial decimation factor > 4 clearly resulted in large departures from the ideal slope (i.e., slope of reference scan) indicating incorrect gated reconstructions (Fig. 7, bottom left and right). The limit on spatial decimation was therefore deduced to be equal to 4. This analysis revealed how quickly we could process the data (through the use of a decimation factor) without causing errors in resulting reconstructions.
Next, we determined the effect of temporal decimation on cardiac period estimates at all slices. We have found that using a temporal decimation factor of 2 (i.e., skipping every other frame) leads to negligible errors in estimated period (the average error in the estimates over the 80 slices was found to be 0.005 seconds, a negligible number compared to the ∼0.45 s cardiac cycle time). However, a temporal decimation factor of 4 and above caused larger errors (>0.025 seconds). With slower systems these larger errors can be reduced by recording over more heartbeats at each slice location. We have typically recorded 1.5 heartbeats at each slice in our standard imaging setup. However, there could be small variations across successive heartbeats even for the same embryo. In a separate imaging experiment with a slightly older quail heart, we recorded just over 2 heartbeats at each slice (120 slices/volume) and quantitatively evaluated the effect of multiple heartbeats on gating accuracy. First, we estimated cardiac cycle time using the entire data (∼2 heartbeats). We then dropped the last ½ heartbeat from each slice position, thereby using ∼1.5 heartbeats/slice to estimate cardiac cycle time. The mean cardiac cycle time across 120 slices was 0.3056 s in the former case, while it was 0.3079 s for the latter case. The difference between the two estimates was 2.3 ms, or, just over half the duration of a frame, a clearly acceptable number. We concluded that our gating algorithm was equally accurate with 2 heartbeats/slice as it was with 1.5 heartbeats/slice. We have further determined that having a shorter temporal sequence at each slice (i.e., < 1.5 heartbeats) sometimes caused bimodal distributions in the correlation-based objective function during relative shift estimation. This in turn leads to false maxima and erroneous relative shifts as verified visually. In the limiting case where we have less than a single heartbeat of data at a slice, the relative shift computation algorithm fails for lack of enough image evidence to determine an accurate relative shift value.
We investigated the errors obtained in absolute shift estimation when we increased the spacing between B-scan frames (or slices) used for computing image correlations. We observed (through visual examination of gated reconstructions) that relative shift estimates obtained by correlating immediately adjacent slices (i.e., native sampling rate in volume direction) were accurate. However, when a slice was correlated with another one spaced two sampling steps away (i.e., ½ of native sampling rate), there were large numerical differences in absolute shift estimates, and the resulting error was confirmed through a visual evaluation of the resulting reconstructions. We concluded that correlating slices separated by a distance more than the spot size of the imaging system (10μm) resulted in erroneous gated reconstructions.
We discuss the possible sources of error in our validation. The small discrepancies in slope value for the validation technique can be attributed to a couple of issues. First, the real-time data were sampled spatially at a lower rate as compared to the out-of-order data, which would cause the small correlation errors leading to errors in slope. Second, the reference and the outof-order data sets employed different scanning geometries and small positioning errors due to minor galvanometer inaccuracies would lead to small correlation errors. However, we found that in spite of these possible sources of error, the frame time offset between the starting and ending slices corresponding to the heart region in the gated data very closely matched that of the reference scan, thereby validating our image-based retrospective gating scheme. We computed errors between actual points on the validation plot and those on a straight line with the ideal slope (i.e., slope of reference scan) for the slices corresponding to the heart region. The following errors are actually worst-case estimates because of the mechanical limitations of the scanning described above. The maximum error was found to be 3.83 frames or ∼18.7ms ms, and the standard deviation of the error for slices in the heart region was found to be 0.9669 frames or 4.7ms (time interval between volumes in gated data = 4.9 ms, see Section 2.2). Systolic dynamics in the heart occur in a time interval of the order of 50 ms; our gating algorithm has an error standard deviation which is an order of magnitude smaller and therefore has sufficient accuracy for visualizing systolic dynamics.
We comment on computational complexity of our algorithm. As compared to our previous (prospective) LDV gating study, the additional benefits gained from retrospective gating far outweigh the longer post-acquisition computation time needed for reconstruction ( Table 1).
As mentioned before, a spatial decimation factor of 4 caused negligible errors in reconstruction. We present results of computation time estimates (Table 1) for some important steps of the gating algorithm on our test data set (80 slices, 120 images of 512 × 500 pixels per slice) with spatial decimation set equal to 4. Without any code optimization, MATLAB® code for cardiac cycle time estimation, data interpolation, and shift estimation took 70 seconds altogether on a Dell Precision 690 workstation with a 3.0 GHZ dual-core Intel Xeon processor and 32 GB of system RAM. Data reassembly, wrap-around, interpolation, and rearrangement on full resolution data to produce the final gated data took about 45 minutes. We are currently performing some code optimizations for reducing the running times of each of the steps above.
In an earlier paper, we have determined using heart wall velocity estimates obtained from realtime volume scans that there could be large tissue displacement errors, as high as 205 μm [3]. With image-based retrospective gating, the reconstructed data set has a maximum displacement error of only 8.5 μm, which would now allow us to see differences in wall velocity at the inflow and outflow of the heart tube (Fig. 4). We had earlier observed with our lower temporal resolution scans [2][3][4] that the contraction appeared to be the same at the inflow and outflow. The increased temporal resolution afforded by our gating algorithm will facilitate precise measurements of many cardiac parameters (wall velocity, contraction propagation, etc.) in the developing heart tube and will allow us to develop more accurate models of cardiac dynamics at these early stages.

Conclusion
We have developed an algorithm for image-based retrospective gating of OCT image data and applied it successfully to reorder out-of-order (2D + time) OCT image data collected from the early embryonic quail heart (stage 14) into a series of heart volumes (4D data). Our scheme works even with decimated spatial image data and does not require any prior denoising, compression, or transformation of data. The proposed robust validation scheme provided quantitative evaluation of retrospective gating quality and helped us establish limits on data decimation in order to guarantee accurate gated reconstructions. Validation of image-based retrospective gating has revealed that errors in reconstruction are very small. Standard deviation of error = 4.7 ms, which is an order of magnitude smaller than the time interval during which systolic dynamics (∼50 ms) occur in the developing embryo. Our gating algorithm has helped us visualize 4D cardiac systolic dynamics for the first time. The extremely high effective temporal resolution afforded as a result of our gating algorithm will allow us to track blood cells and minute, rapid tissue displacements thus allowing us to study mechanical forces regulating tissue homeostasis in the embryonic looping heart.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. A schematic of our image-based gating algorithm. As shown, the major steps are: (i) decimate images, (ii) determine cardiac period at each slice position, (iii) reassemble or wraparound data to produce one correctly ordered, unevenly spaced set of images from a cardiac cycle at each slice, (iv) interpolate linearly to create even temporal spacing, (v) determine relative shift from one slice to the next (vi) determine absolute time shift at each slice, and (vii) reorder image data to create volumes over the cardiac cycle.    4D Visualization of retrospective gated OCT data from stage 14 embryonic quail heart. Movies of 4D gated data (90 volumes in a single cardiac cycle, or ∼270 volumes/s) are shown in sagittal or en face (a) (2.7 MB), transverse (b) (2.9 MB), and coronal (c) (3.2 MB) views. In each case, a cutaway view of the beating heart with a moving orthogonal 2D slicer in the plane of interest is shown. (3.4 MB) Digitally reconstructed radiograph (DRR) volume rendering movie created using gated data from stage 14 embryonic quail heart. A summing of voxel intensities along rays cast through the volume perpendicular to the viewing plane has been used to create a 3D effect with the rotating and beating heart volume. Validation of retrospective gating algorithm using two different reference volume scans A (left) and B (right). A linear fit is obtained through matching volume numbers in gated data as per the method described in Section 2.3. The slope of the fit obtained on gated data is 1.364 ms/ slice for scan A and 1.120 ms/slice for scan B. The slope of reference scan is 1.275 ms/slice. Validation of image-based retrospective gating using input data with different integer decimation factors. In each case, the B-scan frame was decimated by an integer number n in both x and y directions. Following decimation, the gating algorithm (Section 2.2) and the proposed validation scheme (Section 2.3) were applied. We observe that, up to a decimation factor of 8, the slope of the validation plot (slice number versus matching volume number in gated data) closely matches the ideal slope from the reference volume. Computation times for different steps of the image-based retrospective gating algorithm.