Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC4909372

Formats

Article sections

- Abstract
- 1. INTRODUCTION
- 2. MIND DEMONS REGISTRATION
- 3. REGISTRATION PERFORMANCE
- 4. EXPERIMENTAL METHODS
- 5. RESULTS
- 6. DISCUSSION AND CONCLUSION
- References

Authors

Related links

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.2208621PMCID: PMC4909372

NIHMSID: NIHMS780672

S. Reaungamornrat,^{a} T. De Silva,^{b} A. Uneri,^{a} J.-P. Wolinsky,^{c} A. J. Khanna,^{b,}^{d} G. Kleinszig,^{e} S. Vogt,^{e} J. L. Prince,^{b,}^{f} and J. H. Siewerdsen^{a,}^{b}

See other articles in PMC that cite the published article.

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.

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.

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.

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.

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 SyN^{10} approach using a Demons-like optimization framework^{11} 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 coordinates^{12} which allows simple and efficient optimization of the diffeomorphisms; iv) use of modality-independent neighborhood descriptors (MIND)^{13} and a robust and smooth Huber metric^{14} 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.

Registration seeks diffeomorphisms ** ψ**: Ω ×

$${\mathit{\varphi}}_{i}(\mathit{x},0.5-{(-1)}^{i}\tau )={\mathit{\varphi}}_{i}(\mathit{x},0.5)+\underset{0}{\overset{\tau}{\int}}{\mathit{v}}_{i}({\mathit{\varphi}}_{i}(\mathit{x},0.5-{(-1)}^{i}u),0.5-{(-1)}^{i}u)\phantom{\rule{0.16667em}{0ex}}du$$

(1)

where *τ* [0,0.5] denotes elapsed time, *t* = (0.5 − (−1)* ^{i}τ*) [0,1] denotes a time point, and

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

$$\rho ({\mathit{\varphi}}_{i}(\mathit{x}),{\mathit{\varphi}}_{i}(\mathit{x},0.5))=\underset{{\mathit{v}}_{i}}{inf}\{\sqrt{{E}_{v}}\}$$

(2)

with a boundary *ϕ** _{i}*(

A fundamental principle of mechanics states that the Lagrangian (measurement-based) momentum *M** _{i}*(

The method merges the SyN^{10} approach to the Demons^{11} approach by optimizing a pair of time-dependent diffeomorphisms *ϕ** _{i}* for

We constrain the estimated deformations *ϕ** _{i}* to be invertible and smooth using an invertibility constraint
${\mathit{\varphi}}_{i}\circ {\mathit{\varphi}}_{i}^{-1}=Id$ and minimization of their harmonic energy, respectively. The deformations are, additionally, subject to a geodesic length constraint

$$\begin{array}{l}E({\mathit{\varphi}}_{i},{\mathit{\eta}}_{i})=\frac{1}{2}\left\{{\alpha}_{S}^{2}\underset{\mathit{x}\in \mathrm{\Omega}}{\int}S{({I}_{0}\circ {\mathit{\eta}}_{0},{I}_{1}\circ {\mathit{\eta}}_{1},\mathit{x})}^{2}d\mathit{x}+{\alpha}_{U}^{2}(\rho {({\mathit{\varphi}}_{1}^{-1},{\mathit{\eta}}_{0})}^{2}+\rho {({\mathit{\varphi}}_{0}^{-1},{\mathit{\eta}}_{1})}^{2})+{\alpha}_{P}^{2}({\Vert \nabla {\mathit{\varphi}}_{1}^{-1}\Vert}_{2}^{2}+{\Vert \nabla {\mathit{\varphi}}_{0}^{-1}\Vert}_{2}^{2})\right\}\\ \text{subject}\phantom{\rule{0.16667em}{0ex}}\text{to}\phantom{\rule{0.16667em}{0ex}}{\mathit{\varphi}}_{i}(\mathit{x},0.5)={\mathit{\eta}}_{i}(\mathit{x},0.5)=Id(\mathit{x})=\mathit{x},\rho ({\mathit{\varphi}}_{0}(\mathit{x}),{\mathit{\varphi}}_{1}(\mathit{x},0.5))=\rho ({\mathit{\varphi}}_{1}(\mathit{x}),{\mathit{\varphi}}_{1}(\mathit{x},0.5))\phantom{\rule{0.16667em}{0ex}}\text{and}\phantom{\rule{0.16667em}{0ex}}{\mathit{\varphi}}_{1}\circ {\mathit{\varphi}}_{i}^{-1}=Id.\end{array}$$

(3)

where *α _{S}*,

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.

In the first step, given an estimate of
${\mathit{\varphi}}_{j}^{-1}$, we seek *η** _{i}* for

$$\begin{array}{l}{E}_{1}({\mathit{v}}_{i})=\frac{1}{2}\left\{{\alpha}_{S}^{2}\underset{\mathit{x}\in \mathrm{\Omega}}{\int}S{({I}_{0}\circ {\mathit{\varphi}}_{0}(Id+{\mathit{v}}_{0}),{I}_{1}\circ {\mathit{\varphi}}_{1}(Id+{\mathit{v}}_{1}),\mathit{x})}^{2}d\mathit{x}+{\alpha}_{U}^{2}({\Vert L{\mathit{v}}_{0}\Vert}_{2}^{2}+{\Vert L{\mathit{v}}_{1}\Vert}_{2}^{2})\right\}\\ \text{subject}\phantom{\rule{0.16667em}{0ex}}\text{to}\phantom{\rule{0.16667em}{0ex}}{\mathit{\varphi}}_{i}(\mathit{x},0.5)={\mathit{\eta}}_{i}(\mathit{x},0.5)=\mathit{x}\phantom{\rule{0.16667em}{0ex}}\text{and}\phantom{\rule{0.16667em}{0ex}}\rho ({\mathit{\varphi}}_{0}(\mathit{x}),{\mathit{\varphi}}_{0}(\mathit{x},0.5))=\rho ({\mathit{\varphi}}_{1}(\mathit{x}),{\mathit{\varphi}}_{1}(\mathit{x},0.5))\end{array}$$

(4)

where
${\mathit{\eta}}_{i}={\mathit{\varphi}}_{i}\circ (Id+{\mathit{v}}_{i})={\mathit{\varphi}}_{j}^{-1}\circ (Id+{\mathit{v}}_{i})$ from
${\mathit{\varphi}}_{i}={\mathit{\varphi}}_{j}^{-1}$ (i.e., imposition of the geodesic length constraint), *I _{i}*

$$\mathit{u}=\frac{-S{\nabla}_{{\mathit{\varphi}}_{i}}S}{{\alpha}_{U}^{2}/{\alpha}_{S}^{2}+{\Vert {\nabla}_{{\mathit{\varphi}}_{i}}S\Vert}_{2}^{2}}$$

(5)

where _{ϕi}*S* is the gradient of *S* with respect to *ϕ** _{i}* and we use

$$\begin{array}{l}{\mathit{\eta}}_{i}^{k}={\left({\mathit{\varphi}}_{j}^{k}\right)}^{-1}+K\mathit{u}\circ {\left({\mathit{\varphi}}_{j}^{k}\right)}^{-1}\\ ={\left({\mathit{\varphi}}_{j}^{k}\right)}^{-1}+\left({G}_{{\sigma}_{U}}\ast \mathit{u}\right){\left({\mathit{\varphi}}_{j}^{k}\right)}^{-1}\end{array}$$

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

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

$$\begin{array}{l}{E}_{2}({\mathit{\varphi}}_{i})=\frac{1}{2}\left\{({\Vert L({\mathit{\varphi}}_{1}^{-1}-{\mathit{\eta}}_{0})\Vert}_{2}^{2}-{\Vert L({\mathit{\varphi}}_{0}^{-1}-{\mathit{\eta}}_{1})\Vert}_{2}^{2})+\frac{{\alpha}_{P}^{2}}{{\alpha}_{U}^{2}}({\Vert \nabla {\mathit{\varphi}}_{1}^{-1}\Vert}_{2}^{2}+{\Vert \nabla {\mathit{\varphi}}_{0}^{-1}\Vert}_{2}^{2})\right\}\\ \text{subject}\phantom{\rule{0.16667em}{0ex}}\text{to}\phantom{\rule{0.16667em}{0ex}}{\mathit{\varphi}}_{i}\circ {\mathit{\varphi}}_{i}^{-1}=Id.\end{array}$$

(7)

where
$\rho {\left({\mathit{\varphi}}_{j}^{-1},{\mathit{\eta}}_{i}\right)}^{2}=\underset{{\mathit{v}}_{i}}{inf}\{{\Vert L{\mathit{v}}_{i}\Vert}_{2}^{2}\}=\underset{{\mathit{v}}_{i}}{inf}\{{\Vert L({\mathit{\varphi}}_{j}^{-1}-{\mathit{\eta}}_{i})\Vert}_{2}^{2}\}$ from Δ*t* =0.5 and
${\mathit{\eta}}_{i}={\mathit{\varphi}}_{j}^{-1}\circ (Id+{\mathit{v}}_{i})$. By omitting the invertibility constraint, (7) can be optimized using Tikhonov regularization^{22} which leads to a partial differential equation:

$$({L}^{\u2020}L)({\mathit{\varphi}}_{j}^{-1}-{\mathit{\eta}}_{i})+\frac{{\alpha}_{P}^{2}}{{\alpha}_{U}^{2}}{\nabla}^{2}{\mathit{\varphi}}_{j}^{-1}=0$$

(8)

By letting *a* =0 as in previous work^{12,17,23} we have *L* = *Id* and the solution of (8) can be approximated by
${\mathit{\varphi}}_{j}^{-1}\approx {G}_{{\sigma}_{D}}\ast {\mathit{\eta}}_{i}$ (i.e., a solution to an isotropic heat equation) where *G*_{σD} is a Gaussin kernel with width
${\sigma}_{D}=\sqrt{2}{\alpha}_{P}/{\alpha}_{U}$. The invertibility constraint is sequentially imposed by
$\text{minimizing}{\Vert {\mathit{\varphi}}_{i}\circ {\mathit{\varphi}}_{i}^{-1}-Id\Vert}_{2}^{2}$ 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

A MIND descriptor^{13} is constructed using a non-local mean operator^{24,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 * _{S}*.

A MIND descriptor for a voxel ** x** in an image

$${m}_{I,i}=cexp\phantom{\rule{0.16667em}{0ex}}\left(-\frac{d(I,\mathit{x},{\mathit{r}}_{i})}{V(I,\mathit{x})}\right)$$

(9)

where *c* is a normalization factor such that
$\underset{i\in [1,\mid {\mathcal{N}}_{S}\mid ]}{max}\{{m}_{I,i}(\mathit{x})\}=1$ (i.e., for intensity and contrast invariance) and *d*(*I*, ** x**,

$$d(I,\mathit{x},{\mathit{r}}_{i})=\sum _{\mathit{z}\in {\mathcal{N}}_{p}}{G}_{{\sigma}_{p}}(\mathit{z}){(I(\mathit{x}+\mathit{z})-I({\mathit{r}}_{i}+\mathit{z}))}^{2}$$

(10)

where * _{p}* denotes a neighborhood configuration of a patch (e.g., a cube),

Registration is generally ill-posed^{27} 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 *L*^{1} norm tends to be relatively insensitive to noise and outliers^{28} and could yield more reliable estimation (i.e., less confounded by noise) than quadratic norms (e.g., *L*^{2} and square norm); however, numerical optimization involving the *L*^{1} criterion is difficult due to its singularity. To exploit the robustness of *L*^{1} with the twice differentiability of *L*^{2}, we use a Huber distance^{14} as the similarity metric between a MIND descriptor of *I*_{0} *ϕ*_{0}, *m*_{I0}*ϕ*_{0}(** x**), and that of

$$S\left({\mathit{m}}_{{I}_{0}\circ {\mathit{\varphi}}_{0}},{\mathit{m}}_{{I}_{1}\circ {\mathit{\varphi}}_{1}},\mathit{x}\right)=\sum _{i=1}^{\mid {\mathcal{N}}_{S}\mid}\{\begin{array}{ll}{\phi}_{i}^{2}(\mathit{x})/2\epsilon ,\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}0\le \phantom{\rule{0.16667em}{0ex}}\mid {\phi}_{i}(\mathit{x})\mid \phantom{\rule{0.16667em}{0ex}}\le \epsilon \hfill \\ \mid {\phi}_{i}(\mathit{x})\mid -(\epsilon /2),\hfill & \text{otherwise}\hfill \end{array}$$

(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 *L*^{2} norm

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 *I*_{0.5} since the method uses the Lagrangian push-forward of the *I*_{0.5} domain to the domain of *I*_{0} and *I*_{1} 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}

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

Geometric accuracy of registration was measured in terms of TRE – the distance between corresponding anatomical points in *I*_{1} and *I*_{0} after registration:

$$\text{TRE}({\mathit{x}}_{0},{\mathit{x}}_{1};\mathit{\psi})={\Vert {\mathit{x}}_{0}-\mathit{\psi}({\mathit{x}}_{1})\Vert}_{2}$$

(12)

where *x** _{i}* denotes an unambiguous anatomical point (“target”) in

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

$$\text{SSIM}(\mathit{y},\mathit{z})=\frac{\left(2{\mu}_{\mathit{y}}{\mu}_{\mathit{z}}+{C}_{1}\right)\left(2{\sigma}_{\mathit{yz}}+{C}_{2}\right)}{\left({\mu}_{\mathit{y}}^{2}+{\mu}_{\mathit{z}}^{2}+{C}_{1}\right)\left({\sigma}_{\mathit{y}}^{2}+{\sigma}_{\mathit{z}}^{2}+{C}_{2}\right)}$$

(13)

where ** y** denotes a voxel in

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

$$\mathcal{J}(\mathit{\psi},{\mathit{\psi}}^{-1};\mathit{y},\mathit{z})=\frac{1}{2}\{{\Vert (\mathit{\psi}\circ {\mathit{\psi}}^{-1})(\mathit{y})-\mathit{y}\Vert}_{2}+{\Vert ({\mathit{\psi}}^{-1}\circ \mathit{\psi})(\mathit{z})-\mathit{z}\Vert}_{2}\}$$

(14)

where ** y** is a point in

The preservation of topology and lack of folding/tearing^{31} in a given deformation can be characterized in terms of its Jacobian determinant, (** x**) = det(

$$\mathcal{D}(\mathit{\psi})=\underset{\mathit{x}\in {\mathrm{\Omega}}_{V}}{min}\{\mathcal{J}(\mathit{x})\}$$

(15)

where the size of the sliding window Ω* _{V}* is set to 11

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 Demons^{34,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.

Figures 4(a–b) depict a simulated T_{2}-weighted MR moving image (*I*_{0}) and a CT fixed image (*I*_{1}) (1500×1500 pixels, 0.25 mm^{2}). Each 2D simulated image contains 5 simulated (rigid) vertebrae within (deformable) soft-tissue surroundings. The scoliotic curvature in *I*_{0} (Fig. 4(a)) was generated from a radial basis interpolation of a rigid motion^{36} associated with each vertebrae in *I*_{1}. The MR pixel value was computed using T_{2} 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.

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.

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 T_{2}-weighted MR and CT moving images (*I*_{0}), followed by T_{2}-weighted MR and intraoperative CT fixed images (*I*_{1}) 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 T_{E} = 125 ms. An average of two acquisitions was used to reduce noise in the images. The MR *I*_{0} was reconstructed at 0.9×0.9×0.9 mm^{3} with a size of 192×384×128 voxels, and the MR *I*_{1} was reconstructed at 0.5×0.5×0.9 mm^{3} 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 mm^{3} 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 *I*_{0} and CT *I*_{1} 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 (*I*_{0} and *I*_{1}) 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 *I*_{0} were translated to those defined in the MR *I*_{0} using NMI rigid registration. The MR *I*_{1} was visually compared to the MR *I*_{0} after nonlinear transformation.

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.

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 T_{2}-weighted MR and CT acquired for three patients undergoing intervention of cervical, thoracic, and lumbar disorders at our institution.

The T_{2}-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 T_{E} 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 mm^{3}. The MR images of the cervical, thoracic, and lumbar spines were used as moving images *I*_{0} (Figs. 6(a, c, e)) and their corresponding CT images were used as fixed images *I*_{1} (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.

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.

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 > 0. MIND Demons and NMI Demons yielded invertible deformations with sub-voxel ; however, NMI Demons achieved smaller than MIND Demons due to stronger smoothing applied to the update and displacement fields. The invertibility was not measured for the FFD-based methods since the methods only provided forward deformations (i.e., a map from *I*_{0} to *I*_{1}). 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 *I*_{0} after registration. The cyan spheres represent the target points in *I*_{1} 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.

Simulation study results. Transformed *I*_{0} 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.

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 (mean = 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 > 0; however, MI FFD yielded min close to 0 (Table 3 and Fig. 8(c)). The large ranges in 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 *I*_{0} and the cyan CT *I*_{1} after registration. The cyan spheres represent the target points in *I*_{1} 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 *I*_{0} after registration superimposed on the gray MR *I*_{1}. MIND Demons achieved statistically significant improvement in registration accuracy with mean TRE of 1.5 mm.

MR-to-CT registration of the ovine spine phantom. (a–c) TRE, , and as a function of registration methods.

MR-to-CT registration of the ovine spine phantom. (Top) Semi-opaque surface rendering of the pink MR moving image *I*_{0} and the cyan fixed CT image *I*_{1} after registration. Cyan spheres represent the target points in *I*_{1} and yellow lines mark distances between **...**

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 and >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 *I*_{0} and the cyan CT *I*_{1} 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 *I*_{1} 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 *I*_{1} superimposed on the gray MR *I*_{0} after registration, and the bottom row shows checkerboard images between the CT *I*_{1} and the MR *I*_{0} 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.

MR-to-CT registration in clinical studies. (a–c) TRE, , and as a function of the spinal sections: cervical, thoracic, and lumbar spines.

A deformable registration method merging the Demons^{11} and SyN^{10} approaches for symmetric time-dependent diffeomorphisms has been developed. The algorithm incorporates MIND descriptors^{13} and the Huber metric^{14} 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

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 information^{40,41} and/or image features^{42,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 term^{45} or extra-dimension^{46} 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.

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.

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |