1 Introduction

Gliomas are primary brain tumors of central nervous system (CNS) [1]. The gliomas arise from the glia which supports the CNS and can permeate to the neighboring areas. They can be categorized in grade from I to IV based on their histological characteristics defined by the World Health Organization (WHO) [2]. The grade I and II gliomas are classified as low-grade gliomas (LGG) and grade III and IV gliomas are classified as high-grade gliomas (HGG). The LGG grow comparatively slowly but due to their infiltrative attribute and threatening behavior, the mean 10-year survival is 30% [1]. It is generally accepted that the resection of the LGG may increase the survival rate [3].

Intra-operative ultrasound (iUS) was first proposed as a potential tool for guiding resection of intracranial tumors in 1980 [4]. The iUS is still generally used because it enables the surgeons to track the brain tissues and surgical tools in a fast, inexpensive, and real-time way. In addition, the gliomas can often be detected in iUS images even when they are not detectable under the microscope. This can promote accurate resection and helps to obtain better surgical results. However, it is difficult to design effective surgical plans without the high-quality image-guidance. One of the principal reason is that the surgical target and the other tissues can be shifted by intra-operative factors such as tissue removal, change of intracranial pressure and drug administration. However, these shifts may not be easily observed in the surgeon’s field of view.

To estimate and rectify for spatial errors resulting from intra-operative brain shifts, registration of pre-operative magnetic resonance imaging (MRI) to iUS image has been suggested [5,6,7]. This approach helps updating the surgical plans under the continuous brain tissue shift in contrast to comparing directly between pre- and intra-operative images. Many algorithms for registration have been proposed in the past years. However, it is technically demanding because of its intrinsic limits such as the differences between modalities and image qualities. Because of these problems, most MRI to iUS registration methods were conducted by using manually selected fiducial points.

In this paper, we proposed an automatic non-linear MRI-iUS registration algorithm using structural skeletons. First, we conducted several pre-processing steps on the MRI and iUS. Then we calculated the structural skeletons of both modalities. Finally, we calculated the displacement fields using the pairs of skeleton for the MRI-iUS registration.

2 Materials and Methods

2.1 Dataset

As this proposal is submitted to the CURIOUS 2018 challenge, we used the Retrospective Evaluation of Cerebral Tumors (RESECT) dataset [8]. The dataset includes pre-operative 3T MRI images including Gadolinium-enhanced T1w and T2 FLAIR scans, iUS images as a 3D volume covering the entire tumor region after craniotomy but before dura opening and the expert-labeled homologous anatomical landmarks, defined on all image modalities. All reconstructed images were acquired from the same subject and were spatially aligned under the same world coordinate space.

The MR protocol included T1w Gadolinium-enhanced sequence (TE = 2.96 ms, TR = 2000 ms, 192 slices, slice thickness = 1 mm, acquisition matrix = 256 × 256, in-plane resolution = 1.0 × 1.0 mm2) and FLAIR sequence (TE = 388 ms, TR = 5000 ms, 192 slices, slice thickness = 1 mm, acquisition matrix = 256 × 256, in-plane resolution = 1.0 × 1.0 mm2), acquired on a 3T MRI scanner with a 20-channel head coil. For subject 2, 14, 15 (Case 2, 14, 15), the pre-operative MRI included T1w sequence (TE = 2.3 ms, TR = 2500 ms, 176 slices, slice thickness = 1 mm, acquisition matrix = 512 × 496, in-plane resolution = 1.0 × 1.0 mm2) and FLAIR sequence (TE = 333 ms, TR = 6000 ms, 176 slices, slice thickness = 1 mm, acquisition matrix = 256 × 224, in-plane resolution = 1.0 × 1.0 mm2) acquired on a 1.5T MRI scanner with a 12-channel head coil.

The iUS images were acquired using the Sonowand Invite neuronavigation system. In most cases, the 12FLA-L linear probe with a frequency range of 6 to 12 MHz and a footprint of 48 × 13 mm2 was used. For smaller tumors, the 12FLA flat linear array probe with a frequency range of 6 to 12 MHz and a footprint of 32 × 11 mm2 was used. The resolution of reconstructed 3D volume varied from 0.14 × 0.14 × 0.14 mm3 to 0.24 × 0.24 × 0.24 mm3 depending on the probe types and imaging depth.

2.2 MRI Pre-processing

First, we corrected non-uniformity to remove field bias from the MRI image. After the bias removal, we obtained the masks of grey matter (GM), white matter (WM) and cerebrospinal fluid (CSF) from the T1w image. After that, we cropped the T1w image and its masks to match with the location and size of iUS image. Then we resampled the cropped image and masks into 0.2 × 0.2 × 0.2 mm3 voxel resolution. Because the intensity distribution of iUS image was well matched with the GM and CSF mask from the MRI image, we used the inversed WM mask instead of using the summation of GM and CSF mask. The small hole of the inversed mask was removed using flood-fill algorithm with connectivity parameter 4 per each z-slice. The output mask was used to calculate the structural skeleton of MRI image.

2.3 iUS Pre-processing

We resampled the iUS image using the cropped and resampled T1w image as reference. As the voxel resolution of reference image was 0.2 × 0.2 × 0.2 mm3, the iUS image resampled to the same resolution. The order of axis in iUS image was matched with the reference image. Then, the resampled iUS image was blurred with 0.5 mm full width at half maximum (FWHM) Gaussian kernel and filtered with 3 × 3 × 3 median kernel under voxel coordinate space. A binary mask image was obtained by applying threshold to the filtered image and was dilated by the ball shaped structure elements with voxel radius 15. The dilated mask was used for calculating the skeleton. We used a semi-automatic intensity-based segmentation algorithm to compute the binary masks for both MRI and iUS. The masks were structurally enhanced using preprocessing steps such as full with half maximum smoothing and morphological operations, which was performed with typical hyperparameters.

2.4 Structural Skeleton

After the refinement of the binary masks, the Euclidean skeletons were calculated for both masks [9, 10]. For details, we inverted the masks and computed the Euclidean distance map. The ridges of the distance map were computed by watershed algorithm. The final skeletons of each masks were obtained by thresholding the gradient of the ridges.

2.5 Registration Using Deformation Fields

We applied the Demon’s algorithm for non-linear registration [11]. This algorithm performed over the entire space of displacement field. At first, a spatial transformation field was initialized. Then, we iterated the following steps until the error converged:

  1. (1)

    Given field s, compute a correspondence update field u by minimizing the error E which can be defined as following Eq. (1).

    $$ \varvec{E}\left( u \right) = \left| {\left| {\varvec{F} - \varvec{M} \cdot \left( {s + u} \right)} \right|} \right|^{2} + \left( {\frac{{\sigma_{i} }}{{\sigma_{x} }}} \right)\left| {\left| u \right|} \right|^{2} $$
    (1)
  2. (2)

    Let c as s + u.

  3. (3)

    Use diffusion-like regularization by conducting the Gaussian smoothing to the accumulated transformation field c

  4. (4)

    Substitute s with the filtered c.

As mentioned above, the variable s accounts for the given spatial transformation field and u for the corresponding update field. The variable F accounts for the fixed iUS image, M for the moving MRI image, \( \sigma_{i} \) for the noise on the image intensity and \( \sigma_{x} \) for a spatial uncertainty on the correspondences.

2.6 Mean Target Registration Error

We measured the Euclidean distances between the MRI landmarks after the registration and the corresponding iUS landmarks to calculate the target registration errors (TREs). All landmarks were averaged per each case.

3 Results

3.1 Crop and Resampling

Figure 1 shows the results images after cropping and resampling. Figure 1(a) and (b) represent the T1w image and the iUS image respectively. Both images were well aligned to each other so that they could be used for the following processes.

Fig. 1.
figure 1

The MRI image and iUS image obtained after crop and resampling of case 1. (a) is the MRI image and (b) is the iUS image.

Full size image

3.2 Estimated Skeleton

Cropped and resampled MRI and iUS images were used to calculate their structural skeletons. Figure 2 shows the skeletons obtained from the MRI and iUS images. Figure 2(a) shows the MRI image and its skeleton, and Fig. 2(b) shows the iUS image and its skeleton.

Fig. 2.
figure 2

The MRI image and iUS image of case 1 with their structural skeleton. (a) The structural skeleton of the MRI was indicated with green line and (b) the structural skeleton of the iUS was indicated with red line. (Color figure online)

Full size image

3.3 Skeleton Registration Using Displacement Fields

Figure 3 shows the results of the registration using the Demon’s deformation algorithm. Figure 3(a) shows the original skeleton of the MRI image. Figure 3(b) and (c) represents the deformation fields along x-axis, top-to-bottom axis in this image, and y-axis, left-to-right axis in this image respectively. The deformation field along z-axis was also calculated but it was not shown in this figure for better visualization. Figure 3(d) is the skeleton of the MRI image after the registration to iUS image. We compared these results with the original iUS image and its skeleton (Fig. 4). The skeleton of the iUS image was shown as the structural line with red color on both Fig. 4(a) and (b). The structural line with green color on Fig. 4(a) represents the skeleton of MRI image before the registration and the same line on Fig. 4(b) represents the skeleton after the registration.

Fig. 3.
figure 3

The results of the registration using the Demon’s deformation algorithm. (a) is the original skeleton of MRI. (b) and (c) shows the displacement field estimated after the registration on the x-axis and the y-axis respectively. (d) is the skeleton of MRI after the MRI-iUS registration.

Full size image
Fig. 4.
figure 4

The demonstration of the registration using skeletons overlapped on the same iUS image of case 1. (a) shows the skeleton of original MRI with red line and the skeleton of iUS image with green line. (b) shows the skeleton of registered MRI with red line and the skeleton of iUS image with green line. (Color figure online)

Full size image

3.4 Image Registration Using Estimated Fields

As we calculated the deformation fields at each location on every axis, we could obtain the final registered MRI image (Fig. 5). Figure 5(a) shows the original MRI image overlapped with the iUS image. The MRI image was visualized with the green color and the iUS image was visualized with red color. Figure 5(b) shows the registered MRI image overlapped with the iUS image.

Fig. 5.
figure 5

The demonstration of the registration. (a) The MRI image before the registration was shown as green and the original iUS image was shown as red. (b) The MRI image after the registration was shown as green and the original iUS was shown as red. (Color figure online)

Full size image

3.5 Landmark Evaluation

The Euclidean distances between two corresponding landmarks were measured. Then, the mTREs of all cases were measured to evaluate the results (Table 1).

Table 1. The mean target registration error for the all 22 cases.
Full size table