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 three-dimensional 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 ). 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 . 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* ≥ 2

*T*_{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 (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 ), 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 𝓜 = {

*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

*w*_{θ} {

*w* 𝓜|

*Q*{

*w*} = min

_{w}_{′𝓜}
*Q*_{θ} {

*w*′}}. 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 )

Although the synchronized series could be visualized per se (see and ), it is usually desirable to have a volumetric representation of the data in Cartesian coordinates (see ). 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

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* = 2

*r*_{M}/(

*N*_{r} − 1) are the angular and radial spacings, respectively), we can define a cubic spline interpolating function

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 (

*θ*_{i},

*r*_{j})

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* + i

*y*) ≤

*π* is the principal argument of the complex number

*x* + i

*y*, which is restricted to the interval (−

*π*,

*π*].