3D Rigid Registration of Intraoperative Ultrasound and Preoperative MR Brain Images Based on Hyperechogenic Structures

The registration of intraoperative ultrasound (US) images with preoperative magnetic resonance (MR) images is a challenging problem due to the difference of information contained in each image modality. To overcome this difficulty, we introduce a new probabilistic function based on the matching of cerebral hyperechogenic structures. In brain imaging, these structures are the liquid interfaces such as the cerebral falx and the sulci, and the lesions when the corresponding tissue is hyperechogenic. The registration procedure is achieved by maximizing the joint probability for a voxel to be included in hyperechogenic structures in both modalities. Experiments were carried out on real datasets acquired during neurosurgical procedures. The proposed validation framework is based on (i) visual assessment, (ii) manual expert estimations , and (iii) a robustness study. Results show that the proposed method (i) is visually efficient, (ii) produces no statistically different registration accuracy compared to manual-based expert registration, and (iii) converges robustly. Finally, the computation time required by our method is compatible with intraoperative use.


Introduction
Due to its low cost, its real-time imaging capabilities, and its noninvasive nature, ultrasound (US) imaging has become a popular modality. These attributes have been used for a large number of clinical applications. In neurosurgery, ultrasound imaging has been employed in many cases of brain examinations over the last two decades [1]. Several studies demonstrated that ultrasonography can be used for locating tumors, defining their margins, differentiating their internal characteristics, and for detecting of brain shift and residual tumoral tissues [2]. At present, 3D US imaging is integrated within the neuronavigation systems to provide a useful and efficient intraoperative tool [3]. Ultrasound imaging has also been shown to be a promising method for quantifying and for correcting brain shift in Image-Guided Neurosurgery (IGNS) [4][5][6][7][8][9][10][11][12][13][14].
During a neurosurgical procedure, the ultrasound probe is tracked by the neuronavigation system which computes the 3D positions and orientations of the B-scans. Matching between the intraoperative US images and the preoperative MR image is ensured by a rigid registration. In phantom [6] and animal studies [11], the matching accuracy between intraoperative B-scans and preoperative images has been quantified between 1.5 mm and 3 mm. Nevertheless, in clinical context, the matching error can reach 10 mm (see Table 1). This error includes tool calibration errors (the position localizer and the US probe), tool localization errors (tracking system error), and registration errors from the neuronavigation system. International Journal of Biomedical Imaging Table 1: Manual estimation of the initial error in mm (i.e., error of the registration given by the neuronavigation system).
A priori estimation of the registration: initial error in mm Mean (std) Expert  Registration approaches based on classical image similarity measures such as the Sum Square Difference (SSD), Mutual Information (MI), or Correlation Ratio (CR) are known to fail to robustly register MR and US images [15]. Therefore, other options have been investigated.
(a) Landmark-based registration represents the majority of the approaches [6,7,9,13,14,16]. The motivation is bound to the difficulty of finding a function matching US image intensities with MR image intensities. These methods are based on the matching of manually defined points [7], lines representing the vascular system [6,13,14,16], or cortical surface [9].
(b) Intensity-based approaches using histogram-based similarity measures tend to overcome the problem by preprocessing the images in order to register more similar images [4,17].
(c) By introducing the Bivariate Correlation Ratio (BCR), Roche et al. [15] incorporated the transformation of MR to pseudo-US image as a function into the similarity measure.
In this paper, we propose a new objective function based on the matching of the cerebral hyperechogenic structures such as sulci and the cerebral falx, and the lesion when the corresponding tissue is hyperechogenic. The registration is achieved by maximizing the correlation value between the US image and the probabilistic map of hyperechogenic structures estimated from MR image. The proposed method is thus a compromise between landmark and intensity-based approaches.
(i) As with landmark-based approaches, only regions considered as relevant are used to drive the registration procedure. In our method, these regions are the hyperechogenic structures of the brain.
(ii) As with intensity-based methods, the proposed approach does not require segmentation of the US image which is a challenging problem.

Method
Overview. The scheme of the overall workflow is presented in Figure 1. First, the "hyperechogenic" structures present in MR image (i.e., the structures visible in MR image expected to be hyperechogenic in intraoperative US) are detected with the MLvv operator [18,19]. In brain imaging, these structures are the liquid interfaces such as the cerebral falx and the sulci, in addition to the lesions when the corresponding tissue is hyperechogenic (e.g., cavernoma or glioma). The curvature-based MLvv operator was first introduced in [18,19] before being used to detect the sulci and the cerebral falx in [20][21][22]. The US image and the probability map of the hyperechogenic structures extracted from MR image are then registered by maximizing the probability for a voxel to be included in hyperechogenic structures in both modalities. Contrary to histogram-based approaches that match all the information in both images, the proposed approach consists of matching only hyperechogenic structures [23], which makes it more robust to artefacts such as acoustic shadows. Indeed, in US imaging, the bright areas provide information on the underlying structures whereas the dark areas can correspond to the underlying anatomical structure or acoustic shadows [24]. Moreover, the accuracy of sulci matching is an important issue since these structures are used by the neurosurgeon during the neurosurgical procedure [25]. Finally, by using the natural property of US imaging to detect the hyperechogenic structures, the method does not require segmentation of the US image. This way, the method is less sensitive to error of US image segmentation and is less time consuming during the intraoperative stage.

Probabilistic Objective Function.
The proposed registration process is based on the estimation of the transformation T maximizing the joint probability for a voxel X = (x, y, z) to be included in hyperechogenic structures in both modalities: where p(X ∈ Φ US ) is the probability for X to be included in an hyperechogenic structure from the US image and p(X ∈ Φ MR ) is the probability for X to be included in an hyperechogenic structure (in the sense of the ultrasound image) from the T1-w MR image. Assuming that the probabilities are independent, we can write Our objective function can be viewed as the maximization of the correlation value between the two probability maps of hyperechogenic structures extracted from both modalities.

Construction of the Probability Maps.
In order to construct the probability maps, we define a function f matching the intensity of both the US image and the MR image with the probability for X to be included in hyperechogenic structures: where u : Ω → R is an image defined on Ω.  Figure 1: Illustration of the performed workflow to achieve the registration. The skull stripping, the denoising, the MLvv computation, and the segmentation of lesion are performed before the neurosurgical procedure. The 3D reconstruction of intraoperative volume, the reslicing of the MR map, and the estimation of the transformation are then performed during the neurosurgical procedure.

Intraoperative US Image.
For the intraoperative US image U, by definition f is the identity function: The intensity of U is only scaled between 0 and 1 during surgery to fit with our probabilistic framework.

Preoperative MR Image.
For the preoperative MR image V , the evaluation of f is done prior to surgery and is based on both the detection of the liquid interfaces with the MLvv operator and the segmentation of the pathological tissues.
The Lvv operator is a robust intensity-based curvature detector [18] based on the first and second derivatives of the image intensities. The first and second derivatives are combined to obtain an operator less sensitive to flat areas with low gradients. This kind of operator is used to detect ridge-like features in images, with negative value for crests in the intensity domain and positive value for valleys in the intensity domain. In [19], the MLvv has been proposed for multimodal registration of CT and MR images. In our case, as in [20][21][22], the MLvv is used to extract the hyperechogenic structures (sulci, cerebral falx) from T1-w MR image. Overall, curvature information has been used by several other authors to characterize cortical features [26][27][28]. Most of these methods are based on geodesic curvature computed on cortical surfaces.
In T1-w MR images, the sulci are valleys (negative ridges) in the intensity domain. By using the positive values of MLvv, the sulci and the cerebral falx can be efficiently detected [20][21][22]. Finally, our function f is defined as where I M is the indicator function for the set M: such that X belongs to the lesional tissue}.
As for the US intensities, the positive values of the MLvv are scaled between 0 and 1. The MLvv operator is defined in 3D as where is the probability given to X in the segmentation of pathological tissue M 2 . Ψ is used to incorporate a priori on pathology. For pathological tissue such as cavernoma or low-grade glioma, Ψ(X) is high since these tissues are hyperechogenic.

Preprocessing of the MR Data before Surgery.
First, skull stripping is performed from the T1-w MRI sequence [29]. We choose to remove the skull prior to MLvv computation because this structure does not appear in the area of the craniotomy. The raw MR images are then denoised using an optimized Non-Local Means filter (https://www.irisa.fr/visages/benchmarks/) [30] before International Journal of Biomedical Imaging In our experiments, the segmentation of pathology was manually performed by the neuroanatomist before the surgical procedure (see Figure 1). The computational time required by preprocessing steps performed during preoperative stage was 4 minutes for skull stripping and 3 minutes for denoising on a Pentium M 2 GHz. In addition, International Journal of Biomedical Imaging B-scans and the probe positions. From the 2D B-scans and their positions, a 3D volume was reconstructed with the Probe Trajectory method [31]. The experiments were carried out on 3 patients. For each patient, a sequence of images was acquired before opening the dura. Some studies have considered quantitative measurement of brain shift during surgical procedures and showed that nonsignificant displacement occurred before dura opening [32,33]. Thus, we assumed that the transformation between intraoperative US and preoperative MR was rigid.

MR-US Registration of the Neuronavigation System.
During all the neurosurgical procedure, the coordinate system of the preoperative MR image and the coordinate system of the intraoperative field are related by a rigid registration. The rigid registration of the neuronavigation system is based on surface matching between the preoperative MR image and the position of points acquired on the patient's head with the position localizer. First, the skin is extracted from the MR image by manual thresholding. A cloud of points is then continuously acquired on the patient's head close to the eyes region by moving the position localizer. Following this, one point is acquired on each ear with another point on the extremity of the patient's nose. Finally, the neuronavigation system performs a points to surface matching. According to phantom and animal studies, the errors in probe calibration, 3D localization of the probe, and rigid registration performed by the neuronavigation system lead to a global error less than 3 mm [6,10,11]. The error due to the 3D localization of the probe is estimated to 0.35 mm for each marker on a tool from the manufacturer [34]. The error due to the calibration is generally estimated around 1.5 mm [6,11]. In our case, the probe was calibrated with a Z-wire phantom by the manufacturer. Finally, the error due to rigid registration performed by the neuronavigation system has been estimated to be around 1.5 mm in [11].

Pathology of the Patients.
In this study, hyperechogenic pathologies such as cavernoma (patient 1, see Figure 5 and patient 2, see Figure 6) and low-grade glioma (patient 3, see Figure 7) were considered. In T1-w MR images, the central part of cavernoma is usually heterogeneous (hyper-and hyposignal) and the outlying area appears in hyposignal. The low-grade gliomas are more homogeneous and appear in hyposignal in T1-w MR images. In US images, numerous studies showed that all solid brain tumors, metastatic brain lesions, and cavernomas exhibited echogenicity [35][36][37][38][39]. For brain gliomas, the higher its grade (more malignant), the more echogenic it is in US and the less homogeneous it appears. In our study, the corresponding lesional tissues were considered both homogeneous and hyperechogenic in US images. As such, Ψ(X) was set to 1 for all segmentation of pathological tissue M 2 (see (5)). Typical examples of intraoperative images and probability maps are presented in Figures 2, 3, and 4.

Parameter
Settings. The maximization of the joint probability (see (2)) is performed within a multiresolution procedure using the simplex algorithm [40]. During the experiments, the parameters of the simplex algorithm were tolerance = 0.1, stepsize = 1.5, and maximum number of iterations = 100. The coarsest resolution corresponded to the original volumes downsampled by a factor 3 and the finest resolution was that of the original volumes. The registration procedures take less than two minutes on Intel Pentium M at 2 GHz.
As with most derivative-based operators, the MLvv operator uses Gaussian kernel to compute the image derivatives. In [18], the authors showed that the convolution of the image with a derivative Gaussian kernel provides a wellposed approach of the differentiation problem. The standard deviation σ of the Gaussian kernel is called the image scale. This parameter has been shown as very stable for MR image sulci segmentation on numerous works [20][21][22], and thus, no tuning has been done for this parameter throughout our study. In our experiments, an image scale of 2 voxels has been used to compute the MLvv values. This value is consistent with other works [20][21][22] conducted on brain cortical segmentation where the scale parameter was always kept in this range.

Evaluation Framework.
In order to evaluate our method, a validation framework with different approaches is proposed.
(i) First, a visual assessment is proposed.
(ii) Second, a manual validation by experts is presented.
This validation is divided in two parts: a point-based estimation of the rigid registration by 3 experts for the 3 patients and an evaluation of the residual error by all experts for 1 patient (postregistration error).
(iii) Third, a study on convergence robustness was carried out.
The expert manual validation was difficult due to the time required. For each expert, 4 hours were required to perform the a priori estimation of the transformation for 3 patients.

Visual Assessment.
The visual assessment remains a valuable indicator of the registration accuracy. In [41], the observer discernibility of registration errors has been estimated around 0.2 mm. A study on visual inspection for image registration assessment can be found in [42]. In our paper, we propose an overlay of US and MR images before and after registration to assess the registration accuracy.

Validation by Experts.
First, the experts manually evaluate the rigid transformation between the intraoperative International Journal of Biomedical Imaging 7 x-y plane x-z plane y-z plane US and the MR image resliced with the rigid transformation given by the neuronavigation system. This estimation is denoted as a priori estimation of the registration. From this a priori estimation, the initial error (i.e., after the registration performed by the neuronavigation system) and the Target Registration Error (TRE) can be computed. The a priori estimation of the registration is used to show that there are no statistical differences between the expert-based transformations and the transformation estimated with our method in terms of the TRE. The experts estimate the residual error after rigid registration based on a given transformation (either by our method or the point-based expert registrations). This estimation is called a posteriori evaluation of the residual error and is designed to show that experts do not detect significant differences when they inspect the registered volumes with our method or with their own manually defined transformations.

A Priori Estimation of the Registration
Point Picking. The a priori estimation of the registration is based on the location of ten points in the US image and the ten corresponding points in the MR image: each expert defines a set of point in the 3D reconstruction of the intraoperative ultrasound and its corresponding landmark in the resliced MR image. The resliced MR image is obtained with the rigid registration given by the neuronavigation system and has the same resolutions, dimensions, and field of view as the reconstructed US image. During the experiments, the experts used three orthogonal 2D views to define homologous points in the 3D volumes. For each volume, the visualization software was run independently, with the cursors in the two volumes unlinked. Each expert was allowed to choose their set of homologous points.
Initial Error. The initial error is computed by using the mean Euclidean distance between the homologous points defined by the experts in both modalities. The three samples (one per expert) containing the ten error values (one per point) are compared by using a Kruskal-Wallis test.

Target Registration Error.
A leave-one-out procedure is used to compute the TRE of each point (i.e., Euclidean distance between homologous point after rigid transformation). First, one of the ten homologous points is removed from the set of points. The nine remaining homologous points are then used to compute a rigid transformation in the least squares sense. Finally, this rigid transformation is used to compute the TRE of the initially removed point. This procedure is repeated for all the ten points. The final TRE is the mean TRE over all the points. For each patient, the expert-based TRE and the TRE obtained with our method are compared by using a nonparametric Kruskal-Wallis test.  The low intensities of US images are in green and the high in red. In this case, the acoustic shadow artifact was present on the US image. The signal below the lesion was totally dark as shown in Figure 3. The proposed approach allowed to overcome these artifacts without specific detection of the shadows.

A Posteriori Evaluation of the Residual Error
Point Picking. First, the patient images are registered using several transformations. These transformations are (i) the three expert-based transformations ( T 1 , T 2 , T 3 ), (ii) the rigid transformation obtained with our method ( T), and (iii) the transformation computed using all the thirty points defined by the three experts T all . The experts then define ten homologous points on the registered volumes. This procedure is performed for the five studied registrations on the patient 2 dataset. The positions of the points are fixed for all the experts.
Residual Error. As for initial error, the final error or residual error is simply obtained by computing the mean Euclidean distance between the homologous points defined by the experts in both modalities. The statistical comparison of the residual errors is performed on the five samples (one per transformation) containing ten errors values (one per point) with a Kruskal-Wallis test.

Robustness Study.
First, the US and resliced MR images of patient 2 are registered with the transformation T all . Then, 100 rigid transformations are randomly generated with a translation along each axis uniformly distributed between 0 and 5 mm and with a rotation around each axis uniformly distributed between 0 and 5 degrees. Finally, each transformation is applied to the resliced MR image before performing registrations with the proposed method. The warping index ω [43] is used to compute the distance between the estimated transformation by the registration process T and the true transformation T: where · 2 is the L 2 -norm. The success rate is estimated by considering a success as a registration with a warping index inferior to 3.5 mm. This threshold has been chosen close to the upper bound of the TRE estimated by the experts (see distribution for patient 2 in Figure 9). Contrary to TRE estimated over selected points, the warping index is computed as the average error between the volumes over all the voxels.

Visual Assessment.
The registration results are first displayed for visual assessment. The results obtained with our method are presented in Figures 5, 6, and 7. For patient 1 (see Figure 5), even if the lesion was not entirely included in the US volume, the proposed registration procedure converged efficiently. For patient 2 (see Figure 6), acoustic shadows are present on the US image. The signal below the lesion tends International Journal of Biomedical Imaging 9 x-y plane x-z plane y-z plane to zero. The proposed approach overcomes these artifacts without specific detection of the shadows. For patient 3 (see Figure 7), despite the large size of the low-grade glioma and the limited field of view, our approach performed well. Table 1 presents the estimated initial error for the three patients by the three experts. The P value of the Kruskal-Wallis test showed that there was no significant difference between the expert estimations. Table 1 also shows the interindividual variability for the same measure between the experts. Figure 8 summarizes the distribution of the error. The estimated initial errors are significantly higher than values given in [6,11] (<3 mm) or by the manufacturer (<1.5 mm). It is important to note that the ultrasound images used in our experiments were acquired in clinical context during a neurosurgical operation. The real neurosurgery context is likely more difficult than phantom and animal studies. Table 2 shows the TRE estimated by each expert, for each patient dataset. In all the cases, there were no statistically significant differences between the TRE obtained with expert-based estimations and the TRE obtained with our method. Figure 9 shows the result of the Kruskal-Wallis test. In all the cases, the experts and our method provided consistent results. Table 3 shows the expert-based estimation of the a posteriori residual error of the different registrations (manual-based T and automatic T) proposed for patient 2. The Kruskal-Wallis test shows that the errors associated with the transformations ( T 1 , T 2 , T 3 , T all ) and T are not significantly different. Figure 10 shows the statistical distribution of the residual error for each transformation compared. Finally, the experts failed to detect significant differences between the manualbased registrations and our automatic registration. The residual error estimated by experts is around 1-1.5 mm for all the transformations. Table 4 shows the robustness and the warping index results obtained during the experiment. The proposed method obtained 92% of success rate with a mean     Table 2. Figure 11 shows the distribution of the warping index.

Discussion and Conclusion
This paper presents a new framework for the 3D rigid registration of US and T1-w MR brain images. In order to address this challenging problem, we propose an innovative probabilistic objective function that maximizes the joint probability of the (i) a priori most probable locations of hyperechogenic structure in the preoperative MR image and (ii) the highest intensities of the intraoperative US images. We show that the proposed method enables a robust registration of MR and US images in a computational time compatible with clinical use. All our experiments were carried out on real intraoperative data. The expert-based quantitative study shows that our method produces no statistically different registration compared to the a priori estimation of the registration by the experts. Moreover, the a posteriori estimation of the residual registration error shows that the experts failed to detect differences between manual registration and our automatic registration. During our experiments, manual segmentation has been used to build the probability map. This segmentation is always available, since the neurosurgeon performed it before the surgery. In this paper, the used segmentations were the segmentations dedicated to the neurosurgery. However, the segmentation of the lesion could be automated [44,45], and the different parts of pathologies (lesion, coagulated blood, cyst, necrotic tissue, etc.) could be defined. Through this, the simple model of homogeneous hyperechogenic lesion used in our experiment could be improved by using different hyperechogenic levels to the different pathological tissues. To evaluate the robustness of our method to heterogeneous lesion, more datasets are needed although this situation was present in case of patient 1.
The proposed method is related to the segmentation accuracy of the tumor in preoperative MR images. Although the segmentation of the MR image is not considered as difficult, this step may introduce some errors. Our experiments showed that the proposed method produced consistent results with manual segmentation used in clinical routine.
The presented clinical datasets showed that our method is robust to some discrepancies between the features present in both US and MRI probability maps. Based on the correlation of maps where only regions considered as relevant are used to drive the registration procedure, our method is able to deal with partially missing information resulting from a limited field of view or acoustic shadows. In case of patient 1 (see Figures 2 and 5), only a subpart of the lesion was visible in the reconstructed US image. In case of patient 2 (see Figures 3 and 6) the acoustic shadow below the lesion reduced information around ventricle in US. Moreover, the information derived from sulci was much more present in MR maps than in US map. Finally, in case of patient 3 (see Figures 4 and 7) the limited field of view and the large size of glioma reduced the importance of sulcal information. However, experiments using only the segmentation of the lesion or only sulcal information derived from Mlvv operator failed to provide satisfactory registration. This seems to indicate that a certain amount of homologous features has to be present in both probability maps to enable the method working.
Finally, in our opinion, the proposed approach relies on a similar and complementary idea to the vessel-based method proposed by Reinertsen et al. [13,14]. Indeed, in both cases,  an implicit segmentation of salient features in US images (hyperechogenic structures in B-mode or vessels in Doppler) is matched with corresponding structures detected in MR images. Only the selected salient features differ between the methods. In [13,14], the method utilizes vessels extracted from Doppler US images and their segmentation from MR images. Therefore, both methods have the advantage of not requiring segmentation of the US image and also being robust to US artefacts. However, the extraction of the vessel centerlines from MR images is a challenging problem and requires extensive processing. Our method is dedicated to brain US imaging since MLvv operator is relevant for sulci and cerebral falx detection. As such, the application of the proposed framework to another body part requires adaptation of the hyperechogenic structure detection. Moreover, if T2-w MR image or another sequence is used as preoperative MR image, the selected values of the MLvv need to be adapted. Since the final aim of this US/MR registration method is to compensate for the brainshift, further works will investigate extension of our probabilistic objective function to non-rigid deformations.