{'fraction_non_alphanumeric': 0.03652448103705001, 'fraction_numerical': 0.009985109923797844, 'mean_word_length': 4.853627275057677, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sequential Turning Acquisition and Reconstruction (STAR) method for four-dimensional imaging of cyclically moving structures

Optical coherence tomography allows for dynamic, three-dimensional (3D+T) imaging of the heart within animal embryos. However, direct 3D+T imaging frame rates remain insufficient for cardiodynamic analysis. Previously, this limitation has been addressed by reconstructing 3D+T representations of the beating heart based on sets of two-dimensional image sequences (2D+T) acquired sequentially at high frame rate and in fixed (and parallel) planes throughout the heart. These methods either require additional hardware to trigger the acquisition of each 2D+T series to the same phase of the cardiac cycle or accumulate registration errors as the slices are synchronized retrospectively by pairs, without a gating signal. Here, we present a sequential turning acquisition and reconstruction (STAR) method for 3D+T imaging of periodically moving structures, which does not require any additional gating signal and is not prone to registration error accumulation. Similarly to other sequential cardiac imaging methods, multiple fast image series are consecutively acquired for different sections but in between acquisitions, the imaging plane is rotated around the center line instead of shifted along the direction perpendicular to the slices. As the central lines of all image-sequences coincide and represent measurements of the same spatial position, they can be used to accurately synchronize all the slices to a single inherent reference signal. We characterized the accuracy of our method on a simulated dynamic phantom and successfully imaged a beating embryonic rat heart. Potentially, this method can be applied for structural or Doppler imaging approaches with any direct space imaging modality such as computed tomography, ultrasound, or light microscopy.


Introduction
Studies of cardiac development and function increasingly rely on the ability to image the beating embryonic heart in three dimensions at high frame rate. Optical coherence tomography (OCT) allows for non-destructive imaging through several millimeters of biological tissue with single cell resolution. It is being adopted as an accurate imaging method to study in vivo systems, as it provides both information on structure and dynamics. OCT has been used to image whole mammalian and avian embryos [1,2], measure flow in the vasculature [3,4], and image the mammalian and avian embryonic hearts in 4D [5][6][7].
Direct acquisition of volumetric (3D) OCT data, which involves scanning a beam along two directions and remains too slow for imaging fast moving structures such as cardio-vascular structures or the beating heart. Insufficient frame rates lead to image artifacts, such as motion blur, which prevent fine characterization of morphology and function. However, since cardiac structures move periodically, indirect imaging methods can be used to produce 3D+T renderings at high frame rates. One approach consists in acquiring 2D+T image bursts in fixed planes throughout the heart, each lasting the duration of a few heartbeat cycles. When the imaging planes (all parallel) are separated by a small distance, each sequence can be synchronized to its neighbors to finally produce a 3D+T data set of the beating heart with a frame rate identical to that of the 2D+T sequences (Fig. 1). This approach has been used in confocal microscopy [8,9] and OCT [5,6]. This technique is particularly attractive because reconstruction is possible based only on the imaging data itself. However, since the reconstruction proceeds by synchronizing image sequences recursively from one imaging plane to the next, this method is susceptible to error propagation.
Alternatives to this image-intrinsic reconstruction is to use a separate gating signal so imaging can be triggered at a precise time in the cardiac cycle. Known as prospective gating in magnetic resonance imaging (MRI) and X-ray computed tomography (CT) [10][11][12][13][14], gated cardiac imaging does not rely on post-acquisition synchronization but usually requires additional hardware to collect the gating signal. In OCT, gated imaging has been implemented by at least two research groups [15,16], yet it is particularly challenging since two independent OCT scanning heads have to be aligned around the tiny embryo.
Techniques that can take advantage of a single gating signal yet do not require additional hardware are therefore desirable. In CT, Dinkel et al. have proposed to use an intrinsic gating signal obtained from the projections restricted to an ROI that includes the heart [17] and socalled kymograms [11]. ECG-free algorithms have also been applied to ultrasound imaging [18,19]. ECG-free algorithms have also been applied to ultrasound imaging [18,19]. In MRI, Spraggins [20] describes a method in which two image series are acquired alternatively, one in a fixed common area, serving as the gating signal. Alternative k-space sampling procedures have been proposed, which use an interleaved radial k-space sampling trajectory such that each readout samples the k-space center, providing an intrinsic gating signal [21][22][23].
In OCT, Liu et al. proposed to adjust the phase lags between set of parallel planes, by acquiring, in addition to the series of parallel slices, one orthogonal slice [7]. The parallel planes are first aligned with each other then the alignment of each parallel plane is adjusted based on the line in the orthogonal plane with which each plane intersects. Since all parallel planes share a  Fig. 1. Sequential acquisition using parallel slice geometry [5,8,9]. (a) Direct 3D volume is obtained by collecting multiple B-scans (y-planes) (b) Sequential imaging of fast B-scans at a fixed y-positions (c) If triggered at the same time in the heart cycle, or subsequently synchronized, the frames can be re-assembled to form a 3D+time series at high frame rate. different line with the orthogonal reference plane, the orthogonal slice must be chosen carefully and may not intercept the beating heart in its entirety.
Here, we present a method-sequential turning acquisition and reconstruction (STAR)-for imaging periodically moving structures that does not require measuring an additional gating signal yet still overcomes the error accumulation limitation introduced by gating-free recursive synchronization procedures. Our STAR approach proceeds similarly to nongated approaches where multiple fast images series are consecutively acquired, but instead of considering a set of parallel imaging planes, it involves 2D+T series acquired in imaging planes that share a common axis ( Fig. 3(a)). In STAR, all planes share the same common line. This line is chosen, during acquisition, to be of good quality (e.g. in the center of the heart). This geometry is made possible by the fact that traditional OCT data are acquired as a collection of one-dimensional in-depth scans (axial scans or A-scans) serving as single units of information. Two-dimensional images (or B-scans) are obtained by laterally moving the OCT beam with two galvanometercontrolled mirrors (see Fig. 2). Since the direction of the A-scans is fixed along the z-direction the B-scan planes are always vertical. However, the user can freely adjust the orientation and lateral position of the plane by controlling the mirrors' oscillations.
Also, since the central line of all 2D+T time-sequences represents measurements taken over the same physical line (the rotation axis), the central 1D+T signal from each sequence can be used for synchronization. The central lines serve as intrinsic temporal gating signals for postacquisition synchronization. After the synchronization, the data are transformed from cylindrical into orthogonal coordinates to reveal the 3D+T representation of the periodically moving structure. Our method bears some similarity with a recently presented method by Happel et al. [24] for OCT of chicken embryos. There, images are acquired in a geometry identical to ours (with a common A-scan) after which each frame in a sequence at a given angle is assigned a new time stamp (based on an algorithm that assumes perfect, yet unknown, periodicity in the sequence). Although this algorithm can handle variations in the heart frequency that occur in between angles, it does so by reshuffling the images, producing a new sequence in which images from multiple heartbeats are interleaved. Our reconstruction algorithm differs in that it applies a nonuniform, elastic registration that preserves the order of the frames and allows reconstruction even in the presence of slight, local variations in the heartbeat that occur during acquisition at a single angle. By preserving the sequence order, our algorithm produces reconstructions in which aperiodic motion features within a sequence at a given angle, such as moving red blood cells, could be analyzed to infer lateral blood flow velocities within one plane (e.g. via block matching or optical flow techniques).
In Section 2 we formally introduce the acquisition and reconstruction method. Then, in Section 3, we present validation reconstructions obtained with the parallel-slice and proposed STAR method. We also present experimental results from data collected in the beating heart of an embryonic rat.

Method
We consider a dynamic, three-dimensional object, whose measured intensity I(x, y, z,t) is periodic, translating to where T is the period of the cardiac cycle. Since, in practice, signals are never perfectly periodic, we relax this constraint by assuming that, given a reference period T , the length of other periods may only vary by an estimated percentage α of the reference period such that the minimal and maximal lengths of the cycles are T min = (1 − α)T and T max = (1 + α)T , respectively. For cardiac imaging in live embryos, intensity is varying rapidly with time (the embryonic heart rate is approximately 2-3 beats per second) and instantaneous acquisition of threedimensional data in the spatial dimension (corresponding to ideal sampling at a fixed point in time) is technologically challenging. Rather, data along the spatial dimensions are acquired sequentially to build a multi-dimensional image. Images are built by sequentially scanning individual voxels, lines, or, at best, full planes. Here we assume that a 2D plane (B-scan) can be acquired in a single shot with high repetition rate. Specifically, we assume that the imaging system allows for instantaneous acquisition of a full two-dimensional plane parametrized by the coordinates (r, z), defined by a rotation of angle θ ∈ [0, π) around the z-axis of the plane y = 0 (see Fig. 3 (a)). We consider the set of measurements for rotation angles θ ∈ [0, π), radial coordinates r ∈ [−r M , r M ], and axial coordinates z ∈ [0, z m ], over time t ∈ [0, L), each sequence being triggered at a time offset t θ . This set of measurements corresponds to the sequential acquisition of 2D+T image sequences of duration L and is schematically depicted in Fig. 3(b). The acquisition of these sequences is not triggered at any particular time in the cardiac cycle (non-gated imaging), but to ensure that at least one complete cycle of the heartbeat (starting at any arbitrary time in the cycle) can be extracted from each sequence, the acquired sequences must be of length L ≥ 2T max , that is, they must be twice the duration of the longest heartbeat period estimate. These sequences are acquired in a radial, star-like pattern, with all planes intersecting at and sharing the z-axis (see Fig 3 (a)). The first period of the sequence-I 0 (r, z,t) for t = [0, T ) (in the plane y = 0)-is considered the reference and the other sequences are synchronized to it. The nonrigid synchronization problem is equivalent to finding a set of angle-dependent warping functions w θ (t) for all rotation angles θ ∈ [0, π) such as to minimize, along the shared z-axis (that is, the axis defined by r = 0, see Fig. 3(c)), the functional The first double integral over time and the axial direction z compares a temporally warped sequence (a sequence whose temporal axis is deformed) I θ (0, z, w(t)) with the template I 0 (0, z,t), while the second integral keeps the extent of this deformation in check by penalizing warping functions w(t) 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 (λ = 0) or the temporal integrity of the warped sequence (λ → 1). We note that the warping functions w are in the set M = {w ∈ C 1 ([0, T )) | 0 ≤ w(t) < L and 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 ). We then determine a warping function that minimizes the cost function. This minimization is carried out using a previously described dynamic programming algorithm, which traces a minimum-cost path in a two-dimensional similarity graph [25]. For each angle, the optimal warping function produced by this algorithm is For example, if λ = 1, the matching is rigid and warping functions are of the form w(t) = t 0 + t, a uniform offset in time. As λ < 1 decreases, the resulting warping functions can become non-linear to better match the data.
Once the warping functions are determined, we compute time-interpolated series (according to the temporal warping functions w θ ) for the measured intensity series I θ for each angle θ (see Fig. 3 Although the synchronized series could be visualized per se (see Figs. 3(e) and 6(a)), it is usually desirable to have a volumetric representation of the data in Cartesian coordinates (see Fig. 3(f)). The final step of the reconstruction therefore consists in converting, at each time point t ∈ [0, T ) from a cylindrical geometry (θ , r, z) to a Cartesian grid (x, y, z). Given samples I θ i (r j , z,t) for uniformly spaced angles θ i = i∆θ , 0 ≤ i < N θ and uniformly-spaced radial coordinates r j = −r M + j∆r, 0 ≤ j < N r (where ∆θ = π/N θ and ∆r = 2r M /(N r − 1) are the angular and radial spacings, respectively), we can define a cubic spline interpolating functioñ where the centered cubic B-spline function β 3 is defined as the 4-fold convolution of β 0 , the indicator function over the interval [−1/2, 1/2], and where the B-spline coefficients c i, j (z,t) in the expansion (5) are obtained by recursive filtering [26] such that the interpolating function interpolates the data at the sample points Finally, the recovered intensity I r on a Cartesian grid is given by where the Cartesian to polar coordinate transform is given by where −π < atan2(y, x) = Prarg(x + iy) ≤ π is the principal argument of the complex number x + iy, which is restricted to the interval (−π, π].

Results
To evaluate the reconstruction accuracy of STAR and illustrate the advantages over sequential parallel-slice acquisition and reconstruction, we carried out a simulation experiment with a synthetic, periodically deformed heart tube phantom, presented in [8] and shown in Fig4(a). Briefly, the intensity of a static tube is defined for any real value of x, y, and z in three-space by a parametric analytic equation. To generate intensity volumes over time, an affine transform of the spatial coordinates at which the intensity is to be evaluated is carried out, and the heart-tube function is evaluated at the transformed coordinates. The affine transform is defined by a matrix whose coefficients are given by a weighted sum of sinusoidal and cosinusoidal functions at harmonic frequencies. Here, acquisition of data was simulated either by evaluating the phantom values sequentially for different y-sections (see Fig.4(d)) or by acquiring sections in a radial geometry (Fig. 4b) to match our proposed acquisition conditions. For each time-series, we introduced a random time-shift (uniform distribution in [−T, T ]), similarly to the approach in [8].
We next reconstructed the data in the STAR and conventional geometries (Fig. 4(c) and (f), respectively). We compared the resulting data volumes with the reference data volumes and computed the absolute difference voxel by voxel resulting in a difference data volume (see Fig. 4(d) and (g) for the difference images corresponding to reconstructions from sequential parallel-section imaging and STAR imaging, respectively).
By randomly selecting (following a zero-mean normal distribution with σ = 0.1) the weights of sine and cosine functions up to the third harmonic (higher harmonics were set to zero), we generated 100 dynamic reference phantoms as well as corresponding sequential parallel-slice and radial measurement series (50 slices in each case), introducing random time-shifts (following a uniform distribution in [−T, T ]) for each section. We then ran our algorithms to synchronize the slices and recover volumes for both acquisition geometries. Since the ground-truth shifts were known, we could accurately estimate reconstruction errors: for each synchronized slice time series, we compared the local time-shifts of the frames produced by the warping functions, and subtracted them from the actual shift introduced during the simulation. We then calculated the mean of the absolute value of the errors of all frames in a given section. Figure 5  shows the average synchronization error for image series at a given depth or orientation (computed over the 100 simulated experiments) as a function of distance from the reference slice (in the case where volumes were reconstructed from parallel slices) and angular distance from the reference slice (in the case of STAR).
In the case of recursive reconstruction, the average time-synchronization error increases as the distance to the reference slice increases. This is visible in Figure 5. Visually, this translates to high difference images toward the top and bottom in Fig. 4(g)). This degradation can be explained by the fact that the parallel section algorithm proceeds recursively: with 50 slices to align, the error propagates (the central slice was used as a reference). Conversely, for the STAR reconstruction the synchronization proceeds non-recursively (only the central portion of each slice is aligned to the central portion of the reference slice) the error remains in check (Fig. 5), even for large number of slices and for angles very different from that of the reference. Large differences visible in the STAR reconstruction error image near the top and bottom of the tube (see Fig. 4(d)) are due to the fact that the reconstruction is limited to a cylinder. Within the cylinder, the difference image does not show a systematic error. Rather, remaining visible errors are mostly due to the lower sampling step away from the center in the STAR geometry.
To evaluate the feasibility of our method in practice, we acquired image sequences of a live 10.5 day rat embryo using a swept-source OCT (SS-OCT) system as described previously in [27] (procedures were approved by the Animal Care and Use Committee of the University of Houston). Briefly, the experimental system was customized from the Thorlabs SL1325-P16. It utilizes a broadband swept-source laser with a central wavelength of λ 0 = 1325 nm and a bandwidth of ∆λ = 110 nm; output power is P = 12 mW; the A-line scanning frequency was 16 kHz. An interferogram is detected by a balanced-receiver configuration that reduces source intensity noise as well as autocorrelation noise from the sample Thorlabs, PDB140C and is digitized using a 14-bit digitizer. A Mach-Zehnder interferometer MZI-based optical frequency clock is used to calibrate the OCT interference signals from the optical time-delay domain in the frequency domain before application of fast Fourier transform FFT algorithms. An FFT is used to reconstruct an in-depth profile A-scan of OCT intensity from a single scan over the operating wavelength range. By 2-D mechanical scanning over the area of up to 10 × 10 mm, the 3-D structure of the embryo can be revealed.
Each time lapse consisted of 200 frames of 256 × 512 pixels (512 pixels per A-scan) acquired at 51 frames per second. The whole data set consisted of 60 image sequences; the angle between the neighboring time lapses was set to 3 • . The central line (axis of the imaging plane rotation) was positioned over the highly moving structure to ensure strong reference signal for synchronization. We were able to successfully reconstruct the dynamics of the beating heart in the embryo using the proposed algorithm ( Fig. 6 and Media 2, Media 3, and Media 4).

Conclusion
In this paper, we have provided an acquisition and reconstruction algorithm for imaging of dynamic cardio-vascular structures whose motion is cyclical. The STAR algorithm, relies on the intrinsic signal provided by an axis of data common to all slices. Our algorithm has several advantages, in particular over previously proposed non-gated imaging and reconstruction algorithms in which data are acquired over parallel sections instead of radial sections. Our method achieves accurate reconstructions without requiring any additional hardware for acquiring a gating signal. Unlike existing nongated synchronization approaches where slices are synchronized to neighboring slices and errors may accumulate with each round of synchronization, in the STAR method, there is no accumulation of the error since all slices are synchronized to the same signal. If the sample drifts during acquisition, this becomes obvious after reconstruction, since after 180 • rotation, the imaging plane is a mirror-image of the initial position. We demonstrated the method's improved accuracy on simulated data and its practical applicability for in vivo imaging of embryonic rats.
Since sampling occurs in polar coordinates, the spatial resolution is non-uniform: it is finer towards the axis of rotation and coarser as the radial distance increases. In our experiment, where the extent of the OCT scan is 400 µm (and represents the size of a 10.5 dpc embryonic rat heart) and the imaging plane is rotated by 3 • for each time series, the distance between neighboring slices at the radially most distant border is about 10 µm, which is similar to the step used for parallel-slice imaging approaches.
Use of an elastic synchronization method that does not assume the interval between heartbeats to remain constant ensures robustness to slight variations in heartbeat period, which we previously showed to be beneficial in the case of parallel-slice imaging [9].
Our method is general and could be applied for structural or any other (e.g. Doppler) imaging of periodically moving structures with optical coherence tomography, or other modalities that allow acquiring images according to the proposed geometry.