Health-e-Child - IST-2004-027749 - Deliverable D.11.4

Heart Diseases

iLogDemons, a Consistent Incompressible Image Registration Algorithm

1. Introduction: Medical Background and Challenges

Tissue tracking in sequences of medical images is an important task in many applications, either for therapy guidance or diagnosis. However there is no easy way to achieve it, even interactively. A standard approach is now to use non-linear image registration to estimate the spatial transformation between different images, for instance brain shift (Ferrant et al., 2001). However, for such clinical applications one-to-one mapping must often be ensured and the model of transformations that is used must be consistent with the properties of the tissue to track. In particular, constraining the registration to be volume-preserving showed great improvement when tracking incompressible tissues such as the heart (Bistoquet et al., 2008; Mansi et al., 2009).

Mathematical frameworks based on diffeomorphic deformations (Beg et al., 2005; Bossa et al., 2007) have been developed to get one-to-one mappings between the images to register. Among them, logDemons (Vercauteren et al., 2008) is an efficient demons-based registration algorithm that provides diffeomorphic transformations parameterised by stationary velocity fields. However, although mathematical justifications of demons optimisation have been provided (Vercauteren et al., 2008), theoretical foundations of the Gaussian regularisation still has to be consolidated (Pennec et al., 1999), which makes the algorithm difficult to adapt to other deformation models.

Several approaches for incompressible image registration have been proposed. A first method is to constrain the Jacobian determinant of the transformations to equal one. However, this constraint is computationally demanding due to its non-linearity (Rohlfing et al., 2003). Linear approximations have been proposed (Bistoquet et al., 2008) but volume drifts may appear at large deformations. Velocity fields can be made incompressible by constraining them to be divergence-free. Thus, incompressible fluid registration is achieved by projecting the update velocity onto the space of divergence-free vector fields (Saddi et al., 2007). Nevertheless, the fluid model might not be appropriate for tracking biological soft tissues like myocardium.

In Health-e-Child we developed an efficient and consistent framework for demons-based incompressible registration. We first propose a mathematical justification of the Gaussian regularisation, which enables to integrate the incompressibility constraint seamlessly by working on the space of divergence-free velocity fields. The main advantages of this are:

1. Introduction: Medical Background and Challenges

Background: Log-Domain Diffeomorphic Demons

LogDemons algorithm estimates a dense non-linear transformation φ that best aligns a template image T to a reference image R (Vercauteren et al., 2008). φ belongs to the space generated by the one-parameter subgroups of diffeomorphisms. They are parameterised by stationary velocity fields v through the exponential map φ=exp(v) (Arsigny et al., 2006). The images R and T are registered by minimising in the space of velocities, called log-domain, the energy functional:

where σi2 relates to the noise in the images and σd2 controls the regularisation strength. This equation is function of two variables: The velocity field v parameterises the transformation φ to recover whereas vc parameterises an intermediate transformation φc = exp(vc) that models the correspondences between the voxels of the two images. During the optimisation step, E(v,vc) is minimised with respect to vc. This amounts to finding the optimal matching between R and T without considering the regularisation. Under the diffeomorphic update rule φc ← φ o exp(Δv), the optimal update velocity writes

Δv(x) = -(R(x) - T o φ(x)) / (||J(x)||2+σi2/σx2) J(x)

J(x) is the symmetric gradient J(x) = (∇R(x) + ∇(T o φ)(x)) / 2. The correspondence velocity vc is then updated using the first order approximation of Baker-Campbell-Hausdorff (BCH) formula:

vc = Z(v,Δv) = v + Δv + 1/2[v, Δv] + 1/12[v, [v, Δv]] + O(||Δv||2)

where the Lie bracket [·,·] is defined by [v,Δv] = (∇v) Δv - (∇Δv)v.

This approach has experimentally shown promising results in terms of image registration and statistics on diffeomorphisms (Bossa et al., 2007; Vercauteren et al., 2008). Finally, the regularisation step estimates the optimal regularised transformation φ by minimising E(v,vc) with respect to v, which is approximated by smoothing the correspondence velocity vc with a Gaussian kernel Gσ. However, how Gσ relates to E(v,vc) remains to be consolidated (Pennec et al., 1999).

Revisiting Demons Gaussian Regularisation

A consistent mathematical formulation of demons regularisation is required to integrate incompressibility in logDemons. In (Nielsen et al., 1994), the authors demonstrate that Gaussian filtering solves the Tikhonov estimation problem with equal weighting of the spatial derivatives in the Taylor series sense. We thus replace the logDemons regulariser ||∇v||2 by the Tikhonov regulariser to get:

In this equation, Ω is the image domain and ∂ik..il denotes the composition of spatial derivatives ∂ik..∂il. The parameter σx2 has been introduced into the regulariser to simplify calculations. More importantly, the regularisation weight σd2 is now function of the derivative orders to preserve the shape of the impulse response related to the regulariser (Nielsen et al., 1994). The previous equation is minimised by linearising its first term using the zeroth order approximation of BCH formula, log(φ-1 o φc) = vc - v. The resulting equation, which is exactly a Tikhonov estimation problem, is minimised in the Fourier domain. The optimal velocity field v verifies

which is exactly demons Gaussian smoothing Gσ, with σ2 = 2/σd2. The width of the Gaussian kernel corresponds to the strength of the regularisation.

This formulation is built up on two key elements. First, the correspondence field φc decouples the regularisation from the optimisation, making the minimisation of the registration energy independent of the optimisation. Second, the coupling term ||log(φ-1 o φc)||2 is approximated as a least-square problem to get the Gaussian smoothing. These two elements must be ensured in any demons algorithm to justify the Gaussian regularisation.

Incompressible LogDemons

The proposed regulariser enables to integrate the incompressibility constraint into the algorithm. A transformation φ is locally incompressible if its Jacobian determinant det(∇φ)=1. In fluid dynamics, one uses the infinitesimal version: a fluid is said incompressible if the divergence of its velocity v is null. For diffeomorphic transformations one can show that the converse is true: integrating divergence-free velocities over time yields incompressible deformations (Evans, 1998). Making LogDemons incompressible thus consists in constraining the velocity fields v to be divergence-free. Demons optimisation step is not modified, as it optimises vc only, but demons regularisation energy is now optimised under the divergence-free constraint, which amounts to minimising the Lagrange function:

In this equation, the Lagrange multiplier p is a scalar function of the Sobolev space H01(Ω) that vanishes at infinity. Optima of the previous equation are found by solving ∂vP(v,p)=0:

with p=0 at the domain boundaries ∂Ω. The divergence of the previous equation, under the optimal condition ∇·v=0, yields Δp = ∇· vc with 0-Dirichlet boundary conditions, which can be solved independently of v to get p. g = vc - ∇p is thus the L2 projection of vc onto the space of divergence-free vector fields. From the previous section we deduce the optimal incompressible velocity field: v ←Gσ⊗g.

Particular care must be taken when incompressibility is required within a subdomain Γ ⊂ Ω only, for tracking localised incompressible organs like the heart. This is achieved by defining p ∈ H01(Γ), p=0 on Ω/Γ. Although Gaussian smoothing theoretically preserves the divergence, in practice unconstrained velocities close to Γ boundaries (∂Γ) may leak inside the incompressible domain due to the Gaussian convolution, ultimately resulting in volume drifts. Yet, vector derivatives and well-designed Gaussian filters, like Deriche recursive filters, commute. We therefore replace the theoretical "project-and-smooth" strategy by a "smooth-and-project" approach that preserves divergence. To limit numerical instabilities, a smooth domain transition is implemented in a narrow band around Γ by diffusing the pressure field p using heat-transfer equation (Evans, 1998). The following algorithm summarises the main steps of our method, henceforth termed iLogDemons, which was implemented using ITK (Dru et al., 2009) and PETSc.

iLogDemons: Incompressible LogDemons Registration
  • REQUIRE Initial stationary velocity field v0. Usually v0 = 0 (i.e. φ0=Id)
  • LOOP over n until convergence
    • Compute update velocity Δvn, given vn-1.
    • Fluid-like regularisation: Δvn ← Gσf ⊗ Δvn.
    • Update the correspondence velocity: vn ← Z(vn-1,Δvn).
    • Diffusion-like regularisation: vn ← Gσ ⊗ vn.
    • Solve: Δp = ∇·vn, p=0 on the incompressible domain boundaries.
    • Project the velocity field: vn ← vn - ∇p.
  • v, φ = exp(v) and φ-1 = exp(-v).

References

  1. Arsigny, V., Commowick, O., Pennec, X., Ayache, N.: A log-euclidean framework for statistics on diffeomorphisms. In MICCAI'06, LNCS, v. 4190, pp. 924 Springer (2006)
  2. Beg, M., Miller, M., Trouvé, A., Younes, L.: Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. IJCV 61(2), 139-157 Springer (2005)
  3. Bistoquet, A., Oshinski, J., Skrinjar, O.: Myocardial deformation recovery from cine MRI using a nearly incompressible biventricular model. Medical Image Analysis 12(1), 69-85 Elsevier (2008)
  4. Bossa, M., Hernandez, M., Olmos, S.: Contributions to 3D diffeomorphic atlas estimation: application to brain images. In MICCAI'07, LNCS, vol. 10, p. 667 Springer (2007)
  5. Dru, F., Vercauteren, T.: An ITK implementation of the symmetric log-domain diffeomorphic demons algorithm. Insight Journal (May 2009)Download paper
  6. Evans, L.C.: Partial Differential Equations. (June 1998)
  7. Ferrant, M., Nabavi, A., Macq, B., Jolesz, F., Kikinis, R., Warfield, S.: Registration of 3D intraoperative MR images of the brain using a finite element biomechanical model. IEEE TMI 20(12), 1384-1397 (2001)
  8. Hinkle, J., Fletcher, P., Wang, B., Salter, B., Joshi, S.: 4D map image reconstruction incorporating organ motion. In IPMI'07, LNCS, p. 687. Springer (2009)
  9. Mansi, T., Peyrat, J.M., Sermesant, M., Delingette, H., Blanc, J., Boudjemline, Y., Ayache, N.: Physically-constrained diffeomorphic demons for the estimation of 3D myocardium strain from cine-MRI. In FIMH'09, LNCS, v. 5528, p. 201-210 Springer (2009)Download paper
  10. Nielsen, M., Florack, L., Deriche, R.: Regularization, scale-space, and edge detection filters. In JMIV, 291-307 Springer (1997)
  11. Pennec, X., Cachier, P., Ayache, N.: Understanding the "demon's algorithm": 3D non-rigid registration by gradient descent. In MICCAI'99. LNCS, vol. 1679, pp. 597-605. Springer (1999)Download paper
  12. Rohlfing, T., Maurer Jr, C., Bluemke, D., Jacobs, M.: Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint. IEEE TMI 22(6), 730-741 (2003)
  13. Saddi, K.A., Chefd'hotel, C., Cheriet, F.: Large deformation registration of contrast-enhanced images with volume-preserving constraint. In: SPIE Medical Imaging. vol. 6512. (2007)
  14. Vercauteren, T., Pennec, X., Perchant, A., Ayache, N.: Symmetric log-domain diffeomorphic registration: A demons-based approach. In MICCAI'08, LNCS, vol. 5241, pp. 754-761. Springer (2008)Download paper
  15. Moore,C.,Lugo-Olivieri,C.,McVeigh,E.,Zerhouni,E.: Three-dimensional systolic strain patterns in the normal human left ventricle: Characterization with tagged MR imaging. Radiology 214(2), 453-466 (2000)