Keywords

1 Introduction

Brain shift combined with registration and tracking errors reduces the accuracy of image-guided neurosurgery based on neuronavigation systems [1,2,3]. Intraoperative ultrasound, being a real-time imaging modality, has the potential to enable the surgeon to accurately localize the instrument trajectories in the operative field and thus facilitate accurate resection to promote better surgical outcomes. However, registration of intraoperative US with preoperative MR images is a challenging problem due to the different information captured by each image modality. We present a deformable MR-US registration algorithm that uses attribute matching and mutual-saliency weighting and apply it to image-guided neurosurgery.

2 Methods

A registration framework usually consists of three parts: (1) the similarity measure, which defines the criterion to align the two images; (2) the deformation model, which defines the mechanism to transform one image to the other; and (3) the optimization strategy, which is used to determine the best parameters of the deformation model. An open question is how to define the similarity measure for MR-US registration. A deformable registration algorithm known as DRAMMS [4] has shown promise in defining similarity based on optimal Gabor attributes modulated by quantified matching reliabilities. It shows potential in handling large deformations and missing correspondences [5]. However, the original DRAMMS defines similarity by the Sum of Squared Differences (SSD) of attributes which limits its application in MR-US multi-modal registration. We propose the Correlation Coefficient (CC) of attributes to better adapt to MR-US multi-modal registration.

2.1 Problem Formulation

In the original DRAMMS formulation, given two images \( I_{1} :\Omega _{1} \mapsto {\mathbb{R}} \) and \( I_{2} :\Omega _{2} \mapsto {\mathbb{R}} \) in 3D image domains \( \Omega _{i} \) \( (i = 1,2) \)\( {\mathbb{R}}^{3} \), we seek a transformation \( T \) that maps every voxel \( \varvec{u} \) \( \in \) \( \Omega _{1} \) to its correspondence T \( (\varvec{u}) \) \( \in \) \( \Omega _{2} \), by minimizing an overall cost function \( E\left( T \right) \),

$$ \mathop {\hbox{min} }\nolimits_{T} E\left( T \right) = \int_{{\varvec{u} \in\Omega _{1} }} {\underbrace {{ms\left( {\varvec{u},T\left( \varvec{u} \right)} \right)}}_{{{\text{Mutual}} - {\text{Saliency}}}}. \underbrace {{sim \left( {A_{1}^{{ \star }} \left( \varvec{u} \right), A_{2}^{{ \star }} \left( {T\left( \varvec{u} \right)} \right)} \right)d\varvec{u}}}_{{{\text{Attribute}}\,{\text{Matching}}}} + \lambda R\left( T \right)} $$
(1)

where \( A_{1}^{{ \star }} \left( \varvec{u} \right) \left( {i = 1,2} \right) \) is the optimal attribute vector that reflects the geometric and anatomical contexts around voxel \( \varvec{u} \), and \( d \) is its dimension. The term \( ms\left( {\varvec{u}, T\left( \varvec{u} \right)} \right) \) is a continuously-valued mutual-saliency weight between two voxels \( \varvec{u} \) \( \in \) \( \Omega _{1} \) and T\( (\varvec{u}) \) \( \in \) \( \Omega _{2} \). This way, registration is driving by reliably matched voxel pairs which are not necessarily less deformed voxels. The term \( R\left( T \right) \) is a smoothness/regularization term usually corresponding to the Laplacian operator, or the bending energy, of the deformation field \( T \), whereas \( \lambda \) is a parameter that controls the extent of smoothness. The proposed framework is sketched in Fig. 1.

Fig. 1.
figure 1

Non-rigid deformation framework.

2.2 Attribute Matching

The aim of attribute matching is to extract and select optimal attributes that reflect the geometric and anatomic contexts of each voxel. It consists of two parts: attribute extraction and attribute selection.

Attribute Extraction.

A set of multi-scale and multi-orientation Gabor attributes is extracted at each voxel by convolving the images with a set of Gabor filter banks. The parameter settings developed by [6] were adopted: the number of scales, M, is set to 4 and the number of orientations, N, is set to 6, the highest frequency is set at 0.4 Hz and the lowest frequency at 0.05 Hz. Figure 2 shows an example of multi-scale and multi-orientation Gabor attributes extracted from (a) intraoperative US and (b) preoperative MR images. After attribute extraction, each voxel \( {\mathbf{u}} = \left( {x,y,z} \right) \) is characterized by a Gabor attribute vector \( \tilde{A}_{ 1} \left( u \right) \) with dimension \( D = M \times N \times 4 \).

Fig. 2.
figure 2

Multi-scale and multi-orientation Gabor attributes extracted from (a) intraoperative US image and (b) preoperative MR images.

Attribute Selection.

The aim is to select components of attributes to increase the reliability of matching between two images. An expectation-maximization (EM) framework is used. Given the full-length attributes, the E-step finds spatially-scattered (and hence spatially representative) voxel pairs with not just high similarity but more importantly, reliably high similarity. Then on the selected voxel pairs, M-step uses the iterative forward inclusion and backward elimination (iFIBE) feature selection algorithm to find a subset of attribute components that maximize the similarity and matching reliability (defined as mutual-saliency). A major difference from [4] is that the similarity \( sim\left( {\varvec{p},\varvec{q}} \right) \) between a pair of voxels \( \left( {\varvec{p} \in\Omega _{1} ,\varvec{q} \in\Omega _{2} } \right) \) is defined based on the correlation coefficient of their attribute vectors:

$$ sim\left( {\tilde{A}_{1} \left( \varvec{p} \right),\tilde{A}_{2} \left( \varvec{q} \right)} \right) = CC\left( {\tilde{A}_{1} \left( \varvec{p} \right), \tilde{A}_{2} \left( \varvec{q} \right)} \right) \in \left[ {0,1} \right] $$
(2)

where higher correlation coefficient between to attribute vectors indicates higher similarity.

2.3 Mutual-Saliency Map to Modulate Registration

We assign a continuously-valued weight to each voxel, based on the capability of each voxel to establish reliable correspondences across images. This idea is formulated in Eq. (3) and in the associated Fig. 3. Mutual-saliency value, \( ms\left( {\varvec{u}, T\left( \varvec{u} \right)} \right), \) is calculated by dividing the mean similarity between \( \varvec{u} \) and all voxels in the core neighborhood (CN) of \( T\left( u \right) \), with the mean similarity between \( \varvec{u} \) and all voxels in the peripheral neighborhood (PN) of \( T\left( u \right) \).

Fig. 3.
figure 3

The idea of mutual-saliency measure: (a) The matching between a pair of voxels u and T(u) is unique if they are similar to each other and not to anything else in the neighborhood. A delta function in the similarity map indicates a unique matching, and hence high mutual-saliency value. (b) Mutual-saliency function measures the uniqueness of the matching between a voxel pair in a neighborhood. Different colors represent different layers of neighborhoods.

$$ ms\left( {\varvec{u},T\left( \varvec{u} \right)} \right) = \frac{{MEAN_{{\varvec{v} \in CN\left( {T\left( \varvec{u} \right)} \right)}} \left[ {sim\left( {A_{1} \left( \varvec{u} \right), A_{2} \left( \varvec{v} \right)} \right)} \right]}}{{MEAN_{{\varvec{v} \in PN\left( {T\left( \varvec{u} \right)} \right)}} \left[ {sim\left( {A_{1} \left( \varvec{u} \right), A_{2} \left( \varvec{v} \right)} \right)} \right]}} $$
(3)

where sim (∙, ∙) is the attribute-based similarity between two voxels (Eq. 2). The radius of each neighborhood is adaptive to the scale in which Gabor attributes are extracted. For a typical isotropic 3D brain image, the radius of core, transitional and peripheral neighborhood are 2, 5, 8 voxels, respectively [4].

2.4 Deformation Model and Optimization Strategy

The diffeomorphic FFD [7, 8] is chosen because of its flexibility to handle a smooth and diffeomorphic deformation field. We have chosen discrete optimization, a state-of-the-art optimization strategy known for computational efficiency and robustness regarding local optima [9, 10].

3 Results

Preoperative T2-FLAIR MR and intraoperative predurotomy US images were acquired from 22 patients with low-grade gliomas. A set of 15 to 16 homologous landmarks were identified across images (pre-operative MRI vs. US before resection) and used to validate the deformable registration algorithm [11]. The Transforms Module (General Registration Brain) in 3D Slicer was used to determine an initial rigid transformation between each pair of images before applying our method. Figure 4 shows an example of one pair of preoperative T2-FLAIR and intraoperative US images and their alignment before and after deformable registration. Table 1 presents the mean target registration error (mTRE) in mm and for each clinical case, before and after deformable registration.

Fig. 4.
figure 4

Preoperative MR to intraoperative US registration. (a) preoperative T2-FLAIR MR and (b) intraoperative US images. Superimposed preoperative MR to intraoperative US (c) before registration and (d) after non-linear correction. Arrows indicate (c) the misalignment between the tumor boundaries in the different modalities and (d) the alignment after registration.

Table 1. Details of inter-modality landmarks for each clinical case. The number of landmarks and the mean initial Euclidean distances between landmark pairs are shown, and the range (min–max) of the distances is shown in parenthesis after the mean value. The last column shows the results after registration.

4 Conclusions

The proposed registration algorithm reduces the mean target registration error from an initial value of 5.37 ± 4.27 mm to 3.35 ± 1.19 mm for 22 clinical cases. Our future work includes: (i) further comparison of different similarity measures in MR-US registration, (ii) exploring different linear registration methods to initialize deformable registration, (iii) using blurring and gradient information to reduce the negative influence of speckles in the ultrasound image and (iv) investigating the potential to combine landmark-based to voxel-wise registration.