PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ijbiInternational Journal of Biomedical Imaging
 
Int J Biomed Imaging. 2011; 2011: 891585.
Published online 2010 December 8. doi:  10.1155/2011/891585
PMCID: PMC3005125

Diffeomorphic Registration of Images with Variable Contrast Enhancement

Abstract

Nonrigid image registration is widely used to estimate tissue deformations in highly deformable anatomies. Among the existing methods, nonparametric registration algorithms such as optical flow, or Demons, usually have the advantage of being fast and easy to use. Recently, a diffeomorphic version of the Demons algorithm was proposed. This provides the advantage of producing invertible displacement fields, which is a necessary condition for these to be physical. However, such methods are based on the matching of intensities and are not suitable for registering images with different contrast enhancement. In such cases, a registration method based on the local phase like the Morphons has to be used. In this paper, a diffeomorphic version of the Morphons registration method is proposed and compared to conventional Morphons, Demons, and diffeomorphic Demons. The method is validated in the context of radiotherapy for lung cancer patients on several 4D respiratory-correlated CT scans of the thorax with and without variable contrast enhancement.

1. Introduction

In the context of image-based medical diagnostics and treatment, highly deformable anatomies are a problem for multiple-time imaging analysis along the course of treatment. Indeed, a precise tracking of organs is made difficult because of shape and position variations. Nonrigid registration may be used to compute a displacement vector for each voxel of an image [1], enabling the estimation of the spatial variations of the anatomy. The displacement vectors are computed as pointing to the best corresponding location of the voxels in another image according to a metric which is a measure of the image matching and under some constraints on global properties of the resulting deformation, such as invertibility and smoothness.

Several registration methods have been used in the past years to estimate deformations in highly deformable anatomies [25]. Many efforts have been made to improve the quality of displacement estimates and also to reduce the amount of required preprocessing or modeling and improve registration speed [68]. Besides, the choice of a registration method for medical application depends on the characteristics (e.g., modality) of the images to be registered [1]. The existing methods [9] can be divided into parametric, or model-based, methods (B-splines [10], thin-plate splines [11], radial basis functions [12], linear elastic FEM [13], etc.) and nonparametric methods (viscous fluid [14], optical flow [15], etc.). In this second category, the algorithm called Demons [16, 17] is fast, efficient, and easy to use, as it requires no particular preprocessing nor patient-specific modeling. This method aims at calculating a regular displacement field which produces a good matching of the intensities in both images by minimizing a metric, such as the sum of squared differences (SSDs) [18] or the mutual information (MI) [8] between images along with a measure of the field regularity.

In a growing number of applications, the displacement fields resulting from registration are used to deform images from other modalities or other spatial distribution maps (e.g., the dose map associated to CT scans in radiotherapy [19, 20]). Therefore, the matching of structures in images based on their intensities is not a sufficient constraint for producing realistic anatomical deformation estimations [21]. This is the reason why a priori information on the physical characteristics of anatomical deformations has to be included in the registration process. Diffeomorphism is a necessary condition for displacement fields to be physical [22]. Indeed, organs can be compressed and deformed, but cannot undergo noninvertible spatial transformations, for example, showing mirror effects. A method has been proposed in [23] to limit the displacement fields computed by the Demons to a set of diffeomorphic transformations, using diffeomorphic flows and Lie algebra.

In several medical protocols, contrast agents are used in order to facilitate interpretation. This makes the registration problem incompatible with the hypothesis of intensity conservation. Furthermore, an histogram equalization is often not able to correct for contrast agent variability, as different regions will be enhanced in different ways inside the image. Therefore, simple metrics, such as SSD or cross-correlation, are not suitable for matching those images, and methods that are suitable for registering variable contrast images have to be investigated [24, 25].

A method similar to Demons but using a phase-based approach was first proposed in [26] and was called Morphons. The principle of the method is to match transitions (between dark and bright zones) rather than intensities, by looking locally at the spatial oscillations in intensities. This method uses Gaussian smoothing as regularization of the displacement field and additive accumulation during the iterative process. This is nevertheless not sufficient to ensure the invertibility of the deformation [22, 27].

In this paper, a Morphons registration using a diffeomorphic accumulation step is proposed and its accuracy is assessed in the case of thorax image registration, also in presence of different contrast enhancements, and compared to the Demons. The paper is organized as follows. In Section 2, the main mathematical concepts and definitions are presented. Then, in Section 3 a generic nonparametric registration process is presented, and its particularization to Morphons and to diffeomorphisms is proposed in Section 4. In Section 5, different registrations are applied on images of the thorax, without contrast enhancement in the first experiment and with contrast enhancement in the second. The results of these experiments are eventually discussed in Section 6.

2. Mathematical Framework

For the sake of clarity, let us introduce some key mathematical concepts used throughout this paper.

2.1. Images and Deformation Fields

In this paper, we always denote 3D images by lower case letters. For instance, in the process of estimating a displacement field, the fixed and the moving images are written f and m, respectively. We consider them as real valued functions on the volume R3 of points x = (x1, x2, x3), that is, f, m [set membership] F = {g : R3R : x [mapsto] g(x)}. Most of the time, these functions, and also the continuous operations performed on them, such as convolutions or integrals, must be understood as approximated on the discrete voxel grid G = {(x1, x2, x3) [set membership] Z3}, omitting the treatment of volume boundaries. In this study, image convolutions were performed using zero-padding outside the boundaries.

A displacement field on R3 is a vectorial field D [set membership] V = {V : R3R3, x [mapsto] V(x)}. It is associated to the “deformation” operation Δ [triangle, equals] Id + D, that is, Δ(x) [triangle, equals] x + D(x), with Id the identity deformation: Id(x) [triangle, equals] x. The operation Δ, and by extension its vector field D, is said to be diffeomorphic if it is invertible, differentiable, and its inverse is differentiable. For the transformation Δ to be invertible, its Jacobian must not vanish in any point x, that is, if det (J)(x) ≠ 0  for  all  x, with Jij = [partial differential]Δi/[partial differential]xj. Moreover, it has to be positive (det (J)(x) > 0). Indeed, a transformation Δ with negative Jacobians does not correspond to physical deformations (as the mirror operation).

Mathematically, given the images f and m, we will see that our global objective of our study is to estimate D such that the warping of m by D is “close” to f, that is, f[similar, equals]m[composite function (small circle)]Δ with [composite function (small circle)] the common function composition. We will use sometimes the notation

mD  =  mΔ,
(1)

to insist on the warping action of D on m. By extension, this warping symbol can also be used on vector fields themselves, for example, for two displacement fields D1 and D2, D1[open diamond]D2 = D1[composite function (small circle)]Δ2.

In practice, the warping is applied on discrete images. The transformation might therefore need to be truncated (on the volume boundaries) to the closest point inside the volume in order to avoid extrapolation of the images to be warped.

2.2. Compositive Accumulation

In this paper, we promote a particular way to combine, or accumulate, properly two displacement fields D1 and D2. Adding them to form D1 + D2 (as performed by many nonparametric registration methods; see Section 3) is of course computationally efficient, but it breaks the consistency with the composition of the corresponding spatial transformations, as illustrated in Figure 1.

Figure 1
Comparison between additive and compositive field accumulations. Warping is implemented using linear interpolation. In (a) and (b), two different displacement fields are defined on the plane (for visual clarity). In (c), the field D1 warped by D2, that ...

Indeed, one can clearly see in Figure 1(h) that the warping of the image in Figure 1(g) by the sum of two diffeomorphic fields, D1 and D2 does not correspond to the successive warping of this image by D1 and then by D2, which is represented in Figure 1(i).

However, the compositive operation, denoted by [plus sign in circle], solves this issue. It is simply defined as

D1D2Δ1Δ2Id.
(2)

By construction, the deformation operation linked to the deformation field D1 [plus sign in circle] D2 is therefore Δ1[composite function (small circle)]Δ2. If both displacement fields are diffeomorphic, their composition is also diffeomorphic [28].

The operation [plus sign in circle] has some interesting and useful properties. First, the neutral accumulation is of course obtained with the null displacement field, that is, D [plus sign in circle] 0 = 0 [plus sign in circle] D = D. Second, it is easy to prove the associative relations (D1 [plus sign in circle] D2) [plus sign in circle] D3 = D1 [plus sign in circle] (D2 [plus sign in circle] D3) = D1 [plus sign in circle] D2 [plus sign in circle] D3 for three displacement fields D1, D2, and D3. And finally, [plus sign in circle] and [open diamond] are linked through the simple relation

D1D2=D2+D1D2,
(3)

meaning that the displacement field D1 [plus sign in circle] D2 is equivalent to summing the field D2 with the field D1 warped by D2. This is illustrated in Figure 1: the vectors in Figure 1(i), corresponding to the successive warping by D1 and then D2, are the sum of the vectors in Figures 1(a) and 1(c), as shown in Figure 1(f).

2.3. Diffeomorphic Flow and Exponentiation

An important notion used in Section 4.2 is the concept of (continuous) diffeomorphic flow [27, 29, 30]. Given a point x [set membership] R3 and a smooth vector field D [set membership] V, the flow [var phi]D(x, t) is the dynamic solution u(t) [set membership] R3 of the following (autonomous) ordinary differential equation:

ddtu(t)=D(u),u(0)=x.
(4)

At a given “time” t > 0, the position [var phi]D(x, t) is simply a point on the trajectory following D tangentially from the initialization on x (see Figure 2). Following [27], the exponential of a vector field D, that is, exp (D) [set membership] V, is the nonlinear deformation operation obtained by the flow of D at time t = 1, that is, exp (D)(x) = [var phi]D(x, 1). Interestingly, this exponential map acts as the common scalar-valued exponential, that is, exp (αD)[composite function (small circle)]exp (βD) = exp ((α + β)D) for α, β [set membership] R, and it is invertible by simply considering the inverted vector field, that is, exp (−D)[composite function (small circle)]exp (D) = Id. In addition, for differentiable D, exp (D) is also a diffeomorphism on R3. In other words, exp (D) modifies the 3D coordinates with no intersection between the motions of points. Indeed, such a possibility would induce a point x with two different motion vectors, a situation that is forbidden by (4) since D(x) is uniquely defined.

Figure 2
The diffeomorphic flow exp (D) associated to the vector field D is the solution at time t = 1 on the trajectory tangent to D at each point (here represented in 2D). We see that the motion of x induced by exp (D)(x) is more compatible with ...

2.4. Scaling and Squaring

A numerical scheme exists to compute approximately but efficiently exp (D)(x) when x belongs to a regular grid of voxels G. Indeed, when the field D is close enough to zero (i.e., Δ ≈ Id), the exponential of the field can be approximated using the first-order Taylor expansion exp (D) ≈ Id + D = Δ, that is, by the transformation itself. On the other hand, the solution of the flow equation (4) in t = 1 can be approximated by “discretizing” t between 0 and 1. Indeed, as exp (D) = exp (2kD)2k (where the exponent 2k expresses the number of times the deformation operation is combined with itself), one can use the scaling and squaring strategy for computing the exponential [31]. If one chooses k such that the field 2kD is close enough to zero, the first-order approximation can be used to estimate exp (2kD) (based on the Padé approximant near the origin). Then, the solution of the flow equation is computed by performing k recursive compositions of the field by itself, given that such compositions are computationally affordable. Notice that taking k = 0 is equivalent to the simple first-order approximation. The scaling and squaring steps for field exponentiation [22] are depicted hereafter.

  1. Scaling. Divide D by a factor 2k such that 2kD is small enough, for example, when ||2kD|| = maxx||2kD(x)|| < 0.5 voxels.
  2. Exponentiation. Compute first-order explicit integration of the flow: Δ(0)(x) = [var phi]D(x, 2k) ≈ Id(x) + 2kD(x).
  3. Squaring. Perform k recursive squarings (using field composition) of the flow at time 2k in order to obtain the flow at time 1, which is the field exponential. In other words, starting with Δ* = Δ(0), do k times the computation Δ* ← Δ*[composite function (small circle)]Δ*, in order to get Δ*[similar, equals]exp (D).

We see that using this method, only k compositions (and therefore k interpolations) are needed for estimating the exponential. Compared to standard estimation of the flow over a regular discretization of the time interval [0,1] in 2k steps, the scaling and squaring method limits the numerical errors due to composition of vector fields, but it does not decrease the amplification of the error due to the field estimation at time t = 2k.

3. Generic Registration Pipeline

Nonrigid registration methods can be divided into parametric and nonparametric methods. Parametric (or model-based) methods aim at calculating the parameters of a deformation model in a high-dimensional space in order to optimize a global objective function that takes into account image similarity and transformation regularity [10]. In this case, the a priori information is included in the modelization and regularity criteria of the nonrigid transformation. For example, the harmonic energy of transformation can be explicitely included in the objective function [32].

On the other hand, nonparametric methods make it possible to decouple similarity optimization from regularization by directly acting on the displacement field. The a priori information has then to be included in the optimization process by using proper regularization techniques. Decoupled optimization makes the registration computationally efficient [8], mainly because the computation of each displacement vector is independent from others, but it prevents us from easily including more complex regularization constraints in the process, for example, such as in volume preserving registrations [33, 34].

3.1. Multiscale Nonparametric Registration

Most nonparametric registrations are based on an iterative process which is composed of 3 steps: (i) field computation, (ii) field accumulation, and (iii) field regularization. The idea is to progressively build a proper displacement field by iteratively improving the matching between the fixed image and the moving image warped by this displacement field, according to a certain metric. Note that, depending on the nature of the displacement one tries to model, the regularization is applied either on the increment field or on the accumulated field. Regularizing the field increment corresponds to a viscous fluid modeling, while regularizing the global transformation corresponds to an elastic solid modeling [14]. Only the second is considered in this study.

In this paper, our general nonparametric registration framework (e.g., valid for Demons and Morphons) adopts a multiscale approach; that is, the displacement field estimation is stabilized by decomposing the fixed and the moving images in several scales, for example, using a simple smoothing and downsampling procedure [35].

The three steps mentioned above are then applied a certain number of times (until the algorithm reaches a certain stopping criterion) to each scale separately from coarse to fine scales (Figure 3). The general explanation of these three basic blocks is given hereafter. The way they are iteratively applied at each scale is described in Section 3.2.

Figure 3
The nonparametric registration pipeline is composed of 3 main operations (Θ, Φ, and Ψ) and the warping of the moving image. Those operations are performed from coarse to fine scales. At each scale, the process is applied iteratively, ...

3.1.1. Field Computation

At each iteration of the registration process, an update displacement field (Du) is first computed as a function (Θ) of the fixed image (f) and the moving image (m) warped by the displacement field resulting from previous iterations (Da):

DuΘ(f,mΔa),
(5)

where Δa and Δu denote the deformation operations linked to Da and Du, respectively.

Depending on the nature of the images to be registered, this local displacement estimation can be based on different local image metrics, such as SSD [17], mutual information computed on blocks of voxels [8, 36], and local phase [26].

3.1.2. Field Accumulation

After the field computation, the total displacement Da must be increased by the update field

DaΦ(Da,Du).
(6)

This accumulation operation Φ is sometimes implemented as a simple addition of accumulated and update fields (as in [18, 37, 38]). However, as explained in Section 2.2, this accumulation is perhaps computationally efficient but is not consistent with the composition of the corresponding spatial transformations. The solution is, therefore, to replace it by the compositive accumulation [plus sign in circle] introduced earlier. The accumulation Da [plus sign in circle] Du of the displacement fields Da with Du is then compatible with the way Du is estimated. Indeed, since Du is computed from m[composite function (small circle)]Δa, the accumulation of Da and Du must modify Da by Du, a process intrinsically integrated by the operation [plus sign in circle]. Moreover, the associativity of [plus sign in circle] validates the compositive accumulation of displacement fields over several iterations, as illustrated in Figure 3.

3.1.3. Field Regularization

Eventually, the field is regularized in order to get a smoother transformation and reduce the impact of image noise on the registration output:

DaΨ(Da).
(7)

This operation Ψ is achieved by applying a low-pass filter on each component of the displacement field. We assume it to be a Gaussian smoothing with a size σΨ2 of a few voxels, which tends to reduce the harmonic energy of the transformation [32].

It is always possible to produce invertible fields by performing a very strong Gaussian smoothing. This, however, may reduce significantly the accuracy of the estimated displacement by limiting the solution to excessively smooth displacement fields. On the other hand, by preventing the displacement field from being noninvertible, the diffeomorphic accumulation acts in some way as a regularization, allowing the estimation of invertible fields while performing only moderate smoothing.

3.2. Registration Algorithm

Let us explain now the whole multiscale nonparametric registration algorithm relying on the three specific procedures {Θ, Φ, Ψ} defined in Section 3.1.

The algorithm takes as inputs the fixed and the moving images f and m, some parameters described below, and outputs the estimated transformation Δa = Id + Da such that f[similar, equals]m[composite function (small circle)]Δa. The whole procedure described in Table 1 and depicted in Figure 3 involves computations on different scales j [set membership] [0, J], from coarse (j = J) to fine (j = 0). Each scale is associated to a subsampled grid of voxels Gj = κjG, where κ is the subsampling factor (e.g., κ=2 in this study) between scale j and scale j + 1. The functions f and g, defined on the initial grid G = G0 = {(x1, x2, x3) [set membership] Z3}, are downsampled (after antialiasing smoothing) at any scale j by the operation Downj(). An upsampling operator Up(), implemented as a simple linear interpolation, is used to transfer any displacement field defined on a grid Gj+1 to the finer grid Gj using κ as upsampling factor. For each scale j [set membership] [0, J], the accumulated displacement field is iteratively updated until one reaches a particular stopping criterion S (e.g., based on the convergence of Da or on the SSD, as precised in Section 5).

Table 1
Multiscale nonparametric procedure.

4. Diffeomorphic Morphons

Our paper adapts the global registration method explained in the previous section to Morphons [26, 39] by taking care of the invertibility of the accumulated displacement field, that is, by introducing diffeomorphic field accumulations.

As already mentioned above, the particularity of Morphons, compared to other nonparametric methods, is that the field computation (function Θ in (5)) is based on the local phase rather than intensity difference. In other words, knowing the phase difference between periodic signals of the same frequency allows the estimation of the spatial shift between them. Therefore, under the assumption that images can locally be considered as a sum of periodic signals, the computation of the local phase difference is equivalent to the estimation of the local displacement between images. This procedure is stabilized by the multiscale approach described in Section 3.2. Besides, Morphons combines the estimation of displacement vectors with a measure of the confidence we have in these estimations, resulting in a certainty map. Therefore, for Morphons, given two images f and w = m[composite function (small circle)]Δ, the displacement field estimation Θ is actually split into two quantities

DuΘD(f,w),cuΘc(f,w),
(8)

that is, respectively, an update of the displacement field along with an update of the certainty map. A similar split is also performed on subsequent operations Φ and Ψ.

Here are the details about the three steps {Θ, Φ, Ψ} of the pipeline of Section 3 for this specific registration, including our contribution to the field accumulation step.

4.1. Displacement Field Calculation

In Morphons, a displacement field is estimated thanks to the dephasing between the local phases of the fixed and the moving images. This local phase can be probed at a certain frequency and in a particular direction using quadrature filters [40]. More precisely, Morphons method uses a quadrature filter hη of direction η [set membership] R3 (also called loglets [40]) defined in frequency by the polar separable function

Hη(ω)=χ+(ηTω)(ηTω^)2R(||ω||),
(9)

where ω [set membership] R3 is the frequency vector, χ+(λ) = 1 if λ > 0 and 0 else, ||ω||2 = ωTω, ω^=ω/ω is the unit vector supporting ω, and R is a radial function centered on ρ > 0 and defined as R(r) = exp   [−ln 2(r/ρ)/ln 2] for r > 0.

Since their support corresponds to the half volume {ω [set membership] R3 : ηTω > 0} and since (ηTω^)2=cos2ϕ (with ϕ the angle separating ω and η), loglets can be seen as the analytic counterparts of the steerable filters introduced by Freeman and Adelson [41]. As a matter of fact, only a limited number of orientations η are necessary to cover the whole frequency plane. Typically, in 2D, these directions are taken as ηk = (cos ϕk, sin ϕk) with ϕk = /4 for 0 ≤ k ≤ 3, and in 3D, η is taken as the 6 normal vectors {ηk : 0 ≤ k ≤ 5} to the faces of a hemi-icosahedron [42, 43]. Notice also that each filter hk(x) in the spatial domain is centered around the origin with a typical width given by 1/ρ.

Morphons take advantage of the following behavior. Given an image f, defining the filtering

qf(x;k)  =  (fhk)(x),
(10)

with * the common convolution operation and the shorthand hk = hηk, we can write qf(x; k) = Af(x; k)ef(x;k) since qf [set membership] C. Therefore, by processing the warped image w similarly, the local phase difference can be computed as

Δϕk(x)=arg(qf(x;k)qw(x;k)),
(11)

with (·)* the complex conjugation and Δϕk(x) = ϕf(x; k) − ϕw(x; k) the local dephasing between f and w in direction ηk.

An important observation is that the nonnegative value

Af(x;k)=|(fhk)(x)|=|3f(x)Txh¯k(x)dx|
(12)

represents also the correlation between f(x′) and the translated filter Txh¯k(x)=h¯k(x-x); that is, the filter h¯k(x)=hk(-x) translated on x. If the image f was perfectly represented by the latter, that is, if we had locally f(x′) = chk(x′ − x) for any x[set membership] R3 and some constant c [set membership] R, a displacement of f by a displacement field D(x) approximately constant over the support of Txh¯k would induce a dephasing Δϕk(x) = ρηkTD(x) since the frequency vector of Txh¯k is −ρηk. An important implicit assumption is nevertheless that |ρηkTD(x)| < π since the dephasing is known up to modulo 2π. Moreover, only ηkTD and not D can be determined, as another manifestation of the blank wall problem [44].

In practice, for most of x, f(x) is not perfectly represented by one filter but by a linear combination of them where the amplitude Af(x; k) measures the adequacy of the fit between f(x) and Txh¯k. Consequently, the local update displacement Du(x) linking f(x) and w(x) = f(x + Du(x)) in each x [set membership] R3 is estimated by solving the weighted least square optimization

ΘD(f,w)=argmind3k[ck(ρηkTdΔϕk)]2,
(13)

where the ck(x) = Af(x; k)Am(x; k) are the certainty map of the filter hk. As explained above, ck reflects for each voxel how reliable the field estimation is; that is, how contrasted the bandpass-filtered images are.

Numerically, the optimization in (13) is a standard weighted least square minimization; that is, it corresponds the minimization of the energy E(d) = ||C(Nd − Γ)||22, using the diagonal matrix C = diag (c1,…, c6), the matrix N = (η1,…, η6)T, and the vector Γ = (Δϕ1,…, Δϕ6)T. An easy computation shows that the solution of (13) is then given by the Moore-Penrose pseudoinverse (CN) = (NTC2N)−1NTC, that is,

ΘD(f,w)=(CN)CΓ,
(14)

with ΘD(f, w) arbitrary set to 0 when (NTC2N) is not invertible.

Jointly to the estimation (13), a global certainty map associated to the quality of the estimation of ΘD is defined as [43]

Θc(f,w)=kck(x),
(15)

that is, the sum of all certainty measures for each quadrature filter. This update of the certainty map must then be combined with an accumulated certainty computed from previous iterations (see Section 4.2).

In the multiscale approach described in Section 3.2, using the same quadrature filters at decreasing scales Hηk is equivalent to estimating the phase of the bandpass-filtered image around increasing cutoff frequencies, that is, with ρ ← 2ρ each time jj + 1. This sustains the coarse-to-fine displacement estimation, that is, the computation of ΘD and Θc on different scale bands fj and mj of f and m.

Convolutions with quadrature filters can be implemented efficiently in the Fourier domain thanks to the FFT and the convolution theorem. However, since the spatial extent of those filters is small, it is also possible to use efficient spatial convolutions with truncated kernels, as done in this study. As the local phase is invariant to local intensity scaling, the Morphons procedure is suitable for registering images with various contrast enhancements. Besides, some studies indicate that the phase extraction allows a fast convergence and a subvoxel precision in displacement estimation (e.g., see [39]).

4.2. Field Accumulation

In the original Morphons method, the accumulated field is computed as a weighted sum of the update field and the previous accumulated field, as used in damped optimization schemes. The weights are given by the certainty on the update field (cu, as computed from Θc) and the accumulated certainty map (ca). As the certainty map must also be accumulated in order to reflect the confidence in all previous displacement computations, the accumulation step Φ must be divided into two operations ΦD (field accumulation) and Φc (certainty accumulation):

ΦD(Da,Du,ca,cu)=Da+cuca+cuDu,
(16)

Φc(ca,cu)=ca2+cu2ca+cu,
(17)

where in the last formula, similar to the field accumulation, the certainty map is updated by its own certainty [43].

However, as it was explained before, the addition of displacement fields is not really appropriate for accumulating spatial transformations, in contrast to composition. The compositive accumulation may also be damped using the certainty as a weighting factor

ΦD(Da,Du,ca,cu)=Dacuca+cuDu.
(18)

The (SSD-based) Demons registration is a nonparametric algorithm which performs the optimization of the SSD between images. In [27], a diffeomorphic field accumulation is proposed as improvement of the Demons method. The idea is to use an adaptation of the optimization method to Lie groups [45] in order to limit the possible solutions to diffeomorphic transformations. In practice, this is done by replacing the accumulation step of the Demons by an accumulation using the diffeomorphic flow exp () introduced in Section 2. This accumulation reads then

ΦD(Da,Du)=Da(exp(Du)Id),
(19)

where the field exponential exp (Du) can be efficiently estimated using a small number of recursive compositions of the field Du by itself. Consequently, the displacement field ΦD(Da, Du) is linked to the deformation operation Δa[composite function (small circle)]exp (Du).

In the case of the Morphons, the accumulation step can be achieved in the same way. This will produce smoother fields than the traditional addition or composition. However, the accumulation step in the Morphons method involves a damping based on the certainty. Therefore, we propose the following accumulation step for diffeomorphic Morphons:

ΦD(Da,Du,ca,cu)=Da(exp(cuca+cuDu)Id).
(20)

Since exp (0  D) = Id for any vector field D, the accumulation fades away when ca [dbl greater-than sign] cu. The accumulation of the certainty map remains as explained previously in (17).

Notice that, because the field is discretized on a grid of voxels, interpolation is needed for computing the composition of two diffeomorphisms. Therefore, errors due to successive interpolations could potentially lead to noninvertible transformations. However, such problems were not observed in practical experiments using reasonable smoothing of the field.

4.3. Field Regularization

During the displacement estimation step, the relevance of local phase computation is estimated and used as weight for the accumulation. This certainty map may also be used for a smart regularization of the displacement field. Regularization is performed using a normalized convolution [46] of the field by a Gaussian kernel, taking into account the certainty map in order to put greater importance to high certainty locations. The certainty is also regularized in the same way as displacement field components in order to preserve the correspondence between the displacement vectors and their corresponding certainty.

Mathematically, given a positive function h and a filter g (typically a Gaussian kernel of variance σΨ > 0), the normalized convolution of a (scalar) function s by g as involved by the normalization h is

shg(hs)ghg.
(21)

This operation does not increase the maximum amplitude of the filtered function. Indeed, for a nonnegative kernel g, we show easily that ||s*hg|| ≤ ||s||, with ||s|| = maxx |s(x)|. The accumulated displacement field Da and subsequently the certainty map are therefore regularized thanks to this operation using for normalization the certainty map ca, that is,

ΨD(Da,ca)=Dacag,Ψc(ca,ca)=cacag.
(22)

Notice that, for computing ΨD, the normalized convolution is performed separately on all components of the vector field.

This operation tends to propagate the displacement field from high certainty areas to areas which show less significant filter responses. Besides, by setting to zero the certainty outside the volume boundaries, normalized convolution cancels the influence of the padding strategy. This step produces a smooth version of the accumulated field that may reduce the accuracy of image matching resulting from the displacement estimation step, as it limits the possible solutions to smooth displacement fields.

However, if the iterative algorithm is to converge, the solution will be regular and invertible (except for large numerical errors), thanks to accumulation and regularization constraints, but it will also be (at least locally) optimal in terms of local phase difference. Indeed, as the phase is monotonic and smooth, a mismatch between local structures will automatically lead to nonzero field update with a high certainty value, which will tend to improve the displacement estimate and fit the structures together.

The Jacobian of the displacement field may be used as a criterion for validating the physical behavior of the deformation. Indeed, the Jacobian gives for each voxel the change in volume this voxel encounters during deformation. Jacobian indicates expansion when it is greater than 1, and compression when it is smaller than 1. A negative Jacobian means that the voxel is “inverted” (getting a negative volume), which is incompatible with the mass-preservation principle.

In the following, the diffeomorphic version of Demons and Morphons are denoted respectively D-Demons and D-Morphons.

5. Experiments and Results

The methods were first compared for several simple 2D virtual situations in order to demonstrate the interest in chosing the accurate registration method with respect to the images to be registered.

For the clinical validation, Morphons and D-Morphons registrations were first validated on a 10-phase point-validated pixel-based breathing thorax model (POPI-model) from the Léon Bérard Cancer Center, Lyon, France [4], in order to compare the D-Morphons to Morphons, Demons, and D-Demons in the case of intensity conservation between images. Then, it was applied to lung images with different contrast enhancements, in order to illustrate the benefit of a phase-based approach compared to traditional SSD-based registration methods in the case where intensities are not conserved between the images to be registered.

All simulations were performed using Linux, on a single processor Intel Core 2 (2.4 GHz). Our MATLAB implementation used for the prototyping of the methods was also used for simulation. Notice that no efforts were made for achieving good performances in terms of computational cost and memory requirements in the implementations used in this study. The local phase estimation was performed using convolutions with 9 × 9 × 9 quadrature filters. Less than 1 GB of RAM was required for registering two volumes of 256 × 256 × 100 voxels using all registrations. The time required for registering such images, using the parameters presented hereafter, was around 6 minutes for Demons, 42 for Morphons, 7 for D-Demons and 43 for D-Morphons. However, preliminary results based on a C++ implementation of the Morphons, which uses operations in the Fourier domain instead of convolutions (as done in our matlab implementation) and using 4 threads on a quad-core CPU, allowed a division of the computation time by 50, leading to Morphons registrations taking about one minute for such a typical image size.

5.1. Illustrative Virtual Experiments

Two 2D virtual experiments were performed. The first experiment, illustrated in Figure 4, is based on a virtual disk image after blurring. Two images of the same disk were created, the only difference being the scale of intensities (multiplication by 0.75). This experiment shows the interest in using a phase-based method (conventional Morphons in this example) while registering identical shapes with different contrasts, compared to an intensity-based method (conventional Demons).

Figure 4
Results of the registration between 2 identical-sized blurred disks with different contrasts, using Demons and Morphons. In yellow: the contour of the disk. In red: the vector field resulting from the registration. The displacement field resulting from ...

The second virtual experiment is based on two images of a disk (see Figure 5). In the fixed image, a disk of radius r1 + r2 was created, and a hole (disk of radius r2) was added in its center. In the moving image, a disk of radius r1 was created with the same intensity scaling as in the fixed image. This example illustrates the case where a structure is missing in one image compared to the other, as it may occur in practice (e.g., the problem of bowel gas in CT images of the abdomen). This experiment illustrates how the diffeomorphic version of the Morphons algorithm can prevent from producing negative volumes after registration, without increasing the smoothing by using a larger Gaussian regularization kernel.

Figure 5
Results of the registration between 2 images using Morphons and D-Morphons registrations, illustrating the case where a structure (i.e., the bright hole at the center of the fixed image) is missing in the moving image. Both methods lead to deformed images ...

5.2. Accuracy Assessment on a Breathing Thorax Model

The POPI model [4] is composed of 10 volumes reconstructed from a 4D respiration-correlated CT scan (RCCT) of the thorax, each volume corresponding to a particular phase of an average breathing cycle. 41 landmarks were identified by medical experts in each of the 10 images for registration validation.

Conventional Morphons, D-Demons, and D-Morphons were applied between a reference phase and the 9 others. For all methods, the number of scales was set to J = 8, with final resolution of 2 mm × 2 mm × 2 mm. The variance of the Gaussian kernel used for regularization was empirically set to twice the voxel size (σΨ2 = 2 voxels). For this experiment, a minimum of 10 and a maximum of 20 iterations was used at each scale. In between, the iterative process was stopped if the changes, measured in terms of SSD, were inferior to 0.01%. Such a convergence criterion was usually reached before the 20th iteration, supporting the fact that both Demons and Morphons behave like optimization methods.

The results were then compared with each other and with the results from a conventional Demons algorithm as used in [4]. The comparisons were achieved in terms of error in landmark position, SSD between images, harmonic energy, and minimum Jacobian.

  1. The landmark position error evaluates the ability of the registration in finding the physical motion of organs.
  2. The SSD between fixed and deformed images is a measure of the image matching according to the assumption of intensity conservation. It is computed as ∑x(fm[composite function (small circle)]Δ)2.
  3. The harmonic energy [29, 32] of the displacement field D indicates how regular the field is and is computed as (1/2)∑x(||[nabla]D1||2 + ||[nabla]D2||2 + ||[nabla]D3||2).
  4. The Jacobian of the field indicates the volume change of each voxel. Recall that negative values of the Jacobian correspond to inverted volumes, which is not acceptable in a physical point of view. The Jacobian is computed as det (J), with Jij = [partial differential]Δi/[partial differential]xj = δij + [partial differential]di/[partial differential]xj, where δij is Kronecker's delta (δij = 1 if i = j, 0 else) and di is the ith component of the displacement field. In practice, the partial derivatives [partial differential]di/[partial differential]xj can be computed using centered finite difference approximations.

The comparisons of landmark position errors (expressed in mm) resulting from the different registrations can be seen in Table 2 with, from left to right, the error in landmark position (norm of the difference) before registration, using Demons (values from the POPI website), Morphons, D-Demons, and D-Morphons. Position errors are noted as follows: mean/std (max). On average, for Morphons, D-Demons, and D-Morphons, the error in landmark position was equal or inferior to 1 mm, which is half the size of the voxels at the finest scale of the registration process.

Table 2
Results for the POPI experiment: error in landmark position.

Results showed that all registrations greatly improved the matching of intensities. The SSD between fixed and deformed image was similar for Morphons, D-Demons, and D-Morphons (see Figure 6). The harmonic energy of the fields resulting from these registrations was also comparable (see Figure 6).

Figure 6
Results for the 9 registered phases of the POPI model. (a) Boxplots of the SSD before registration (in yellow) and after all 4 registrations. (b) Boxplots of the energy of deformation after all 4 registrations. (c) Boxplots of the minimum Jacobian after ...

The matching and the harmonic energy obtained by Demons (as presented by the authors of [4] on the POPI website) was slightly less good than for the 3 other methods. However, this is most likely due to the parameters used for registration (e.g., the number of scales, the variance for smoothing, etc.). In particular, for very similar images (first 2 phases of the RCCT), the algorithm was not able to find a smooth displacement field that reduced the SSD.

The minimum Jacobian of the displacement fields resulting from conventional methods gets down to −0.5 for both Demons and Morphons (see Figure 6), as, respectively, 67 and 460 voxels were inverted for the corresponding phase when applying the field on the moving image (which is composed of almost 6 mega voxels). However, when using diffeomorphic accumulation, the minimum Jacobian was raised to 0.2 for the Demons and 0.1 for the Morphons, showing that the diffeomorphic accumulation step prevented the field from inverting voxels.

5.3. Application to Images of the Thorax with and without Lodine Contrast Agent

The breathing-correlated motion of tumor is a typical feature of lung cancer that has to be dealt with in radiotherapy planning. RCCT images provide information about the tumor motion throughout the breathing cycle. From the different respiratory phases, an adequate margin around the tumor (the ITV, i.e., the Internal Target Volume) can be estimated, integrating thus all tumor positions through the respiratory cycle [20].

However, the lack of contrast enhancement, as well as the high noise level and the presence of artifacts that characterize 4D RCCT, may significantly impair the accurate delineation of the target volumes on these images. More particularly, the iodine contrast agent is of prime importance to help at differentiating tumor extents from vascular structures in the centrally located lung tumors. In this context, the acquisition of a conventional contrast-enhanced CT (CE-CT) acquired during free breathing should be considered for the delineation task, while the 4D RCCT is used to estimate the motion range of the tumor during breathing. To automatize this process, the delineated tumor volume at the CE-CT can be deformed on the various respiratory phase images from the 4D RCCT using nonrigid registration to finally get the ITV, as illustrated on Figure 7.

Figure 7
Schematic representation of the ITV creation (with only 4 phases). The CTV delineated on a reference image with contrast enhancement (on the left) is deformed towards every phases (middle) using displacement fields estimated by registration, and their ...

The purpose of this experiment is to compare Demons and Morphons algorithms (conventional and diffeomorphic versions) for the registration between images with and without contrast enhancement, while keeping the same setting as for the POPI experiment.

A CE-CT scan of 3 lung cancer patients was acquired as well as a 4D RCCT scan at another time point. The first CT scan was taken in free breathing using an iodine contrast agent. The 4D RCCT scan was acquired without any contrast agent and was reconstructed into 10 phases. Histogram equalization was not able to correct for localized contrast differences between the CE-CT and RCCT phase images. For all 3 patients, Demons, Morphons, D-Demons, and D-Morphons were applied between each of the 10 RCCT images and the CE-CT, with the same registration parameters as for the POPI simulation.

The displacement fields resulting from these registrations were compared in terms of harmonic energy and minimum Jacobian (see Figure 8). The resulting images were compared in terms of SSD and mutual information.

Figure 8
Results for the variable contrast experiment on 30 phases (3 patients with 10 phases each). (a) Boxplots of the energy of deformation after all 4 registrations. (b) Boxplots of the minimum Jacobian after all 4 registrations. From (a) to (b), these registrations ...

The harmonic energy of displacement fields resulting from Demons and D-Demons was quite higher than that with the Morphons and D-Morphons, and the minimum Jacobian of the displacement fields was positive only for registrations using the diffeomorphic accumulation. In the worst case, 7455 and 1114 voxels were inverted using respectively Demons and Morphons without diffeomorphic accumulation (on an image of 5 mega voxels). An example of area leading to bad transformations (with negative Jacobians) using conventional methods is depicted in Figure 9. D-Morphons lead to the smoothest transformation, with minimum Jacobian values around 0.2. These quite low values, however, were very sporadic within the image volume.

Figure 9
Illustration of negative Jacobians resulting from nondiffeomorphic registrations. Left: moving and fixed images. Right: fields resulting from registrations (red arrows) and their Jacobian (grayscale images). The negative Jacobians regions are contoured ...

We noticed that, unlike the results obtained with the POPI simulation, the SSD resulting from the Morphons and D-Morphons was a bit higher than the SSD resulting from Demons and D-Demons. However, as illustrated in the example of Figure 4, the SSD does not reflect the matching in variable contrast areas. On the other hand, no significant differences in terms of mutual information were observed between images resulting from the different registrations. This is likely due to the very low contrasts in the noncontrasted images within the regions corresponding to contrast-enhanced tissues in the other image whereas the main differences in terms of displacement field were located in these regions, as illustrated in Figure 10.

Figure 10
Illustration of the results for the variable contrast experiment. (a) and (b): Fixed image (left) and moving image (right). (c) and (d): Deformed image with deformed contours and displacement field resulting from D-Demons (left) and from D-Morphons (right). ...

In order to illustrate the effect of the registration on contrast-enhanced tissues, one phase of the RCCT scan of one of the 3 patients was chosen as example. For this patient, the tumor was located close to contrasted tissues. The tumor and the blood vessels were delineated by a physician, on the contrast-enhanced scan and on one phase of the RCCT scan. The delineations on the phase image were deformed according to the fields resulting from the different registrations. The results are illustrated in Figure 10.

The change in volume due to warping was computed, as well as the harmonic energy inside the delineated stuctures and the difference between the center of mass of the tumor with and without registration.

The change in volume was very small when using a phase-based field computation for both the vessels (around 0%) and for the tumor delineations (around 1%), while it rose up to 23% for the vessels and to 6% for the tumor while using the Demons. In the same way, the harmonic energy and the error on the center of mass of the tumor were much smaller for the phase-based registration methods. These results are summarized in Table 3. One can notice that the diffeomorphic accumulation of the field in the Morphons did not change the results in terms of harmonic energy and volume changes compared to conventional Morphons. This is due to the fact that the displacement of the considered organs is small and smooth.

Table 3
Comparison of volume change, harmonic energy, and errors in center of mass of the delineations of the vessels and tumor on a single phase.

6. Discussion

The first medical experiment showed that D-Morphons and D-Demons lead to similar matching of both image intensities and anatomical landmarks. This shows that for monomodal registration of lung CT scans, the phase difference has an efficiency comparable to the efficiency of the SSD metric. Furthermore, the D-Morphons produced displacement fields as smooth as those obtained with D-Demons. In opposition to conventional Demons and Morphons, both diffeomporphic methods produced invertible displacement fields which are physically meaningful.

The second medical experiment illustrates the limitations in registering images with various levels of contrast enhancement with the Demons method. Indeed, the intensity matching resulting from Demons was better than that from Morphons, but the field was obviously wrong, as the Demons results in a global shrinking of the contrasted tissues (arteries) that does not reflect a proper anatomical behavior, but that is due to the fact that the Demons registration is based on the minimization of the SSD, which produces an improper displacement estimation when the intensities of identical tissues are different in the fixed and moving images. This mismatch between registered anatomical structures is clearly visible on Figure 10. As illustrated in the example of Figure 4, the field produced by Demons tries to match structures of same intensity, which do not correspond to identical anatomical structures because of the difference in contrast agent concentration. Therefore, the field resulting from Demons (see the field on the left part of Figure 10) is far less smooth than it should be and can lead to wrong deformation estimations as it illustrated in the example (see Table 3). In this case, the difference in intensity between the images with and without contrast enhancement lead to important volume changes for vessels and tumor by using Demons or D-Demons, while almost no changes in volume were observed for these tissues when using a phase-based approach. Besides, the harmonic energy inside these tissues shows that the field is much more smooth using the phase-based registration. It is important to notice that these effects are mostly limited by the regularization of the displacement field during the Demons and D-Demons registrations, and that they will still be worse if less regularization is used (smaller variance of the Gaussian kernel used for smoothing the displacement field). This is not the case for the fields produced by the Morphons and diffeomorphic Morphons, which are much smoother and preserve the anatomical topology even with contrast variations between images (see Figure 10). Notice that the reduction of the smallest segmentation that can be observed in the Morphons results is mostly due to interslices motion, as confirmed by the Jacobian close to 1 in this area that shows that there is no important volume changes within this segmented region. Finally, one can see that the invertibility of the displacement field is observed with both diffeomorphic registrations.

These results can be summarized by classifying the different registration strategies according to the smoothness (harmonic energy) and the invertibility (minimum Jacobian) of the resulting displacement fields (see Table 4) for the variable contrast experiment.

Table 4
Classification of the registration algorithms for variable contrast enhancement.

One can notice that the D-Morphons algorithm combines both advantages: the field is invertible and smooth, which suggests that it is likely a better estimation of the real transformation which is known to be smooth in this area.

7. Conclusion

The D-Morphons is a multiresolution registration algorithm which computes a diffeomorphic displacement field based on the minimization of the local intensity phase. The method managed to estimate the deformations in a breathing thorax, with an accuracy comparable to the accuracy of the D-Demons, and leads to the same requisite property of invertibility of the field. Moreover, the D-Morphons managed to accurately estimate the deformations between images with variable contrast, while the conventional SSD-based methods led to misalignment of anatomical structures affected by the contrast variation.

Acknowledgments

The authors thank the Medical Informatics Laboratory of Linköping University for sharing their implementation of the Morphons. G. Janssens and J. O. de Xivry are funded by an FRIA grant. L. Jacques is a Postdoctoral Researcher funded by the Belgian FRS-FNRS.

References

1. Maintz JBA, Viergever MA. A survey of medical image registration. Medical Image Analysis. 1998;2(1):1–36. [PubMed]
2. Wang H, Dong L, O’Daniel J, et al. Validation of an accelerated ’demons’ algorithm for deformable image registration in radiation therapy. Physics in Medicine and Biology. 2005;50(12):2887–2905. [PubMed]
3. Pevsner A, Davis B, Joshi S, et al. Evaluation of an automated deformable image matching method for quantifying lung motion in respiration-correlated CT images. Medical Physics. 2006;33(2):369–376. [PubMed]
4. Vandemeulebroucke J, Sarrut D, Clarysse P. The popi-model,a point-validated pixel-based breathing thorax model. In: Proceedings of the 15th International Conference on the Use of Computers in Radiation Therapy (ICCR ’07); 2007.
5. Rietzel E, Chen GTY. Deformable registration of 4D computed tomography data. Medical Physics. 2006;33(11):4423–4430. [PubMed]
6. Zitová B, Flusser J. Image registration methods: a survey. Image and Vision Computing. 2003;21(11):977–1000.
7. Klein A, Andersson J, Ardekani BA, et al. Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration. NeuroImage. 2009;46(3):786–802. [PMC free article] [PubMed]
8. Modat M, Vercauteren T, Ridgway GR, Hawkes DJ, Fox NC, Ourselin S. Diffeomorphic demons using normalized mutual information, evaluation on multimodal brain MR images. In: Medical Imaging 2010: Image Processing, vol. 7623; 2010; Proceedings of SPIE.
9. Holden M. A review of geometric transformations for nonrigid body registration. IEEE Transactions on Medical Imaging. 2008;27(1):111–128. [PubMed]
10. Rueckert D, Sonoda LI, Hayes C, Hill DLG, Leach MO, Hawkes DJ. Nonrigid registration using free-form deformations: application to breast mr images. IEEE Transactions on Medical Imaging. 1999;18(8):712–721. [PubMed]
11. Rohr K, Stiehl H, Sprengel R, et al. Point-based elastic registration of medical image data using approximating thin-plate splines. Visualization in Biomedical Computing. 1996;1131:297–306.
12. Fornefett M, Rohr K, Stiehl HS. Radial basis functions with compact support for elastic registration of medical images. Image and Vision Computing. 2001;19(1-2):87–96.
13. Ferrant M, Macq B, Nabavi A, Warfield SK. Deformable modeling for characterizing biomedical shape changes. Discrete Geometry for Computer Imagery. 2000:235–248.
14. Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Transactions on Image Processing. 1996;5(10):1435–1447. [PubMed]
15. Horn B, Schunck B. Determining optical dlow. MIT Tech. Rep. 1981
16. Thirion JP. Fast Non-Rigid Matching of 3D Medical Images. HAL—CCSd—CNRS, 1995.
17. Thirion J-P. Image matching as a diffusion process: an analogy with Maxwell’s demons. Medical Image Analysis. 1998;2(3):243–260. [PubMed]
18. Pennec X, Cachier P, Ayache N. Medical Image Computing and Computer-Assisted Intervention. New York, NY, USA: Springer; 1999. Understanding the “demon’s algorithm”: 3d non-rigid registration by gradient descent; pp. 597–605.
19. Rosu M, Chetty IJ, Balter JM, Kessler ML, McShan DL, Ten Haken RK. Dose reconstruction in deforming lung anatomy: dose grid size effects and clinical implications. Medical Physics. 2005;32(8):2487–2495. [PubMed]
20. Orban de Xivry J, Janssens G, Bosmans G, et al. Tumour delineation and cumulative dose computation in radiotherapy based on deformable registration of respiratory correlated CT images of lung cancer patients. Radiotherapy and Oncology. 2007;85(2):232–238. [PubMed]
21. Rohlfing T. Transformation model and constraints cause bias in statistics on deformation fields. Medical Image Computing and Computer-Assisted Intervention. 2006;9(1):207–214. [PubMed]
22. Arsigny V, Commowick O, Pennec X, Ayache N. A log-Euclidean framework for statistics on diffeomorphisms. Medical Image Computing and Computer-Assisted Intervention. 2006;9(1):924–931. [PubMed]
23. Vercauteren T, Pennec X, Perchant A, Ayache N. Non-parametric diffeomorphic image registration with the demons algorithm. Medical Image Computing and Computer-Assisted Intervention. 2007;10(2):319–326. [PubMed]
24. Rodríguez-Vila B, Pettersson J, Borga M, García-Vicente F, Gómez EJ, Knutsson H. 3D deformable registration for monitoring radiotherapy treatment in prostate cancer. In: Proceedings of the 15th Scandinavian conference on image analysis (SCIA ’07); June 2007.
25. Melbourne A, Atkinson D, Hawkes D. Influence of organ motion and contrast enhancement on image registration. Medical Image Computing and Computer-Assisted Intervention. 2008;11(2):948–955. [PubMed]
26. Knutsson H, Andersson M. Morphons: segmentation using elastic canvas and paint on priors. In: Proceedings of the IEEE International Conference on Image Processing (ICIP ’05); September 2005; pp. 1226–1229.
27. Vercauteren T, Pennec X, Perchant A, Ayache N. Diffeomorphic demons: efficient non-parametric image registration. NeuroImage. 2009;45(1):S61–S72. [PubMed]
28. Ashburner J. A fast diffeomorphic image registration algorithm. NeuroImage. 2007;38(1):95–113. [PubMed]
29. Jost J. Riemannian Geometry and Geometric Analysis. New York, NY, USA: Springer; 2005.
30. Pennec X. Probabilities and statistics on riemannian manifolds: a geometric approach. In: IEEE Workshop on Nonlinear Signal and Image Processing; 2004; pp. 194–198.
31. Higham NJ. The scaling and squaring method for the matrix exponential revisited. SIAM Journal on Matrix Analysis and Applications. 2005;26(4):1179–1193.
32. Fischer B, Modersitzki J. Fast diffusion registration. (AMS Contemporary Mathematics).Inverse Problems, Image Analysis, and Medical Imaging. 2002;313:117–129.
33. Rohlfing T, Maurer CR, Jr., Bluemke DA, Jacobs MA. Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint. IEEE Transactions on Medical Imaging. 2003;22(6):730–741. [PubMed]
34. Haber E, Modersitzki J. Image registration with guaranteed displacement regularity. International Journal of Computer Vision. 2007;71(3):361–372.
35. Mallat S. A Wavelet Tour of Signal Processing: The Sparse Way. 3rd edition. Academic Press; 2008.
36. Suárez E, Santana JA, Rovaris E, Westin C-F, Ruiz-Alzola J. Fast entropy-based nonrigid registration. In: Proceedings of the International Conference On Computer Aided Systems Theory (EUROCAST ’04), vol. 2809; 2004; pp. 607–615. Lecture Notes in Computer Science.
37. Ibanez L, Schroeder W, Ng L, Cates J. The ITK Software Guide Second Edition Updated for ITK Version 2.4. Kitware; 2003.
38. Plumat J, Andersson M, Janssens G, de Xivry JO, Knutsson H, Macq B. Image registration using the morphon algorithm: an itk implementation. Insight Journal. 2009
39. Wrangsjo A, Pettersson J, Knutsson H. Non-rigid registration using morphons. In: Proceedings of the 14th Scandinavian Conference on Image Analysis (SCIA ’05); 2005; Springer;
40. Knutsson H, Andersson M. Implications of invariance and uncertainty for local structure analysis filter sets. Signal Processing: Image Communication. 2005;20(6):569–581.
41. Freeman WT, Adelson EH. The design and use of steerable filters. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1991;13(9):891–906.
42. Knutsson H. Representing local structure using tensors. In: Proceedings of the 6th Scandinavian Conference on Image Analysis (SCIA ’89); June 1989; Oulu, Finland. pp. 244–251.
43. Bergnéhr L. Segmentation and alignment of 3-D transaxial myocardial perfusion images and automatic dopamin transporter quantification. Linköping, Sweden: Institutionen för medicinsk teknik. Linköpings Universitet; 2008. Ph.D. thesis.
44. Lucas BD, Kanade T. An iterative image registration techniquewith an application to stereo vision. In: Proceedings of the International Joint Conference on Artificial Intelligence, vol. 3; 1981; pp. 674–679.
45. Mahony R, Manton J. The geometry of the Newton method on non-compact lie groups. Journal of Global Optimization. 2002;23(3):309–327.
46. Pettersson J, Knutsson H, Borga M. Normalised convolution and morphons for non-rigid registration. In: Proceedings of the SSBA Symposium on Image Analysis; March 2006; Umea, Sweden. SSBA;

Articles from International Journal of Biomedical Imaging are provided here courtesy of Hindawi Publishing Corporation