PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Proc SPIE Int Soc Opt Eng. Author manuscript; available in PMC 2016 June 15.
Published in final edited form as:
Proc SPIE Int Soc Opt Eng. 2016 February 27; 9786: 97860H.
Published online 2016 March 18. doi:  10.1117/12.2208621
PMCID: PMC4909372
NIHMSID: NIHMS780672

MIND Demons for MR-to-CT Deformable Image Registration In Image-Guided Spine Surgery

Abstract

Purpose

Localization of target anatomy and critical structures defined in preoperative MR images can be achieved by means of multi-modality deformable registration to intraoperative CT. We propose a symmetric diffeomorphic deformable registration algorithm incorporating a modality independent neighborhood descriptor (MIND) and a robust Huber metric for MR-to-CT registration.

Method

The method, called MIND Demons, solves for the deformation field between two images by optimizing an energy functional that incorporates both the forward and inverse deformations, smoothness on the velocity fields and the diffeomorphisms, a modality-insensitive similarity function suitable to multi-modality images, and constraints on geodesics in Lagrangian coordinates. Direct optimization (without relying on an exponential map of stationary velocity fields used in conventional diffeomorphic Demons) is carried out using a Gauss-Newton method for fast convergence. Registration performance and sensitivity to registration parameters were analyzed in simulation, in phantom experiments, and clinical studies emulating application in image-guided spine surgery, and results were compared to conventional mutual information (MI) free-form deformation (FFD), local MI (LMI) FFD, and normalized MI (NMI) Demons.

Result

The method yielded sub-voxel invertibility (0.006 mm) and nonsingular spatial Jacobians with capability to preserve local orientation and topology. It demonstrated improved registration accuracy in comparison to the reference methods, with mean target registration error (TRE) of 1.5 mm compared to 10.9, 2.3, and 4.6 mm for MI FFD, LMI FFD, and NMI Demons methods, respectively. Validation in clinical studies demonstrated realistic deformation with sub-voxel TRE in cases of cervical, thoracic, and lumbar spine.

Conclusions

A modality-independent deformable registration method has been developed to estimate a viscoelastic diffeomorphic map between preoperative MR and intraoperative CT. The method yields registration accuracy suitable to application in image-guided spine surgery across a broad range of anatomical sites and modes of deformation.

Keywords: deformable image registration, Demons algorithm, symmetric diffeomorphism, multimodality image registration, MIND, CT, MRI, image-guided surgery

1. INTRODUCTION

Spinal disorders are the main cause of disability limiting locomotion, daily activities, and ability to work.1 They cover a broad spectrum of pathologies, such as spinal injury in 25% of trauma patients,2 spine metastases identified in 10% of patients with cancer,3 and scoliosis presenting in 2–32% of adults and 68% of elderly patients.4 Such spinal diseases are treatable by surgery; however, the complexity of spinal structure and function can challenge safe and accurate intervention. Image-guided spine surgery (IGSS) has been shown to improve surgical accuracy and outcomes in pedicle screw placement,5,6 correction of spinal deformities,5,6 trauma surgery,7 percutaneous vertebroplasty,8 and resection of tumors.9 In IGSS, preoperative MR often provides a basis for definition of target anatomy (e.g., vertebral levels and tumors) and critical structures (e.g., nervous and vascular systems), and localization in intraoperative CT requires multimodality registration capable of resolving significant deformation (e.g., supine MR to prone, kyphosed intraoperative CT). We propose a new symmetric diffeomorphic deformable registration method called MIND Demons, which estimates a pair of time-dependent diffeomorphisms as in the SyN10 approach using a Demons-like optimization framework11 under constraints on geodesics, invertibility, and smoothness of the diffeomorphisms. The method demonstrates several advances including: i) a single energy functional that incorporates constraints on smoothness of both the velocity fields and the incremental diffeomorphisms; ii) imposition of a geodesic length constraint satisfying an inverse condition (i.e., the estimated diffeomorphisms are inverse of each other); iii) conservation of momentum in Lagrangian coordinates12 which allows simple and efficient optimization of the diffeomorphisms; iv) use of modality-independent neighborhood descriptors (MIND)13 and a robust and smooth Huber metric14 for non-linear MR-CT image intensity non-correspondence; v) a Gauss-Newton method for fast convergence;15 and vi) estimation of diffeomorphic deformations without relying on the exponential map of stationary velocity fields as in conventional diffeomorphic Demons.11 The registration method is described in Section 2. The figures of merit used to evaluate registration performance are given in Section 3. Analysis of parameter sensitivity and registration performance in comparison to well-established methods is described in Sections 4 and 5, including simulation, physical phantom experiments, and clinical studies. A discussion on advantages and limitations as well as future work is provided in Section 6.

2. MIND DEMONS REGISTRATION

2.1 Notation and Basic Derivation

Registration seeks diffeomorphisms ψ: Ω × t [set membership] [0,1] → Ω on a closed domain Ω [subset or is implied by] Rn (n = 2 or 3 for 2D and 3D, respectively) that nonlinearly map a moving image I0 to a fixed image I1 as I0 [composite function (small circle)] ψ = I1. Since a group of diffeomorphisms is closed under composition,10 one can decompose ψ into a pair of diffeomorphisms ϕi for i [set membership]{0,1} as shown in Fig. 1. The diffeomorphism ϕi(x, t) is a flow of time-dependent velocity fields:

ϕi(x,0.5-(-1)iτ)=ϕi(x,0.5)+0τvi(ϕi(x,0.5-(-1)iu),0.5-(-1)iu)du
(1)

where τ [set membership] [0,0.5] denotes elapsed time, t = (0.5 − (−1)iτ) [set membership] [0,1] denotes a time point, and νi [set membership] L2([0,1], V) is a velocity field in a Sobolev space V.16 The diffeomorphisms ϕi and velocity fields νi are defined with respect to the Lagrangian frame of a virtual image I0.5 domain defined at time 0.5 such that ϕi(0.5) = Id and νi(0.5) = 0. The diffeomorphisms ϕi(t) push the Lagrangian frame (defined at time 0.5) forward to time t. As shown in Fig. 1, ϕ0(x, 0) maps a point x defined at time 0.5 to a point y defined at time 0, and ϕ1(x, 1) maps x to a point z defined at time 1, implying I0 [composite function (small circle)] ϕi(x, 0) = I0.5(x) = I1 [composite function (small circle)] ϕ1(x, 1). For brevity, the time arguments are dropped for the diffeomorphisms at the endpoints as ϕi(x) = ϕi(x, i).

Figure 1
Lagrangian description of the diffeomorphisms (ϕ0, ϕ1, and ψ) and their associated time-dependent velocity fields.17

The energy of the flow ϕi is computed as the time-integrated square norm of νi, Ev=00.5Lvi(0.5-(-1)iτ)22dτ,, where the differential operator L = (Ida2[nabla]2) for a [set membership] R is used to defined a norm on V as ||νi(t)||V = ||Lνi(t)||2.16,17 The energy minimizing path is the optimal geodesic12 of ϕi and the geodesic shortest length – measured using the left-invariance metric (i.e., inverse invariance)18 – is defined in terms of the minimizing energy as:

ρ(ϕi(x),ϕi(x,0.5))=infvi{Ev}
(2)

with a boundary ϕi(x, 0.5) = x and ϕ0(x) = y and ϕ1(x) = z.

A fundamental principle of mechanics states that the Lagrangian (measurement-based) momentum Mi(t) = Lνi(t) is constant over time along optimal geodesics (i.e., energy minimizing paths) in the absence of external forces.12 Since ϕi and νi are defined in the Lagrangian frame, we have Mi(t) = Lνi(t) = Mi(0.5), and since the energy of the flow can be computed in terms of Mi as E(ϕi)=00.5Mi(0.5)-(-1)iτ22,12 we can impose a constraint on the conservation of momentum by using a time step Δt = 0.5. In this work, we use Δt = 0.5 to impose such a constraint and to reduce computational complexity.

2.2 MIND Demons

The method merges the SyN10 approach to the Demons11 approach by optimizing a pair of time-dependent diffeomorphisms ϕi for i [set membership] {0,1} that yields ψ=ϕ0ϕ1-1 under the constraint on the smoothness of ϕi in an alternating optimization enabled through introduction of intermediate diffeomorphisms ηi The alternating optimization allows separation between maximization of image alignment and imposition of a smoothness prior. This separation simplifies the optimization process, allows estimation of various deformation models (e.g., fluid, elastic, and viscoelastic models), and could allow simple integration of additional prior information (e.g., a rigidity constraint on vertebral bodies19); however, it could yield a local optimum that is different from the true optimum of the energy functional. The optimization is performed in the Lagrangian frame with the time step of Δt=0.5 to impose the conservation of momentum. We use Lagrangian push-forward (i.e, diffeomorphisms at each time point are defined with respect to the Lagrangian frame at time 0.5) to bypass the computation of Jacobian change of variables.10,16 Since the method can be used with various similarity metrics and image representations, we present a similarity metric in an abstract form as S(I0, I1, x). The Huber metric and the MIND descriptor representation, used in this work to define S(I0, I1, x), are described in the next section.

(a) Constrained Energy Functional

We constrain the estimated deformations ϕi to be invertible and smooth using an invertibility constraint ϕiϕi-1=Id and minimization of their harmonic energy, respectively. The deformations are, additionally, subject to a geodesic length constraint ρ(ϕ0(x), ϕ1(x, 0.5)) = ρ(ϕ1(x), ϕ1(x, 0.5)) which yields a relation ϕi=ϕj-1 for ij [set membership] {0,1} from the inverse invariance property of the geodesic shortest length18 and the uniqueness of the ordinary differential equation [partial differential]tϕi(x, t) = νi(ϕi(x, t), t) with initial condition ϕi(x, 0.5) = x.20 We use this relation ϕi=ϕj-1 to impose the geodesic length constraint in the energy functional as

E(ϕi,ηi)=12{αS2xΩS(I0η0,I1η1,x)2dx+αU2(ρ(ϕ1-1,η0)2+ρ(ϕ0-1,η1)2)+αP2(ϕ1-122+ϕ0-122)}subjecttoϕi(x,0.5)=ηi(x,0.5)=Id(x)=x,ρ(ϕ0(x),ϕ1(x,0.5))=ρ(ϕ1(x),ϕ1(x,0.5))andϕ1ϕi-1=Id.
(3)

where αS, αU, and αP control regularization strengths, and ηi = ϕi [composite function (small circle)] (Id + νi) for Δt =0.5. The relation ϕi=ϕj-1 and the invertibility constraint ϕiϕi-1=Idlead to the inverse condition ϕi [composite function (small circle)] ϕj = Id (i.e., ϕ0 and ϕ1 are inverse of each other).

(b) Alternating Optimization

The functional (3) is optimized by alternating between two simple steps: (i) maximization of image alignment with imposition of the geodesic length constraint; and (ii) imposition of the smoothness prior and invertibility constraint.

Inexact Image Matching using a Gauss-Newton Method

In the first step, given an estimate of ϕj-1, we seek ηi for ij [set membership] {0,1} that maximizes alignment of I0 and I1. The optimization of ηi is equivalent to minimization of an energy functional defined in terms of νi (i.e., ηi is a function of νi) as

E1(vi)=12{αS2xΩS(I0ϕ0(Id+v0),I1ϕ1(Id+v1),x)2dx+αU2(Lv022+Lv122)}subjecttoϕi(x,0.5)=ηi(x,0.5)=xandρ(ϕ0(x),ϕ0(x,0.5))=ρ(ϕ1(x),ϕ1(x,0.5))
(4)

where ηi=ϕi(Id+vi)=ϕj-1(Id+vi) from ϕi=ϕj-1 (i.e., imposition of the geodesic length constraint), Ii [composite function (small circle)] ηi Ii [composite function (small circle)] ϕi(Id + νi), and ρ(ϕj-1,ηi)2=infvi{00.5Lvi22dt}=infvi{Lvi22} for Δt = 0.5. The functional (4) consists of the first two terms in (3). The first term measures similarity between I0 [composite function (small circle)] ϕ0 and I1 [composite function (small circle)] ϕ1 in the Lagrangian frame, and the second term measures the energy of the diffeomorphisms. We optimize the functional (4) using a Gauss-Newton (GN) method that uses approximations to both first and second–order derivatives of (4) and can generally achieve convergence in fewer iterations than a method using only a first derivative.15 As in the classic Demons method,11 we use the Sherman-Morrison formula to compute inverse of the second-order term in GN, resulting in simplified momenta:

u=-SϕiSαU2/αS2+ϕiS22
(5)

where [nabla]ϕiS is the gradient of S with respect to ϕi and we use αS(x) = 1/S(x) to penalize noise in a spatially varying manner. Using the fact that (αU|S(x)| − ||[nabla]ϕiS(x)||2)2 ≥ 0, we have a constraint on the length of the momenta as ||u(x)||2 ≤ 1/2αU.21 The momentum field equation in (5) is similar to the update equation of the classic Demons method11 except that the regularization strengths in (3) are defined differently. We can retrieve a velocity field from a momentum field using a Green kernel K for L as KMi = KLνi = νi.12 We approximate the Green kernel K for L with a = 1 by a Gaussian kernel GσU with width σU.10 The intermediate diffeomorphisms are updated using the momenta in (5) and the Green kernel K under the geodesic length constraint as

ηik=(ϕjk)-1+Ku(ϕjk)-1=(ϕjk)-1+(GσUu)(ϕjk)-1
(6)

where * is a convolution operator, and k is an iteration number. After both intermediate diffeomorphisms are estimated, the method continues to the second step.

Tikhonov Regularization

The diffeomorphisms ϕi are regularized under smoothness and invertibility constraints by minimizing the energy functional consisting of the last two terms in (3) as

E2(ϕi)=12{(L(ϕ1-1-η0)22-L(ϕ0-1-η1)22)+αP2αU2(ϕ1-122+ϕ0-122)}subjecttoϕiϕi-1=Id.
(7)

where ρ(ϕj-1,ηi)2=infvi{Lvi22}=infvi{L(ϕj-1-ηi)22} from Δt =0.5 and ηi=ϕj-1(Id+vi). By omitting the invertibility constraint, (7) can be optimized using Tikhonov regularization22 which leads to a partial differential equation:

(LL)(ϕj-1-ηi)+αP2αU22ϕj-1=0
(8)

By letting a =0 as in previous work12,17,23 we have L = Id and the solution of (8) can be approximated by ϕj-1GσDηi (i.e., a solution to an isotropic heat equation) where GσD is a Gaussin kernel with width σD=2αP/αU. The invertibility constraint is sequentially imposed by minimizingϕiϕi-1-Id22 using a gradient descent method. We do not use GN here, since it requires inversion of the second-order term (which cannot be simplified using the Sherman-Morrison formula), and the increase in computation time from the required matrix inversion could adversely compromise the convergence rate of GN. After both diffeomorphisms ϕi have been estimated, the diffeomorphic map between I0 and I1 is ψ=ϕ0ϕ1-1.

2.3 Modality Independent Neighborhood Descriptor (MIND)

A MIND descriptor13 is constructed using a non-local mean operator24,25 similar to a self-similarity descriptor,26 and is capable of capturing local structures in MR and CT images. It is a vector representation of each voxel, and its computation involves other voxels in its neighborhood. The configuration of neighboring voxels used in the calculation is called a stencil and is given the symbol [mathematical script N]S.13 Stencils can be arranged in a variety of patterns—e.g., the 2D and 3D examples shown in Fig. 2.

Figure 2
MIND stencil configurations used in this work.

A MIND descriptor for a voxel x in an image I is mI(x) = [mI,i(x)] for i = 1 − |[mathematical script N]S|. Each element mI,i in the descriptor corresponds to a voxel ri [set membership] [mathematical script N]S in the stencil, and its value is computed as:

mI,i=cexp(-d(I,x,ri)V(I,x))
(9)

where c is a normalization factor such that maxi[1,NS]{mI,i(x)}=1 (i.e., for intensity and contrast invariance) and d(I, x, ri) measures the distance between a patch of x and that of ri as

d(I,x,ri)=zNpGσp(z)(I(x+z)-I(ri+z))2
(10)

where [mathematical script N]p denotes a neighborhood configuration of a patch (e.g., a cube), z is an offset from the center voxel in [mathematical script N]p, and Gσp is a discrete Gaussian kernel— with width σp mm and truncation (tail cut-off) errors tp = 1 − Σz[set membership][mathematical script N]p Gσp(z)—used to increase the importance of the central voxel. The term V(I, x) in (9) approximates the local variance at x by an average of patch distances between the patch of x and patches of the nearest neighboring voxels of x (i.e., 4 and 6 nearest neighbors for 2D and 3D images, respectively).

2.4 Huber Distance Metric

Registration is generally ill-posed27 implying that small changes in images (e.g., due to noise and artifacts) can lead to large changes in the estimated deformations. A metric such as the L1 norm tends to be relatively insensitive to noise and outliers28 and could yield more reliable estimation (i.e., less confounded by noise) than quadratic norms (e.g., L2 and square norm); however, numerical optimization involving the L1 criterion is difficult due to its singularity. To exploit the robustness of L1 with the twice differentiability of L2, we use a Huber distance14 as the similarity metric between a MIND descriptor of I0 [composite function (small circle)] ϕ0, mI0[composite function (small circle)]ϕ0(x), and that of I1 [composite function (small circle)] ϕ1, mI1[composite function (small circle)]ϕ1(x). Note that mIi[composite function (small circle)]ϕi is computed on images after transformation Ii [composite function (small circle)] ϕi and mIi[composite function (small circle)]ϕimIi[composite function (small circle)]ϕi. We denote the difference between an element i of the MIND descriptors as [var phi]i(x) = mI0[composite function (small circle)]ϕ0i(x) − mI1[composite function (small circle)]ϕ1, i. The Huber distance between the descriptors is:

S(mI0ϕ0,mI1ϕ1,x)=i=1NS{φi2(x)/2ε,if0φi(x)εφi(x)-(ε/2),otherwise
(11)

where ε is the threshold between the quadratic and linear parts.

The Huber metric yields a more reliable deformation estimation, since the influence of outliers on gradients of the metric is less than that on gradients of quadratic norms. Moreover, edges in images can be preserved since it penalizes large differences less than the quadratic penalty in the L2 norm

2.5. MIND Demons Registration Method

We incorporate a multiresolution strategy to improve robustness against local minima. As shown in Fig. 3, a multiresolution pyramid (defining coarse-to-fine evolution in each pyramid level) is constructed only for the virtual image I0.5 since the method uses the Lagrangian push-forward of the I0.5 domain to the domain of I0 and I1 using ϕ0 and ϕ1, respectively. In each level of the pyramid, the optimization described in Section 2.2 is performed until it reaches either convergence or a maximum number of iterations. Since a MIND descriptor is not deformation invariant, it is recomputed in every optimization iteration. The parameters intrinsic to the MIND Demons method along with their nominal ranges and values (see Section 4.1 and 5.1) are given in Table 1. The algorithm was implemented using the Insight Segmentation and Registration Toolkit (ITK).29

Figure 3
Flowchart for the MIND Demons algorithm.
Table 1
MIND Demons Parameters and Nominal Settings

3. REGISTRATION PERFORMANCE

Registration performance was quantified in terms of geometric accuracy [target registration error (TRE) and structural similarity metric (SSIM)] and diffeomorphic properties [invertibility (J) and preservation of topology] of the estimated deformations.

3.1. Target Registration Error (TRE)

Geometric accuracy of registration was measured in terms of TRE – the distance between corresponding anatomical points in I1 and I0 after registration:

TRE(x0,x1;ψ)=x0-ψ(x1)2
(12)

where xi denotes an unambiguous anatomical point (“target”) in Ii, and ψ is the estimated deformation.

3.2. Structural Similarity Metric (SSIM)

The difference between images from a similar modality can be quantified in terms of SSIM,30 computed using local statistical features as:

SSIM(y,z)=(2μyμz+C1)(2σyz+C2)(μy2+μz2+C1)(σy2+σz2+C2)
(13)

where y denotes a voxel in I0 after registration, z denotes a corresponding voxel in I1, μi is a local mean computed within a sliding window of size Nw centered on i, σi2 denotes a local variance within the window, and σyz is a local covariance in the window centered on y and z. Prior to computation of local statistical features, a Gaussian kernel of size Nw was applied to images to avoid blocking artifacts from the rectangular sliding window. In this work, we set Nw = 11×11 pixels, C1 = (0.01R)2, and C2 = (0.03R)2, where R is the dynamic range of the pixel values.30 In the simulations (below), SSIM was computed with respect to the deformed moving image (I0 after registration) and a known true image (a simulated fixed MR image) of the same modality (i.e., image intensity).

3.3. Invertibility (J)

A desirable characteristic of diffeomorphisms as described above is their invertibility, which can be characterized in terms of the residual:

J(ψ,ψ-1;y,z)=12{(ψψ-1)(y)-y2+(ψ-1ψ)(z)-z2}
(14)

where y is a point in I0, z is a point in I1, ψ is a diffeomorphism, and ψ−1 is its inverse. The residual J = 0 for invertible or singular transformations.

3.4. Minimum of Jacobian Determinant (D)

The preservation of topology and lack of folding/tearing31 in a given deformation can be characterized in terms of its Jacobian determinant, J(x) = det(Dxψ(x)) where Dxψ are spatial Jacobians of ψ, and det(·) denotes a matrix determinant. A deformation is nonsingular if J ≠ 0 and preserves topology if J > 0 (i.e., local reflection occurs for J < 0, which could lead to self-intersection). To quantify invertibility and preservation of orientation as well as topology, we measure the minimum value of J within a sliding window ΩV as

D(ψ)=minxΩV{J(x)}
(15)

where the size of the sliding window ΩV is set to 11n voxels (for n = 2 or 3 for 2D and 3D, respectively).

4. EXPERIMENTAL METHODS

The registration performance of MIND Demons was analyzed in comparison to other well-established reference methods, including elastix free-form deformation (FFD) with MI and local MI (LMI),32,33 and NMI Demons34,35 with a symmetric energy formulation. Each method was evaluated using nominal parameter settings identified through analysis of parameter sensitivity in simulation. Comparison of the overall registration performance was performed in a 3D physical phantom experiment (ovine spine phantom), and validation of the registration performance under realistic imaging conditions was performed in clinical studies using MR and CT images of the cervical, thoracic, and lumbar spine.

4.1. Simulation Studies

(a) Formation of Simulated Images

Figures 4(a–b) depict a simulated T2-weighted MR moving image (I0) and a CT fixed image (I1) (1500×1500 pixels, 0.25 mm2). Each 2D simulated image contains 5 simulated (rigid) vertebrae within (deformable) soft-tissue surroundings. The scoliotic curvature in I0 (Fig. 4(a)) was generated from a radial basis interpolation of a rigid motion36 associated with each vertebrae in I1. The MR pixel value was computed using T2 and proton density of tissue for a spin-echo pulse sequence with a constant main magnetic field at 1.5 T.37 As shown in Fig. 4, 31 target points were defined in each image at the corners of the vertebrae, the centroid of the tumors, and the unambiguous points in soft-tissue for TRE measurement.

Figure 4
2D Simulated images emulating simple coronal curvature (scoliosis) of the spine, each with 31 target points. (a) T2-weighted MR moving image (I0). (b) CT fixed image (I1).

(b) Analysis of Parameter Sensitivity and Registration Performance

The sensitivity to individual parameter settings in each method was investigated in univariate analysis to identify nominal parameter range and value. The nominal parameter values were subsequently used in evaluating the performance of each method in comparison to MIND Demons. A similar morphological pyramid was constructed for the Demons-based methods using Gaussian kernel widths of [16, 8, 4, 2, 1] pixels and downsampling factors of [32, 16, 8, 4, 2] pixels, while a pyramid for the FFD methods was constructed using only Gaussian smoothing.

4.2. Physical Phantom Studies (Ovine Spine)

(a) Image Acquisition

As shown in Fig. 5, an ovine spine was enclosed in a MR-CT compatible and bendable plastic cylinder filled with polyvinyl alcohol (to simulate soft-tissue). The phantom was imaged first in a preoperative setup with scoliotic curvature for T2-weighted MR and CT moving images (I0), followed by T2-weighted MR and intraoperative CT fixed images (I1) with the spine straightened to natural posture. The MR scans were acquired with 3D acquisition type on a 1.5 T Magnetom Avanto (Siemens Healthcare, Malvern PA) and one echo with TE = 125 ms. An average of two acquisitions was used to reduce noise in the images. The MR I0 was reconstructed at 0.9×0.9×0.9 mm3 with a size of 192×384×128 voxels, and the MR I1 was reconstructed at 0.5×0.5×0.9 mm3 with a size of 192×384×144 voxels. The CT images were acquired with a Somatom Definition Flash scanner (Siemens Healthcare, Erlangen, Germany) (100 kVp, 291 mAs) and reconstructed at 0.6×0.6×0.8 mm3 with a size of 256×256×312 voxels. Example sagittal slices of MR and CT images are shown in Figs. 5(b, c). The MR I0 and CT I1 images were used as the moving and fixed images, respectively, in the following phantom experiments. For visualization and target point definition (TRE calculation), both CT images (I0 and I1) were segmented (simple thresholding at 200 HU), and 32 target points were defined on unambiguous anatomical features (tips of the spinous and transverse processes). The target points and segmentation defined in the CT I0 were translated to those defined in the MR I0 using NMI rigid registration. The MR I1 was visually compared to the MR I0 after nonlinear transformation.

Figure 5
Ovine spine phantom. (a) Phantom assembly – spine encased in polyvinyl alcohol (PVA) hydrogel within a flexible cylinder. (b) CT images with the scoliotic spine (I0) and the straight spine (I1). (c) T2-weighted MR images with the scoliotic spine ...

(b) Assessment of MR-to-CT Registration Performance

The registration performance of MIND Demons was evaluated in comparison to that of MI FFD, LMI FFD, and NMI Demons implemented using the nominal parameter settings. The three-level image pyramids for the Demons-based methods were constructed with Gaussian kernel widths of [4, 2, 1] voxels and downsampling factors of [8, 4, 2] voxels, while those for the FFD-based methods were constructed using only Gaussian smoothing without downsampling.

4.3. Clinical Studies

An institutional review board (IRB) approved retrospective study was performed to validate the registration performance of MIND Demons for clinically realistic patient images. The study used three pairs of T2-weighted MR and CT acquired for three patients undergoing intervention of cervical, thoracic, and lumbar disorders at our institution.

(a) Image Acquisition

The T2-weighted MR images and their corresponding CT images (for the cervical, thoracic, and lumbar spines) exhibit realistic variations in imaging protocols and image quality (Fig. 6). The MR scans were acquired with 2D (sagittal-slice) acquisition type on a 3T Signa HDxt (GE Healthcare, Little Chalfont, UK), a 1.5T Aera (Siemens Healthcare, Erlangen, Germany), or a 1.5T Vantage Titan (Toshiba Corporation, Tokyo, Japan) with slice thickness of ~3 mm and TE varied from 100–120 ms. The CT images were acquired using a LightSpeed Ultra scanner (GE Healthcare, Little Chalfont, UK) or a Somatom Definition Flash scanner (Siemens Healthcare, Erlangen, Germany) with scan techniques varied from 120–140 kVp and 80–165 mAs, and reconstructed at approximately 0.3×0.3×0.5 mm3. The MR images of the cervical, thoracic, and lumbar spines were used as moving images I0 (Figs. 6(a, c, e)) and their corresponding CT images were used as fixed images I1 (Figs. 6(b, d, f)). Example sagittal slices of the MR and CT images are shown in Fig. 6. For visualization and target point definition, the vertebrae in MR were manually segmented, and those in the CT images were segmented using simple bone-thresholding after median filtering. Twelve, eight, and eleven target points were identified for the cervical, thoracic, and lumbar spine, respectively.

Figure 6
Clinical MR and CT image data. (a) T2-weighted MR (I0) and (b) CT (I1) images of the cervical spine. (c, d) The same, in the thoracic spine. (e, f) The same, in the lumbar spine.

(b) Validation of MR-to-CT Registration Performance

The registration performance of MIND Demons with clinical data was evaluated using the nominal parameter settings. The studies used a four-level morphological pyramid with Gaussian kernel widths of [8, 4, 2, 1] voxels and downsampling factors of [16, 8, 4, 2] voxels. MIND Demons was initialized using NMI rigid registration.

5. RESULTS

5.1. Simulation Studies

Table 2 and Figure 7 summarize results from the simulation studies using the nominal parameter settings for each method. All methods were initialized with a rigid transformation with mean TRE 13.0 mm and interquartile range (IQR) of (9.0, 15.9) mm, and registration accuracy of each method is summarized in Table 2 and Fig. 7. The p-values (Table 2) measure statistical significance in the differences between the measured mean values of each metric from that of MIND Demons. All methods were found to preserve topology (lack of tissue folding and tearing) with D > 0. MIND Demons and NMI Demons yielded invertible deformations with sub-voxel J; however, NMI Demons achieved smaller J than MIND Demons due to stronger smoothing applied to the update and displacement fields. The invertibility J was not measured for the FFD-based methods since the methods only provided forward deformations (i.e., a map from I0 to I1). MIND Demons yielded mean TRE = 1.3 mm, compared to 2.1 mm for MI FFD, 3.9 mm for LMI FFD, and 2.0 mm for NMI Demons. Additionally, MIND Demons outperformed the others in terms of SSIM. Figure 7 shows I0 after registration. The cyan spheres represent the target points in I1 and the green lines mark the shortest distances between the target points after registration. MIND Demons demonstrated the ability to resolve large deformations with robustness against spurious distortion that is evident in the MI FFD, LMI FFD, and NMI Demons methods.

Figure 7
Simulation study results. Transformed I0 image following (a) MI FFD, (b) LMI FFD, (c) NMI Demons, and (d) MIND Demons. Note the reduced TRE (alignment of cyan target points) and robustness against spurious distortion for the MIND method.
Table 2
Summary of Registration Results in Simulation Studies

5.2. Physical Phantom Studies

For the ovine spine phantom studies, each registration method was initialized with a rigid transformation with mean TRE 6.0 mm and IQR of (3.7, 7.5) mm. Table 3 and Figures (8, ,9)9) summarize the overall registration performance of MIND Demons compared to the reference methods. Figure 8(b) demonstrates sub-voxel J (mean J = 0.008 mm) with a small range and IQR for MIND Demons compared to that for NMI Demons. All of the methods were capable of preserving structural topology with D > 0; however, MI FFD yielded min D close to 0 (Table 3 and Fig. 8(c)). The large ranges in D associated with MI FFD and MIND Demons reveal a large change in local volume (i.e., expansion and compression) from large local motions. Table 3, Figures 8(a) and and99 summarize the registration accuracy of each method. The top row in Fig. 9 depicts semi-opaque overlays of the pink MR I0 and the cyan CT I1 after registration. The cyan spheres represent the target points in I1 and the yellow line segments mark the distance between the corresponding target points after registration. The bottom row shows the yellow Canny-edges of the MR I0 after registration superimposed on the gray MR I1. MIND Demons achieved statistically significant improvement in registration accuracy with mean TRE of 1.5 mm.

Figure 8
MR-to-CT registration of the ovine spine phantom. (a–c) TRE, J, and D as a function of registration methods.
Figure 9
MR-to-CT registration of the ovine spine phantom. (Top) Semi-opaque surface rendering of the pink MR moving image I0 and the cyan fixed CT image I1 after registration. Cyan spheres represent the target points in I1 and yellow lines mark distances between ...
Table 3
Summary of Registration Results in Phantom Studies

5.3. Clinical Studies

Figures 10 and and1111 summarize the registration performance of MIND Demons in clinical studies using MR and CT images of the cervical (C), thoracic (T), and lumbar (L) spines. As illustrated in Figs. 10(b, c), MIND Demons yielded invertible deformations with preservation of topology and sub-voxel values of J and D >0. Figures 10(a) and 11 summarize the registration accuracy of MIND Demons. The top row in Fig. 11 depicts semi-opaque overlays of the pink MR I0 and the cyan CT I1 after registration for each spinal section, each section with results from NMI rigid and MIND Demons registration. The cyan spheres represent the target points in I1 and the yellow line segments mark the distance between the corresponding target points after registration. The middle row depicts the yellow Canny edges of the CT I1 superimposed on the gray MR I0 after registration, and the bottom row shows checkerboard images between the CT I1 and the MR I0 after registration. MIND Demons was able to resolve deformation induced by variation in patient positioning and provided improved registration accuracy, with TRE improved from 3.3±1.2 mm after rigid registration to 1.6±0.6 mm after MIND Demons for the C spine, from 4.4±1.8 mm to 1.7±0.6 mm for the T spine, and from 4.3±1.7 mm to 1.9±0.5 mm for the L spine.

Figure 10
MR-to-CT registration in clinical studies. (a–c) TRE, J, and D as a function of the spinal sections: cervical, thoracic, and lumbar spines.
Figure 11
MR-to-CT registration in clinical studies. (Top) Semi-opaque surface rendering of the pink MR moving image I0 and the cyan fixed CT image I1 after NMI rigid and MIND Demons registration for each spinal section. Cyan spheres represent the target points ...

6. DISCUSSION AND CONCLUSION

A deformable registration method merging the Demons11 and SyN10 approaches for symmetric time-dependent diffeomorphisms has been developed. The algorithm incorporates MIND descriptors13 and the Huber metric14 for robust multimodality registration and the Gauss-Newton approach for fast convergence.15 The sensitivity analysis showed that the Huber metric with a small quadratic region (i.e., ε = 0.001–0.01) was able to reject outliers from local structural differences captured by corresponding MIND descriptors and provide reliable estimation of the deformation. Locality of the MIND descriptor—determined through the configuration of the stencil and the patches (i.e., values of σp and tp)—led to robustness against distortion, and its patch-based computation reduced sensitivity to image noise. Viscoelastic deformations, with adjustable strength of fluid and elastic models, were able to resolve large deformation in realistically noisy data and attain accurate registration (sub-voxel TRE < 2.0 mm in clinical studies – within the range required for spinal intervention38,39). The estimated deformation was diffeomorphic to the extent that topology was preserved with sub-voxel J < 0.008 mm and D > 0.

MI-based methods have been somewhat widely used for MR-to-CT volumetric registration; however, MI-based metrics are sensitive to intensity non-uniformity, and they lose robustness when used as local measures.40 Incorporation of spatial information40,41 and/or image features42,43 could improve their sensitivity to intensity distortion. However, due to the challenge associated with the assumption of tissue-class correspondence, we adopted MIND descriptors. The descriptors are not deformation invariant, but this could increase their discriminative power.44

A diffeomorphism is a bijective map; it therefore assumes consistent anatomical structures to appear in both images. This inhibits applications of the method to resolve deformation involving content mismatch (e.g., due to insertion and removal of surgical tools/implants and/or tissue resection). A method using an asymmetric energy formulation of MIND Demons and a geometrical penalty term45 or extra-dimension46 to handle added and/or missing structures could be applicable in this case. Within the current implementation, we have demonstrated registration accuracy suitable to registration of preoperative MR and intraoperative CT before the delivery of surgical implants – i.e., at the beginning of the case, allowing localization of target and critical structures. Future work will consider registration to intraoperative images featuring a high density of implants (e.g., spine screws).

In this work, the time step used in the integration of time-dependent velocity fields was fixed at 0.5 to impose the conservation of momentum in Lagrangian coordinates.12 The smaller time step might yield improved registration accuracy, but lack of the constraint and more importantly increase in computational complexity and expense is discouraging. The clinical studies used three image pairs for the cervical, thoracic, and lumbar spines. An IRB-approved retrospective study using a database of clinical data to further validate application of the method in realistic clinical setups is a subject of future work. Since our main objective is application in intraoperative use, application of the method to intraoperative CBCT images is also to be investigated. Owing to voxel-wise computation of the algorithm, distributed and/or parallel computing can be used to improve computational time. Future work can additionally include application of deformation-invariant descriptors, sensitivity analysis to image noise and intensity distortion, as well as incorporation of other prior knowledge, such as rigidity of the vertebrae to further constrain the solution space.

Acknowledgments

This work was supported in part by the National Institutes of Health grant number R01 -EB-017226, collaboration with Siemens Healthcare XP, and the Thai Royal Government Scholarship. The authors gratefully acknowledge Dr. Adam Wang (Biomedical Engineering, Johns Hopkins University) and Dr. Amir Pourmorteza (Biomedical Engineering, Johns Hopkins University) for assistance with the physical phantom assembly and the MR scans.

References

1. Haldeman S, Kopansky-Giles D, Hurwitz EL, Hoy D, Mark Erwin W, Dagenais S, Kawchuk G, Strömqvist B, Walsh N. Advancements in the Management of Spine Disorders. Best Pract Res Clin Rheumatol. 2012;26(2):263–280. [PubMed]
2. Whang PG, Patel AA, Vaccaro AR. The development and evaluation of the subaxial injury classification scoring system for cervical spine trauma. Clin Orthop Relat Res. 2011;469(3):723–731. 2010/09/22 ed. [PMC free article] [PubMed]
3. Delank KS, Wendtner C, Eich HT, Eysel P. The treatment of spinal metastases. Dtsch Arztebl Int. 2011;108(5):71–79. 2011/02/12 ed. [PubMed]
4. Smith JS, Shaffrey CI, Glassman SD, Berven SH, Schwab FJ, Hamill CL, Horton WC, Ondra SL, Sansur CA, et al. Risk-benefit assessment of surgery for adult scoliosis: an analysis based on patient age. Spine (Phila Pa 1976) 2011;36(10):817–824. 2010/08/05 ed. [PubMed]
5. Larson AN, Polly DWJ, Guidera KJ, Mielke CH, Santos ERG, Ledonio CGT, Sembrano JN. The Accuracy of Navigation and 3D Image-Guided Placement for the Placement of Pedicle Screws in Congenital Spine Deformity. J Pediatr Orthop. 2012;32(6):e23–e29. [PubMed]
6. Flynn J, Sakai D. Eur Spine J. 2. Vol. 22. Springer-Verlag; 2013. Improving safety in spinal deformity surgery: advances in navigation and neurologic monitoring; pp. 131–137. [PMC free article] [PubMed]
7. Schouten R, Lee R, Boyd M, Paquette S, Dvorak M, Kwon BK, Fisher C, Street J. Intra-operative cone-beam CT (O-arm) and stereotactic navigation in acute spinal trauma surgery. J Clin Neurosci. 2012;19(8):1137–1143. 2012/06/23 ed. [PubMed]
8. Yang C-D, Chen Y-W, Tseng C-S, Ho H-J, Wu C-C, Wang K-W. Non-invasive, fluoroscopy-based, image-guided surgery reduces radiation exposure for vertebral compression fractures: A preliminary survey. Formos J Surg. 2012;45(1):12–19.
9. Bandiera S, Ghermandi R, Gasbarrini A, Barbanti Bròdano G, Colangeli S, Boriani S. Eur Spine J. 6. Vol. 22. Springer; Berlin Heidelberg: 2013. Navigation-assisted surgery for tumors of the spine; pp. 919–924. [PMC free article] [PubMed]
10. Avants BB, Epstein CL, Grossman M, Gee JC. Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Med Image Anal. 2008;12(1):26–41. 2007/07/31 ed. [PMC free article] [PubMed]
11. Vercauteren T, Pennec X, Perchant A, Ayache N. Diffeomorphic demons: efficient non-parametric image registration. Neuroimage. 2009;45(1 Suppl):S61–S72. 2008/12/02 ed. [PubMed]
12. Miller MI, Trouvé A, Younes L. Geodesic Shooting for Computational Anatomy. J Math Imaging Vis. 2006;24(2):209–228. [PMC free article] [PubMed]
13. Heinrich MP, Jenkinson M, Bhushan M, Matin T, Gleeson FV, Brady SM, Schnabel JA. MIND: Modality independent neighbourhood descriptor for multi-modal deformable registration. Med Image Anal. 2012;16(7):1423–1435. [PubMed]
14. Huber PJ. Ann Stat. 5. Vol. 1. Institute of Mathematical Statistics; 1973. Robust Regression: Asymptotics, Conjectures and Monte Carlo; pp. 799–821.
15. Hernandez M, Olmos S. Gauss-Newton optimization in Diffeomorphic registration. 5th IEEE Int Symp Biomed Imaging From Nano to Macro. 2008:1083–1086.
16. Beg MF, Miller M, Trouvé A, Younes L. Int J Comput Vis. 2. Vol. 61. Kluwer Academic Publishers; 2005. Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms; pp. 139–157.
17. Miller MI, Trouve A, Younes L. On the metrics and euler-lagrange equations of computational anatomy. Annu Rev Biomed Eng. 2002;4:375–405. 2002/07/16 ed. [PubMed]
18. Miller MI, Younes L. Int J Comput Vis. 1–2. Vol. 41. Kluwer Academic Publishers; 2001. Group Actions, Homeomorphisms, and Matching: A General Framework; pp. 61–84.
19. Reaungamornrat S, Wang AS, Uneri A, Otake Y, Khanna AJ, Siewerdsen JH. Deformable image registration with local rigidity constraints for cone-beam CT-guided spine surgery. Phys Med Biol. 2014;59(14):3761–3787. 2014/06/18 ed. [PMC free article] [PubMed]
20. Younes L. Shapes and Diffeomorphisms. Springer Berlin Heidelberg; Berlin, Heidelberg: 2010.
21. Cachier P, Pennec X, Ayache N. Algorithm. 1999. Fast Non Rigid Matching by Gradient Descent: Study and Improvements of the ‘Demons’
22. Mansi T, Pennec X, Sermesant M, Delingette H, Ayache N. Int J Comput Vis. 1. Vol. 92. Springer US; 2011. iLogDemons: A Demons-Based Registration Algorithm for Tracking Incompressible Elastic Biological Tissues; pp. 92–111.
23. Arnold V. Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications àl’hydrodynamique des fluides parfaits. Ann l’institut Fourier. 1966
24. Buades A, Coll B, Morel JM. A non-local algorithm for image denoising. IEEE Comput Soc Conf Comput Vis Pattern Recognit (CVPR) 2005;2:60–65.
25. Gilboa G, Osher S. Nonlocal Operators with Applications to Image Processing. Multiscale Model Simul. 2009;7(3):1005–1028.
26. Shechtman E, Irani M. Matching Local Self-Similarities across Images and Videos. IEEE Comput Soc Conf Comput Vis Pattern Recognit (CVPR) 2007:1–8.
27. Bernd F, Jan M. Ill-posed medicine—an introduction to image registration. Inverse Probl. 2008;24(3):34008.
28. Wedel A, Pock T, Zach C, Bischof H, Cremers D. An Improved Algorithm for TV-L 1 Optical Flow. In: Cremers D, Rosenhahn B, Yuille A, Schmidt F, editors. Statistical and Geometrical Approaches to Visual Motion Analysis. Springer; Berlin Heidelberg: 2009. pp. 23–45.
29. Yoo TS, Ackerman MJ, Lorensen WE, Schroeder W, Chalana V, Aylward S, Metaxas D, Whitaker R. Engineering and Algorithm Design for an Image Processing API: A Technical Report on ITK - the Insight Toolkit. In: Westwood JD, editor. Proc Med Meets Virtual Real. Vol. 85. IOS Press; Amsterdam: 2002. pp. 586–592. [PubMed]
30. Zhou W, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process. 2004;13(4):600–612. [PubMed]
31. Yang X, Xue Z, Liu X, Xiong D. Topology preservation evaluation of compact-support radial basis functions for image registration. Pattern Recogn Lett. 2011;32(8):1162–1177.
32. Klein S, van der Heide UA, Lips IM, van Vulpen M, Staring M, Pluim JPW. Automatic segmentation of the prostate in 3D MR images by atlas matching using localized mutual information. Med Phys. 2008;35(4):1407–1417. [PubMed]
33. Klein S, Staring M, Murphy K, Viergever MA, Pluim JPW. elastix: A Toolbox for Intensity-Based Medical Image Registration. IEEE Trans Med Imaging. 2010;29(1):196–205. [PubMed]
34. Gholipour A, Kehtarnavaz N, Yousefi S, Gopinath K, Briggs R. Symmetric deformable image registration via optimization of information theoretic measures. Image Vis Comput. 2010;28(6):965–975.
35. Ying L, Yonggang S, Fa J, Zhiwen L, Yong Y. Multi-modal diffeomorphic demons registration based on mutual information. 4th Int Conf Biomed Eng Informatics (BMEI) 2011;2:800–804.
36. Moreno A, Chambon S, APS, JPR, Angelini E, Bloch I. Combining a breathing model and tumor-specific rigidity constraints for registration of CT-PET thoracic data. Comput Aided Surg. 2008;13(5):281–298. 2008/09/30 ed. [PubMed]
37. Kato H, Kuroda M, Yoshimura K, Yoshida A, Hanamoto K, Kawasaki S, Shibuya K, Kanazawa S. Composition of MRI phantom equivalent to human tissues. Med Phys. 2005;32(10):3199–3208. [PubMed]
38. Shahidi R, Clarke L, Bucholz RD, Fuchs H, Kikinis R, Robb RA, Vannier MW. Comput Aided Surg. 3. Vol. 6. John Wiley & Sons, Inc; 2001. White paper: Challenges and opportunities in computer-assisted interventions January 2001; pp. 176–181. [PubMed]
39. Cleary K, Anderson J, Brazaitis M, Devey G, DiGioia A, Freedman M, Grönemeyer D, Lathan C, Lemke H, et al. Comput Aided Surg. 3. Vol. 5. John Wiley & Sons, Inc; 2000. Final report of the Technical Requirements for Image-Guided Spine Procedures Workshop*, April 17–20, 1999, Ellicott City, Maryland, USA; pp. 180–215. [PubMed]
40. Xiahai Z, Arridge S, Hawkes DJ, Ourselin S. A Nonrigid Registration Framework Using Spatially Encoded Mutual Information and Free-Form Deformations. IEEE Trans Med Imaging. 2011;30(10):1819–1828. [PubMed]
41. Loeckx D, Slagmolen P, Maes F, Vandermeulen D, Suetens P. Nonrigid Image Registration Using Conditional Mutual Information. IEEE Trans Med Imaging. 2010;29(1):19–29. [PubMed]
42. Jonghye W, Stone M, Prince JL. Multimodal Registration via Mutual Information Incorporating Geometric and Spatial Context. IEEE Trans Image Process. 2015;24(2):757–769. [PMC free article] [PubMed]
43. Rivaz H, Karimaghaloo Z, Collins DL. Self-similarity weighted mutual information: A new nonrigid image registration metric. Med Image Anal. 2014;18(2):343–358. [PubMed]
44. Bay H, Ess A, Tuytelaars T, Van Gool L. Speeded-Up Robust Features (SURF) Comput Vis Image Underst. 2008;110(3):346–359.
45. Berendsen FF, Kotte AN, de Leeuw AA, Jurgenliemk-Schulz IM, Viergever MA, Pluim JP. Registration of structurally dissimilar images in MRI-based brachytherapy. Phys Med Biol. 2014;59(15):4033–4045. 2014/07/06 ed. [PubMed]
46. Nithiananthan S, Schafer S, Mirota DJ, Stayman JW, Zbijewski W, Reh DD, Gallia GL, Siewerdsen JH. Extra-dimensional Demons: A method for incorporating missing tissue in deformable image registration. Med Phys. 2012;39(9):5718–5731. [PubMed]