International Journal of Biomedical Imaging

International Journal of Biomedical Imaging / 2011 / Article
Special Issue

Modern Mathematics in Biomedical Imaging

View this Special Issue

Research Article | Open Access

Volume 2011 |Article ID 891585 | 16 pages |

Diffeomorphic Registration of Images with Variable Contrast Enhancement

Academic Editor: James G. Nagy
Received30 Apr 2010
Revised28 Jul 2010
Accepted24 Sep 2010
Published08 Dec 2010


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 [2–5]. 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 [6–8]. 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 𝑓 and 𝑚, respectively. We consider them as real valued functions on the volume ℝ3 of points 𝑥=(𝑥1,𝑥2,𝑥3), that is, 𝑓,𝑚∈ℱ={𝑔∶ℝ3→ℝ∶𝑥↦𝑔(𝑥)}. 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 𝒢={(𝑥1,𝑥2,𝑥3)∈ℤ3}, omitting the treatment of volume boundaries. In this study, image convolutions were performed using zero-padding outside the boundaries.

A displacement field on ℝ3 is a vectorial field 𝐷∈𝒱={𝑉∶ℝ3→ℝ3,𝑥↦𝑉(𝑥)}. It is associated to the “deformation” operation Δ≜Id+𝐷, that is, Δ(𝑥)≜𝑥+𝐷(𝑥), with Id the identity deformation: Id(𝑥)≜𝑥. The operation Δ, and by extension its vector field 𝐷, 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 𝑥, that is, if det(𝒥)(𝑥)≠0forall𝑥, with 𝒥𝑖𝑗=𝜕Δ𝑖/𝜕𝑥𝑗. Moreover, it has to be positive (det(𝒥)(𝑥)>0). Indeed, a transformation Δ with negative Jacobians does not correspond to physical deformations (as the mirror operation).

Mathematically, given the images 𝑓 and 𝑚, we will see that our global objective of our study is to estimate 𝐷 such that the warping of 𝑚 by 𝐷 is “close” to 𝑓, that is, 𝑓≃𝑚∘Δ with ∘ the common function composition. We will use sometimes the notation 𝑚⋄𝐷=𝑚∘Δ,(1) to insist on the warping action of 𝐷 on 𝑚. By extension, this warping symbol can also be used on vector fields themselves, for example, for two displacement fields 𝐷1 and 𝐷2, 𝐷1⋄𝐷2=𝐷1∘Δ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 𝐷1 and 𝐷2. Adding them to form 𝐷1+𝐷2 (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.

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, 𝐷1 and 𝐷2 does not correspond to the successive warping of this image by 𝐷1 and then by 𝐷2, which is represented in Figure 1(i).

However, the compositive operation, denoted by ⊕, solves this issue. It is simply defined as 𝐷1⊕𝐷2≜Δ1∘Δ2−Id.(2) By construction, the deformation operation linked to the deformation field 𝐷1⊕𝐷2 is therefore Δ1∘Δ2. If both displacement fields are diffeomorphic, their composition is also diffeomorphic [28].

The operation ⊕ has some interesting and useful properties. First, the neutral accumulation is of course obtained with the null displacement field, that is, 𝐷⊕0=0⊕𝐷=𝐷. Second, it is easy to prove the associative relations (𝐷1⊕𝐷2)⊕𝐷3=𝐷1⊕(𝐷2⊕𝐷3)=𝐷1⊕𝐷2⊕𝐷3 for three displacement fields 𝐷1, 𝐷2, and 𝐷3. And finally, ⊕ and ⋄ are linked through the simple relation𝐷1⊕𝐷2=𝐷2+𝐷1⋄𝐷2,(3) meaning that the displacement field 𝐷1⊕𝐷2 is equivalent to summing the field 𝐷2 with the field 𝐷1 warped by 𝐷2. This is illustrated in Figure 1: the vectors in Figure 1(i), corresponding to the successive warping by 𝐷1 and then 𝐷2, 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 𝑥∈ℝ3 and a smooth vector field 𝐷∈𝒱, the flow 𝜑𝐷(𝑥,𝑡) is the dynamic solution 𝑢(𝑡)∈ℝ3 of the following (autonomous) ordinary differential equation:dd𝑡𝑢(𝑡)=𝐷(𝑢),𝑢(0)=𝑥.(4) At a given “time” 𝑡>0, the position 𝜑𝐷(𝑥,𝑡) is simply a point on the trajectory following 𝐷 tangentially from the initialization on 𝑥 (see Figure 2). Following [27], the exponential of a vector field 𝐷, that is, exp(𝐷)∈𝒱, is the nonlinear deformation operation obtained by the flow of 𝐷 at time 𝑡=1, that is, exp(𝐷)(𝑥)=𝜑𝐷(𝑥,1). Interestingly, this exponential map acts as the common scalar-valued exponential, that is, exp(𝛼𝐷)∘exp(𝛽𝐷)=exp((𝛼+𝛽)𝐷) for 𝛼,𝛽∈ℝ, and it is invertible by simply considering the inverted vector field, that is, exp(−𝐷)∘exp(𝐷)=Id. In addition, for differentiable 𝐷, exp(𝐷) is also a diffeomorphism on ℝ3. In other words, exp(𝐷) modifies the 3D coordinates with no intersection between the motions of points. Indeed, such a possibility would induce a point 𝑥 with two different motion vectors, a situation that is forbidden by (4) since 𝐷(𝑥) is uniquely defined.

2.4. Scaling and Squaring

A numerical scheme exists to compute approximately but efficiently exp(𝐷)(𝑥) when 𝑥 belongs to a regular grid of voxels 𝒢. Indeed, when the field 𝐷 is close enough to zero (i.e., Δ≈Id), the exponential of the field can be approximated using the first-order Taylor expansion exp(𝐷)≈Id+𝐷=Δ, that is, by the transformation itself. On the other hand, the solution of the flow equation (4) in 𝑡=1 can be approximated by “discretizing” 𝑡 between 0 and 1. Indeed, as exp(𝐷)=exp(2−𝑘𝐷)2𝑘 (where the exponent 2𝑘 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 𝑘 such that the field 2−𝑘𝐷 is close enough to zero, the first-order approximation can be used to estimate exp(2−𝑘𝐷) (based on the Padé approximant near the origin). Then, the solution of the flow equation is computed by performing 𝑘 recursive compositions of the field by itself, given that such compositions are computationally affordable. Notice that taking 𝑘=0 is equivalent to the simple first-order approximation. The scaling and squaring steps for field exponentiation [22] are depicted hereafter.(i)Scaling. Divide 𝐷 by a factor 2𝑘 such that 2−𝑘𝐷 is small enough, for example, when ‖2âˆ’ğ‘˜ğ·â€–âˆž = max𝑥‖2−𝑘𝐷(𝑥)‖<0.5 voxels.(ii)Exponentiation. Compute first-order explicit integration of the flow: Δ(0)(𝑥)=𝜑𝐷(𝑥,2−𝑘)≈Id(𝑥)+2−𝑘𝐷(𝑥).(iii)Squaring. Perform 𝑘 recursive squarings (using field composition) of the flow at time 2−𝑘 in order to obtain the flow at time 1, which is the field exponential. In other words, starting with Δ∗=Δ(0), do 𝑘 times the computation Δ∗←Δ∗∘Δ∗, in order to get Δ∗≃exp(𝐷).

We see that using this method, only 𝑘 compositions (and therefore 𝑘 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 2𝑘 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 𝑡=2−𝑘.

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.

3.1.1. Field Computation

At each iteration of the registration process, an update displacement field (𝐷𝑢) is first computed as a function (Θ) of the fixed image (𝑓) and the moving image (𝑚) warped by the displacement field resulting from previous iterations (ğ·ğ‘Ž):𝐷𝑢⟵Θ𝑓,ğ‘šâˆ˜Î”ğ‘Žî€¸,(5) where Î”ğ‘Ž and Δ𝑢 denote the deformation operations linked to ğ·ğ‘Ž and 𝐷𝑢, 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 ğ·ğ‘Ž must be increased by the update fieldğ·ğ‘Žî€·ğ·âŸµÎ¦ğ‘Ž,𝐷𝑢.(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 ⊕ introduced earlier. The accumulation ğ·ğ‘ŽâŠ•ğ·ğ‘¢ of the displacement fields ğ·ğ‘Ž with 𝐷𝑢 is then compatible with the way 𝐷𝑢 is estimated. Indeed, since 𝐷𝑢 is computed from ğ‘šâˆ˜Î”ğ‘Ž, the accumulation of ğ·ğ‘Ž and 𝐷𝑢 must modify ğ·ğ‘Ž by 𝐷𝑢, a process intrinsically integrated by the operation ⊕. Moreover, the associativity of ⊕ 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:ğ·ğ‘Žî€·ğ·âŸµÎ¨ğ‘Žî€¸.(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 𝑓 and 𝑚, some parameters described below, and outputs the estimated transformation Î”ğ‘Ž=Id+ğ·ğ‘Ž such that ğ‘“â‰ƒğ‘šâˆ˜Î”ğ‘Ž. The whole procedure described in Table 1 and depicted in Figure 3 involves computations on different scales 𝑗∈[0,𝐽], from coarse (𝑗=𝐽) to fine (𝑗=0). Each scale is associated to a subsampled grid of voxels 𝒢𝑗=𝜅𝑗𝒢, where 𝜅 is the subsampling factor (e.g., √𝜅=2 in this study) between scale 𝑗 and scale 𝑗+1. The functions 𝑓 and 𝑔, defined on the initial grid 𝒢=𝒢0={(𝑥1,𝑥2,𝑥3)∈ℤ3}, are downsampled (after antialiasing smoothing) at any scale 𝑗 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 𝒢𝑗+1 to the finer grid 𝒢𝑗 using 𝜅 as upsampling factor. For each scale 𝑗∈[0,𝐽], the accumulated displacement field is iteratively updated until one reaches a particular stopping criterion 𝒮 (e.g., based on the convergence of ğ·ğ‘Ž or on the SSD, as precised in Section 5).

Inputs and parameters:
 (i) Images 𝑓 and 𝑚 defined on 1 D 4 A 2 .
 (ii) Number of scales 𝐽 .
 (iii) A stopping criterion 𝒮 .
 (iv) Gaussian kernel variance ğœŽ 2 Ψ of Ψ .
Output: The displacement field 𝐷 ğ‘Ž .
 (4)  Initialization:
  Set scale to 𝑗 = 𝐽 and initialize 𝐷 ğ‘Ž = 0 on 𝒢 𝐽 + 1 .
  (5)  Transfer on grid 𝒢 𝑗 :
  Compute 𝑚 𝑗 = Downj ( 𝑚 ) , 𝑓 𝑗 = Downj ( 𝑓 ) ,
  and assign 𝐷 ğ‘Ž ←   Up ( 𝐷 ğ‘Ž ) .
  (6) While 𝒮 is false, do:
  (i) Warping: 𝑤 = 𝑚 𝑗 ∘ Δ ğ‘Ž
  (ii) Field computation: 𝐷 𝑢 ← Θ ( 𝑓 𝑗 , 𝑤 )
  (iii) Accumulation: 𝐷 ğ‘Ž ← Φ ( 𝐷 ğ‘Ž , 𝐷 𝑢 )
  (iv) Regularization: 𝐷 ğ‘Ž ← Ψ ( 𝐷 ğ‘Ž )
  (7) If 𝑗 = 0 , stop and return 𝐷 ğ‘Ž ,
  else, set 𝑗 ← 𝑗 − 1 and return to step (5).

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 𝑓 and 𝑤=𝑚∘Δ, the displacement field estimation Θ is actually split into two quantities 𝐷𝑢⟵Θ𝐷𝑐(𝑓,𝑤),𝑢⟵Θ𝑐(𝑓,𝑤),(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 â„Žğœ‚ of direction 𝜂∈ℝ3 (also called loglets [40]) defined in frequency by the polar separable function 𝐻𝜂(𝜔)=𝜒+𝜂𝑇𝜔𝜂𝑇𝜔2𝑅(‖𝜔‖),(9) where 𝜔∈ℝ3 is the frequency vector, 𝜒+(𝜆)=1 if 𝜆>0 and 0 else, ‖𝜔‖2=𝜔𝑇𝜔, 𝜔=𝜔/‖𝜔‖ is the unit vector supporting 𝜔, and 𝑅 is a radial function centered on 𝜌>0 and defined as 𝑅(𝑟)=exp[−ln2(𝑟/𝜌)/ln2] for 𝑟>0.

Since their support corresponds to the half volume {𝜔∈ℝ3∶𝜂𝑇𝜔>0} and since (𝜂𝑇𝜔)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 𝜂𝑘=(cos𝜙𝑘,sin𝜙𝑘) with 𝜙𝑘=𝑘𝜋/4 for 0≤𝑘≤3, and in 3D, 𝜂 is taken as the 6 normal vectors {𝜂𝑘∶0≤𝑘≤5} to the faces of a hemi-icosahedron [42, 43]. Notice also that each filter â„Žğ‘˜(𝑥) 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 𝑓, defining the filtering ğ‘žğ‘“î€·(𝑥;𝑘)=ğ‘“âˆ—â„Žğ‘˜î€¸(𝑥),(10) with ∗ the common convolution operation and the shorthand â„Žğ‘˜=â„Žğœ‚ğ‘˜, we can write ğ‘žğ‘“(𝑥;𝑘)=𝐴𝑓(𝑥;𝑘)𝑒𝑖𝜙𝑓(𝑥;𝑘) since ğ‘žğ‘“âˆˆâ„‚. Therefore, by processing the warped image 𝑤 similarly, the local phase difference can be computed as Î”ğœ™ğ‘˜î€·ğ‘ž(𝑥)=arg𝑓(𝑥;𝑘)ğ‘žâˆ—ğ‘¤î€¸,(𝑥;𝑘)(11) with (⋅)∗ the complex conjugation and Δ𝜙𝑘(𝑥)=𝜙𝑓(𝑥;𝑘)−𝜙𝑤(𝑥;𝑘) the local dephasing between 𝑓 and 𝑤 in direction 𝜂𝑘.

An important observation is that the nonnegative value 𝐴𝑓||(𝑥;𝑘)=ğ‘“âˆ—â„Žğ‘˜î€¸||=||||(𝑥)ℝ3ğ‘“î€·ğ‘¥î…žî€¸ğ‘‡ğ‘¥â„Žğ‘˜î€·ğ‘¥î…žî€¸dğ‘¥î…ž||||(12) represents also the correlation between 𝑓(𝑥′) and the translated filter ğ‘‡ğ‘¥â„Žğ‘˜(ğ‘¥î…ž)=â„Žğ‘˜(ğ‘¥î…žâˆ’ğ‘¥); that is, the filter â„Žğ‘˜(ğ‘¥î…ž)=â„Žğ‘˜(âˆ’ğ‘¥î…ž) translated on 𝑥. If the image 𝑓 was perfectly represented by the latter, that is, if we had locally 𝑓(ğ‘¥î…ž)=ğ‘â„Žğ‘˜(ğ‘¥î…žâˆ’ğ‘¥) for any ğ‘¥î…žâˆˆâ„3 and some constant 𝑐∈ℝ, a displacement of 𝑓 by a displacement field 𝐷(𝑥) approximately constant over the support of ğ‘‡ğ‘¥â„Žğ‘˜ would induce a dephasing Δ𝜙𝑘(𝑥)=𝜌𝜂𝑇𝑘𝐷(𝑥) since the frequency vector of ğ‘‡ğ‘¥â„Žğ‘˜ is −𝜌𝜂𝑘. An important implicit assumption is nevertheless that |𝜌𝜂𝑇𝑘𝐷(𝑥)|<𝜋 since the dephasing is known up to modulo 2𝜋. Moreover, only 𝜂𝑇𝑘𝐷 and not 𝐷 can be determined, as another manifestation of the blank wall problem [44].

In practice, for most of 𝑥, 𝑓(𝑥) is not perfectly represented by one filter but by a linear combination of them where the amplitude 𝐴𝑓(𝑥;𝑘) measures the adequacy of the fit between 𝑓(𝑥) and ğ‘‡ğ‘¥â„Žğ‘˜. Consequently, the local update displacement 𝐷𝑢(𝑥) linking 𝑓(𝑥) and 𝑤(𝑥)=𝑓(𝑥+𝐷𝑢(𝑥)) in each 𝑥∈ℝ3 is estimated by solving the weighted least square optimizationΘ𝐷(𝑓,𝑤)=argmin𝑑∈ℝ3𝑘𝑐𝑘𝜌𝜂𝑇𝑘𝑑−Δ𝜙𝑘2,(13) where the 𝑐𝑘(𝑥)=𝐴𝑓(𝑥;𝑘)𝐴𝑚(𝑥;𝑘) are the certainty map of the filter â„Žğ‘˜. As explained above, 𝑐𝑘 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 𝐸(𝑑)=‖𝐶(𝑁𝑑−Γ)‖22, using the diagonal matrix 𝐶=diag(𝑐1,…,𝑐6), the matrix 𝑁=(𝜂1,…,𝜂6)𝑇, and the vector Γ=(Δ𝜙1,…,Δ𝜙6)𝑇. An easy computation shows that the solution of (13) is then given by the Moore-Penrose pseudoinverse (𝐶𝑁)†=(𝑁𝑇𝐶2𝑁)−1𝑁𝑇𝐶, that is, Θ𝐷(𝑓,𝑤)=(𝐶𝑁)†𝐶Γ,(14) with Θ𝐷(𝑓,𝑤) arbitrary set to 0 when (𝑁𝑇𝐶2𝑁) is not invertible.

Jointly to the estimation (13), a global certainty map associated to the quality of the estimation of Θ𝐷 is defined as [43] Θ𝑐(𝑓,𝑤)=𝑘𝑐𝑘(𝑥),(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 𝐻𝜂𝑘 is equivalent to estimating the phase of the bandpass-filtered image around increasing cutoff frequencies, that is, with 𝜌←2𝜌 each time 𝑗←𝑗+1. This sustains the coarse-to-fine displacement estimation, that is, the computation of Θ𝐷 and Θ𝑐 on different scale bands 𝑓𝑗 and 𝑚𝑗 of 𝑓 and 𝑚.

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 (𝑐𝑢, as computed from Θ𝑐) and the accumulated certainty map (ğ‘ğ‘Ž). 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 Φ𝐷 (field accumulation) and Φ𝑐 (certainty accumulation):Î¦ğ·î€·ğ·ğ‘Ž,𝐷𝑢,ğ‘ğ‘Ž,𝑐𝑢=ğ·ğ‘Ž+ğ‘ğ‘¢ğ‘ğ‘Ž+𝑐𝑢𝐷𝑢,Φ(16)ğ‘î€·ğ‘ğ‘Ž,𝑐𝑢=𝑐2ğ‘Ž+𝑐2ğ‘¢ğ‘ğ‘Ž+𝑐𝑢,(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 Î¦ğ·î€·ğ·ğ‘Ž,𝐷𝑢,ğ‘ğ‘Ž,𝑐𝑢=ğ·ğ‘ŽâŠ•ğ‘ğ‘¢ğ‘ğ‘Ž+𝑐𝑢𝐷𝑢.(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Î¦ğ·î€·ğ·ğ‘Ž,𝐷𝑢=ğ·ğ‘ŽâŠ•î€·î€·ğ·exp𝑢−Id,(19) where the field exponential exp(𝐷𝑢) can be efficiently estimated using a small number of recursive compositions of the field 𝐷𝑢 by itself. Consequently, the displacement field Φ𝐷(ğ·ğ‘Ž,𝐷𝑢) is linked to the deformation operation Î”ğ‘Žâˆ˜exp(𝐷𝑢).

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:Î¦ğ·î€·ğ·ğ‘Ž,𝐷𝑢,ğ‘ğ‘Ž,𝑐𝑢=ğ·ğ‘ŽâŠ•î‚µî‚µğ‘expğ‘¢ğ‘ğ‘Ž+𝑐𝑢𝐷𝑢−Id.(20) Since exp(0𝐷)=Id for any vector field 𝐷, the accumulation fades away when ğ‘ğ‘Žâ‰«ğ‘ğ‘¢. 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 ℎ and a filter 𝑔 (typically a Gaussian kernel of variance ğœŽÎ¨>0), the normalized convolution of a (scalar) function 𝑠 by 𝑔 as involved by the normalization ℎ is ğ‘ âˆ—â„Žğ‘”â‰œ(â„Žğ‘ )∗𝑔.â„Žâˆ—ğ‘”(21) This operation does not increase the maximum amplitude of the filtered function. Indeed, for a nonnegative kernel 𝑔, we show easily that â€–ğ‘ âˆ—â„Žğ‘”â€–âˆžâ‰¤â€–ğ‘ â€–âˆž, with â€–ğ‘ â€–âˆž=max𝑥|𝑠(𝑥)|. The accumulated displacement field ğ·ğ‘Ž and subsequently the certainty map are therefore regularized thanks to this operation using for normalization the certainty map ğ‘ğ‘Ž, that is, Î¨ğ·î€·ğ·ğ‘Ž,ğ‘ğ‘Žî€¸=ğ·ğ‘Žâˆ—ğ‘ğ‘ŽÎ¨ğ‘”,ğ‘î€·ğ‘ğ‘Ž,ğ‘ğ‘Žî€¸=ğ‘ğ‘Žâˆ—ğ‘ğ‘Žğ‘”.(22) Notice that, for computing Ψ𝐷, 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).

The second virtual experiment is based on two images of a disk (see Figure 5). In the fixed image, a disk of radius 𝑟1+𝑟2 was created, and a hole (disk of radius 𝑟2) was added in its center. In the moving image, a disk of radius 𝑟1 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.

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 𝐽=8, with final resolution of 2mm×2mm×2mm. 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.(i)The landmark position error evaluates the ability of the registration in finding the physical motion of organs. (ii)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 ∑𝑥(𝑓−𝑚∘Δ)2.(iii)The harmonic energy [29, 32] of the displacement field 𝐷 indicates how regular the field is and is computed as ∑(1/2)𝑥(‖∇𝐷1‖2+‖∇𝐷2‖2+‖∇𝐷3‖2).(iv)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(𝒥), with 𝒥𝑖𝑗=𝜕Δ𝑖/𝜕𝑥𝑗=𝛿𝑖𝑗+𝜕𝑑𝑖/𝜕𝑥𝑗, where 𝛿𝑖𝑗 is Kronecker's delta (𝛿𝑖𝑗=1 if 𝑖=𝑗, 0 else) and 𝑑𝑖 is the 𝑖th component of the displacement field. In practice, the partial derivatives 𝜕𝑑𝑖/𝜕𝑥𝑗 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.

RCCT phases Original Demons [4] Morphons D-Demons D-Morphons

Phase 1 0.5/0.5 (2.4) 1.3/0.3 (1.8) 0.7/0.3 (1.6) 0.7/0.3 (1.6) 0.7/0.3 (1.6)
Phase 2 0.5/0.6 (2.6) 1.4/0.2 (2.1) 0.7/0.4 (2.1) 0.7/0.4 (1.6) 0.7/0.4 (2.1)
Phase 3 2.2/1.8 (6.6) 1.4/0.4 (2.3) 1.2/0.6 (2.5) 1.2/0.6 (2.5) 1.2/0.6 (2.4)
Phase 4 4.3/2.5 (10) 1.2/0.4 (2.3) 1.0/0.4 (2.2) 1.0/0.5 (2.5) 1.0/0.4 (2.2)
Phase 5 5.8/2.6 (12) 1.3/0.5 (2.6) 1.1/0.5 (2.7) 1.1/0.5 (2.5) 1.1/0.5 (2.8)
Phase 6 6.1/2.9 (14) 1.1/0.4 (2.0) 1.0/0.5 (2.1) 1.1/0.6 (2.8) 1.0/0.5 (2.1)
Phase 7 5.0/2.3 (12) 1.3/0.5 (2.4) 1.1/0.6 (2.8) 1.2/0.6 (2.7) 1.1/0.6 (2.8)
Phase 8 3.7/1.6 (6.2) 1.1/0.3 (1.7) 0.8/0.4 (1.9) 0.8/0.4 (1.8) 0.8/0.4 (1.8)
Phase 9 2.1/1.1 (4.5) 1.1/0.3 (1.9) 0.8/0.4 (2.0) 0.8/0.4 (1.7) 0.8/0.4 (2.0)

All phases 3.3/2.0 (14) 1.2/0.4 (2.6) 0.9/0.5 (2.8) 1.0/0.5 (2.8) 0.9/0.5 (2.8)

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).

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.

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.

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.

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.

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.

Original Demons Morphons D-Demons D-Morphons

Volume change (vessels) [in %] 23 0 21 0
Volume change (tumor) [in %] 6 1 6 1
Harmonic energy (vessels) [ × 1 0 − 3 ] 89 8 69 8
Harmonic energy (tumor) [ × 1 0 − 3 ] 39 4 34 4
Error on COM (tumor) [in m m ] 2.1 1.5 1.1 1.5 1.1

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.

Low harmonic energyHigh harmonic energy

Invertible ( 𝐽 m i n > 0 ) D-Morphons D-Demons
Noninvertible ( 𝐽 m i n < 0 ) Morphons Demons

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.


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.


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

Copyright © 2011 Guillaume Janssens et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

2087 Views | 780 Downloads | 49 Citations
 PDF  Download Citation  Citation
 Download other formatsMore
 Order printed copiesOrder