Keywords

1 Introduction

During neurosurgery, Image-Guided Neurosurgical Systems (IGNSs) provide a patient-to-image mapping that relates the preoperative image data to an intraoperative patient coordinate system, allowing surgeons to infer the locations of their surgical instruments relative to preoperative image data and helping them to achieve a radical tumor resection while avoiding damage to surrounding functioning brain tissue.

Commercial IGNSs assume a rigid registration between preoperative imaging and patient coordinates. However, intraoperative deformation of the brain, also known as brain shift, invalidates this assumption. Since brain shift progresses during surgery, the rigid patient-to-image mapping of IGNS becomes less and less accurate. Consequently, most surgeons only use IGNS to make a surgical plan but justifiably do not trust it throughout the entire operation [1, 2].

Related Work. As one of the most important error sources in IGNS, intraoperative brain shift must be compensated in order to increase the accuracy of neurosurgeries. Registration between the intraoperative MRI (iMRI) image and preoperative MRI (preMRI) image (preop-to-intraop registration) has been a successful strategy for brain shift compensation [3,4,5,6]. However, iMRI acquisition is disruptive, expensive and time consuming, making this technology unavailable for most clinical centers worldwide. More recently, 3D intraoperative Ultrasound (iUS) appears to be a promising replacement for iMRI. Although some progress has been made by previous work on preMRI-to-iUS registration [7,8,9,10,11,12,13], yet there are still no clinically accepted solutions and no commercial neuro-navigation systems that provide brain shift compensation. This is due to three reasons: (1) Most non-rigid registration methods can not handle artifacts and missing structures in iUS; (2) The multi-modality of preMRI-to-iUS registration makes the already difficult problem even more challenging; (3) A few methods [14] can achieve a reasonable alignment, yet they take around 50 min for an US pair and are too slow to be clinically applicable. Another shortcoming of existing brain shift compensation approaches is the lack of an uncertainty measure. Brain shift is a complex spatio-temporal phenomenon and, given the state of registration technology and the importance of the result, it seems reasonable to expect an indication (e.g. error bars) of the confidence level in the estimated deformation.

In this paper, we propose a novel feature-driven active framework for brain shift compensation. Here, landmarks and their displacement are first estimated from a pair of US images using corresponding local image features. Subsequently, a Gaussian Process (GP) model [15] is used to interpolate a dense deformation field from the sparse landmarks. Kernels of the GP are estimated by using variograms and a discrete grid search method. If necessary, for areas that are difficult to align, the user can actively add new landmarks based on the image context and visualization of the uncertainty measure provided by the GP to further improve the registration accuracy. We retrospectively demonstrate the efficacy of our method on clinical data.

Contributions and novelties of our work can be summarized as follows:

  1. 1.

    The proposed feature-based registration is robust for aligning iUS image pairs with missing correspondence and is fast.

  2. 2.

    We explore applying the GP model and variograms for image registration.

  3. 3.

    Registration uncertainty in transformation parameters can be naturally obtained from the GP model.

  4. 4.

    To the best of our knowledge, the proposed active registration strategy is the first method to actively combine user expertise in brain shift compensation.

2 Method

2.1 The Role of US-to-US Registration

In order to alleviate the difficulty of preop-to-intraop registration, instead of directly aligning iMRI and iUS images, we choose an iterative compensation approach which is similar to the work in [16].

As shown in Fig. 1 the acquisition processes for pre-duraUS (preUS) and post-resectionUS (postUS) take place before opening the dura and after (partial) tumor resection, respectively. Since most brain-shift occurs after taking the preUS, a standard multi-modal registration may be suffice to achieve a good alignment \(T_{\mathrm {multi}}\) between preMRI and preUS [12]. Next, we register the preUS to postUS using the proposed feature-driven active framework to acquire a deformable mapping \(T_{\mathrm {mono}}\). After propagating \(T_{\mathrm {multi}}\) and \(T_{\mathrm {mono}}\) to the preMRI, surgeons may use it as an updated view of anatomy to compensate for brain shift during the surgery.

Fig. 1.
figure 1

Pipeline of the US-based brain shift compensation.

Full size image

2.2 Feature-Based Registration Strategy

Because of tumor resection, compensating for brain shift requires non-rigid registration algorithms capable of aligning structures in one image that have no correspondences in the other image. In this situation, many image registration methods that take into account the intensity pattern of the entire image will become trapped in incorrect local minima.

We therefore pursue a Feature-Based Registration (FBR) strategy due to its robustness in registering images with missing correspondence [17]. FBR mainly consists of 3 steps: feature-extraction, feature-matching and dense deformation field estimation. An optional “active registration” step can be added depending on the quality of FBR.

Fig. 2.
figure 2

Pipeline of the feature-based active preUS-to-postUS registration.

Full size image

Feature Extraction and Matching. As illustrated in Fig. 2(a) and (b), distinctive local image features are automatically extracted and identified as key-points on preUS and postUS images. An automatic matching algorithm searches for a corresponding postUS key-point for each key-point on the preUS image [17].

For a matched key-point pair, let \(\mathbf {x}_i\) be the coordinates of the preUS key-point and \(\mathbf {x}^{\mathrm {post}}_i\) be the coordinate of its postUS counterpart. We first use all matched PreUS key-points as landmarks, and perform a land-mark based preUS-to-postUS affine registration to obtain a rough alignment. \(\mathbf {x}^{\mathrm {post}}_i\) becomes \(\mathbf {x}^{\mathrm {affine}}_i\) after the affine registration. The displacement vector, which indicates the movement of landmark \(\mathbf {x}_i\) due to the brain shift process, can be calculated as \(\mathbf {d}(\mathbf {x}_i)=\mathbf {x}^{\mathrm {affine}}_i-\mathbf {x}_i\). where \(\mathbf {d}=[d_x,d_y,d_z]\).

Dense Deformation Field. The goal of this step is to obtain a dense deformation field from a set of N sparse landmark and their displacements \(\mathcal {D}=\{ (\mathbf {x}_i,\mathbf {d}_i),i=1:N \}\), where \(\mathbf {d}_i=\mathbf {d}(\mathbf {x}_i)\) is modeled as a observation of displacements.

In the GP model, let \(\mathbf {d}(\mathbf {x})\) be the displacement vector for the voxel at location \(\mathbf {x}\) and define a prior distribution as \(d(\mathbf {x})\sim \mathrm {GP}(\mathrm {m}(\mathbf {x}),\mathrm {k}(\mathbf {x},\mathbf {x}'))\), where \(\mathrm {m}(\mathbf {x})\) is the mean function, which usually is set to 0, and the GP kernel \(\mathbf {k}(\mathbf {x},\mathbf {x}')\) represents the spatial correlation of displacement vectors.

By the modeling assumption, all displacement vectors follow a joint Gaussian distribution \(p(\mathbf {d}\mid \mathbf {X})=\mathcal {N} (\mathbf {d}\mid \mathbf {\mu },\mathbf {K}) \), where \(K_{ij}=\mathbf {k}(\mathbf {x},\mathbf {x}')\) and \(\mathbf {\mu } = (\mathrm {m}(\mathbf {x}_1) ,...,\mathrm {m}(\mathbf {x}_N)) \). As a result, the displacement vectors \(\mathbf {d}\) for known landmarks and \(N_*\) unknown displacement vectors \(\mathbf {d}_*\) at location \(\mathbf {X}_*\), which we want to predict, have the following relationship:

$$\begin{aligned} \begin{pmatrix} \mathbf {d}\\ \mathbf {d}_* \\ \end{pmatrix}\sim \mathrm {GP}\left( \begin{pmatrix} \mathbf {\mu } \\ \mathbf {\mu }_*\\ \end{pmatrix}, \begin{pmatrix} \mathbf {K} &{} \mathbf {K}_*\\ \mathbf {K}^T_* &{} \mathbf {K}_{**}\\ \end{pmatrix}\right) _. \end{aligned}$$
(1)

In Eq. 1, \(\mathbf {K}=\mathrm {k}(\mathbf {X},\mathbf {X})\) is the \(N\times N\) matrix, \(\mathbf {K}_*=\mathrm {k}(\mathbf {X},\mathbf {X_*})\) is a similar \(N \times N_*\) matrix, and \(\mathbf {K_{**}}=\mathrm {k}(\mathbf {X_*},\mathbf {X_*})\) is a \(N_* \times N_*\) matrix. The mean \(\mu _*=[\mu _{*x},\mu _{*y},\mu _{*z}]\) represents values of voxel-wise displacement vectors and can be estimated from the posterior Gaussian distribution \(p(\mathbf {d}_*\mid \mathbf {X_*},\mathbf {X},\mathbf {d})= \mathcal {N}(\mathbf {d}_*\mid \mu _*,\Sigma _*)\) as

$$\begin{aligned} \mu _*= \mu (\mathbf {X_*})+\mathbf {K}^T_*\mathbf {K}^{-1}(\mathbf {d}-\mu (\mathbf {X})). \end{aligned}$$
(2)

Given \(\mu (\mathbf {X})= \mu (\mathbf {X_*})=0\), we can obtain the dense deformation field for the preUS image by assigning \(\mu _{*x}\),\(\mu _{*y}\),\(\mu _{*z}\) to \(d_x\), \(d_y\) and \(d_z\), respectively.

Active Registration. Automatic approaches may have difficulty in the preop-to-intraop image registration, especially for areas near the tumor resection site. Another advantage of the GP framework is the possibility of incorporating user expertise to further improve the registration result.

From Eq. 1, we can also compute the covariance matrix of the posterior Gaussian \(p(\mathbf {d}_*\mid \mathbf {X_*},\mathbf {X},\mathbf {d})\) as

$$\begin{aligned} \Sigma _*= \mathbf {K}_{**}-\mathbf {K}^T_*\mathbf {K}^{-1}\mathbf {K}_*. \end{aligned}$$
(3)

Entries on the diagonal of \(\Sigma _*\) are the marginal variances of predicted values. They can be used as an uncertainty measure to indicates the confidence in the estimated transformation parameters.

Fig. 3.
figure 3

(a) \(\Vert {d_x(\mathbf {x}_2)-d_x(\mathbf {x}_1)}\Vert \) and h; (b) Empirical variogram cloud; (c) Variogram cloud divided into bins with their means marked as blue.

Full size image

If users are not satisfied by the FBR alignment result, they could manually, guided by the image context and visualization of registration uncertainty, add new corresponding pairs of key-points to drive the GP towards better results.

2.3 GP Kernel Estimation

The performance of GP registration depends exclusively on the suitability of the chosen kernels and its parameters. In this study, we explore two schemes for the kernel estimation: Variograms and discrete grid search.

Variograms. The variogram is a powerful geostatistical tool for characterizing the spatial dependence of a stochastic process [18]. While being briefly mentioned in [19], it has not yet received much attention in the medical imaging field.

In the GP registration context, where \(\mathbf {d}(\mathbf {x})\) is modelled as a random quantity, variograms can measure the extent of pairwise spatial correlation between displacement vectors with respect to their distance, and give insight into choosing a suitable GP kernel.

In practice, we estimate the empirical variogram of landmarks’ displacement vector field using

$$\begin{aligned} \begin{aligned} \hat{\gamma }(h\pm \delta )&:=\frac{1}{2|N(h\pm \delta )|}\sum _{(i,j)\in {N(h\pm \delta )}}^{}\Vert {\mathbf {d}(\mathbf {x}_i)-\mathbf {d}(\mathbf {x}_j)}\Vert ^2.&\end{aligned} \end{aligned}$$
(4)

For the norm term \(\Vert {\mathbf {d}(\mathbf {x}_i)-\mathbf {d}(\mathbf {x}_j)}\Vert \), we separate its 3 components \(d_x\) \(d_y\) \(d_z\) and construct 3 variograms respectively. As shown in Fig. 3(a), for displacement vectors \(\mathbf {d}(x_1)\) and \(\mathbf {d}(x_2)\), \(\Vert {d_x(\mathbf {x}_2)-d_x(\mathbf {x}_1)}\Vert \) is the vector difference with respect to the \(\mathbf {x}\) axis, etc. h represents the distance between two key-points.

To construct an empirical variogram, the first step is to make a variogram cloud by plotting \(\Vert {d(\mathbf {x}_2)-d(\mathbf {x}_1)}\Vert ^2\) and \(h_{ij}\) for all displacement pairs. Next, we divide the variogram cloud into bins with a bin width setting to 2\(\delta \). Lastly, the mean of each bin is calculated and further plotted with the mean distance of that bin to form an empirical variogram. Figure 4(a) shows an empirical variogram of a real US image pair that has 71 corresponding landmarks.

In order to obtain the data-driven GP kernel function, we further fit a smooth curve, generated by pre-defined kernel functions, to the empirical variogram. As shown in Fig. 4(b), a fitted curve is commonly described by the following characteristics:

Nugget:

The non-zero value at \(h=0\).

Sill:

The value at which the curve reaches its maximum.

Range:

The value of distance h where the sill is reached.

Fig. 4.
figure 4

(a) X-axis empirical variogram of a US images pair; (b) Sill, range and nugget; (c) Fitting a continuous model to an empirical variogram.

Full size image

Fitting a curve to an empirical variogram is implemented in most geostatistics software. A popular choice is choosing several models that appear to have the right shape and use the one with smallest weighted squared error [18]. In this study, we only test Gaussian curves

$$\begin{aligned} \gamma (h)=c_0+c\{1-exp(- \frac{h^2}{a} )\}. \end{aligned}$$
(5)

Here, \(c_0\) is the nugget, \(c=\mathrm {Sill}-c_0\) and a is the model parameter. Once the fitted curve is found, we can obtain a from the equation (5) and use it as the Gaussian kernel scale in the GP interpolation.

Discrete Grid Search. The variogram scheme often requires many landmarks to work well [18]. For US pairs that have fewer landmarks, we choose predefined Gaussian kernels, and use cross validation to determine the scale parameter in a discrete grid search fashion [15].

3 Experiments

The experimental dataset consists of 6 sets 3D preUS and postUS image pairs. The US signals were acquired on a BK Ultrasound 3000 system that is directly connected to the Brainlab VectorVision Sky neuronavigaton system during surgery. Signals were further reconstructed as 3D volume using the PLUS [20] library in 3D Slicer [21] (Table 1).

Table 1. Registration evaluation results (in mm)
Full size table

We used the mean euclidean distance between the predicted and ground truth of key-points’ coordinates, measured in mm, for the registration evaluation. During the evaluation, we compared: affine, thin-plate kernel FBR, variograms FBR and gaussian kernel FBR. For US pairs with fewer than 50 landmarks, we used leave-one-out cross validation, otherwise we used 5-fold cross validation. All of the compared methods were computed in less than 10 min.

The pre-defined Gaussian kernel with discrete grid search generally yield better result than the variogram scheme. This is reasonable as the machine learning approach stresses the prediction performance, while the geostatistical variogram favours the interpretability of the model. Notice that the cross validation strategy is not an ideal evaluation, this could be improved by using manual landmarks in public datasets, such as RESECT [22] and BITE [23].

In addition, we have performed preliminary tests on active registration as shown in Fig. 5, which illustrate the use of a colour map of registration uncertainty to guide the manual placement of 3 additional landmarks to improve the registration. By visual inspection, we can see the alignment of tumor boundary substantially improved.

Fig. 5.
figure 5

(a) FBR result of the preUS with a tumor boundary outlined in green; (b) Overlaying the visualization of uncertainty on the preUS image. A characteristic of GP is that voxels near landmarks tend to have smaller uncertainty. In this example, all landmarks happen to be located near the large sulcus, hence the incertitude looks high everywhere else except around the sulcus. (c) Active registration result of the preUS with a tumor boundary outlined in blue; (d) Overlaying the green and blue tumor boundary on the target image.

Full size image

4 Discussion

One key point of our framework is the “active registration” idea that aims to overcome the limitation of automatic image registration. Human and machines have complementary abilities; we believe that the element of simple user interaction should be added to the pipeline for some challenging medical imaging applications. Although the proposed method is designed for brain shift compensation, it is also applicable to other navigation systems that require tracking of tissue deformation. The performance of FBR is highly correlated with the quality of feature matching. In future works, we plan to test different matching algorithms [24], and also perform more validation with public datasets.