Open Access
1 September 2005 Four-dimensional cardiac imaging in living embryos via postacquisition synchronization of nongated slice sequences
Michael Liebling, Arian S. Forouhar, Morteza Gharib, Scott E. Fraser, Mary E. Dickinson
Author Affiliations +
Abstract
Being able to acquire, visualize, and analyze 3D time series (4D data) from living embryos makes it possible to understand complex dynamic movements at early stages of embryonic development. Despite recent technological breakthroughs in 2D dynamic imaging, confocal microscopes remain quite slow at capturing optical sections at successive depths. However, when the studied motion is periodic--such as for a beating heart--a way to circumvent this problem is to acquire, successively, sets of 2D+time slice sequences at increasing depths over at least one time period and later rearrange them to recover a 3D+time sequence. In other imaging modalities at macroscopic scales, external gating signals, e.g., an electro-cardiogram, have been used to achieve proper synchronization. Since gating signals are either unavailable or cumbersome to acquire in microscopic organisms, we have developed a procedure to reconstruct volumes based solely on the information contained in the image sequences. The central part of the algorithm is a least-squares minimization of an objective criterion that depends on the similarity between the data from neighboring depths. Owing to a wavelet-based multiresolution approach, our method is robust to common confocal microscopy artifacts. We validate the procedure on both simulated data and in vivo measurements from living zebrafish embryos.

1.

Introduction

Confocal laser scanning microscopy (CLSM) has emerged as a popular method for high-resolution imaging of fluorescent labels, particularly in thick or scattering samples. By placing a pinhole in the conjugate optical plane, before the detector, out-of-focus light from above and below the focal plane is rejected from the image, enhancing the axial resolution.1 By collecting images from defined optical slices at successive depths, the three-dimensional arrangement of fluorescently-labeled structures can be derived. In traditional point scanning confocal systems, images are collected in a pixel-by-pixel manner and acquisition speeds for sequences with frame size 512×512pixels are on the order of only a few frames per second. Recent advances in beam shaping, the availability of fast CCD line detectors, and the implementation of efficient hardware for data transmission have made possible the development of a fast laser scanning microscope, the LSM 5 LIVE.2 A blade-shaped beam focused to a line (instead of a single point for conventional CLSM) permits the parallel acquisition of a whole line of pixels and reduces the scanning dimensionality to one direction. This microscope allows for the acquisition of 2D image sequences (512×512pixels) at frame rates of up to 120 frames per second. This opens new avenues for a variety of fields.

In developmental biology, one major goal is to gain a better understanding of the mechanisms that influence the development of the cardiovascular system. In particular, it is desirable to assess the influence of genetic as well as epigenetic factors such as blood flow, heart wall forces, shear stress, etc.3, 4 While the frame rates of typical confocal microscopes are suitable to study many dynamic processes occurring in living systems (e.g., cell migration, division, etc.), cell motions in the cardiovascular system (e.g., heart-wall motions, blood flow, etc.) typically occur at several millimeters per second, 2 to 3 orders of magnitude faster than cell migration. The significant improvement in frame rate offered by parallel scanning systems now makes it possible to collect image data from single optical sections of fast-moving structures. However, resolving rapid three-dimensional motions in real time still remains a challenge because it is not currently possible to scan the z direction as fast as the xy plane.

Other imaging modalities such as magnetic resonance imaging (MRI), computerized tomography (CT), or ultrasound (US) suffer from similar limitations. However, if the imaged body undergoes the same deformation at regular intervals and the acquisition is always triggered at a particular phase in the cycle, it is possible to assemble the data to recover a whole volume over one full period. For larger organisms (from mice to humans), it is relatively easy to gate the acquisition with respect to electrocardiograms (ECG) or respiratory signals—a technique known as prospective gating or triggering—and reconstruct volumes at a fixed moment in the cycle.5, 6, 7 Remaining motion artifacts may then be reduced by the use of various elastic registration procedures that warp the spatial data.8, 9, 10 In cases where gating is not possible or unreliable, nongated dynamic datasets have been registered by a variety of methods and for various purposes. For instance, in nuclear medicine, noise reduction may be performed through temporal averaging of nongated signals.11 Thompson and McVeigh have used the imaging data from flow-encoded MRI to retrospectively perform the gating.12 Using specific modifications to conventional MRI pulse sequences, it is also possible to generate and extract a signal that varies in synchrony with the cardiac cycle for later reconstruction.13, 14, 15, 16 For CT, various methods have been developed, either to recover an imaged volume of the heart in a defined motion state at a single time point17 or for 4D imaging, by tracking the projection’s center of mass.18 ECG-free algorithms have also been used for US imaging.19, 20

In this paper, we present a technique and the associated image processing that make it possible to reconstruct dynamic 3D volumes of microscopic objects that are periodically moving, using currently available CLSM technology. We sequentially acquire slice sequences at different depths and reassemble them a posteriori to recover dynamic 3D volumes. For smaller organisms, such as the zebrafish embryos we study, reliable triggering signals to gate the acquisition are difficult and cumbersome to acquire. We have therefore devised a method for postacquisition synchronization based upon information within the recorded nongated data itself.

Our synchronization algorithm registers pairs of slice sequences with respect to time by minimizing a least-squares intensity difference criterion. The core element of our method is thus reminiscent of standard methods for image registration,21 the latter being a particular instance of the more general problem of motion estimation;22, 23, 24 however, the nature of the data requires special adjustments in order to achieve stable and repeatable results with minimal operator input. First, since our problem naturally calls for periodic boundary conditions, we have to (automatically) crop the data to cover a whole number of periods. Second, the method must be robust to various acquisition artifacts that are specific to confocal microscopy.25, 26, 27, 28 Last, the large amount of data that is involved imposes a subtle balance between memory and time constraints.

To solve some of the above requirements, our synchronization algorithms rely on the wavelet transform for robustness and rapid execution. For similar purposes, various authors have taken advantage of the favorable wavelet properties to implement affine and elastic registration algorithms of 2D or 3D datasets,29, 30, 31, 32, 33, 34, 35, 36 although their methods are not directly related to ours.

The paper is organized as follows. In Sec. 2, we discuss the measurement process as well as the requirements and limitations of the method. In Sec. 3, we present the synchronization algorithm, along with the tools that are required for period determination and noise reduction. In Sec. 4, we present two experiments, one based on a simulated acquisition procedure and the second, based on experimental in vivo measurements. Finally, we discuss the method’s overall performance and limitations as well as further developments.

2.

Measurement

2.1.

Measurement Strategy

We image a slice of a 3D object, subject to periodic motions and deformations over typically two to four periods. We assume that the object is given by the local intensity I(x,z,t)[0,Imax] , with x=(x,y) and that the periodic deformations are such that at any fixed spatial position (x,z) we have

Eq. 1

I(x,z,t)I(x,z,t+T)Imax,
where T is the deformation period. Between acquisitions, the object is moved axially ( z direction) to be sequentially imaged in its entirety (see Fig. 1 ). Because the acquisition is triggered at a random moment in the heart cycle, the different sequences are not synchronized. The measured intensity can be modeled by

Eq. 2

Im(x,zk,t)=I(x,z,tsk)PSF(xx,zkz)dxdz
where the sk are the unknown time shifts (defined modulo the period) at each measured depth zk=kh , k=0,,Nz1 , with h the axial slice spacing. The ideal point spread function (PSF) can be expressed by the Dirac delta function

Eq. 3

PSF(x,z)=δ(x)δ(z).
In practice, the sampling is not ideal and we assume that the PSF has a spatial extent in the z direction that is larger than the axial slice spacing h .

Fig. 1

(a) Nongated slice-sequence acquisition procedure. (b) Two 2D slice sequences before and (c) after synchronization.

054001_1_006505jbo1.jpg

The subsequent algorithms aim at finding the unknown sequence sk in order to retrieve the original volume I(x,z,t) from the measurements Im(x,zk,t) .

2.2.

Synchronization

The core of the synchronization procedure rests on the registration of slice-sequence pairs with respect to time. We seek solutions that, for a given time shift, maximize the similarity (in some metric to be defined) between two adjacent slices. This similarity hypothesis is reasonable if the axial sampling step h (the distance between two adjacent slices) is smaller than the PSF extent in z or that the object undergoes sufficiently smooth and homogeneous deformations. Indeed, while the axial resolution drops as the axial extent of the PSF increases, the similarity between two adjacent slices improves as both measurements contain information from the same physical region. For the same slice spacing, ideal sampling induces better axial resolution, to the detriment of the similarity hypothesis. While a rigorous investigation about all possible motions that may or may not be imaged using this technique is beyond the scope of this paper, we have heuristically determined that a unique and correct dynamic object can be recovered in the case of periodic, continuous, and homogeneous transforms even in the unfavorable case of ideal sampling. In Sec. 4(A), we present a simulation that supports this observation. Deformations that are nonhomogeneous with respect to the z axis may result in incorrect reconstructions when the axial slice spacing h is too large, that is, larger than the axial extent of the PSF. In practice, such cases may only be dealt with by considering a region of interest where the deformation is known to be homogeneous (see Sec. 4(A)) or by the use of external information (ECG, etc.).

3.

Algorithms

In this section we present the synchronization methodology. The following steps are involved in the data processing:

  • 1. period determination, data interpolation, and cropping,

  • 2. determination of relative shifts between pairs of slices,

  • 3. determination of the slices’ absolute shifts with respect to the first slice,

  • 4. synchronization and postprocessing.

3.1.

Period Determination

In order to ensure proper synchronization, the heartbeat period must be known precisely and be the same for all slice sequences at different depths. The image sequences are acquired at times τi=ihτ , i=0,,Nτ1 , where hτ is the acquisition sampling step and Nτ is the number of acquired frames. We achieved precise and automatic period determination by use of a technique inspired from astronomy.37, 38 For a given slice sequence and a candidate period T , the time positions of every pixel are brought back to the first period (phase locking)

Eq. 4

τi=τiτiTT
and a bijective mapping i=i(j) [respectively j=j(i) ] is defined such that τj1τjτj+1 . An estimate of the phase-locked signal’s dispersion is given by the length of the graph (τj,Im[x,zk,τi(j)])j=0,,Nτ1 that joins the newly ordered samples on a normalized time scale (see Fig. 2 ), cumulated over the whole image, i.e.,

Eq. 5

D(zk,T)=mZ2j=1Nτ1{Im[xm,zk,τi(j)]Im[xm,zk,τi(j1)]2 +τjτj12T2}12,
with xm=mhxy . Here, for simplicity, we consider that the sampling step in the 0xy plane, hxy , is 1. The correct period T(zk) is found by minimizing the above expression, viz.,

Eq. 6

T(zk)=argminTD(zk,T).
Starting from an initial guess of the range T(zk)[Tmin,Tmax] , we solve Eq. 6 iteratively using a combined parabolic and golden section search algorithm,39, 40, 41 which usually converges to a subsampling-step accuracy in less than 10 iterations.

Fig. 2

Dispersion measure for period estimation. (a) Sampled intensity variation over time (with unknown period T ) at one location x . (b) For a candidate period T1 the samples are brought back to the first period (phase locking). The dispersion of the samples is given by the length of the curve that joins the newly ordered samples (bold curve). (c) When the candidate period corresponds to the actual period, i.e., T2=T , the dispersion of the samples, hence the length of this curve, is minimized.

054001_1_006505jbo2.jpg

In order for the periodic boundary conditions in the time direction to be applied during subsequent operations, we crop and resample the data to cover an integer number of periods. We used linear spline interpolation, which offers a fair compromise between the accuracy of higher order interpolation schemes and the time efficiency of nearest neighbor interpolation.42, 43 The samples are taken at times ti=iht , i=0,,Nt1 , with sampling step ht=LNt , where Nt is the number of considered frames over the total time L=NTT and NT the number of considered periods. This also allows for temporal stretching or compression in cases where the periods at different depths are not the same. From this point onward, we consider that the measured signal Im(x,zk,t) is known for xR2 and t[0,L) (possibly via the interpolation of samples that are uniformly distributed over that domain) and that periodic boundary conditions in time apply.

3.2.

Determination of Relative Shifts

Our automatic synchronization algorithm is based on the minimization of an objective criterion to measure the similarity between the data from neighboring depths zk and zk . We have chosen a least-squares criterion that has been shown to be effective for registration algorithms44

Eq. 7

Qk,k(s)=R20LIm(x,zk,t)Im(x,zk,ts)2dtdx
where the shift sR . We can further write

Eq. 8

Qk,k(s)=R20LIm(x,zk,t)2+Im(x,zk,ts)2dtdx2R20LIm(x,zk,t)Im(x,zk,ts)dtdx =C2R20LIm(x,zk,t)Im(x,zk,ts)dtdx,
where the integral of the second quadratic term does not depend on s because of the periodicity with respect to time. Since the above expression has the form of a correlation and periodic boundary conditions apply, we can compute Qk,k(s) (up to the constant C ) for a number of regularly spaced shifts s=iht , with ht=LNt , i=0,,Nt1 , at a limited cost using the fast Fourier transform (FFT). The relative shifts sk,k between any two pairs of z slices are obtained by finding the shifts s that minimizes Qk,k(s) . They may be represented by the antisymmetric matrix S , whose elements are

Eq. 9

sk,k=argmins=kT,k=1,,NtQk,k(s).
Note that this matrix not only includes slice-sequence pairs that are immediate neighbors but also pairs that lie farther apart. We also compute the correlation for such slices in order to reduce synchronization errors that may quickly propagate due to the sequential alignment and the limited discrete steps the FFT imposes.

Before we derive the method for the determination of the shifts relative to the first slice sequence (absolute shifts), we refine the above correlation technique to make it time and memory effective as well as robust. Indeed, computing Eqs. 8, 5 naively would require considerable time and memory resources as the multidimensional data rapidly exceeds the storage capacity of even the latest available desktop computers. Another concern that complicates the equations’ direct implementation is that the images are corrupted by noise. As a consequence, the objective functions are as well. Yet another caveat is the presence of features that are characteristic of the studied structure but do not comply with the similarity hypothesis. For example, red blood cells are confined to the inside of the heart tube and have a movement that is in synchrony with the heart movement, however, the individual cells do not occupy the same positions from slice to slice. The correct extremum determination is thus severely affected. We have chosen to take advantage of the multiresolution45 and noise decorrelation46, 47 properties that the wavelet decomposition offers to solve these issues.

We consider a separable orthogonal wavelet basis of L2(R2) ,

Eq. 10

{ψj,m1(x),ψj,m 2(x),ψj,m 3(x)}jZ,mZ2
where the two-dimensional wavelets

Eq. 11

ψj,m p(x)=12jψ p(x2jm)
are constructed with separable products of the 1D scaling function ϕ(x) and wavelet ψ(x)

Eq. 12

ψ 1(x)=ϕ(x)ψ(y),ψ 2(x)=ψ(x)ϕ(y),ψ 3(x)=ψ(x)ψ(y).
For the sake of brevity, we index the basis functions with a single index k that includes the scale jZ , translation mZ2 , and wavelet type p{1,2,3} :

Eq. 13

ψk(x)=ψj,m p(x),k=(p,j,m).
With this notation, we may expand an image at a fixed depth zk and time point t in the wavelet basis as

Eq. 14

Im(x,zk,t)=kck,k(t)ψk(x)
where the coefficients are given by the inner products (wavelet transform)

Eq. 15

ck,k(t) =Im(,zk,t),ψk

Eq. 16

=Im(x,zk,t)ψk(x)dx.

Since the basis functions are orthogonal, i.e., ψk,ψk=δk,k , we may rewrite Eq. 8 as

Eq. 17

Qk,k(jht)=C20Lkck,k(t)ck,k(tjht)dtC2htki=0Nt1ck,k(iht)ck,k[(ij)ht].

In practice, we only consider a finite number of scales and translations for k (because of the finite resolution and support of the image, and appropriate boundary conditions). Furthermore, we discard the fine resolution coefficients thus downsizing the data’s complexity to a tractable size (see Fig. 3 ). Since wavelet transforms induce concise signal representations we make sure that the most important information is still present. Also, at coarse scales, individual blood cells are not resolved. Since they are confined to the inside of the heart tube, their global position contributes to a useful correlation signal. However, since confocal images are subject to bleaching (whose consequence is the presence of a nonuniform background), we discard the low-pass coefficients (that contain most of the background energy) as well. We then apply a soft threshold to the remaining coefficients to limit the influence of other noise sources.

Fig. 3

Synchronization based on wavelet coefficients. (a) Original slice sequence at a given depth. (b) A 2D wavelet transform is applied to each frame. (c) Fine-scale wavelet coefficients are discarded (data reduction) as well as the low-pass coefficients at the coarsest scale. A threshold is applied to the remaining coefficients to increase robustness to noise. The reduced data is interpolated to increase synchronization accuracy.

054001_1_006505jbo3.jpg

Similarly, we may apply Eq. 5 to the reduced data set of wavelet coefficients instead of the sampled image pixels, i.e.,

Eq. 18

D̃(zk,T)=kj=1Nτ1{ck,k[τi(j)]ck,k[τi(j1)]2 +τjτj12T2}12,
thus gaining robustness, reducing the required memory, and decreasing the computation time. Although Eqs. 5, 18 are not formally equivalent, the latter may be compared to applying the former to a sequence of images whose main features (edges) have been enhanced. Indeed, the wavelet transform essentially acts as an oriented differential operator at multiple scales.

We did not notice significant differences in the overall behavior of the algorithm depending on the choice of the wavelet basis, which must, however, be orthogonal to ensure validity of Eq. 17. We chose to work with the Daubechies 97 wavelets.48 Although they are not orthogonal (but nearly49), they have good approximation properties and are symmetric. The latter property allows the implementation of an algorithm which does not require that the image size be a multiple of a power of two50 and which is thus well suited for region-of-interest processing.

Finally, to increase the synchronization accuracy, we linearly interpolate the processed wavelet coefficients with respect to time in order to obtain a finer time step when computing Eq. 17. This interpolation is fast since the amount of data is reduced.

3.3.

Absolute Shifts Determination

To determine the slice-sequence shifts with respect to the uppermost sequence (absolute shifts) sk , we consider their relation to the relative shifts sk,k :

19.

s1=0
sksk=sk,k,withk,k=1,,Nzandk<k.
Since slice-sequence pairs that are separated by a larger depth are less trustworthy, we assign different weights wkk to equations that involve the estimated shifts sk,k depending on the distance j=kk . We set lower weights wj to equations for slice pairs less likely to exhibit similarities, that is, when the distance kk between them increases. For a system with Nz=5 and wj=0 for j> 2 , we can rewrite Eq. 19 in matrix form

Eq. 20

054001_1_006505jbod1.jpg
along with the diagonal weighting matrix

Eq. 21

W=diag(1,w1,w1,w1,w1,w2,w2,w2).
We determine the weighted least-squares solution of Eq. 20, which is equivalent to solving

Eq. 22

AWWAt=AWWs
where () denotes transposition. Equation 20 may easily be modified to include supplementary information (not image-intrinsic) that may become available in the future, such as ECGs. Depending on the accuracy of the signals, we may then set appropriate weights in Eq. 21.

3.4.

Synchronization and Postprocessing

The original slice sequences are finally circularly shifted by the computed absolute shifts (using linear interpolation and resampling). The synchronized data may then be visualized using 4D-capable software packages.51 Noise reduction steps may be applied. We made use of a rolling-ball background removal algorithm (see Refs. 52, 53) to normalize the background. The 4D data series may also be analyzed to follow individual cell movements. The higher dimensionality of the data should also make it possible to take advantage of more sophisticated noise removal algorithms that have proven to be effective for other high-dimensional modalities.54, 55 Finally, the synchronized data might be suitable for subsequent deconvolution.

4.

Results and Discussion

4.1.

Simulation

We validate our approach by simulating the acquisition procedure on a periodically deformed test body. We have considered the following, much simplified, heart-tube phantom. At time t=0 , the contributing intensity at every location (x,z) is given by

Eq. 23

I(x,z,0)=I0β3((x2+y2)r0(z)w4)[1+γcos(2πfdx)cos(2πfdy)cos(2πfdz)]
where the central wall radius is given by

Eq. 24

r0(t,z)=R0+ΔRsin(2παz)
and where w is the wall thickness, α controls the tube’s geometry, R0 is the average tube radius, ΔR is the radius movement amplitude, γ is the amplitude of a regular pattern of frequency fd , and the cubic B-spline56 is given by

Eq. 25

β3(x)={23x2+x32,0x<1(2x)36,1x<20,2x.}
Typical heart motions include rotation, expansion, contraction, and shear.57 We model the intensity at subsequent times by a general periodic affine transformation of the coordinate system corresponding to an homogeneous deformation of the original body. The intensity at position (x,y,z) and time t is given by

Eq. 26

I(x,z,t)=I(x,z,0)
where, using homogeneous coordinates,58

Eq. 27

(xyz1)=(a11(t)a12(t)a13(t)a14(t)a21(t)a22(t)a23(t)a24(t)a31(t)a32(t)a33(t)a34(t)0001)(xyz1).
The time-periodic affine transformation matrix A(t) can be decomposed as a combination of translation, rotation, scaling, and shear:

Eq. 28

A(t)=T(t)R(t)H(t)C(t)
where T(t) , R(t) , H(t) , and C(t) are the matrices corresponding to the respective transformations. Twelve parameters control the deformation matrix and each of them is a periodic function of time, which we specify through the coefficients of its Fourier series (see Appendix A). The latter may be chosen randomly to cover the full range of possible transformations. In Fig. 4 , we show several time points of such a random, periodic, and continuous deformation cycle that includes shear, rotation, translation, and scaling.

Fig. 4

Five time points of a simulated heart-tube deformation cycle (homogeneous transform). The corresponding movie is available at Ref. 62.

054001_1_006505jbo4.jpg

To assess the performance of our method, we generated a set of 100 deformation cycles using at each time different (normally distributed) random variables for the second and third harmonics of each parameter function, as well as random shifts s¯kUn(T,T) (uniform probability distribution). We considered the simplified PSF of Eq. 3 with Nz=20 , Nt=40 , hr=1 (normalized time units), and a period T=19.5 . From these simulated measurements, we then applied our algorithm (using 80 time points to compute the correlations, that is, after the cropping step, approximately 2× oversampling) to retrieve the shifts sk . Since the true absolute shifts were known, we could compute the absolute error using the following formula (which takes into account periodicity, e.g., comparing shifts s1=δ and s2=Tδ yields an error ϵ=2δ ):

Eq. 29

ϵ=min(WT(s¯k)WT(sk),TWT(s¯k)WT(sk))
where WT(x)=xxTT . The mean error over the 100 experiments was ϵ¯=0.31±0.08 frames. This result confirms that for the vast class of periodic homogeneous transforms our method is highly reliable, even when the considered sampling is ideal, i.e., when there is no axial overlap of the PSF. The error may be reduced by linearly interpolating the wavelet coefficients at a finer sampling rate in time. For different oversampling rates, we obtained the following errors: ϵ¯1×=0.41±0.12 , ϵ¯2×=ϵ¯=0.31±0.08 , ϵ¯4×=0.27±0.06 , and ϵ¯8×=0.25±0.06 . However, visual inspection of the reconstructions from in vivo measurements showed no significant improvement of the results above 2× oversampling. Also, although the accuracy of correlation-based registration methods is known to be inherently limited,59 in practice, the current limiting factors are the irregularities in the heartbeat periodicity of the biological samples themselves.

4.2.

Experimental Measurements

With the aim of a better understanding of the zebrafish cardiac development, we applied our method to slice sequences from an early embryonic, 48hours post fertilization (h.p.f.), beating heart. In Fig. 5 , we show a bright-field microscopy image of a 48h.p.f. zebrafish embryo where the heart position has been indicated. The study focuses on zebrafish for several reasons: they are vertebrates that reproduce externally and rapidly, they are relatively transparent, and it is possible to genetically engineer fish strains that express vital fluorescent markers in specific tissues (for instance, heart wall, or blood cells). Here, we have chosen to study Tg(gata1∷GFP) zebrafish embryos whose endo- and myocardial cells as well as erythrocytes are fluorescent.60 The embryos were anesthetized in order to limit the imaged movements to those of the heart. Images were acquired using a Zeiss LSM 5 LIVE laser scanning microscope prototype2 at a frame rate of 151Hz for the duration of three to four heartbeats. The images had 256×256pixels and a sampling step of 0.9μmperpixel ( 40× AchroPlan water-immersion lens NA=0.8 ). The stage was then moved axially in increments of 5μm before a new sequence was acquired. A total of about 20 positions could be imaged per embryo. A complete description of the experimental aspects of the measurement procedure is reported elsewhere.61

Fig. 5

Bright-field images of a zebrafish embryo at 48h.p.f. ; h=heart ; e=eye ; mb=midbrain ; ot=otocyst ; yolk=yolk mass; A=atrium ; V=ventricle .

054001_1_006505jbo5.jpg

The heartbeat of the studied zebrafish appeared to remain steady over the usual acquisition time for one slice (three to four heartbeats). However, we observed changes in the rate of up to about 10% between the sequence at the first and last z position. Once identified, these variations—mainly due to ambient temperature changes—could subsequently be controlled to limit the period change. We considered three periods per slice sequence. In Fig. 6 , we show the experimental correlation curve for one slice-sequence pair. The curve’s three main maxima correspond to admissible periodic shifts (one peak shift per imaged heartbeat).

Fig. 6

Experimental correlation curve [see Eq. 8] for two adjacent slice sequences. The maximum correlation occurs for a shift of 259 frames.

054001_1_006505jbo6.jpg

In Fig. 7 , we show the rendering over five frames of two reconstructed embryo hearts. Mostly erythrocytes, but also endo- and myocardial cells are fluorescent and visible. The orientation with respect to the z axis is different for the two samples, yet the reconstructions show similar features, which supports the hypothesis that the method is suitable for accurate imaging of the wall deformations. These reconstructed images allow the visualization of complex flow and wall movement patterns that previously could not be studied.

Fig. 7

Reconstructed 4D datasets of beating 48h.p.f. Tg(gata1∷GFP) zebrafish hearts. In the first frame, the red blood cells fill the atrium. In the following frames, the red blood cells are pumped through the atrial ventricular canal into the ventricle. The atrium fills again while the cycle is completed. The grid spacing is 20μm . Movies are available at Ref. 62.

054001_1_006505jbo7.jpg

The computation time on a 2-GHz PowerPC G5, for a set of 20 slice sequences of size 256×256pixels and 220 time frames, is distributed as follows: preprocessing (wavelet transform): 1min ; period retrieval: 10s ; time interpolation, resampling of wavelet coefficients, and FFT: 7s ; shift determination: 7s (absolute shift determination takes less than 0.01s ); original data shifting, interpolation, and sampling: 40s . Finally, our implementation’s memory requirements (RAM) are below 512MB for the above dataset.

5.

Conclusion

We have presented a procedure for the synchronization of nongated confocal slice sequences to build dynamic 3D volumes. We have investigated the ability of our method to achieve this goal and found that it performs well. We have validated the approach both through simulation and in vivo measurements. The described algorithms appear to be robust and lead to coherent results. Provisions are made in the method for the subsequent inclusion of a priori data to relieve the current requirements on the movements that can be studied with this technique.

Appendices

Appendix: Periodic Affine Transformations

The transformation matrix in Eq. 28 can be decomposed using the following matrices for scaling

Eq. A1

C(t)=(Sx(t)0000Sy(t)0000Sz(t)00001),
translation

Eq. A2

T(t)=(100Tx(t)010Ty(t)001Tz(t)0001),
shear

Eq. A3

H(t)=(1Sxy(t)Sxz(t)001Syz(t)000100001),
and rotation

Eq. A4

R(t)=(cos[ψ(t)]sin[ψ(t)]00sin[ψ(t)]sin[ψ(t)]0000100001)

Eq. A5

(10000cos[ϑ(t)]sin[ϑ(t)]00sin[ϑ(t)]sin[ϑ(t)]00001)

Eq. A6

(cos[ϕ(t)]sin[ϕ(t)]00sin[ϕ(t)]sin[ϕ(t)]0000100001),
where the twelve coefficients

Eq. A7

θ(t)=[θ1(t),,θ12(t)] =[ϕ(t),ϑ(t),ψ(t),Sx(t),Sy(t),Sz(t),Sxy(t),Sxz(t),Syz(t),Tx(t),Ty(t),Tz(t)]
are periodic functions that can conveniently be expressed by their Fourier series

Eq. A8

θi(t)=a0i+k=1akicos[2π(kT)t]+k=1bkisin[2π(kT)t] i=1,,12.

Acknowledgments

We thank S. Lin for the transgenic Tg(gata1∷GFP) zebrafish line. We thank Mike Tyszka for his comments on the manuscript. M. Liebling is supported by a fellowship from the Swiss National Science Foundation, PBEL2-104418.

References

1. 

M. Minsky, “Memoir on inventing the confocal scanning microscope,” Scanning, 10 (4), 128 –138 (1988). 0161-0457 Google Scholar

2. 

Carl Zeiss AG, “LSM 5 LIVE,” http://www.zeiss.de/ Google Scholar

3. 

J. R. Hove, R. W. Köster, A. S. Forouhar, G. Acevedo-Bolton, S. E. Fraser, and M. Gharib, “Intracardiac fluid forces are an essential epigenetic factor for embryonic cardiogenesis,” Nature (London), 421 172 –177 (2003). https://doi.org/10.1038/nature01282 0028-0836 Google Scholar

4. 

E. A. V. Jones, M. H. Baron, S. E. Fraser, and M. E. Dickinson, “Measuring hemodynamic changes during mammalian development,” Am. J. Physiol. Heart Circ. Physiol., 287 (4), H1561 –H1569 (2004). 0363-6135 Google Scholar

5. 

R. Jerecic, M. Bock, S. Nielles-Vallespin, C. Wacker, W. Bauer, and L. R. Schad, “ECG-gated Na-23-MRI of the human heart using a 3D-radial projection technique with ultra-short echo times,” Magn. Reson. Mater. Phys., Biol., Med., 16 (6), 297 –302 (2004). 1352-8661 Google Scholar

6. 

M. Markl, F. P. Chan, M. T. Alley, K. L. Wedding, M. T. Draney, C. J. Elkins, D. W. Parker, R. Wicker, C. A. Taylor, R. J. Herfkens, and N. J. Pelc, “Time-resolved three-dimensional phase-contrast MRI,” J. Magn. Reson Imaging, 17 (4), 499 –506 (2003). https://doi.org/10.1067/mva.1993.38251 1053-1807 Google Scholar

7. 

S. Skare and J. L. R. Andersson, “On the effects of gating in diffusion imaging of the brain using single shot EPI,” Magn. Reson. Imaging, 19 (8), 1125 –1128 (2001). 0730-725X Google Scholar

8. 

C. Dornier, M. K. Ivancevic, P. Thévenaz, and J. P. Vallée, “Improvement in the quantification of myocardial perfusion using an automatic spline-based registration algorithm,” J. Magn. Reson Imaging, 18 (2), 160 –168 (2003). 1053-1807 Google Scholar

9. 

Z. P. Liang, H. Jiang, C. P. Hess, and P. C. Lauterbur, “Dynamic imaging by model estimation,” Int. J. Imaging Syst. Technol., 8 (6), 551 –557 (1997). https://doi.org/10.1002/(SICI)1098-1098(1997)8:6<551::AID-IMA7>3.0.CO;2-9 0899-9457 Google Scholar

10. 

T. Pan, T. Y. Lee, E. Rietzel, and G. T. Y. Chen, “4D-CT imaging of a volume influenced by respiratory motion on multi-slice CT,” Med. Phys., 31 (2), 333 –340 (2004). https://doi.org/10.1118/1.1639993 0094-2405 Google Scholar

11. 

C. Le Rest, O. Couturier, A. Turzo, and Y. Bizais, “Post-synchronization of dynamic images of periodically moving organs,” Nucl. Med. Commun., 21 (7), 677 –684 (2000). https://doi.org/10.1097/00006231-200007000-00012 0143-3636 Google Scholar

12. 

R. B. Thompson and E. R. McVeigh, “Flow-gated phase-contrast MRI using radial acquisitions,” Magn. Reson. Med., 52 (3), 598 –604 (2004). https://doi.org/10.1002/mrm.20187 0740-3194 Google Scholar

13. 

M. E. Crowe, A. C. Larson, Q. Zhang, J. Carr, R. D. White, D. B. Li, and O. P. Simonetti, “Automated rectilinear self-gated cardiac cine imaging,” Magn. Reson. Med., 52 (4), 782 –788 (2004). https://doi.org/10.1002/mrm.20212 0740-3194 Google Scholar

14. 

A. C. Larson, R. D. White, G. Laub, E. R. McVeigh, D. B. Li, and O. P. Simonetti, “Self-gated cardiac cine MRI,” Magn. Reson. Med., 51 (1), 93 –102 (2004). https://doi.org/10.1002/mrm.10664 0740-3194 Google Scholar

15. 

T. A. Spraggins, “Wireless retrospective gating—application to cine cardiac imaging,” Magn. Reson. Imaging, 8 (6), 675 –681 (1990). 0730-725X Google Scholar

16. 

R. D. White, C. B. Paschal, M. E. Clampitt, T. A. Spraggins, and G. W. Lenz, “Electrocardiograph-independent, wireless cardiovascular cine MR imaging,” J. Magn. Reson Imaging, 1 (3), 347 –355 (1991). 1053-1807 Google Scholar

17. 

M. Grass, R. Manzke, T. Nielsen, P. Koken, R. Proksa, M. Natanzon, and G. Shechter, “Helical cardiac cone beam reconstruction using retrospective ECG gating,” Phys. Med. Biol., 48 (18), 3069 –3084 (2003). https://doi.org/10.1088/0031-9155/48/18/308 0031-9155 Google Scholar

18. 

M. Kachelrieß, D. A. Sennst, W. Maxlmoser, and W. A. Kalender, “Kymogram detection and kymogram-correlated image reconstruction from subsecond spiral computed tomography scans of the heart,” Med. Phys., 29 (7), 1489 –1503 (2002). https://doi.org/10.1118/1.1487861 0094-2405 Google Scholar

19. 

S. A. de Winter, R. Hamers, M. Degertekin, K. Tanabe, P. A. Lemos, P. W. Serruys, J. R. T. C. Roelandt, and N. Bruining, “Retrospective image-based gating of intracoronary ultrasound images for improved quantitative analysis: the intelligate method,” Catheterization Cardiovascular Interventions, 61 (1), 84 –94 (2004). Google Scholar

20. 

G. M. Treece, R. W. Prager, A. H. Gee, C. J. C. Cash, and L. Berman, “Grey-scale gating for freehand 3D ultrasound,” 993 –996 (2002). Google Scholar

21. 

Q. Tian and M. N. Huhns, “Algorithms for subpixel registration,” Comput. Vis. Graph. Image Process., 35 (2), 220 –233 (1986). 0734-189X Google Scholar

22. 

L. G. Brown, “A survey of image registration techniques,” ACM Comput. Surv., 24 (4), 325 –376 (1992). https://doi.org/10.1145/146370.146374 0360-0300 Google Scholar

23. 

C. Stiller and J. Konrad, “Estimating motion in image sequences—a tutorial on modeling and computation of 2D motion,” IEEE Signal Process. Mag., 16 (4), 70 –91 (1999). https://doi.org/10.1109/79.774934 1053-5888 Google Scholar

24. 

J. L. Barron, D. J. Fleet, and S. S. Beauchemin, “Performance of optical-flow techniques,” Int. J. Comput. Vis., 12 (1), 43 –77 (1994). https://doi.org/10.1007/BF01420984 0920-5691 Google Scholar

25. 

J. Cox and C. J. R. Sheppard, “Digital image processing of confocal images,” Image Vis. Comput., 1 (1), 52 –56 (1983). https://doi.org/10.1016/0262-8856(83)90008-2 0262-8856 Google Scholar

26. 

A. Diaspro, “Introduction: image processing in optical microscopy,” Microsc. Res. Tech., 64 (2), 87 –88 (2004). https://doi.org/10.1002/jemt.20079 1059-910X Google Scholar

27. 

Y. L. Sun, B. Rajwa, and J. P. Robinson, “Adaptive image-processing technique and effective visualization of confocal microscopy images,” Microsc. Res. Tech., 64 (2), 156 –163 (2004). https://doi.org/10.1002/jemt.20064 1059-910X Google Scholar

28. 

M. Kozubek, P. Matula, P. Matula, and S. Kozubek, “Automated acquisition and processing of multidimensional image data in confocal in vivo microscopy,” Microsc. Res. Tech., 64 (2), 164 –175 (2004). https://doi.org/10.1002/jemt.20068 1059-910X Google Scholar

29. 

R. L. Allen, F. A. Kamangar, and E. M. Stokely, “Laplacian and orthogonal wavelet pyramid decompositions in coarse-to-fine registration,” IEEE Trans. Signal Process., 41 (12), 3536 –3541 (1993). https://doi.org/10.1109/78.258092 1053-587X Google Scholar

30. 

M. J. Carroll and K. E. Britton, “Wavelet based image registration and its applications in nuclear medicine,” Eur. J. Nucl. Med., 26 (9), 1008 –1008 (1999). 0340-6997 Google Scholar

31. 

J. Deubler and J. C. Olivo, “A wavelet-based multiresolution method to automatically register images,” J. Math. Imaging Vision, 7 (3), 199 –209 (1997). 0924-9907 Google Scholar

32. 

I. D. Dinov, M. S. Mega, T. Poon, P. Thompson, and A. W. Toga, “Volumetric registration of brain data: a new technique for unbiased assessment and comparison of different registration methods using efficient wavelet representation,” Neurology, 58 (7), A141 –A141 (2002). 0028-3878 Google Scholar

33. 

S. Gefen, O. Tretiak, L. Bertrand, G. D. Rosen, and J. Nissanov, “Surface alignment of an elastic body using a multiresolution wavelet representation,” IEEE Trans. Biomed. Eng., 51 (7), 1230 –1241 (2004). https://doi.org/10.1109/TBME.2004.827258 0018-9294 Google Scholar

34. 

J. Le Moigne, W. J. Campbell, and R. F. Cromp, “An automated parallel image registration technique based on the correlation of wavelet features,” IEEE Trans. Geosci. Remote Sens., GE–40 (8), 1849 –1864 (2002). 0196-2892 Google Scholar

35. 

R. Sharman, J. M. Tyler, and O. S. Pianykh, “A fast and accurate method to register medical images using wavelet modulus maxima,” Pattern Recogn. Lett., 21 (6–7), 447 –462 (2000). https://doi.org/10.1016/S0167-8655(00)00002-7 0167-8655 Google Scholar

36. 

M. Unser, P. Thévenaz, C. H. Lee, and U. E. Ruttimann, “Registration and statistical-analysis of PET images using the wavelet transform,” IEEE Eng. Med. Biol. Mag., 14 (5), 603 –611 (1995). https://doi.org/10.1109/51.464777 0739-5175 Google Scholar

37. 

M. M. Dworetsky, “A period-finding method for sparse randomly spaced observations or how long is a piece of string,” Mon. Not. R. Astron. Soc., 203 (3), 917 –924 (1983). 0035-8711 Google Scholar

38. 

R. F. Stellingwerf, “Period determination using phase dispersion minimization,” Astrophys. J., 224 (3), 953 –960 (1978). https://doi.org/10.1086/156444 0004-637X Google Scholar

39. 

Matlab Function Reference, 2 The MathWorks, Inc., Natick, MA (2002). Google Scholar

40. 

R. P. Brent, Algorithms for Minimization without Derivatives, Prentice Hall, Englewood Cliffs, NJ (1973). Google Scholar

41. 

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed.Cambridge University Press, Cambridge (1992). Google Scholar

42. 

E. H. W. Meijering, W. J. Niessen, and M. A. Viergever, “Quantitative evaluation of convolution-based methods for medical image interpolation,” Med. Image Anal, 5 (2), 111 –126 (2001). 1361-8415 Google Scholar

43. 

P. Thévenaz, T. Blu, and M. Unser, “Image interpolation and resampling,” Handbook of Medical Imaging, Processing and Analysis, 393 –420 Academic Press, San Diego, CA (2000). Google Scholar

44. 

P. Thévenaz, U. E. Ruttimann, and M. Unser, “A pyramid approach to subpixel registration based on intensity,” IEEE Trans. Image Process., 7 (1), 27 –41 (1998). https://doi.org/10.1109/83.650848 1057-7149 Google Scholar

45. 

S. G. Mallat, “A theory for multiresolution signal decomposition—the wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell., 11 (7), 674 –693 (1989). https://doi.org/10.1109/34.192463 0162-8828 Google Scholar

46. 

D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” Biometrika, 81 (3), 425 –455 (1994). 0006-3444 Google Scholar

47. 

J. B. Weaver, X. Yansun, D. M. Healy Jr., and L. D. Cromwell, “Filtering noise from images with wavelet transforms,” Magn. Reson. Med., 21 (2), 288 –295 (1991). 0740-3194 Google Scholar

48. 

M. Antonini, M. Barlaud, M. Mathieu, and I. Daubechies, “Image coding using wavelet transform,” IEEE Trans. Image Process., 1 (2), 205 –220 (1992). https://doi.org/10.1109/83.136597 1057-7149 Google Scholar

49. 

M. Unser and T. Blu, “Mathematical properties of the JPEG2000 wavelet filters,” IEEE Trans. Image Process., 12 (9), 1080 –1090 (2003). https://doi.org/10.1109/TIP.2003.812329 1057-7149 Google Scholar

50. 

C. M. Brislawn, “Classification of nonexpansive symmetric extension transforms for multirate filter banks,” Appl. Comput. Harmon. Anal., 3 (4), 337 –357 (1996). https://doi.org/10.1006/acha.1996.0026 1063-5203 Google Scholar

51. 

Bitplane AG, “Imaris,” http://www.bitplane.ch/ Google Scholar

52. 

S. R. Sternberg, “Biomedical image-processing,” Computer, 16 (1), 22 –34 (1983). 0018-9162 Google Scholar

53. 

W. Rasband, “ImageJ,” http://rsb.info.nih.gov/ij/ Google Scholar

54. 

E. D. Angelini, A. F. Laine, S. Takuma, J. W. Holmes, and S. Homma, “LV volume quantification via spatiotemporal analysis of real-time 3-D echocardiography,” IEEE Trans. Med. Imaging, 20 (6), 457 –469 (2001). https://doi.org/10.1109/42.929612 0278-0062 Google Scholar

55. 

R. M. Willett and R. D. Nowak, “Platelets: a multiscale approach for recovering edges and surfaces in photon-limited medical imaging,” IEEE Trans. Med. Imaging, 22 (3), 332 –350 (2003). https://doi.org/10.1109/TMI.2003.809622 0278-0062 Google Scholar

56. 

M. Unser, “Splines: a perfect fit for signal and image processing,” IEEE Signal Process. Mag., 16 (6), 22 –38 (1999). https://doi.org/10.1109/79.799930 1053-5888 Google Scholar

57. 

M. Sermesant, C. Forest, X. Pennec, H. Delingette, and N. Ayache, “Deformable biomechanical models: application to 4D cardiac image analysis,” Med. Image Anal, 7 (4), 475 –488 (2003). 1361-8415 Google Scholar

58. 

J. D. Foley and A. Van Dam, Fundamentals of Interactive Computer Graphics, Addison-Wesley, Reading, MA (1984). Google Scholar

59. 

D. Robinson and P. Milanfar, “Fundamental performance limits in image registration,” IEEE Trans. Image Process., 13 (9), 1185 –1199 (2004). https://doi.org/10.1109/TIP.2004.832923 1057-7149 Google Scholar

60. 

Q. M. Long, A. M. Meng, H. Wang, J. R. Jessen, M. J. Farrell, and S. Lin, “GATA-1 expression pattern can be recapitulated in living transgenic zebrafish using GFP reporter gene,” Development, 124 (20), 4105 –4111 (1997). 0950-1991 Google Scholar

61. 

A. S. Forouhar, M. Liebling, R. Wolleschensky, R. Ankerhold, B. Zimmermann, S. E. Fraser, M. Gharib, and M. E. Dickinson, “Rapid, confocal microscopy reveals functional changes in embryonic heart development,” Google Scholar
©(2005) Society of Photo-Optical Instrumentation Engineers (SPIE)
Michael Liebling, Arian S. Forouhar, Morteza Gharib, Scott E. Fraser, and Mary E. Dickinson "Four-dimensional cardiac imaging in living embryos via postacquisition synchronization of nongated slice sequences," Journal of Biomedical Optics 10(5), 054001 (1 September 2005). https://doi.org/10.1117/1.2061567
Published: 1 September 2005
Lens.org Logo
CITATIONS
Cited by 160 scholarly publications and 1 patent.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Heart

Wavelets

Cardiac imaging

Point spread functions

Confocal microscopy

Data acquisition

Microscopes

Back to Top