Keywords

1 Introduction

Image-guided neurosurgery involves the display of preoperative anatomy and surgical plans in intraoperative reference frame to increase the accuracy of pathological tissue resection and to reduce damage to the surrounding structures. Preoperative MRI can reveal information such as nerve fiber tracts and brain activation areas. Interventional MRI (iMRI) can image intraoperative deformations due to cerebrospinal fluid (CSF) drainage, gravity and edema (collectively, brain shift) [9]. Non-rigid registration of preoperative MRI to intraoperative iMRI enables surgical guidance using propagated preoperative plans [2].

Correspondences missing due to resection present a challenge to the registration. Daga et al. [2] estimated brain shift intraoperatively by masking out voxels lying outside a brain mask. However, automated brain extraction such as using FSL-BET [11] can be inaccurate near the resection cavity due to fluid accumulation and surgical gauze in the cavity. Another challenge to registration arises from contrast changes due to CSF drainage, bleeding, edema, MRI bias field and low signal to noise ratio (SNR) of iMRI.

We consider registration of a T1-weighted (T1w) image pair. The local normalized correlation coefficient (LNCC, [1]) captures a local affine intensity relationship. LNCC involves smoothing based on convolution with a Gaussian kernel, which includes voxels located outside a tissue or resection boundary in the statistics of a local neighborhood centered inside the boundary. This potentially reduces the matching specificity near the resection margin.

The bilateral filter was introduced for edge-preserving image smoothing and weights the voxels in the local neighborhood based on their spatial distance and intensity difference from the central voxel [13]. Bilateral filtering was used for locally adaptive patch-based similarity cost evaluation in a stereo reconstruction problem [14]: as pixels on a surface tended to have similar colors, the estimated disparity map became more accurate. In a T1w image pair, the voxels from the same tissue tend to have similar intensities and we suggest that bilateral weighting can lead to more accurate brain shift estimation.

We propose to introduce adaptive bilateral weighting into LNCC calculation as illustrated in Fig. 1. We evaluated the modified measure in registration experiments on three datasets and found an improvement in registration accuracy.

2 Methods

2.1 Bilateral Adaptively Weighted LNCC Similarity

LNCC was first used in the context of image registration by [1]. Let R be the reference image and F the floating image in the same coordinate space, then LNCC for the local neighborhood of a point \({\varvec{v}}\) is defined as

$$\begin{aligned} {\mathrm {LNCC}_{{\varvec{v}}}(R,F)}^2 = \frac{ {\langle {R,F}\rangle _{{\varvec{v}}}}^2 }{ \langle {R,R}\rangle _{{\varvec{v}}} \cdot \langle {F,F}\rangle _{{\varvec{v}}}}, \end{aligned}$$
(1)

where the \(\langle {R,R}\rangle _{{\varvec{v}}}\) and \(\langle {F,F}\rangle _{{\varvec{v}}}\) are the local variances and \(\langle {R,F}\rangle _{{\varvec{v}}}\) is the local covariance. The latter is defined as \( \langle {R,F}\rangle _{{\varvec{v}}} = \overline{R \cdot F}_{{\varvec{v}}} - \overline{R}_{{\varvec{v}}} \cdot \overline{F}_{{\varvec{v}}}, \) where \(\overline{R}_{{\varvec{v}}}\) and \(\overline{F}_{{\varvec{v}}}\) are the respective local means. The local variances are defined analogously. The local mean for R is defined as \( \overline{R}_{{\varvec{v}}} = \frac{1}{N} \sum _{{\varvec{x}}} { R({\varvec{v}}-{\varvec{x}}) w_{{\varvec{v}}}({\varvec{x}}) }, \) where N is the number of voxels in the neighborhood of \({\varvec{v}}\), \({\varvec{x}}\) is the offset relative to \({\varvec{v}}\) and \(w_{{\varvec{v}}}({\varvec{x}})\) are the weights, here given by a generic term that depends on \({\varvec{v}}\). The local mean for F is defined analogously.

LNCC uses a Gaussian kernel for the local weights, \( w_{{\varvec{v}}}({\varvec{x}})= G_{\beta }({\varvec{x}})= \frac{1}{\sqrt{2\pi }\beta } \mathrm {exp} \left( {-\frac{{|{\varvec{x}}|}^2}{2\beta ^2}} \right) \!, \) where \(\beta \) controls the neighborhood’s size (negligible for \(|{\varvec{x}}| > 3\beta \)). Since \(G_{\beta }({\varvec{x}})\) does not depend on \({\varvec{v}}\), the local mean can be implemented using convolution \( \overline{I}_{{\varvec{v}}} = \left( G_{\beta } *I \right) ({\varvec{v}}) \).

Fig. 1.
figure 1

(a) T1-weighted reference image and (b) intensity-based weights for a point (blue cross). (c) T1-weighted floating image and (d) intensity-based weights for the point. (e) Final weights based on (b), (d) and distance from the point. (Color figure online)

Full size image

We introduce bilateral adaptive weighting and refer to the modified measure as LNCC-AW. A bilateral filtered smoothing of an arbitrary image I is

$$\begin{aligned} \overline{I}_{{\varvec{v}}}^{\mathrm {bilat.}} = \frac{1}{N} \sum \limits _{{\varvec{x}}} { I({\varvec{v}}-{\varvec{x}}) \cdot G_\beta ({\varvec{x}}) \cdot G_\alpha \left( I({\varvec{v}}-{\varvec{x}}) - I({\varvec{v}}) \right) }, \end{aligned}$$
(2)

where \(G_{\alpha }(d)\) is a range kernel i.e. a kernel for the intensity difference \(d=I({\varvec{v}}-{\varvec{x}})-I({\varvec{v}})\). The edge-preserving property arises as the voxels beyond an intensity rise/drop are excluded. Given images R and F to register, we guide the adaptive weighting by both the images as in [14] by using the composite term

$$\begin{aligned} w_{{\varvec{v}}}({\varvec{x}}) = G_\beta ({\varvec{x}}) \cdot G_\alpha \left( R({\varvec{v}}-{\varvec{x}})-R({\varvec{v}}) \right) \cdot G_\alpha \left( F({\varvec{v}}-{\varvec{x}})-F({\varvec{v}}) \right) \end{aligned}$$
(3)

as illustrated in Fig. 1. Since \(w_{{\varvec{v}}}({\varvec{x}})\) vary spatially, we can no longer implement the local mean using convolution nor take advantage of kernel separability.

The low SNR of iMRI can potentially cause the weights of adjacent neighborhoods in homogeneous areas to vary along with the varying intensity of the central voxels. In order to reduce spatial inconsistencies in similarity values, we replace the Gaussian range kernel, as used in [13, 14], with a kernel shaped as Student’s t-distribution, which down-weights rather than suppresses differing intensities:

$$\begin{aligned} G_\alpha (d)=\frac{ \mathrm {\Gamma }(\frac{\nu +1}{2})}{\sqrt{\nu \pi \alpha ^2} \mathrm {\Gamma }(\frac{\nu }{2})} {\left( 1+\frac{d^2}{\nu \alpha ^2} \right) }^{-\frac{\nu +1}{2} }. \end{aligned}$$
(4)

We selected \(\nu =2\) as it has a gradual drop-off and provides a trade-off between boundary-preservation and robustness to noise. For \(\alpha =\infty \), the weighting reduces to locally non-varying.

2.2 Registration Using a Discrete Optimization Framework

The derivation of analytical gradient of the similarity measure, for instance with respect to a voxel-based deformation field, for use in gradient-based non-rigid registration schemes becomes complicated when using adaptive weighting, because the gradient depends on the local weights which in turn depend on the deformation. However, [4] reformulated non-rigid registration as a discrete Markov Random Field (MRF) optimization problem, for which the similarity measure gradient is not needed. We employ the proposed measure in a related discrete optimization scheme of [6]. A grid \(\mathcal {P}\) of B-spline transformation control points \(p\in \mathcal {P}\) with positions \({\varvec{c}}_p\) is overlaid onto the reference image. The control point displacements in the floating image are \({\varvec{u}}_p=[u_p,v_p,w_p]\) with discrete valued components. For efficiency, a minimum spanning tree \(\mathcal {N}\) of the most relevant edges \((p,q)\in \mathcal {N}\) is optimized rather than a full MRF. Displacements are sought minimizing the energy

$$\begin{aligned} \sum \limits _{p \in \mathcal {P}}{\left( 1-\Vert \text {LNCC}_{{\varvec{c}}_{\varvec{p}}}(R(\varvec{\xi }),F(\varvec{\xi }+{\varvec{u}}_{\varvec{p}}))\Vert \right) } +\alpha \sum \limits _{(p,q) \in \mathcal {N}}{\frac{{\Vert {\varvec{u}}_{\varvec{p}}-{\varvec{u}}_{\varvec{q}}\Vert }^2}{\Vert {\varvec{x}}_{\varvec{p}}-{\varvec{x}}_{\varvec{q}}\Vert }.} \end{aligned}$$
(5)

3 Experiments

3.1 Patch Matching on 2D Synthetic Phantom

We compare matching accuracy for two 2D synthetic phantoms. We place a fixed patch representing a local neighborhood in the reference image and a moving patch in the floating image to plot the similarity profile of LNCC and LNCC-AW, respectively. We assess two phantoms. A contrast-enhanced lesion near a resection phantom is shown in Fig. 2(a–d). The similarity profile for LNCC has a mild maximum at the true zero displacement due to voxels included from the resected area. The similarity profile for LNCC-AW has a clear maximum due voxels down-weighted in the resected area. A phantom of the medial longitudinal fissure is shown in Fig. 2(e–h). The reference and floating image are the same axial slice from the BrainWeb database. The patch is centered next to the medial longitudinal fissure that contains dark voxels in the CSF and the falx cerebri. The similarity profile for LNCC has a band of false matches due to voxels included in the fissure. The similarity profile for LNCC-AW has a unique maximum at the true zero displacement due to these voxels being down-weighted.

Fig. 2.
figure 2

2D numerical phantoms. (a–d) A contrast-enhanced lesion near a resection. (e–h) Medial longitudinal fissure. (a, e) Reference image. The outline shows the fixed patch. (b, f) Floating image. The inner and outer outline show the moving patch at zero and maximum displacement, respectively. (c, g) Similarity profile of LNCC as a function of displacement. (d, h) Same for LNCC-AW.

Full size image

3.2 Recovery of a 3D Synthetic Deformation

We perform a registration experiment on a BrainWeb dataset. The reference image is made by inserting a synthetic resection cavity in the right temporal lobe. The floating image is resampled using B-spline interpolation from the BrainWeb image using a synthetic sinusoidal deforming field (period 100 mm in all directions, displacement amplitude 4 mm). The voxel intensities are normalized to 0–1 range. We use 5 discrete registration grid levels with grid spacing (7, 6, 5, 4 and 3 voxels), search radius (6, 5, 4, 3 and 2 control point grid spacings) and discretization step (5, 4, 3, 2 and 1 voxels). The floating image is updated between levels using B-spline interpolation. We run the scheme for LNCC (\(\beta = 5\,\mathrm {mm}\)) and twice for LNCC-AW (\(\beta =5\,\mathrm {mm},\alpha = 0.30\) and \(\alpha = 0.10\)).

We quantify registration accuracy using landmarks found in the reference using 3D-SIFT [12]. We include 43 landmarks from a 2 cm region from the resection margin. We propagate the landmarks using the true and recovered deformations. The target registration error (TRE) is shown in Fig. 3(c). TRE for LNCC-AW is significantly lower (for both \(\alpha = 0.30\) and \(\alpha = 0.10\)) than for LNCC (paired t-tests, \(p<0.001\)). The log of Jacobian determinant maps for true and recovered deformations are shown in Fig. 3(d–g). The deformations recovered using LNCC-AW follow the true deformation closer than using LNCC.

Fig. 3.
figure 3

Axial view of 3D BrainWeb based phantom. (a) Reference image (inserted resection). (b) Floating image (synthetic deformation). (c) Target registration error. (d) Map of log Jacobian determinant for ground truth deformation (forward field). (e) Same map for fields recovered using LNCC, (f) LNCC-AW with \(\alpha =0.30\) and (g) LNCC-AW with \(\alpha =0.10\).

Full size image

3.3 Evaluation on an iMRI Surgical Dataset

We validate the measure on 12 cases of anterior temporal lobe resection. The dataset is described in [2]. We skull-strip the pre- and the intraoperative image, normalize the 1st–99th intensity percentile linearly to the range 0–1, crop the intraoperative image to only contain the brain, resample the intraoperative image to a resolution \(1.1\times 1.1\times 1.1\) mm and register the preoperative image affinely to the intraoperative reference [8]. We use the bilateral filter as per Eq. 2 on the reference and floating image pair in order to generate a guidance image pair using settings (\(\beta = 2.2\) mm, \(\alpha =0.03\)) that we found to produce a mild smoothing in homogeneous areas whilst preserving edges. We perform a non-rigid registration for LNCC (\(\beta = 5.5\,\mathrm {mm}\)) and two non-rigid registrations for LNCC-AW (\(\beta = 5.5\,\mathrm {mm},\alpha = 0.30\) and \(\alpha = 0.10\)), using the guidance image pair to construct the weights in Eq. 3. The discrete optimization parameters are identical as in Sect. 3.2 (in voxels). The registration takes approx. 10 h per subject using 4 threads on a computing cluster node.

For each case, we annotate 50–60 landmarks pairs in the pre/intraoperative image a few cm from the resection margin. We propagate the landmarks using the recovered deformations. The mean TRE for all cases is shown in Fig. 4(a) and is significantly lower for registrations based on LNCC-AW with \(\alpha = 0.30\) (paired t-test, \(p = 0.0236\)) and LNCC-AW with \(\alpha = 0.10\) (\(p = 0.0054\)), respectively, than for registrations based on LNCC. The effect size is below the image resolution, potentially as few reliably identifiable landmark pairs exist near the resection margin. We evaluate the smoothness of the recovered deformations and assess the absolute log Jacobian determinant maps in a region of interest (ROI) in the brain less than 2 cm from the base of the resection cavity (located in iMRI). The means within the ROI are shown in Fig. 4(b) and are significantly lower for LNCC-AW with \(\alpha = 0.30\) (paired t-test, \(p = 0.0133\)) and for LNCC-AW with \(\alpha = 0.10\) (\(p<0.001\)), respectively, than for LNCC.

Fig. 4.
figure 4

Registration results for 12 iMRI cases. (a) Target registration error. (b) Mean (in vicinity of the resection) of abs. log Jacobian determinant map.

Full size image

3.4 Segmentation Propagation Experiment

We explore how the adaptive weighting affects registration accuracy for brain structures. We use a database of 35 T1w scans with parcellations of 140 key structures provided by Neuromorphometrics for the MICCAI 2012 Grand Challenge and Workshop on Multi-Atlas LabelingFootnote 1. We normalize image intensities and use each image as a reference image and the remaining images as floating images. For each of the 1190 image pairs, we perform affine registration and non-rigid registrations using LNCC (\(\beta =5\,\mathrm {mm}\)) and LNCC-AW (\(\beta =5\,\mathrm {mm}, \alpha =0.10\) only) using discrete registration parameters as above. We propagate the floating image segmentations using nearest-neighbor interpolation and calculate Dice score for each label. The average Dice score for 1190 affine registration image pairs is \(0.422\,\pm \,0.00187\), for LNCC based non-rigid registrations it is \(0.517\,\pm \,0.0101\) and for LNCC-AW based registrations it is \(0.526 \pm 0.00947\). Average Dice score is significantly higher when using LNCC-AW than LNCC (\(p<10^{-6}\)).

4 Discussion and Conclusion

We introduced bilateral adaptive weighting into a local similarity measure (LNCC). The modification facilitated a more accurate landmark localization in several T1w registration experiments. In a study on clinical iMRI data, we recovered a smoother deformation near the resection margin, which is biomechanically more plausible and potentially enables more accurate guidance near the resection margin. The brain shift we assessed arose from CSF leakage and postural drainage, but in principle our approach can improve accuracy near distinct intensity edges at margins of tumors, collapsed cysts or haematomas from bleeding into the brain, which should be confirmed in a future study.

Unoptimized bilateral weighting introduces a time bottleneck that precludes intraoperative application. We note that the discrete optimization steps collectively take approx. one minute. However, powerful options are open toward optimizing the bilateral weighting, such as guided image filtering [5]. The proposed method could be extended to a multi-channel local similarity measure such as LCCA [7]. A related approach to ours is to constrain the deforming field using bilateral filtering [10] and a unified scheme should be investigated. The analytical gradient could potentially be derived using the approach of [3].