The temporal alignment of nongated slice-sequences acquired

at different axial positions in the living embryonic zebrafish

heart permits the reconstruction of dynamic, three-dimensional

data. This approach overcomes the current acquisition-

speed limitation of confocal microscopes for real-time

three-dimensional imaging of fast processes. Current synchronization

methods align and uniformly scale the data in

time, but do not compensate for slight variations in the heart

rhythm that occur within a heartbeat. Therefore, they impose

constraints on the admissible data quality. Here, we derive

a nonuniform registration procedure based on the minimization

of the absolute value of the intensity difference between

adjacent slice-sequence pairs. The method compensates for

temporal intra-sample variations and allows the processing of

a wider range of data to build functional, dynamic models of

the beating embryonic heart. We show reconstructions from

data acquired in living, fluorescent zebrafish embryos

Dickinson, Mary E.

Forouhar, Arian S.

Fraser, Scott E.

Gharib, Morteza

Liebling, Michael

Vermot, Julien

Caltech Authors - Main

NONUNIFORM TEMPORAL ALIGNMENT OF SLICE SEQUENCES FOR FOUR-
DIMENSIONAL IMAGING OF CYCLICALLY DEFORMING EMBRYONIC STRUCTURES
Michael Liebling1, Julien Vermot1, Arian S. Forouhar2,
Morteza Gharib2, Mary E. Dickinson3, Scott E. Fraser1
1Biological Imaging Center, Beckman Institute, California Institute of Technology, Pasadena, CA
2Option of Bioengineering, California Institute of Technology, Pasadena, CA
3Dept. of Molecular Physiology and Biophysics, Baylor College of Medicine, Houston, TX
ABSTRACT
The temporal alignment of nongated slice-sequences acquired
at different axial positions in the living embryonic zebrafish
heart permits the reconstruction of dynamic, three-dimen-
sional data. This approach overcomes the current acquisi-
tion-speed limitation of confocal microscopes for real-time
three-dimensional imaging of fast processes. Current syn-
chronization methods align and uniformly scale the data in
time, but do not compensate for slight variations in the heart
rhythm that occur within a heartbeat. Therefore, they impose
constraints on the admissible data quality. Here, we derive
a nonuniform registration procedure based on the minimiza-
tion of the absolute value of the intensity difference between
adjacent slice-sequence pairs. The method compensates for
temporal intra-sample variations and allows the processing of
a wider range of data to build functional, dynamic models of
the beating embryonic heart. We show reconstructions from
data acquired in living, fluorescent zebrafish embryos.
1. INTRODUCTION
With speeds that reach several millimeters per second, heart-
wall motions and blood flow in the embryonic heart require
fast frame-rates to be imaged under the microscope without
inducing motion artifacts such as blurring or aliasing [1]. Con-
focal microscopes enable selective, three-dimensional imag-
ing of fluorescent cells at depths up to 100-200 µm. Re-
trieving the temporal evolution of fluorescently labeled cells
or blood flow gives invaluable insights for understanding the
driving forces behind the development of the cardio-vascular
system in early embryos.
Although recently developed confocal microscopes can
capture images of two-dimensional slices at frame-rates as
fast as 120 frames per s for 512×512 pixel images, direct
recording of dynamic, three-dimensional data at similar frame-
rates is currently not possible. To overcome this limitation,
Corresponding author is M. Liebling, California Institute of Technology,
Mail Code 139-74, Pasadena, CA 91125, USA, email: liebling@caltech.edu.
z z
(a)
(b) (c) (d)
1 2
t
Fig. 1. Four-dimensional imaging: acquisition and re-
construction procedure. (a) Sequential acquisition of two-
dimensional slices (confocal microscopy) as a time-series at
increasing depths from a cyclically deformed object (cone).
(b) Direct reconstruction is not possible from the nongated
data. (c) Temporal alignment procedure (d) Reconstruction.
an acquisition and reconstruction procedure based on the se-
quential, nongated acquisition of slice-sequences during sev-
eral heart-beats and subsequent synchronization (see Fig. 1)
was recently proposed [2]. In that work, the registration re-
sulted from scaled, temporal shifting of the slice sequences.
This technique is efficient if, during a single two-dimensional
temporal sequence acquisition, the motion is periodic. The
reconstruction precision and quality eventually depend on the
conformity of the measurements to this hypothesis. From our
observations, we noticed that slight irregularities (variation in
the heart beat) of 1% can induce severe visual artifacts that
prevent us from building accurate models and also force us to
discard measurements that do not comply with the regularity
hypothesis.
Here, we propose a nonuniform time-registration method
that recursively aligns pairs of slice-sequences to build a four-
dimensional reconstruction of the beating heart. The method
is based on the minimization of the absolute intensity differ-
ence of the warped sequences’ two-dimensional spatial wave-
let coefficients, an optimization problem that we solve in a
dynamic programming framework.
11560-7803-9577-8/06/$20.00 ©2006 IEEE ISBI 2006
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 
While, to the best of our knowledge, the use of dynamic
programming to solve our image reconstruction problem is
new, similar registration techniques are widespread for differ-
ent purposes, in particular, for speech analysis [3, 4].
In Section 2 we describe the proposed acquisition and re-
construction method. In Section 3 we present reconstructions
that we obtained from in vivo imaging of fluorescent zebrafish
embryos. In Section 4 we discuss our results and future re-
search directions.
2. METHOD
2.1. Acquisition Model and Cost Function Definition
We model the measured intensity Im as follows,
Im(x, zk, t) =
∫∫∫
I[x′, z, τk(t)] h(x − x′, zk − z) dx′dz,
(1)
where τk(t) is an unknown warping function at depth zk =
khz, k = 0, . . . , Nz. The warping functions τk : R+ →
R+, t → τk(t) are strictly increasing, i.e. if t1 < t2, we have
τk(t1) < τk(t2). The optical system’s point spread function
is denoted h(x, z). We assume that the four-dimensional in-
tensity function I(x, z, τ), has a period T , i.e.,
|I[x, z, τ ] − I[x, z, τ + T ]|  Imax. (2)
We assume that at depth zk̄, the corresponding warping func-
tion is identity, i.e. τk̄(t) = t.
In order to recover an estimate of the volume I(x, z, τ),
τ ∈ [0, L], L > T from the measurements Im(x, zk, t),
where t ∈ [tk, tk + L
′] with L′ > σL (σ > 1 is the max-
imal stretching factor), we consider the following objective
criterion to measure the discrepancy between the warped data
from neighboring depths zk and zk′ ,
Qk,k′,wk{w
′
k} =
∫ L
0
∫∫
R2
∣∣Im[x, zk, wk(τ)]
− Im[x, zk′ , wk′(τ)]
∣∣ dx dτ. (3)
We aim at recovering the warping functions wk′ (τ) for k
′ =
{0, . . . , Nz − 1}\{k̄}. We start from the depth zk̄ for which
wk̄ is known (identity) and find the functions wk′ recursively
by aligning neighboring pairs, i.e. (k, k′ = k − 1) for pairs
that are below k̄ and (k, k′ = k + 1) for pairs above the ref-
erence depth k̄. The reconstructed, dynamic volume is com-
puted by warping the measurements, viz.
IR(x, zk, τ) = Im(x, zk, wk(τ)). (4)
To find the appropriate functions wk in practice, we proceed
as follows. We consider discrete measurements (in space and
time): for the image sequence that is to be aligned, zk′ , we
consider measurements taken at times t = iht, with i =
0, . . . , Nt − 1 and, for the reference sequence at depth zk,
0 1 2 3
1
2
3
4
5
6
7
8
2 4 6 8 10 12 14 16
0
1
2
3
Test function
Reference
Fig. 2. Nonuniform temporal alignment via dynamic pro-
gramming. A matrix is generated by taking the absolute-
difference between samples of a test function (top) and a
reference function (right, triangles). Proceeding recursively
from left to right, a cost is associated to each position in the
matrix, based on the value of the matrix at that position, the
cost computed at the most favorable position in the neigh-
boring left column (yielding a vertex), and a possible penalty
if the pick induces stretching or compression of the time axis.
Next, starting from locations on the bottom line, the path with
lowest cost that reaches the top is kept, yielding the optimal
warping function (bold) to interpolate the test function (right,
circles).
measured at times τ = ihτ , with i = 0, . . . , N̄r − 1, we
consider interpolated samples at times τ = ihτ/s, with i =
0, . . . , Nr, Nr = (sN̄r) − 1 and s ∈ N
∗
+ an over-sampling
factor (typically, s = 4). The (discrete) problem boils down to
finding the appropriate path in a two-dimensional graph, that
pairs the test samples with the interpolated reference samples.
2.2. Dynamic Programming for Temporal Alignment
We denote Ir[] the discrete, interpolated reference slice se-
quence and It[k]) the test slice sequence. We compute the dif-
ference matrix ∆[, k] = ‖Ir[]−It[k]‖, with k = 0, . . . , Nt−
1 and for  = 0, . . . , Nr − 1, and where the norm ‖ · ‖ cor-
responds to a discretized version of Eq. (3). Exploring each
of the (Nr)
Nt possible paths across the matrix and finding
the one that minimizes Eq. (3) would be of prohibitive com-
plexity and possibly yield solutions that do not comply with
the monotonicity hypothesis. In order to limit the number of
paths to be considered, we use a dynamic programming [5]
approach, resulting in a complexity O(NrNt).
We define K(, ′) = α(|′−−w0|)
γ+Kmin, α, γ ∈ R∗+.
The limits for stretchability are given by wmin ≤ wmax < 0,
wmin, wmax ∈ Z
−. We recursively compute, for k = 0, . . . , Nt−
1 and for  = 0, . . . , Nr − 1, the following quantities:
L[, k] =
{
L[, k − 1] + 1 if  > 0 and k > 0
1 otherwise,
(5)
1157
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 
80 µm
D
P A
V
D
R L
V
R
P A
L
Eye
Heart
Tube
Heart
Tube
Heart
Tube
RBC
Yolk Sack
(a) (b)
*
*
Fig. 3. Orthogonal slice views of a 38 h.p.f. zebrafish heart reconstruction. (a) Before alignment, temporally misaligned slices
result in black streaks in the heart tube sections (two such slices are indicated by asterisks ∗). (b) After alignment, the heart
tube can be recognized. Red blood cells (RBC) appear as black dots inside the heart tube. Movies are available online [6].
g[, k] = min
′=+wmin,...,+wmax
g′(k, , ′), (6)
̄[, k] = arg min
′=+wmin,...,+wmax
g′(k, , ′), (7)
that is, the length of the path, a cost function, and the index of
the most favorable left neighbor, respectively, and where
g′(k, , ′)
=
⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩
(g[′, k − 1]L[′, k − 1]
+(1 − λ)∆[, k] + λK(, ′))
/(L[′, k − 1] + 1) if ′ ≥ 0 and k > 0
(1 − λ)∆[, k] + λKmin if 
′ ≥ 0 and k = 0
(1 − λ)∆[, k] + λKmin if 
′ < 0
and ′ −  = w0
∞ otherwise.
(8)
The weighting factor λ permits to adjust the warping stiffness.
From the resulting graph (see Fig.2), we determine the lead
of the optimal path by finding the vertex on the lower image
border with the lowest cost and that yields an acceptable path:
k0 = arg min
k′∈B
g[Nr − 1, k
′] (9)
with B = {k′|1 < k′ < Nt, and ̃k′ [0] ≤ 0} and where the
paths
{
̃k′ [k]
}
0≤k<Nt
are defined for 1 < k′ < Nt
̃k′ [k] =
⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩
Nr − 1 if k = k
′
̄
[
̃k′ [k + 1], k
′ + 1
]
if 0 ≤ k < k′
and ̃k′ [k + 1] ≥ 0
2̃k′ [k − 1] − ̃k′ [k − 2] if k
′ < k < Nt
2̃k′ [k + 1] − ̃k′ [k + 2] otherwise.
(10)
The optimal path is then given by
{
(k0 − i, ̃k0 [i])
}
0<i<Nt
unless B ≡ ∅. The warping function k̄(j), j = s, and
 = 0, ..., N̄r is obtained by inverting ̄k0(j) by use of lin-
ear interpolation.
Similarly to the approach taken in [2], we consider a lim-
ited set of wavelet coefficients instead of pixel values to com-
pute the norm (3). We therefore ensure that the method is
robust and that the computation of the difference matrix is
fast.
3. RESULTS
Wild-type zebrafish (danio rerio) eggs were spawned using
standard techniques [7] and the 38 hours post fertilization
(h.p.f) old embryos were soaked in a green fluorescent dye
(BODIPY FL C5-ceramide, Molecular Probes) for 3 hours.
We acquired two-dimensional image sequences at 60 frames
per s for 2 s, corresponding to 3–4 heartbeats, using a fast
slit-scanning confocal microscope (Zeiss LSM 5 LIVE, Carl
Zeiss Jena GmbH, Germany) with 488nm excitation light and
a 500–525nm band-pass filter to filter emission signal. Be-
tween the acquisition of two sequences, the stage was moved
axially by 10 µm or 5 µm, depending on the magnification
of the microscope objective (10 × or 40 ×, respectively). Per
embryo, we could image a total of 20–30 depths. The im-
ages of 512×512 pixels had a sampling step of 1.78 µm ×
1.78µm per pixel (Zeiss Plan-Neofluar 10×/0.3 microscope
objective), or 0.62 µm× 0.62 µm per pixel (Zeiss C-Apochro-
mat 40×/1.2 W microscope objective). The images were re-
aligned by use of a MATLAB implementation of the method
described above. The four-dimensional data was analyzed
and visualized using appropriate software (Imaris, Biplane
1158
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 
0 ms 102 ms 204 ms
408 ms306 ms 510 ms
100 µm
Eye Heart
Yolk Sack
Fig. 4. Four-dimensional reconstruction of a 38 h.p.f living zebrafish. The heart-tube is visible in the center, pumping blood.
The corresponding movie is available online [6].
AG, Zurich, Switzerland). The reconstructions prior and af-
ter alignment are presented in Fig. 3 (a) and (b), respectively.
The reconstructed data can be rendered (see Fig. 4) and fur-
ther analyzed.
4. DISCUSSION AND CONCLUSION
The presented non-uniform temporal alignment procedure al-
lows reconstruction of four-dimensional volumes from slice
sequences acquired in living zebrafish embryos, even in cases
where the time series are not strictly periodical. This con-
stitutes a major improvement over previously proposed tech-
niques that only included uniform scaling and translation in
time. We plan further investigations on the method, includ-
ing extensive evaluation and comparisons on synthetic data
as well as validation on a large body of experimentally ac-
quired data. Possible extensions also include deconvolution,
segmentation, and flow analysis of the reconstructed data.
5. ACKNOWLEDGMENTS
We thank Sean Megason, Le Trinh, and Chris Waters for help
with zebrafish preparation and advice on imaging techniques.
M.L. would like to thank Michael Unser for helpful discus-
sions and for sharing unpublished work on related topics.
6. REFERENCES
[1] J. R Hove, R. W. Kö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, vol. 421, pp. 172–177, Jan. 2003.
[2] M. Liebling, A. S. Forouhar, M. Gharib, S. E. Fraser,
and M. E. Dickinson, “Four-dimensional cardiac imag-
ing in living embryos via postacquisition synchronization
of nongated slice sequences,” J. Biomed. Opt., vol. 10,
no. 5, pp. 054001 1–10, 2005.
[3] H. Sakoe and S. Chiba, “Dynamic-programming algo-
rithm optimization for spoken word recognition,” IEEE
Trans. Acoust., Speech, Signal Processing, vol. 26, no. 1,
pp. 43–49, 1978.
[4] C. Yang and M. Stone, “Dynamic programming method
for temporal registration of three-dimensional tongue sur-
face motion from multiple utterances,” Speech Commun.,
vol. 38, no. 1, pp. 201–209, 2002.
[5] R. Bellman, Dynamic Programming, Princeton Univ.
Press, Princeton NJ, USA, 1957.
[6] http://bioimaging.caltech.edu/publications/isbi2006/.
[7] M. Westerfield, The Zebrafish Book, University of Ore-
gon Press, Eugene, 1995.
1159
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 


English

Nonuniform temporal alignment of slice sequences for four-dimensional imaging of cyclically deforming embryonic structures

NONUNIFORM TEMPORAL ALIGNMENT OF SLICE SEQUENCES FOR FOUR-
DIMENSIONAL IMAGING OF CYCLICALLY DEFORMING EMBRYONIC STRUCTURES
Michael Liebling1, Julien Vermot1, Arian S. Forouhar2,
Morteza Gharib2, Mary E. Dickinson3, Scott E. Fraser1
1Biological Imaging Center, Beckman Institute, California Institute of Technology, Pasadena, CA
2Option of Bioengineering, California Institute of Technology, Pasadena, CA
3Dept. of Molecular Physiology and Biophysics, Baylor College of Medicine, Houston, TX
ABSTRACT
The temporal alignment of nongated slice-sequences acquired
at different axial positions in the living embryonic zebraﬁsh
heart permits the reconstruction of dynamic, three-dimen-
sional data. This approach overcomes the current acquisi-
tion-speed limitation of confocal microscopes for real-time
three-dimensional imaging of fast processes. Current syn-
chronization methods align and uniformly scale the data in
time, but do not compensate for slight variations in the heart
rhythm that occur within a heartbeat. Therefore, they impose
constraints on the admissible data quality. Here, we derive
a nonuniform registration procedure based on the minimiza-
tion of the absolute value of the intensity difference between
adjacent slice-sequence pairs. The method compensates for
temporal intra-sample variations and allows the processing of
a wider range of data to build functional, dynamic models of
the beating embryonic heart. We show reconstructions from
data acquired in living, ﬂuorescent zebraﬁsh embryos.
1. INTRODUCTION
With speeds that reach several millimeters per second, heart-
wall motions and blood ﬂow in the embryonic heart require
fast frame-rates to be imaged under the microscope without
inducingmotion artifacts such as blurring or aliasing [1]. Con-
focal microscopes enable selective, three-dimensional imag-
ing of ﬂuorescent cells at depths up to 100-200 µm. Re-
trieving the temporal evolution of ﬂuorescently labeled cells
or blood ﬂow gives invaluable insights for understanding the
driving forces behind the development of the cardio-vascular
system in early embryos.
Although recently developed confocal microscopes can
capture images of two-dimensional slices at frame-rates as
fast as 120 frames per s for 512×512 pixel images, direct
recording of dynamic, three-dimensional data at similar frame-
rates is currently not possible. To overcome this limitation,
Corresponding author is M. Liebling, California Institute of Technology,
Mail Code 139-74, Pasadena, CA 91125, USA, email: liebling@caltech.edu.
z z
(a)
(b) (c) (d)
1 2
t
Fig. 1. Four-dimensional imaging: acquisition and re-
construction procedure. (a) Sequential acquisition of two-
dimensional slices (confocal microscopy) as a time-series at
increasing depths from a cyclically deformed object (cone).
(b) Direct reconstruction is not possible from the nongated
data. (c) Temporal alignment procedure (d) Reconstruction.
an acquisition and reconstruction procedure based on the se-
quential, nongated acquisition of slice-sequences during sev-
eral heart-beats and subsequent synchronization (see Fig. 1)
was recently proposed [2]. In that work, the registration re-
sulted from scaled, temporal shifting of the slice sequences.
This technique is efﬁcient if, during a single two-dimensional
temporal sequence acquisition, the motion is periodic. The
reconstruction precision and quality eventually depend on the
conformity of the measurements to this hypothesis. From our
observations, we noticed that slight irregularities (variation in
the heart beat) of 1% can induce severe visual artifacts that
prevent us from building accurate models and also force us to
discard measurements that do not comply with the regularity
hypothesis.
Here, we propose a nonuniform time-registration method
that recursively aligns pairs of slice-sequences to build a four-
dimensional reconstruction of the beating heart. The method
is based on the minimization of the absolute intensity differ-
ence of the warped sequences’ two-dimensional spatial wave-
let coefﬁcients, an optimization problem that we solve in a
dynamic programming framework.
11560-7803-9577-8/06/$20.00 ©2006 IEEE ISBI 2006
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 
While, to the best of our knowledge, the use of dynamic
programming to solve our image reconstruction problem is
new, similar registration techniques are widespread for differ-
ent purposes, in particular, for speech analysis [3, 4].
In Section 2 we describe the proposed acquisition and re-
construction method. In Section 3 we present reconstructions
that we obtained from in vivo imaging of ﬂuorescent zebraﬁsh
embryos. In Section 4 we discuss our results and future re-
search directions.
2. METHOD
2.1. Acquisition Model and Cost Function Deﬁnition
We model the measured intensity Im as follows,
Im(x, zk, t) =
∫∫∫
I[x′, z, τk(t)]h(x − x′, zk − z) dx′dz,
(1)
where τk(t) is an unknown warping function at depth zk =
khz, k = 0, . . . , Nz. The warping functions τk : R+ →
R+, t → τk(t) are strictly increasing, i.e. if t1 < t2, we have
τk(t1) < τk(t2). The optical system’s point spread function
is denoted h(x, z). We assume that the four-dimensional in-
tensity function I(x, z, τ), has a period T , i.e.,
|I[x, z, τ ]− I[x, z, τ + T ]|  Imax. (2)
We assume that at depth zk¯, the corresponding warping func-
tion is identity, i.e. τk¯(t) = t.
In order to recover an estimate of the volume I(x, z, τ),
τ ∈ [0, L], L > T from the measurements Im(x, zk, t),
where t ∈ [tk, tk + L
′] with L′ > σL (σ > 1 is the max-
imal stretching factor), we consider the following objective
criterion to measure the discrepancy between the warped data
from neighboring depths zk and zk′ ,
Qk,k′,wk{w
′
k} =
∫ L
0
∫∫
R2
∣∣Im[x, zk, wk(τ)]
− Im[x, zk′ , wk′(τ)]
∣∣ dx dτ. (3)
We aim at recovering the warping functions wk′ (τ) for k
′ =
{0, . . . , Nz − 1}\{k¯}. We start from the depth zk¯ for which
wk¯ is known (identity) and ﬁnd the functions wk′ recursively
by aligning neighboring pairs, i.e. (k, k′ = k − 1) for pairs
that are below k¯ and (k, k′ = k + 1) for pairs above the ref-
erence depth k¯. The reconstructed, dynamic volume is com-
puted by warping the measurements, viz.
IR(x, zk, τ) = Im(x, zk, wk(τ)). (4)
To ﬁnd the appropriate functions wk in practice, we proceed
as follows. We consider discrete measurements (in space and
time): for the image sequence that is to be aligned, zk′ , we
consider measurements taken at times t = iht, with i =
0, . . . , Nt − 1 and, for the reference sequence at depth zk,
0 1 2 31
2
3
4
5
6
7
8
2 4 6 8 10 12 14 16
0
1
2
3
Test function
Reference
Fig. 2. Nonuniform temporal alignment via dynamic pro-
gramming. A matrix is generated by taking the absolute-
difference between samples of a test function (top) and a
reference function (right, triangles). Proceeding recursively
from left to right, a cost is associated to each position in the
matrix, based on the value of the matrix at that position, the
cost computed at the most favorable position in the neigh-
boring left column (yielding a vertex), and a possible penalty
if the pick induces stretching or compression of the time axis.
Next, starting from locations on the bottom line, the path with
lowest cost that reaches the top is kept, yielding the optimal
warping function (bold) to interpolate the test function (right,
circles).
measured at times τ = ihτ , with i = 0, . . . , N¯r − 1, we
consider interpolated samples at times τ = ihτ/s, with i =
0, . . . , Nr, Nr = (sN¯r) − 1 and s ∈ N
∗
+ an over-sampling
factor (typically, s = 4). The (discrete) problem boils down to
ﬁnding the appropriate path in a two-dimensional graph, that
pairs the test samples with the interpolated reference samples.
2.2. Dynamic Programming for Temporal Alignment
We denote Ir[] the discrete, interpolated reference slice se-
quence and It[k]) the test slice sequence. We compute the dif-
ferencematrix∆[, k] = ‖Ir[]−It[k]‖, with k = 0, . . . , Nt−
1 and for  = 0, . . . , Nr − 1, and where the norm ‖ · ‖ cor-
responds to a discretized version of Eq. (3). Exploring each
of the (Nr)
Nt possible paths across the matrix and ﬁnding
the one that minimizes Eq. (3) would be of prohibitive com-
plexity and possibly yield solutions that do not comply with
the monotonicity hypothesis. In order to limit the number of
paths to be considered, we use a dynamic programming [5]
approach, resulting in a complexityO(NrNt).
We deﬁneK(, ′) = α(|′−−w0|)
γ+Kmin, α, γ ∈ R∗+.
The limits for stretchability are given by wmin ≤ wmax < 0,
wmin, wmax ∈ Z
−. We recursively compute, for k = 0, . . . , Nt−
1 and for  = 0, . . . , Nr − 1, the following quantities:
L[, k] =
{
L[, k − 1] + 1 if  > 0 and k > 0
1 otherwise,
(5)
1157
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 
80 µm
D
P A
V
D
R L
V
R
P A
L
Eye
Heart
Tube
Heart
Tube
Heart
Tube
RBC
Yolk Sack
(a) (b)
*
*
Fig. 3. Orthogonal slice views of a 38 h.p.f. zebraﬁsh heart reconstruction. (a) Before alignment, temporally misaligned slices
result in black streaks in the heart tube sections (two such slices are indicated by asterisks ∗). (b) After alignment, the heart
tube can be recognized. Red blood cells (RBC) appear as black dots inside the heart tube. Movies are available online [6].
g[, k] = min
′=+wmin,...,+wmax
g′(k, , ′), (6)
¯[, k] = arg min
′=+wmin,...,+wmax
g′(k, , ′), (7)
that is, the length of the path, a cost function, and the index of
the most favorable left neighbor, respectively, and where
g′(k, , ′)
=
⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩
(g[′, k − 1]L[′, k − 1]
+(1− λ)∆[, k] + λK(, ′))
/(L[′, k − 1] + 1) if ′ ≥ 0 and k > 0
(1− λ)∆[, k] + λKmin if 
′ ≥ 0 and k = 0
(1− λ)∆[, k] + λKmin if 
′ < 0
and ′ −  = w0
∞ otherwise.
(8)
The weighting factor λ permits to adjust the warping stiffness.
From the resulting graph (see Fig.2), we determine the lead
of the optimal path by ﬁnding the vertex on the lower image
border with the lowest cost and that yields an acceptable path:
k0 = arg min
k′∈B
g[Nr − 1, k
′] (9)
with B = {k′|1 < k′ < Nt, and ˜k′ [0] ≤ 0} and where the
paths
{
˜k′ [k]
}
0≤k<Nt
are deﬁned for 1 < k′ < Nt
˜k′ [k] =
⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩
Nr − 1 if k = k
′
¯
[
˜k′ [k + 1], k
′ + 1
]
if 0 ≤ k < k′
and ˜k′ [k + 1] ≥ 0
2˜k′ [k − 1]− ˜k′ [k − 2] if k
′ < k < Nt
2˜k′ [k + 1]− ˜k′ [k + 2] otherwise.
(10)
The optimal path is then given by
{
(k0 − i, ˜k0 [i])
}
0<i<Nt
unless B ≡ ∅. The warping function k¯(j), j = s, and
 = 0, ..., N¯r is obtained by inverting ¯k0(j) by use of lin-
ear interpolation.
Similarly to the approach taken in [2], we consider a lim-
ited set of wavelet coefﬁcients instead of pixel values to com-
pute the norm (3). We therefore ensure that the method is
robust and that the computation of the difference matrix is
fast.
3. RESULTS
Wild-type zebraﬁsh (danio rerio) eggs were spawned using
standard techniques [7] and the 38 hours post fertilization
(h.p.f) old embryos were soaked in a green ﬂuorescent dye
(BODIPY FL C5-ceramide, Molecular Probes) for 3 hours.
We acquired two-dimensional image sequences at 60 frames
per s for 2 s, corresponding to 3–4 heartbeats, using a fast
slit-scanning confocal microscope (Zeiss LSM 5 LIVE, Carl
Zeiss Jena GmbH, Germany) with 488nm excitation light and
a 500–525nm band-pass ﬁlter to ﬁlter emission signal. Be-
tween the acquisition of two sequences, the stage was moved
axially by 10 µm or 5 µm, depending on the magniﬁcation
of the microscope objective (10 × or 40 ×, respectively). Per
embryo, we could image a total of 20–30 depths. The im-
ages of 512×512 pixels had a sampling step of 1.78 µm ×
1.78µm per pixel (Zeiss Plan-Neoﬂuar 10×/0.3 microscope
objective), or 0.62µm× 0.62µmper pixel (Zeiss C-Apochro-
mat 40×/1.2 W microscope objective). The images were re-
aligned by use of a MATLAB implementation of the method
described above. The four-dimensional data was analyzed
and visualized using appropriate software (Imaris, Biplane
1158
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 
0 ms 102 ms 204 ms
408 ms306 ms 510 ms
100 µm
Eye Heart
Yolk Sack
Fig. 4. Four-dimensional reconstruction of a 38 h.p.f living zebraﬁsh. The heart-tube is visible in the center, pumping blood.
The corresponding movie is available online [6].
AG, Zurich, Switzerland). The reconstructions prior and af-
ter alignment are presented in Fig. 3 (a) and (b), respectively.
The reconstructed data can be rendered (see Fig. 4) and fur-
ther analyzed.
4. DISCUSSION AND CONCLUSION
The presented non-uniform temporal alignment procedure al-
lows reconstruction of four-dimensional volumes from slice
sequences acquired in living zebraﬁsh embryos, even in cases
where the time series are not strictly periodical. This con-
stitutes a major improvement over previously proposed tech-
niques that only included uniform scaling and translation in
time. We plan further investigations on the method, includ-
ing extensive evaluation and comparisons on synthetic data
as well as validation on a large body of experimentally ac-
quired data. Possible extensions also include deconvolution,
segmentation, and ﬂow analysis of the reconstructed data.
5. ACKNOWLEDGMENTS
We thank Sean Megason, Le Trinh, and Chris Waters for help
with zebraﬁsh preparation and advice on imaging techniques.
M.L. would like to thank Michael Unser for helpful discus-
sions and for sharing unpublished work on related topics.
6. REFERENCES
[1] J. R Hove, R. W. Ko¨ster, A. S. Forouhar, G. Acevedo-
Bolton, S. E. Fraser, and M. Gharib, “Intracardiac ﬂuid
forces are an essential epigenetic factor for embryonic
cardiogenesis,” Nature, vol. 421, pp. 172–177, Jan. 2003.
[2] M. Liebling, A. S. Forouhar, M. Gharib, S. E. Fraser,
and M. E. Dickinson, “Four-dimensional cardiac imag-
ing in living embryos via postacquisition synchronization
of nongated slice sequences,” J. Biomed. Opt., vol. 10,
no. 5, pp. 054001 1–10, 2005.
[3] H. Sakoe and S. Chiba, “Dynamic-programming algo-
rithm optimization for spoken word recognition,” IEEE
Trans. Acoust., Speech, Signal Processing, vol. 26, no. 1,
pp. 43–49, 1978.
[4] C. Yang and M. Stone, “Dynamic programming method
for temporal registration of three-dimensional tongue sur-
face motion from multiple utterances,” Speech Commun.,
vol. 38, no. 1, pp. 201–209, 2002.
[5] R. Bellman, Dynamic Programming, Princeton Univ.
Press, Princeton NJ, USA, 1957.
[6] http://bioimaging.caltech.edu/publications/isbi2006/.
[7] M. Westerﬁeld, The Zebraﬁsh Book, University of Ore-
gon Press, Eugene, 1995.
1159
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on April 14,2010 at 18:51:07 UTC from IEEE Xplore.  Restrictions apply. 


https://authors.library.caltech.edu/22120/1/Liebling2006p89842007_4Th_Ieee_International_Symposium_On_Biomedical_Imaging_Macro_To_Nano_Vols_1-3.pdf

Nonuniform temporal alignment of slice sequences for four-dimensional imaging of cyclically deforming embryonic structures

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main

Caltech Authors - Main