![TABLE I WARPING INDEX (GEOMETRIC ERROR) IN PIXELS FOR THE FOUR TEST SEQUENCES DESCRIBED IN SECTION III-A FOR THE NEW ALGORITHM (ST-Reg), AND THE PREVIOUS METHOD (C-Reg) [37]](/figures/table-i-warping-index-geometric-error-in-pixels-for-the-four-1s4nf5ev.webp){kind=tracking}
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 24, NO. 9, SEPTEMBER 2005 1113
Spatio-Temporal Nonrigid Registration for
Ultrasound Cardiac Motion Estimation
María J. Ledesma-Carbayo*, Member, IEEE, Jan Kybic, Member, IEEE, Manuel Desco,
Andrés Santos, Senior Member, IEEE, Michael Sühling, Student Member, IEEE, Patrick Hunziker, and
Michael Unser, Fellow, IEEE
Abstract—We propose a new spatio-temporal elastic registration
algorithm for motion reconstruction from a series of images. The
specific application is to estimate displacement fields from two-di-
mensional ultrasound sequences of the heart. The basic idea is to
find a spatio-temporal deformation field that effectively compen-
sates for the motion by minimizing a difference with respect to
a reference frame. The key feature of our method is the use of a
semi-local spatio-temporal parametric model for the deformation
using splines, and the reformulation of the registration task as a
global optimization problem. The scale of the spline model con-
trols the smoothness of the displacement field. Our algorithm uses
a multiresolution optimization strategy to obtain a higher speed
and robustness.
We evaluated the accuracy of our algorithm using a synthetic
sequence generated with an ultrasound simulation package, to-
gether with a realistic cardiac motion model. We compared our
new global multiframe approach with a previous method based on
pairwise registration of consecutive frames to demonstrate the ben-
efits of introducing temporal consistency. Finally, we applied the al-
gorithm to the regional analysis of the left ventricle. Displacement
and strain parameters were evaluated showing significant differ-
ences between the normal and pathological segments, thereby il-
lustrating the clinical applicability of our method.
Index Terms—Cardiac motion, elastic registration, parametric
models, splines, temporal models.
I. INTRODUCTION
T
HE ESTIMATION of cardiac motion constitutes an
important aid for the quantification of the elasticity and
contractility of the myocardium. Localized regions exhibiting
Manuscript received May 9, 2005; revised May 17, 2005. This work was sup-
ported in part by the Swiss Science Foundation, the Spanish Health Ministry
under Grant 3200-059517.99/1 (research project PI041495 and Red Temática
IM3 G03/185) and in part by the Czech Ministry of Education under Project
MSM6840770012. The Associate Editor responsible for coordination the re-
view of this paper and recommending its publication was M. Viergever.
Asterisk
indicates corresponding author.
*M. J. Ledesma-Carbayo is with ETSI Telecomunicación, Universidad
Politécnica de Madrid, Ciudad Universitaria s/n, E-28040 Madrid, Spain
(e-mail: mledesma@die.upm.es).
J. Kybic was with Biomedical Imaging Group, EPFL, Switzerland. He is now
with Center for Machine Perception, Czech Technical University, Prague, Czech
Republic (e-mail: kybic@fel.cvut.cz).
M. Desco is with Medicina y Cirugía Experimental, Hospital G. U.
Gregorio Marañón, Dr. Esquerdo 46, E-28007 Madrid, Spain (e-mail:
desco@mce.hggm.es).
A. Santos is with ETSI Telecomunicación, Universidad Politécnica de
Madrid, Ciudad Universitaria s/n, E-28040 Madrid, Spain (e-mail: an-
dres@die.upm.es).
M. Sühling and M. Unser are with the Biomedical Imaging Group, Swiss Fed-
eral Institute of Technology Lausanne (EPFL), CH-1015 Lausanne, Switzerland
(e-mail: michael.suehling@epfl.ch; michael.unser@epfl.ch).
P. Hunziker is with the University Hospital Basel, CH-4031 Basel, Switzer-
land (e-mail: PHunziker@uhbs.ch).
Digital Object Identifier 10.1109/TMI.2005.852050
movement abnormalities are indicative of the existence of
ischemic segments, which are caused by insufficient tissue
microcirculation. Currently, the reference modality for motion
estimation is tagged magnetic resonance (MR) imaging, which
allows us to obtain cardiac displacement fields and derived
parameters, such as the myocardial strain, with high accuracy
[1]–[4]. Most approaches using magnetic resonance imaging
(MRI), single photon emission comuted tomography (SPECT),
and computed tomography (CT) are based on deformable and
mechanical models, and they require a presegmentation step
[2], [3], [5]–[7]. Other methods use energy-based registration
[4], [8], [9] and optical flow techniques [10] to compute the
displacement of the myocardium.
Sequence alignment and registration methods for motion de-
tection have been investigated in computervision[11]–[13].Reg-
istration methods have been used in cardiac imaging; they are
usually applied to data acquired at the same time point in the car-
diac cycle, with the aim of achieving either multimodal integra-
tion [14] or to compensate for small misalignments [15]. Image
registration has also been successfully applied for estimating car-
diac motion in tagged MR data [4], [6], [16]–[19]. Some of these
methods impose spline temporal models to assure temporal con-
sistency and better motion tracking [4], [11], [16], [17].
Our work concentrates on two-dimensional (2-D) echocar-
diography, as it is ubiquitous, and is the most widely used
imaging method to assess cardiac function. The techniques
proposed for cardiac motion recovery in other modalities
cannot be applied directly, because of the especific features of
echocardiographic data.
1) Thesignal-to-noiseratio(SNR)is relativelylow anddepen-
dent on the angle of incidence and depth. Signal dropouts
may appear because of “shadowing,” even though the de-
velopment of recent image acquisition techniques, such as
second harmonic imaging, allow for better performance in
this respect.
2) The complex three-dimensional motion of the heart results
in a partially decorrelated speckle when 2-D sequences are
analyzed [20], [21], therefore making interframe relation
weaker. Similar effects are observed for out-of-plane mo-
tion, which may cause intracardiac structures (for example
papillary muscles or valve cordae) entering and leaving the
view plane, a limitation shared by other 2-D modalities.
Different approaches have been proposed for motion re-
covery from 2-D echocardiagraphic sequences. The most
popular is myocardial border segmentation using deformable
models [22]–[25]. Some of these methods try to overcome
0278-0062/$20.00 © 2005 IEEE
1114 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 24, NO. 9, SEPTEMBER 2005
the intrinsic complexity of the data by introducing a priori
knowledge by modeling a statistical representation of the pos-
sible motions and shapes [24], [25]. However these techniques
estimate the cardiac motion considering the myocardial borders
only; this can lead to inaccurate motion estimations when the
movement is parallel to the border. An additional problem is
that the borders are usually not well defined in echocardio-
graphic data. Another approach is to use optical flow methods
to compute local myocardial movements [26]–[31]. Specific
designs for echocardiographic data consider the Rayleigh
statistics of the signal in the process [29]. Both differential and
block matching techniques add mechanical, spatial or temporal
constraints to overcome the well-known aperture problem with
solutions adapted for echocardiography [26], [27], [30], [31].
The third approach with promising results is to obtain myocar-
dial motion and deformation using speckle tracking [32] and
elastographic techniques [33]–[35]. These methods are based
on the processing of the RF signal to obtain the displacement of
one or several consecutive lines of response, using correlation
and phase shift techniques.
In this paper, we propose using a nonrigid parametric motion
estimation algorithm developed to overcome some of the under-
lying problems inherent to echocardiographic image tracking.
Our approach is global, in the sense that it considers all the
frames in the sequence together, and that it tries to find the most
globally plausible dense spatio-temporal motion field. The de-
formation field is represented using a parametric model based
on B-spline basis functions.
The algorithm does not require any preliminary segmenta-
tion, which would be particularly difficult in the case of car-
diac ultrasound images. The spatio-temporal parametric model
together with a multiresolution optimization strategy provide a
good framework for tracking both global shape and texture. The
multiresolution approach increases speed and tolerance to noise.
The underlying assumption for our approach is that the echo
signal is due to the presence of strong scatterers in the tissue
that produce large and bright speckles [36]. While these scat-
terers move in-plane, they produce a signal component with
stable and visible texture pattern whose displacement is directly
linked with the in-plane cardiac motion. On the other hand,
the out-of-plane movement of the scatterers produces a speckle
component decorrelated with time [20], [21]. Our algorithm is
designed to lock on the temporally coherent part of the signal,
while suppressing the second component as much as possible.
This is achieved by imposing temporal and spatial smoothness
constraints on the deformation field.
The paper is organized as follows. In Section II, we present
our method in detail, covering all the methodological aspects.
In Section III, we evaluate the algorithm using simulated se-
quences derived from a realistic cardiac motion model. In Sec-
tion IV, we present the results obtained from a clinical trial in
which regional cardiac analysis of the left ventricle was per-
formed in a population of patients and in healthy controls.
II. S
PATIO-TEMPORAL REGISTRATION
A. Problem Definition
Let us consider an image sequence
with
and , where is the
intensity at time
and position . Our goal is to find a dense
displacement field over the whole sequence; to this end, we
introduce the deformation function,
, which represents
the position at time
of a point that was at position at time
, i.e., the so-called Lagrangian representation. In other
words, we are using the first frame as a spatial reference,
implying
.
B. Consecutive Registration
This registration method is described in [37], and is based on
the registration of consecutive pairs of images obtained from the
sequence, using an algorithm derived from [38]. This approach
calculates the interframe displacement fields
. The total
deformation field,
, is then obtained from the contribution
of the partial fields.
Registration is performed twice, in the forward and backward
directions, to minimize any error accumulation, and the mean
of the two displacements is used as the final result. We have
also imposed a periodicity on the measurements, as the sequence
encompasses a complete cycle.
In the remainder of the paper, we shall denote this method as
“consecutive elastic registration,” or C-Reg.
C. Spatio-Temporal Registration
In contrast to the consecutive registration method, the
new algorithm presented in this article works globally on all
the images of the sequence simultaneously. It searches for
a spatio-temporal deformation field,
, expressed by a
parametric B-spline model. By applying this field to warp the
original sequence,
, we obtain a motion-corrected sequence,
, that should resemble the reference
frame as much as possible. In other words,
should
appear stationary. The key features of the algorithm are the
similarity criterion (Section II-D), and the spatio-temporal
deformation model (Section II-F).
Fig. 1 shows an example of the spatio-temporal registration
process and the results obtained. The upper part of Fig. 1 shows
three images from the original sequence covering the entire car-
diac cycle. The lower part shows the corresponding images from
the warped (i.e., motion-corrected) sequence.
D. Optimization Criterion
Our registration procedure seeks a minimum value,
, for a criterion, , which is defined
as the mean value obtained from the entire sequence of an
image similarity criterion,
(1)
(2)
where
is the total number of images in the sequence, is the
set of coordinates specifying the spatial region of interest, and
is the corresponding number of pixels. The image criterion,
, is the average of the square of the differences with respect
LEDESMA-CARBAYO et al.: SPATIO-TEMPORAL NON-RIGID REGISTRATION FOR ULTRASOUND CARDIAC MOTION ESTIMATION 1115
Fig. 1. The spatio-temporal registration process, showing the original sequence images (top), and the corresponding images in the warped sequence (bottom),
which tend to appear stationary.
to the reference image at time (and is equivalent to the
sum of squared differences, SSD).
We chose to use the SSD criterion because of its simplicity,
fast computation time, and smoothness of the resulting crite-
rion space. We extended this criterion to the temporal dimen-
sion, and observed that this criterion performed well, even in
the presence of noise and partially decorrelated speckle. We se-
lected the end-diastolic frame as the reference frame, because it
is easily identified in ultrasound sequences from the R-wave of
the electrocardiogram (ECG).
E. Interpolation
It is necessary to have a continuous version of
to be able
to calculate the warped sequence,
, by interpolation, as
well as to be able to evaluate the criterion derivatives. To this
end, we chose to represent
using a 2-D spline interpolation
(3)
where is the tensor product of centered uniform B-splines
of degree
. (Note, was used in all our experiments).
The coefficients,
, were obtained from the pixel values, ,
using filtering [39]. The spline model has the advantage of good
accuracy, low computational complexity, and allows for the pos-
sibility of evaluating spatial derivatives analytically.
F. Spatio-Temporal Model
The deformation function,
, is represented by a linear
model, which is separable in time and space, with parameters,
(4)
where
and define the set of spatial and temporal
parameter indices. The parameter
defines the basis func-
tions in the spatial direction, and is responsible for the spatial
smoothness, and
are the time-axis basis functions that im-
pose the temporal coherence of the deformation. As shown in
[11], [38], and [40], B-splines constitute a good choice for the
spatial basis functions,
. We also used B-splines for the tem-
poral basis functions,
, [4], [11], [16], [17], because of their
computational simplicity, good approximation properties, and
implicit smoothness (minimum curvature property). We found
that temporal B-splines performed at least as well as harmonic
functions (as used in [1], [2], and [5]) in terms of registration ac-
curacy, with the advantage that the criterion minimization was
easier, thanks to their compact support [41]. Specifically, we
used the following basis functions:
(5)
where
and
(6)
1116 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 24, NO. 9, SEPTEMBER 2005
Fig. 2. The B-splines temporal model. The axial and longitudinal displacement of a point in the myocardium is de
fined by the temporal B-spline basis functions.
The individual scaled basis functions are shown, along with their sum, the total trajectory.
The basis functions, , were placed on a uniform rect-
angular spatial grid, and the
were placed at regularly
spaced time intervals. The scale parameters
(space) and
(time) govern the knot spacing and, therefore, the total number
of parameters. These parameters also control the rigidity of
the solution. Typically, we used quadratic B-splines,
,as
, and cubic B-splines, ,as , with and .In
Section III-B2 we analyze the influence of the knot spacings in
more detail. We found that, using cubic B-splines in the spatial
direction did not improve the accuracy significantly and was
not worth the additional computational effort.
G. Motion Field Constraints
The motion model of (5) can be further constrained by using
a priori knowledge of the motion field. This increases the ro-
bustness of the registration process by taking out superfluous
degrees of freedom. First, we know that the displacement at
the reference frame
must be zero. This removes one
degree of freedom from our problem, and leads to a modified
basis function set that only generates displacements satisfying
this constraint
(7)
Similarly, if our sequence contains a full cycle, then we set
the reference frame to , and we impose , leading
to a cyclic set of basis functions defined by
(8)
Fig. 2 shows how these modified basis functions define the
axial and longitudinal displacement for a point in the my-
ocardium. It depicts the individual basis functions scaled using
proper coefficients, as well as the overall trajectory
.
H. Multiresolution and Optimization Strategy
The solution to our registration problem is a deformation
field,
, that minimizes the criterion, . This is found by
using a multidimensional optimization algorithm acting on the
parameters
. The required partial derivatives of can be
calculated explicitly
(9)
while the partial derivatives of
as given by (4) and
as in (3) are
(10)
(11)
We used a gradient descent method with an automatic
step-size update [38]. We applied a multiresolution optimiza-
tion strategy that ensured a robust and efficient approach.
A pyramid of progressively reduced versions of the original
sequence was created by fitting the data using splines with
coarser levels of resolution (a spatio-temporal wavelet-like
pyramid) [42]. This pyramid was compatible with our sequence
model (3), and was optimal in the
-sense. We also used mul-
tiresolution for the motion model, beginning with a coarsely
defined deformation function,
, with a few parameters, ,
and then increasing the number of parameters until the finest
representation of the model was achieved. After converging at a
given level, the result was then used as an initial estimation for
the ensuing, finer level. The projection onto the finer space was
achieved using no approximations, thanks to the embedding
properties of the underlying B-spline spaces.
To summarize, the optimization process proceeded in a
coarse-to-fine fashion for both the image sequence and the mo-
tion field model. The optimization stopped when the changes in
were below a given a priori threshold, . The convergence
speed depends on the number of parameters and the sequence
size. In a typical 350
230 image size, with 35 frames, and
parameters
, with 8 multiresolution levels and
LEDESMA-CARBAYO et al.: SPATIO-TEMPORAL NON-RIGID REGISTRATION FOR ULTRASOUND CARDIAC MOTION ESTIMATION 1117
Fig. 3. First and tenth frames of simulated sequences
S
(top) and
S
(bottom).
a threshold of , the current version of the algorithm
coded in Python needed about half an hour using a standard
PC (
iterations for each multiresolution level). We expect
a fivefold-to-tenfold reduction in computation time when the
algorithm is completely recoded in C. We observed that the
algorithm always converged to a sensible solution.
III. E
XPERIMENTS WITH SIMULATED DATA
This section discusses evaluation experiments on simulated
data. We analyzed the benefits of the temporal model and the
influence of the different algorithm parameters. The use of sim-
ulated sequences allows us to quantify the accuracy of the re-
constructed motion, which would not be possible to obtain with
real data.
As the true cardiac motion was not available, we generated
a realistic cardiac motion field, as described in Appendix A.
The corresponding model was separable, and consisted of two
components: an affine spatial component that simulated radial
myocardial contraction or expansion, and a temporal component
that modulated this movement in a realistic fashion throughout
the cardiac cycle.
A. Simulated Sequences
We generated two different sets of simulated sequences using
the model mentioned above. The first set was defined to explore
the behavior of the algorithm under controlled noise. To gen-
erate this first set, we took one real, end-diastole apical view
image, and deformed it according to the motion model (12).
We corrupted the deformed images using different levels of ad-
ditive Gaussian noise. We generated three sequences with in-
creasing noise levels:
(noiseless), ( dB), and
( dB). Fig. 3 shows the first and tenth frames of
the
and simulated sequences.
Second, we generated a sequence
using the FIELD II ul-
trasound simulation package [43], [44]. This package provides
an excellent framework to simulate ultrasound fields. It incor-
porates realistic transducer features, even though the latest ul-
trasound imaging acquisition technologies, such as second har-
monic imaging or fusion imaging, are not included.
The main purpose of generating this sequence was to include
more realistic ultrasonic features, while keeping known motion.
The simulation of the ultrasound field was based on the compu-
tation of the spatial impulse response, including the excitation
scheme (dynamic focusing and apodization). The images were
generated from a map of independent scatterers with determined
positions and amplitudes [45]. To generate our test sequence,
we specified a typical cardiac transducer (3-MHz central fre-
quency, 64-element phased array with Hanning apodization for
both transmission and reception, and with single focus in trans-
mission and multiple focusing in reception mode).
The phantom scatterer amplitudes and positions were gener-
ated from a sequence of scattering strength maps. These maps
modeled the different densities and speed of sound in the dif-
ferent tissues [45]. We designed the first frame of the scattering
map using a real end-diastole image as a template. The entire se-
quence of maps was then generated by deforming the first map
according to the motion model (12). For each image, 200 000
scatterers were generated using random positions, simulating a
1-cm-thick slice of the heart. The amplitude of the scatterers fol-
lowed a Gaussian distribution that was determined by the scat-
tering map value at each position. The final image was calcu-
lated by summing the responses of all the scatterers, which were
specified by their positions and amplitudes [45]. We used 128
scanning lines to define the image sector, and the resultant image
size was 359
256 pixels. Alone, this process would generate
an unrealistically decorrelated sequence in time as we run inde-
pendent random processes for each frame to select the scatterer
positions. Therefore, we imposed some interframe correlation
through the introduction of stable scatterers in the myocardium.
These composed about 5% of the total number of scatterers in
the myocardium, and their position in the tissue did not change.
Therefore, the sequence,
, was weakly correlated in time,