# Symmetric and Transitive Registration of Image Sequences

Authors: Oskar Škrinjar; Arnaud Bistoquet; Hemant Tagare
Journal: International Journal of Biomedical Imaging
 (2008)

---

## Abstract

This paper presents a method for constructing symmetric and transitive algorithms for registration of image sequences from
image registration algorithms that do not have these two properties. The method is applicable to both rigid and nonrigid registration
and it can be used with linear or periodic image sequences. The symmetry and transitivity properties are satisfied exactly (up to
the machine precision), that is, they always hold regardless of the image type, quality, and the registration algorithm as long as
the computed transformations are invertable. These two properties are especially important in motion tracking applications since
physically incorrect deformations might be obtained if the registration algorithm is not symmetric and transitive. The method was tested on two sequences of cardiac magnetic resonance images using two different nonrigid image registration
algorithms. It was demonstrated that the transitivity and symmetry errors of the symmetric and transitive modification of the
algorithms could be made arbitrary small when the computed transformations are invertable, whereas the corresponding errors
for the nonmodified algorithms were on the order of the pixel size. Furthermore, the symmetric and transitive modification of the
algorithms had higher registration accuracy than the nonmodified algorithms for both image sequences.


---

## Body

## 1. Introduction

The process of aligning images so that the corresponding features can be related is called image registration [1]. Image registration methods have been discussed and
classified in books [1–4] and surveys [5–10].
Most registration methods are ad hoc with assumptions often violated in
practical applications. This results in a behavior that is often not
predictable. A way to reduce the ad hoc nature of registration methods is to
require them to satisfy certain properties. Researchers have realized the
importance of symmetry and transitivity of registration methods [11–20]. In [11], Ashburner et al. proposed an approximately symmetric
image registration method that uses symmetric priors. In [12], Christensen and Johnson
proposed a registration algorithm that approximately satisfies the symmetry
property. (Christensen and Johnson used the term “inverse consistency” for
what we refer to as “symmetry.”) Their idea is to estimate the forward and
reverse transformation simultaneously by minimizing an objective function
composed of terms that measure the similarity between the two images and the
consistency of the forward and reverse transformations. This approach requires
maintaining two transformations, computing their inverses and it has a tradeoff
among the terms in the objective function. In [13], Rogelj and Kovačič proposed a
registration method that uses symmetrically designed forces that deform the two
images. The method is approximately symmetric, it requires maintaining forward
and backward transformation, but it does not use their inverses. In [14], Škrinjar and Tagare proposed
an exactly symmetric registration method that is based on a symmetrically
designed objective function, but it requires the computation of the inverse
transformation. In [15], Lorenzen et al. proposed an exactly symmetric
registration method that is based on a symmetrically designed fluid model. The
method uses two transformations whose compositions define the forward and
backward transformations in such a way that they are inverses of each other.
Beg and Khan in [18]
and Avants et al. in [19] used an exactly symmetric registration method that
maintains two functions, which when composed appropriately give  forward and backward transformations that are exact inverses of each other. In [16], Cachier and Rey analyzed
the reasons behind the asymmetry in registration, proposed symmetrized
similarity and smoothing energies, but their implementation of the method was
not exactly symmetric. In [17], Tagare et al. proposed an exactly symmetric
registration method that does not require to maintain both forward and reverse
transformations and compute their inverses. Instead, the objective function,
which can be based on intensity differences (e.g., mean squared difference,
normalized cross-correlation) or distributions (e.g., mutual information,
normalized mutual information) is modified such that the method is symmetric
and only the forward transformation is needed. Consequently, the objective
function has only one term, which avoids the problem of tradeoff among multiple
terms. Their implementation is symmetric up to the machine precision. In [20], Christensen and Johnson
realized the importance of transitivity of image registration but did not provide a way to satisfy it.The above methods are either approximately or exactly
symmetric but none of them is transitive. In this paper, we propose a method to
modify any image registration algorithm such that it is provably symmetric and
transitive on an image sequence. Symmetry and transitivity are especially
important in motion tracking applications; they insure that a physical point is
tracked in the same way regardless of the order in which the images are
registered. If there are topological changes present in the image sequence, the
two properties can hold only over corresponding regions.Registration of image sequences has a wide
applicability in medical imaging problems. Motion within the body is present at
the system level, organ level, tissue level, cellular level, subcellular level,
and molecular level. In addition to the normal motion, pathology-induced motion
or changes can occur (e.g., osteoporosis, multiple sclerosis, and tumor
growth). In both normal and pathology-induced motion or changes, it is often
useful to compute the motion, that is, to register image sequences. Such
information can improve our understanding of the normal function and diseases
as well as help develop better treatments. The presented approach is
illustrated on sequences of cardiac MR images, which if accurately registered
can provide clinically useful myocardial displacement and strain information.
However, the same or similar approach can be used for the registration of any
other image sequence.


## 2. Methods

### 2.1. Notation

Letℝ denote a set of
real numbers and 𝕊 an N-dimensional
metric space [21]. An N-dimensional
intensity image is a function I:𝕊↦ℝ. Intensity images will be referred to as just images.
Without loss of generality, it is assumed that all the images have the same
domain. The set of all images is denoted as ℐ. An N-dimensional
geometric transformation is a function T:𝕊↦𝕊. Geometric transformations will be referred to as
just transformations. The set of all transformations is denoted as 𝒯. Let Tid denote the
identity transformation, that is, Tid(r)=r,∀r∈𝕊. Let ∘ denote the
composition of transformations.


### 2.2. Image Registration Operator and Its Properties

An image  registration operator is a functionΓ:ℐ2↦𝒯. Ideally, any image registration operator should have
the following three properties ∀I,J,K∈ℐ.(i)
Identity.
An image registration operator, when applied to two identical images should
generate the identity transformation. Formally,(1)Γ(I,I)=Tid.(ii)
Symmetry.
The result of the registration should not depend on the order of images, that
is, when an image registration operator is applied to two images, the obtained
transformation should be the inverse of the transformation obtained when the
order of images is reversed. Formally,(2)Γ(J,I)∘Γ(I,J)=Tid.This is illustrated in Figure 1(a).(iii)
Transitivity.
For any three images, the generated transformation from the second to the third
image composed with the generated transformation from the first to the second
image should be equal to the generated transformation from the first to the
third image. Formally,(3)Γ(J,K)∘Γ(I,J)=Γ(I,K).This is illustrated in Figure 1(b).If two points (A1 and A2) in two images
(I1 and I2) correspond,
as sketched in (a), then the registration algorithm should associate the two
points regardless of the order of images. This is the symmetry property. Let I1, I2, and I3 represent
images of the same deforming object taken at three-time points, and let A1, A2, and A3 represent the
location of the same physical point in the three images, as sketched in (b). If
the registration algorithm, when applied to images I1 and I2, associates points A1 and A2, and, when applied to images I2 and I3, associates points A2 and A3, then it should, when applied to images I1 and I3, associate points A1 and A3. This is the transitivity property.
(a)(b)


### 2.3. Reference-Based Registration

The proposed
approach is simple; select a reference image and then perform the registration
of any two images from the sequence of images through the reference. The
reference can be an image from the sequence of images or an image similar to
the images in the sequence. Let the reference be denoted asR and let Γ represent an
image registration operator that does not necessarily have any of the
properties from Section 2.2. A new image registration operator Ψ is defined
as(4)Ψ(I,J)=[Γ(J,R)]−1∘Γ(I,R),where I and J are any two
images from the sequence. It is assumed that the transformations generated by Γ from images in
the sequence to the reference image are invertable, that is, that [Γ(J,R)]−1 always exists.Theorem 1.
Ψ satisfies
Identity, Symmetry, and Transitivity.Proof.
Identity holds since(5)Ψ(I,I)=[Γ(I,R)]−1∘Γ(I,R)by(4)=Tid. Symmetry holds
since(6)Ψ(J,I)∘Ψ(I,J)=[Γ(I,R)]−1∘Γ(J,R)∘[Γ(J,R)]−1∘Γ(I,R)by(4)=[Γ(I,R)]−1∘Γ(I,R)=Tidby(4). Transitivity holds
since(7)Ψ(J,K)∘Ψ(I,J)=[Γ(K,R)]−1∘Γ(J,R)∘[Γ(J,R)]−1∘Γ(I,R)by(4)=[Γ(K,R)]−1∘Γ(I,R)=Ψ(I,K)by(4).It should be noted that the only requirement for the
above result to hold is that[Γ(J,R)]−1 exists. This
means that the obtained transformations and their inverses do not need to be
differentiable. However, a number of registration methods involve
regularization terms that use transformation derivatives, in which case the
registration operator needs to generate diffeomorphic transformations.If the Jacobian of the transformationΓ(J,R) is positive
then the inverse transformation [Γ(J,R)]−1 exists. If the
Jacobian of the transformation
Γ(J,R) is zero or
negative, the inverse transformation does not exist and the
reference-based registration operator given by (4) cannot be used or it can be
used only over the part of the domain where the Jacobian is positive. Many
registration methods control the Jacobian either directly [22–26] or indirectly [12, 27–30] by using a smoothness term that penalizes extreme
warps to prevent singularities (zero Jacobian) or folding of the space
(negative Jacobian), in which case the inverse transformation exists and the
reference-based registration can be used.


## 2.1. Notation

Letℝ denote a set of
real numbers and 𝕊 an N-dimensional
metric space [21]. An N-dimensional
intensity image is a function I:𝕊↦ℝ. Intensity images will be referred to as just images.
Without loss of generality, it is assumed that all the images have the same
domain. The set of all images is denoted as ℐ. An N-dimensional
geometric transformation is a function T:𝕊↦𝕊. Geometric transformations will be referred to as
just transformations. The set of all transformations is denoted as 𝒯. Let Tid denote the
identity transformation, that is, Tid(r)=r,∀r∈𝕊. Let ∘ denote the
composition of transformations.


## 2.2. Image Registration Operator and Its Properties

An image  registration operator is a functionΓ:ℐ2↦𝒯. Ideally, any image registration operator should have
the following three properties ∀I,J,K∈ℐ.(i)
Identity.
An image registration operator, when applied to two identical images should
generate the identity transformation. Formally,(1)Γ(I,I)=Tid.(ii)
Symmetry.
The result of the registration should not depend on the order of images, that
is, when an image registration operator is applied to two images, the obtained
transformation should be the inverse of the transformation obtained when the
order of images is reversed. Formally,(2)Γ(J,I)∘Γ(I,J)=Tid.This is illustrated in Figure 1(a).(iii)
Transitivity.
For any three images, the generated transformation from the second to the third
image composed with the generated transformation from the first to the second
image should be equal to the generated transformation from the first to the
third image. Formally,(3)Γ(J,K)∘Γ(I,J)=Γ(I,K).This is illustrated in Figure 1(b).If two points (A1 and A2) in two images
(I1 and I2) correspond,
as sketched in (a), then the registration algorithm should associate the two
points regardless of the order of images. This is the symmetry property. Let I1, I2, and I3 represent
images of the same deforming object taken at three-time points, and let A1, A2, and A3 represent the
location of the same physical point in the three images, as sketched in (b). If
the registration algorithm, when applied to images I1 and I2, associates points A1 and A2, and, when applied to images I2 and I3, associates points A2 and A3, then it should, when applied to images I1 and I3, associate points A1 and A3. This is the transitivity property.
(a)(b)


## 2.3. Reference-Based Registration

The proposed
approach is simple; select a reference image and then perform the registration
of any two images from the sequence of images through the reference. The
reference can be an image from the sequence of images or an image similar to
the images in the sequence. Let the reference be denoted asR and let Γ represent an
image registration operator that does not necessarily have any of the
properties from Section 2.2. A new image registration operator Ψ is defined
as(4)Ψ(I,J)=[Γ(J,R)]−1∘Γ(I,R),where I and J are any two
images from the sequence. It is assumed that the transformations generated by Γ from images in
the sequence to the reference image are invertable, that is, that [Γ(J,R)]−1 always exists.Theorem 1.
Ψ satisfies
Identity, Symmetry, and Transitivity.Proof.
Identity holds since(5)Ψ(I,I)=[Γ(I,R)]−1∘Γ(I,R)by(4)=Tid. Symmetry holds
since(6)Ψ(J,I)∘Ψ(I,J)=[Γ(I,R)]−1∘Γ(J,R)∘[Γ(J,R)]−1∘Γ(I,R)by(4)=[Γ(I,R)]−1∘Γ(I,R)=Tidby(4). Transitivity holds
since(7)Ψ(J,K)∘Ψ(I,J)=[Γ(K,R)]−1∘Γ(J,R)∘[Γ(J,R)]−1∘Γ(I,R)by(4)=[Γ(K,R)]−1∘Γ(I,R)=Ψ(I,K)by(4).It should be noted that the only requirement for the
above result to hold is that[Γ(J,R)]−1 exists. This
means that the obtained transformations and their inverses do not need to be
differentiable. However, a number of registration methods involve
regularization terms that use transformation derivatives, in which case the
registration operator needs to generate diffeomorphic transformations.If the Jacobian of the transformationΓ(J,R) is positive
then the inverse transformation [Γ(J,R)]−1 exists. If the
Jacobian of the transformation
Γ(J,R) is zero or
negative, the inverse transformation does not exist and the
reference-based registration operator given by (4) cannot be used or it can be
used only over the part of the domain where the Jacobian is positive. Many
registration methods control the Jacobian either directly [22–26] or indirectly [12, 27–30] by using a smoothness term that penalizes extreme
warps to prevent singularities (zero Jacobian) or folding of the space
(negative Jacobian), in which case the inverse transformation exists and the
reference-based registration can be used.


## 3. Results

While the
result of the previous section holds for any registration operator that
generates invertable transformations, here we illustrate the approach on two
sequences of cardiac magnetic resonance images using two nonrigid image
registration algorithms.


### 3.1. MR Protocols

We acquired a
3D anatomical cine MRI scan together with a 3D tagged cine MRI scan of a
healthy volunteer and then repeated the acquisitions four months later. The
volunteer was a 27-year-old male subject with no history of heart disease. The
purpose of the tagged scan was to validate the myocardial deformation recovered
from the anatomical scan. The scans were acquired using steady-state
free-precession short axis cine imaging (flip angle  =  65°, TR  =  3.4 ms, TE  =  1.7 ms) covering the entire heart
on a 1.5 T clinical MRI scanner (Intera, Philips Medical Systems, Best, The
Netherlands). All the scans had contiguous short-axis slices with similar field
of view and phases covering the entire cardiac cycle and their parameters are
given in Table 1. The tags were applied immediately after the detection of the
R-wave from the EKG signal and the first frame was acquired at a delay of 15
milliseconds after the R-wave. Two orthogonal sets of parallel planar tags with
about 9 mm plane separation were oriented orthogonal to the image planes.Table 1
Parameters of the anatomical and tagged MR scans.
ParameterScan 1Scan 2AnatomicalTaggedAnatomicalTaggedIn-plane resolution [mm]1.411.241.441.44Number of slices12151717Slice thickness [mm]101088Temporal resolution [ms]38303530For both scans, for each acquired slice the scanner
recorded the rigid body transformation from the scanner coordinate system to
the slice. This allowed us to map all the slices to a common coordinate system,
that is, to spatially align the anatomical and tagged scans. Similarly, the
scanner recorded the start time for each phase (frame) relative to the peak of
the R wave, which allowed us to temporally align the anatomical and tagged
scans. Since the heart rate, that is, the duration of the cardiac cycle, was
not the same for the anatomical and tagged scans, we used relative time (as a
percentage of the cardiac cycle) for the temporal alignment.


### 3.2. Myocardial Deformation Recovery

To recover the
myocardial deformation, we use thin plate splines (TPS) [31] to represent the
transformation between any two frames and then maximize the normalized mutual
information [32] to determine
the transformation parameters (TPS node positions). We use normalized mutual
information since it was shown to outperform several other images similarity
measures [33]. Since
myocardium is nearly incompressible and its volume does not change by more than
4% over the cardiac cycle [26], we constrain the optimization of the TPS node
positions such that the Jacobian of the transformation never deviates from 1
(which corresponds to exact incompressibility) by more than 4%. The details of
the method are given in [26]. (The purpose of this section is to illustrate the
approach of Section 2.3 (construction of symmetric and
transitive registration algorithms from nonsymmetric and nontransitive
registration algorithms), and instead of the method given in [26] we could have used any
other registration method. For this reason we did not present here all the
details of the used registration method and the interested reader is referred
to [26].)   This registration
algorithm we denote as Γ1, while Γ2 represents its
unconstrained version. Given that near incompressibility is a physical property
of the myocardium, Γ1 is expected to
be more accurate than Γ2. Each of the two operators was used to recover
myocardial deformation from the two cardiac MR image sequences in two ways:
sequential and reference-based. In the sequential approach, the deformation was
recovered from the first to the second frame, then from the second to the third
frame, and so on. In the reference-based approach, the deformation was
recovered directly from the reference frame to any given frame. Figure 2 shows
a short-axis section through the 3D image and 3D LV wall model surface for the
sequential and referenced-based recoveries using the two registration
operators.The recovered myocardial deformation for a normal subject over the
cardiac cycle (first row: end diastole; third row: end systole) is shown by
means of the endocardial and endocardial surface model contours overlaid over a
midventricular short-axis MRI slice. The myocardium was segmented in the first
frame (shown in the first row), a surface model was generated around the
segmented myocardium and the recovered deformation for the rest of the sequence
was applied to the surface model. The two red contours represent a
cross-section through the surface model. The four columns correspond to (a)
sequential recovery byΓ1, (b) sequential recovery by Γ2, (c) reference-based recovery by Γ1, and (d) reference-based recovery by Γ2. The registration algorithms Γ1 and Γ2 and the
difference between sequential and reference-based recovery are explained in
Section 3.2. Note that the best deformation recovery,
that is, the best agreement of the red contours and the edges of the left
ventricular wall, was achieved for the reference-based recovery by Γ1, shown in column (c).
(a)(b)(c)(d)To quantitatively evaluate the deformation recovery
accuracy we compared the cardiac deformation recovered from the anatomical cine
MRI against the corresponding information from the tagged cine MRI. We
developed a tool for interactive positioning of virtual tag planes in tagged MRI
scans. The tag planes are modeled as splines that the user can position and
deform by moving control points. This allows the user to fit the virtual tag
planes to the tagged image as well as to compute tag plane intersections. Once
the cardiac deformation is recovered from the anatomical cine MRI using the
proposed method, it is applied to the virtual tag planes at ED and then
compared to the interactively positioned tag planes in other frames. In each
image slice, we compute the distances between the corresponding intersections
of orthogonal virtual tag lines (in-slice cross-sections of the virtual tag
planes) generated interactively and by the model. This allows for in-plane
(short-axis) deformation recovery validation. The out-of-plane (long-axis) deformation
is not evaluated with this procedure since the tag planes, being perpendicular
to the short-axis image slices, do not contain information about the
out-of-plane motion. Virtual tag lines for the sequential and referenced-based
recoveries for the two operators are shown in Figure3. Table 2 contains the
distances between corresponding intersections of virtual and real tag lines at
end systole, which is the most deformed state relative to end diastole.Table 2
The mean (± std) distance [mm]  between corresponding
intersections of virtual and real tag lines at end systole for the four
algorithms for both scans.
SequentialΓ1SequentialΓ2Reference-basedΓ1Reference-basedΓ2Scan 12.3 ± 0.52.6 ± 0.61.9 ± 0.32.2 ± 0.4Scan 22.2 ± 0.42.4 ± 0.41.7 ± 0.22.1 ± 0.3The virtual tag lines and the corresponding short-axis slices of the tagged MRI scan are shown over the cardiac cycle (first row: end diastole; third row: end
systole) for (a) sequential recovery byΓ1, (b) sequential recovery by Γ2, (c) reference-based recovery by Γ1, and (d) reference-based recovery by Γ2. The registration algorithms Γ1 and Γ2 and the
difference between sequential and reference-based recovery are explained in
Section 3.2. The virtual tag lines were manually
positioned over the tagged MR image in the first frame (shown in the first row),
and then the deformation recovered from the anatomical image sequence was applied
to the virtual tag lines and they were overlaid over the tagged MR images in
the corresponding frames.
(a)(b)(c)(d)


### 3.3. Identity, Symmetry, and Transitivity Errors

LetΨ1 and Ψ2 represent the
symmetric and transitive modifications of Γ1 and Γ2, respectively. Operators Ψ1 and Ψ2 are defined by
(4) (the end diastole frame is used for R), which
involves transformation inversion. Since we use TPS for transformation representation
and the inverse of a TPS transformation cannot be obtained analytically, we
invert the transformation numerically using the Newton-Raphson method for
solving nonlinear systems of equations [34]. The numerical error of the computation of the
inverse transform is denoted as ϵ. If the Jacobian of the transformation is positive
then the inverse transformation exists and ϵ can be
specified to be an arbitrary small positive number, that is, the inverse
transformation can be computed with an arbitrary small error. While ϵ can be set to
an arbitrary small positive value, in practical applications little is gained
if ϵ is set to a
value smaller than two orders of magnitude smaller than the pixel size. The
reason for that is that the registration error is on the order of the pixel
size, and by setting ϵ to one
hundredth the pixel size the error of the computation of the inverse
transformation is already negligible compared to the registration error, and
further reducing ϵ does not
improve the registration accuracy.For a given image registration operatorΓ, we define the following three errors.(i)
Identity error of image I is(8)Eiden=maxr∈𝕊∥Γ(I,I)(r)−r∥.(ii)
Symmetry
error of images I and J is(9)Esym=maxr∈𝕊∥Γ(J,I)(Γ(I,J)(r))−r∥.(iii)
Transitivity
error of images I, J, and K is(10)Etran=maxr∈𝕊∥Γ(J,K)(Γ(I,J)(r))−Γ(I,K)(r)∥.Identity, symmetry and transitivity errors forΓ1, Γ2, Ψ1, and Ψ2 are given in
Tables 3, 4, and 5, respectively. The errors for Ψ1 and Ψ2 in the three
tables were computed using ϵ  =  0.001 mm. The
dependence of Eiden, Esym, and Etran on ϵ is depicted in
Figure 4 for Ψ1.Table 3
Identity errors
[mm]  forΓ1, Γ2, Ψ1, and Ψ2 are given for
random frames (F) of the two cardiac cine MRI scans.
Scan 1Scan 2FΓ1Γ2Ψ1Ψ2Γ1Γ2Ψ1Ψ23000.00120.0014000.00110.00147000.00130.0011000.00120.001112000.00120.0012000.00140.001314000.00110.0013000.00130.001119000.00140.0011000.00110.0012Table 4
Symmetry errors
[mm]  forΓ1, Γ2, Ψ1, and Ψ2 are given for
random frames (F1 and F2) of the two
cardiac cine MRI scans.
Scan 1Scan 2F1F2Γ1Γ2Ψ1Ψ2Γ1Γ2Ψ1Ψ24121.11.30.00230.00181.10.90.00170.0021790.60.80.00190.00220.60.70.00200.001914190.80.90.00210.00200.81.10.00150.00231861.11.50.00160.00151.21.50.00170.0016250.60.50.00220.00260.70.80.00190.0025Table 5
Transitivity
errors [mm]  forΓ1, Γ2, Ψ1, and Ψ2 are given for
random frames (F1, F2, and F3) of the two cardiac cine MRI scans.
Scan 1Scan 2F1F2F3Γ1Γ2Ψ1Ψ2Γ1Γ2Ψ1Ψ21371.61.80.00260.00192.22.50.00220.001549151.72.20.00220.00211.81.70.00180.0023612182.21.90.00160.00232.52.60.00230.002082142.22.30.00170.00201.92.30.00170.001656171.92.10.00190.00222.52.70.00210.0022Figure 4
The dependence
ofEiden (dotted), Esym (solid), and Etran (dashed) on ϵ for Ψ1 is shown in the
log-log axes for a representative image, image pair, and image triple,
respectively. The three curves represent interpolations of the errors
corresponding to ϵ of 0.000001 mm,
0.00001 mm, 0.0001 mm, 0.001 mm, 0.01 mm, 0.1 mm, and 1.0 mm. The squared
correlation coefficient between Eiden and ϵ is .9998,
between Esym and ϵ is .99997, and
between Etran and ϵ is .9998,
indicating a strong dependence of the three errors on ϵ.


## 3.1. MR Protocols

We acquired a
3D anatomical cine MRI scan together with a 3D tagged cine MRI scan of a
healthy volunteer and then repeated the acquisitions four months later. The
volunteer was a 27-year-old male subject with no history of heart disease. The
purpose of the tagged scan was to validate the myocardial deformation recovered
from the anatomical scan. The scans were acquired using steady-state
free-precession short axis cine imaging (flip angle  =  65°, TR  =  3.4 ms, TE  =  1.7 ms) covering the entire heart
on a 1.5 T clinical MRI scanner (Intera, Philips Medical Systems, Best, The
Netherlands). All the scans had contiguous short-axis slices with similar field
of view and phases covering the entire cardiac cycle and their parameters are
given in Table 1. The tags were applied immediately after the detection of the
R-wave from the EKG signal and the first frame was acquired at a delay of 15
milliseconds after the R-wave. Two orthogonal sets of parallel planar tags with
about 9 mm plane separation were oriented orthogonal to the image planes.Table 1
Parameters of the anatomical and tagged MR scans.
ParameterScan 1Scan 2AnatomicalTaggedAnatomicalTaggedIn-plane resolution [mm]1.411.241.441.44Number of slices12151717Slice thickness [mm]101088Temporal resolution [ms]38303530For both scans, for each acquired slice the scanner
recorded the rigid body transformation from the scanner coordinate system to
the slice. This allowed us to map all the slices to a common coordinate system,
that is, to spatially align the anatomical and tagged scans. Similarly, the
scanner recorded the start time for each phase (frame) relative to the peak of
the R wave, which allowed us to temporally align the anatomical and tagged
scans. Since the heart rate, that is, the duration of the cardiac cycle, was
not the same for the anatomical and tagged scans, we used relative time (as a
percentage of the cardiac cycle) for the temporal alignment.


## 3.2. Myocardial Deformation Recovery

To recover the
myocardial deformation, we use thin plate splines (TPS) [31] to represent the
transformation between any two frames and then maximize the normalized mutual
information [32] to determine
the transformation parameters (TPS node positions). We use normalized mutual
information since it was shown to outperform several other images similarity
measures [33]. Since
myocardium is nearly incompressible and its volume does not change by more than
4% over the cardiac cycle [26], we constrain the optimization of the TPS node
positions such that the Jacobian of the transformation never deviates from 1
(which corresponds to exact incompressibility) by more than 4%. The details of
the method are given in [26]. (The purpose of this section is to illustrate the
approach of Section 2.3 (construction of symmetric and
transitive registration algorithms from nonsymmetric and nontransitive
registration algorithms), and instead of the method given in [26] we could have used any
other registration method. For this reason we did not present here all the
details of the used registration method and the interested reader is referred
to [26].)   This registration
algorithm we denote as Γ1, while Γ2 represents its
unconstrained version. Given that near incompressibility is a physical property
of the myocardium, Γ1 is expected to
be more accurate than Γ2. Each of the two operators was used to recover
myocardial deformation from the two cardiac MR image sequences in two ways:
sequential and reference-based. In the sequential approach, the deformation was
recovered from the first to the second frame, then from the second to the third
frame, and so on. In the reference-based approach, the deformation was
recovered directly from the reference frame to any given frame. Figure 2 shows
a short-axis section through the 3D image and 3D LV wall model surface for the
sequential and referenced-based recoveries using the two registration
operators.The recovered myocardial deformation for a normal subject over the
cardiac cycle (first row: end diastole; third row: end systole) is shown by
means of the endocardial and endocardial surface model contours overlaid over a
midventricular short-axis MRI slice. The myocardium was segmented in the first
frame (shown in the first row), a surface model was generated around the
segmented myocardium and the recovered deformation for the rest of the sequence
was applied to the surface model. The two red contours represent a
cross-section through the surface model. The four columns correspond to (a)
sequential recovery byΓ1, (b) sequential recovery by Γ2, (c) reference-based recovery by Γ1, and (d) reference-based recovery by Γ2. The registration algorithms Γ1 and Γ2 and the
difference between sequential and reference-based recovery are explained in
Section 3.2. Note that the best deformation recovery,
that is, the best agreement of the red contours and the edges of the left
ventricular wall, was achieved for the reference-based recovery by Γ1, shown in column (c).
(a)(b)(c)(d)To quantitatively evaluate the deformation recovery
accuracy we compared the cardiac deformation recovered from the anatomical cine
MRI against the corresponding information from the tagged cine MRI. We
developed a tool for interactive positioning of virtual tag planes in tagged MRI
scans. The tag planes are modeled as splines that the user can position and
deform by moving control points. This allows the user to fit the virtual tag
planes to the tagged image as well as to compute tag plane intersections. Once
the cardiac deformation is recovered from the anatomical cine MRI using the
proposed method, it is applied to the virtual tag planes at ED and then
compared to the interactively positioned tag planes in other frames. In each
image slice, we compute the distances between the corresponding intersections
of orthogonal virtual tag lines (in-slice cross-sections of the virtual tag
planes) generated interactively and by the model. This allows for in-plane
(short-axis) deformation recovery validation. The out-of-plane (long-axis) deformation
is not evaluated with this procedure since the tag planes, being perpendicular
to the short-axis image slices, do not contain information about the
out-of-plane motion. Virtual tag lines for the sequential and referenced-based
recoveries for the two operators are shown in Figure3. Table 2 contains the
distances between corresponding intersections of virtual and real tag lines at
end systole, which is the most deformed state relative to end diastole.Table 2
The mean (± std) distance [mm]  between corresponding
intersections of virtual and real tag lines at end systole for the four
algorithms for both scans.
SequentialΓ1SequentialΓ2Reference-basedΓ1Reference-basedΓ2Scan 12.3 ± 0.52.6 ± 0.61.9 ± 0.32.2 ± 0.4Scan 22.2 ± 0.42.4 ± 0.41.7 ± 0.22.1 ± 0.3The virtual tag lines and the corresponding short-axis slices of the tagged MRI scan are shown over the cardiac cycle (first row: end diastole; third row: end
systole) for (a) sequential recovery byΓ1, (b) sequential recovery by Γ2, (c) reference-based recovery by Γ1, and (d) reference-based recovery by Γ2. The registration algorithms Γ1 and Γ2 and the
difference between sequential and reference-based recovery are explained in
Section 3.2. The virtual tag lines were manually
positioned over the tagged MR image in the first frame (shown in the first row),
and then the deformation recovered from the anatomical image sequence was applied
to the virtual tag lines and they were overlaid over the tagged MR images in
the corresponding frames.
(a)(b)(c)(d)


## 3.3. Identity, Symmetry, and Transitivity Errors

LetΨ1 and Ψ2 represent the
symmetric and transitive modifications of Γ1 and Γ2, respectively. Operators Ψ1 and Ψ2 are defined by
(4) (the end diastole frame is used for R), which
involves transformation inversion. Since we use TPS for transformation representation
and the inverse of a TPS transformation cannot be obtained analytically, we
invert the transformation numerically using the Newton-Raphson method for
solving nonlinear systems of equations [34]. The numerical error of the computation of the
inverse transform is denoted as ϵ. If the Jacobian of the transformation is positive
then the inverse transformation exists and ϵ can be
specified to be an arbitrary small positive number, that is, the inverse
transformation can be computed with an arbitrary small error. While ϵ can be set to
an arbitrary small positive value, in practical applications little is gained
if ϵ is set to a
value smaller than two orders of magnitude smaller than the pixel size. The
reason for that is that the registration error is on the order of the pixel
size, and by setting ϵ to one
hundredth the pixel size the error of the computation of the inverse
transformation is already negligible compared to the registration error, and
further reducing ϵ does not
improve the registration accuracy.For a given image registration operatorΓ, we define the following three errors.(i)
Identity error of image I is(8)Eiden=maxr∈𝕊∥Γ(I,I)(r)−r∥.(ii)
Symmetry
error of images I and J is(9)Esym=maxr∈𝕊∥Γ(J,I)(Γ(I,J)(r))−r∥.(iii)
Transitivity
error of images I, J, and K is(10)Etran=maxr∈𝕊∥Γ(J,K)(Γ(I,J)(r))−Γ(I,K)(r)∥.Identity, symmetry and transitivity errors forΓ1, Γ2, Ψ1, and Ψ2 are given in
Tables 3, 4, and 5, respectively. The errors for Ψ1 and Ψ2 in the three
tables were computed using ϵ  =  0.001 mm. The
dependence of Eiden, Esym, and Etran on ϵ is depicted in
Figure 4 for Ψ1.Table 3
Identity errors
[mm]  forΓ1, Γ2, Ψ1, and Ψ2 are given for
random frames (F) of the two cardiac cine MRI scans.
Scan 1Scan 2FΓ1Γ2Ψ1Ψ2Γ1Γ2Ψ1Ψ23000.00120.0014000.00110.00147000.00130.0011000.00120.001112000.00120.0012000.00140.001314000.00110.0013000.00130.001119000.00140.0011000.00110.0012Table 4
Symmetry errors
[mm]  forΓ1, Γ2, Ψ1, and Ψ2 are given for
random frames (F1 and F2) of the two
cardiac cine MRI scans.
Scan 1Scan 2F1F2Γ1Γ2Ψ1Ψ2Γ1Γ2Ψ1Ψ24121.11.30.00230.00181.10.90.00170.0021790.60.80.00190.00220.60.70.00200.001914190.80.90.00210.00200.81.10.00150.00231861.11.50.00160.00151.21.50.00170.0016250.60.50.00220.00260.70.80.00190.0025Table 5
Transitivity
errors [mm]  forΓ1, Γ2, Ψ1, and Ψ2 are given for
random frames (F1, F2, and F3) of the two cardiac cine MRI scans.
Scan 1Scan 2F1F2F3Γ1Γ2Ψ1Ψ2Γ1Γ2Ψ1Ψ21371.61.80.00260.00192.22.50.00220.001549151.72.20.00220.00211.81.70.00180.0023612182.21.90.00160.00232.52.60.00230.002082142.22.30.00170.00201.92.30.00170.001656171.92.10.00190.00222.52.70.00210.0022Figure 4
The dependence
ofEiden (dotted), Esym (solid), and Etran (dashed) on ϵ for Ψ1 is shown in the
log-log axes for a representative image, image pair, and image triple,
respectively. The three curves represent interpolations of the errors
corresponding to ϵ of 0.000001 mm,
0.00001 mm, 0.0001 mm, 0.001 mm, 0.01 mm, 0.1 mm, and 1.0 mm. The squared
correlation coefficient between Eiden and ϵ is .9998,
between Esym and ϵ is .99997, and
between Etran and ϵ is .9998,
indicating a strong dependence of the three errors on ϵ.


## 4. Discussion

Theorem1 says
that reference-based modification of any registration operator satisfies
identity, symmetry, and transitivity on an image sequence. While the theorem
always holds and the modified registration operator satisfies the three
properties exactly, Section 3 demonstrates a practical application of the
reference-based registration using a relatively accurate registration algorithm
(Γ1) and its less accurate version (Γ2). The method was applied to two cardiac cine MRI
scans, which are periodic image sequences, but the same conclusions hold in the
case of linear image sequences.As expected, the constrained registration (Γ1) outperformed its unconstrained version (Γ2), which can be seen in Figures 2 and 3 and in Table
2 for both sequential and reference-based approaches. The two figures and the
table also show that in this case the reference-based registration was more
accurate than the sequential registration. The advantage of the sequential
registration is that the difference between any two consecutive frames is
relatively small, while the reference-based registration deals with larger
deformations (e.g., from end diastole to end systole). However, the
disadvantage of the sequential registration is that there is accumulation of
error from frame to frame, which seems to overweigh the advantage of small
frame-to-frame difference. There is no accumulation of error for
reference-based registration. The reference-based constrained registration
(column (c) in Figures 2 and 3) was more accurate than the other three algorithms
(Table 2). The difference in accuracy of the four algorithms is most
pronounced at end systole, and it can be seen as the different level of
agreement between the model contours and the underlying image in Figure 2 and
between the virtual tag lines and the underlying image in Figure 3.SinceΓ1 and Γ2 use Tid as the initial
transformation in searching for the transformation that maximizes the
normalized mutual information, both operators satisfy (1) exactly and
consequently have Eiden=0 for any image,
as it can be seen in Table 3. Operators Ψ1 and Ψ2 involve
transformation inversion as defined in (4), which is done numerically with an
accuracy of ϵ. This is why the identity errors for Ψ1 and Ψ2 in Table 3
are approximately equal to ϵ.The symmetry errors (Table4) for Γ1 and Γ2 are
approximately 1 mm, while for Ψ1 and Ψ2 they are
approximately 2ϵ. The reason for this is that the evaluation of (9)
involves a composition of the operators, each contributing approximately ϵ to Esym due to the
numerical inversion of transformation in (4).Similarly, the transitivity errors (Table5) for Γ1 and Γ2 are
approximately 2 mm, while for Ψ1 and Ψ2 they are
approximately 2ϵ. The reason for this is that the evaluation of (10)
involves a composition of the operators, each contributing approximately ϵ to Etran due to the
numerical inversion of transformation in (4).It should be noted thatΨ1 and Ψ2 have similar
values for Eiden, Esym, and Etran. The reason for this is that the three errors depend
only on the accuracy of the computation of the inverse transformation, and not
on the registration accuracy (Ψ1 is more accurate than Ψ2). In fact, if
the inverse transformation could be computed exactly, the three errors would be
zero regardless of the registration operator. The three errors scale with ϵ, as shown in Figure 4, and they can be made as small
as the machine precision. Figure 4 also shows that Eiden≈ϵ, Esym≈2ϵ, and Etran≈2ϵ. Thus, the reference-based registration slightly
worsens the identity error and it significantly improves the symmetry and
transitivity errors over the sequential registration.For the two tested image sequences, the symmetric and
transitive registration methodsΨ1 and Ψ2 had smaller
registration errors than their nonsymmetric and nontransitive counterparts Γ1 and Γ2, respectively (Table 2). While all four registration
methods had either zero or nearly zero identity errors (Table 3), Ψ1 and Ψ2 had very small
(nearly zero) symmetry and transitivity errors and Γ1 and Γ2 had these
errors on the order of the pixel size or even larger (Tables 4 and 5). Thus,
in this limited study, very small (nearly zero) symmetry and transitivity
errors were accompanied by reduced registration errors. However, to determine
if this is always the case one would need to prove it mathematically or at
least repeat the experiment on a large number of image sequences.To simplify the notation, it was assumed that all the
images had the same domain, but the same conclusions hold when the domains are
different. Furthermore, the method involves the transformation inversion, which
is done only once (after both images are registered to the reference), as
opposed to the methods proposed in [12, 14–16] that require computing the
inverse transformation in each iteration of the optimization.We used the normalized mutual information as the image
similarity measure for the registration algorithmsΓ1 and Γ2. We repeated all the experiments by using the mutual
information, mean square difference, and normalized cross-correlation as
alternative image similarity measures and the obtained results were nearly
identical to those reported in Section 3 and for this reason they are not
included in the paper. These repeated experiments confirm that the conclusions
of the paper are not specific to the normalized mutual information.


## 5. Conclusion

The
reference-based registration of image sequences is provably symmetric and
transitive. This conclusion is independent of the images and registration
algorithm used. Furthermore, a limited study showed that the reference-based
registration was more accurate than the sequential registration, although this
cannot be guaranteed to always hold.



---

*Source: 686875-2009-03-24.xml*