We propose an iterative two-step method to compute a diffeomorphic non-rigid transformation between images of anatomical structures with rigid parts, without any user intervention or prior knowledge on the image intensities. First we compute spatially sparse, locally optimal rigid transformations between the two images using a new block matching strategy and an efficient numerical optimiser (BOBYQA). Then we derive a dense, regularised velocity field based on these local transformations using matrix logarithms and M-smoothing. These two steps are iterated until convergence and the final diffeomorphic transformation is defined as the exponential of the accumulated velocity field. We show our algorithm to outperform the state-of-the-art log-domain diffeomorphic demons method on dynamic cervical MRI data.

In medical image analysis, one is often confronted with the problem of registering anatomical structures containing both hard, rigid (typically, bones) and soft, non-rigid (most other tissues) parts. Such problems are met for instance when following-up spinal cord lesions in MRI for the diagnosis of multiple sclerosis [1], or when assessing cervical injuries using dynamic/kinematic MR imaging with positional changes [2]. Many methods have been developed for both fully global rigid registration and fully local non-rigid registration separately [3], but the literature on hybrid methods, allowing for adequate registration of the structures depending on the stiffness of their components, is still quite sparse.

The earliest work we know of is that of Little et al. [4], who showed how to incorporate rigid structures into a deformation field, using radial basis functions; this was later improved by others to make the field invertible and even diffeomorphic [5–7]. However, these methods require the user to specify which structures are rigid, which led to the development of semi-automated methods in which rigidity can be locally favoured/enforced through a regularisation term in the criterion to be minimised [8]. This idea was later improved to allow for this term to be adaptively tuned to the structures to register, through prior segmentation of the rigid parts or design of a stiffness map (typically computed from the image intensities; e.g. bones have high intensities in CT) [9–11]. Instead of segmenting rigid parts, it was also proposed to define several anchors, to which is attached an unknown polyaffine transformation, which can be subsequently estimated using a modified EM-ICP algorithm [12].

In this paper, we propose an iterative two-step method to compute a diffeomorphic non-rigid transformation between images of structures with rigid parts, without any user intervention or prior knowledge on the image intensities (to compute rigid parts or anchors). First we compute spatially sparse, locally optimal rigid transformations between the two images by adopting a new (as opposed to classical, translation-based) block matching strategy, made possible by the use of an efficient numerical optimiser (BOBYQA) (Sec. 2.1). The rationale behind this original strategy is our hope to recover both large rotations and subvoxel displacements. Then we derive a dense, regularised velocity field based on these local transformations using matrix logarithms and M-smoothing (Sec. 2.2). The floating image is then resampled and the two steps are iterated until convergence; the final diffeomorphic transformation is defined as the exponential of the accumulated velocity field. We finally compare our algorithm with the state-of-the-art log-domain diffeomorphic demons method [13] on dynamic cervical and multiple sclerosis MR images (Sec. 3).

To compute a diffeomorphism T between a reference image I and a floating image J, we iterate between two steps: computation of a sparse set of locally optimal rigid transformations using block matching between I and J ∘ T^l (Sec. 2.1) and computation of a dense velocity field δLT^l computed from these locally estimated transformations (Sec. 2.2). Given that the transformation T is initially set to the identity (T⁰ = Id), and that the initial velocity field is set to LT⁰ = log T⁰ = 0, the velocity field is then updated as LT^l+1 = LT^l + δLT^l. This two-step algorithm stops at the iteration l when δLT^l is close to 0, and the final diffeomorphism is computed as T = T^l = exp(LT^l). The complete algorithm is outlined in Sec. 2.3. For the sake of clarity, we detail the two steps using the simpler notations I, J and δLT (Sec. 2.1 and 2.2).

2.2. Estimating a dense velocity field

The set of estimated transformations (R₁, …, R_m) is spatially sparse, but is also noisy and likely to contain outliers (due to the noise in the images to be registered and the potential errors in local registrations). How to estimate a dense (n = card(I)) and smooth velocity field δLT from (R₁, …, R_m)? We propose to use the logarithms of these m transformations, defined in the space of 4 × 4 real matrices ( (ℝ)) restricted to those whose last row contains only zeros, and to estimate n intermediate matrices in the same space (that we name log S₁, …, log S_n by analogy) as the minimisers of a criterion C: where:

• ρ: ℝ → ℝ⁺ is a robust error norm,
• ||.|| is the Frobenius norm in ( ℝ), and |.| is the Euclidean norm in ℝ³,
• v_j is the coordinate of the central voxel of the block where R_j was estimated,
• v_i is the coordinate of the voxel where S_i is to be computed,
• V_i is a neighbourhood around the position v_i; note that the sum over j ∈ V_i must read: “the sum over all the points v_j where a rigid transformation R_j was estimated, and which are inside V_i”,
• w_j is the weight defined in Sec. 2.1,
• d(.) : ℝ³ → ℝ⁺ is a (spatial) error norm.

It must be clear that we do not estimate S_i and then its logarithm; we do estimate log S_i directly; we use this notation here only as a convention for the sake of simplicity. Solving this minimisation problem is known as local M-smoothing [15], due to the use of a robust error norm and of spatial neighbourhoods to design C; C can be minimised using gradient descent, where each transformation can be estimated independently of the others. Using a particular adaptive, data-dependent step size leads to an easy-to-interpret update formula for each log S_i [15]:

It can be seen from this formula that ρ′ acts as a tonal kernel, while d acts as a spatial kernel. After convergence, each finally estimated log S_i is a linear combination of the logarithms of the rigid transformations R_j; following Arsigny et al. [5], we define the final dense velocity field δLT as δLT (v_i) = log(S_i).v_i, ∀i = 1, …, n. In practice, we define ρ as the Welsch function, which leads to ρ′(a²) = exp(−a²/2λ²), and we define d as d(b²) = exp(−b²/2θ²); this leads to two similar expressions for the two kernels (with different bandwidths). V_i is spherical with radius 2θ (to achieve an approximate 95% confidence interval for a Gaussian law). To initialise the gradient descent algorithm, is computed as the solution of the update formula by setting ρ′(a²) = 1 (i.e. no tonal kernel).

We propose to assess our algorithm quantitatively on ten patients with traumatic cervical cord injury, who got dynamic cervical MRI (T2-w, size 384 × 384 × 14, voxel size 0.8 × 0.8 × 3 mm³) with two different positions each: either flexion/neutral, or extension/neutral [2]. For each patient, we manually defined landmarks on the cervical/thoracic vertebrae C1-C3-C6-T1 (and T4 when visible), the pontomedullary junction, and the gnathion (lower border of the mandible) on each of the two MRI. We considered the neutral position as the reference image in the extension/neutral setting, and the flexion as the reference image in the flexion/neutral setting. For a given patient, the registration accuracy was evaluated as the root mean square error (RMSE) computed over the homologous landmarks after registration using four different methods: global rigid registration (M₁) [16], log-domain diffeomorphic demons (M₂) [13], our algorithm (M₃), and M₂ initialised using M₃ (M₄); both M₂ and M₃ are initialised using M₁. We also assessed our algorithm visually on two patients with MS lesions in the spinal cord, and two other patients with tumours in the spinal cord, with two time points each (T1-w, size 256 × 256 × 64, voxel size 1 × 1 × 1 mm³).

The box-and-whisker plot in Fig. 1 computed from the ten patients shows our algorithm (M₃) to significantly outperform both M₁ (paired t-test: p = 3×10⁻⁴) and M₂ (paired t-test: p = 2 × 10⁻³), with much smaller error and much smaller error dispersal. M₃ is also slightly better than M₄ (paired t-test: p = 5 × 10⁻²). This suggests that the log-domain diffeomorphic demons performs worse than our algorithm even when properly initialised (using M₃ instead of M₁). The results of M₁, M₂ and M₃ on one of the ten patients are shown in Fig. 2, and that of M₁ and M₃ are shown on one of the MS patients in Fig. 3.

It appears that our strategy for non-rigid registration, based on the computation of locally optimal rigid transformations in the first place, allows us to recover displacements and deformations of piecewise rigid structures (as seen e.g. in dynamic cervical MRI) much better than standard methods which are implicitly based on locally optimal translations, such as the log-domain diffeomorphic demons algorithm. This original strategy was made possible by (i) the use of an up-to-date very efficient optimiser and (ii) the design of a specific regularising procedure on the (sparse) set of locally estimated rigid transformations, based on robust estimation techniques. A future line of research could be to combine our regularisation technique with those previously proposed in this context [9–11]. Our intuition is also that our algorithm could perform very well in more general problems, without necessarily rigid structures involved, and on other image modalities; we will evaluate this in a near future.