Alignment of Cortical Vessels viewed through the Surgical Microscope with Preoperative Imaging to Compensate for Brain Shift

2.1. Method overview

Our approach updates the pre-operative neurosurgical plan (particularly the tumor position) by taking into account the intra-operative brain shift. As shown in the pipeline of figure 1, at any time t, inputs to our registration system include a single monocular image, $\mathcal{I}$, acquired intra-operatively via the surgical microscope and a 3D model of the cortical vessels, $\mathcal{V}$, derived from the pre-operative MRI. A 3D/2D non-rigid registration method is applied to align the vessels to the image. This registration deforms the vessels while preserving its physical consistency. This deformation estimates a shift of the whole brain (global) from a sub-part of the brain surface (local). The local deformation is converted to a global deformation using a biomechanical model that propagates the location deformation of cortical veins to deep structures including vessels and the tumor.

2.2. Data processing and problem formulation

The first stage of our pipeline is processing the input intra-operative image $\mathcal{I}$ and pre-operative vessels $\mathcal{V}$. We chose to rely on vascular network graphs (or centerlines) to perform our non-rigid registration (see figure 2). Indeed, centerlines represent an effective shape abstraction over using the complete vascular network and can be manipulated more easily.

For this preliminary study, extracting centerlines from $\mathcal{I}$ is done manually, by interactively selecting a set of points representing the starting and ending points of vessel segments. More advanced methods exist, either supervised or unsupervised^22,23 and could be easily plugged in our pipeline. We denote $U=\left\{u_j\in {ℝ}^2\right\}$, the set of 2D points corresponding to the intra-operative centerlines. Extracting the centerlines from $\mathcal{V}$ is done by first segementing the MRI scans and building a triangulated surface mesh of the vessels. This mesh is then iteratively thinned using mean curvature skeletonization²⁴ to obtain the final centerlines that we denote $V=\left\{v_i\in {ℝ}^3\right\}$.

Our goal is to find the transformation function Φ so that U = Φ_(П,R,t,δu)(V) where П is the camera matrix that transforms camera coordinates to world coordinates, R and t are respectively the rotation and translation matrices representing rigid transformation components and δu is the non-rigid transformation component.

2.3. Modelling the vessels and its surrounding tissues

Let S be a 3D triangulated surface mesh of the brain generated from pre-operative MRI scans that includes the parenchyma, internal structures such as vessels or tumors, and constraints from surrounding anatomical information such as the skull.²¹ The mesh S allows us to build a physical model incorporating tissue properties and biomechanical behavior (see Figure 3). This physical model is characterized by a geometry S and a stiffness matrix K that describe physical properties such as tissue elasticity, damping or viscosity. In order to account for tissue heterogeneity, the stiffness matrix of the whole brain is a mechanical composition of paranchyma and vessel stiffness properties. This mechanical coupling permits the propagation of vessel deformation to the paranchyma and the tumor. In practice, the global stiffness matrix is built so that each element of the parenchyma model that intersects a vessel centerline has a higher stiffness.

The deformation is specified by its stress-strain relationship. The brain deformation is computed using a non-linear geometric finite element method, with a linear constitutive law,²⁵ where the finite element model is meshed with linear hexahedrons. This allows for rotation in the model, while relying on a linear expression of the stress-strain relationship. The equation relating the external forces f to the nodal displacements x can be written as f = K(x)δx.

2.4. Energy-based Non-rigid Registration

As most of the registration algorithms, our non-registration has an intialization phase. This consists of a rigid alignment, done manually, that permits to initialize the rotation matrix R and translation vector t. With R and t set, a correspondences vector $C=\left\{c_i\in \mathbb{N}\right\}$ can be built. This vector is built so that if a source point v_i corresponds to a target point u_j then c_i = j for each point of the source set. Having the biomechanical model on one hand and the set of correspondences C on the other hand, we can update the pre-operative 3D mesh S by registering the pre-operative centerlines V with the intra-operative centerlines U that minimizes the physical energy of the biomechanical model. This leads to the following energy minimization expression:

where κ is the stiffness coefficient that ties the source and target centerlines to compute the bending energy. Note the subscript c_i that denotes the correspondence indices between the two point sets. The image fitting term can be considered as sightlines coming from the camera position to the 2D centerlines that enforces the physical model to fit the image.

3.1. Simulated Data

We tested our method on synthetic data of a complete human brain. The vessel centerlines were extracted from the vessels and linearly down-sampled to permit a point-to-point vessel matching. The number of centerlines nodes after down-sampling is 102 nodes. The biomechanical brain model is derived from the triangular mesh using Delaunay volume triangulation. In our tests, both hemispheres are meshed with approximately 1000 tetrahedrons. The stiffness matrix is built with a Young’s modulus E set to 3000 Pa and Poisson’s ratio to 0.45. We used the SOFA framework to build the biomechanical model.²⁶

We quantitatively evaluated our method by measuring the registration error on tumors. We considered 3 tumors at different locations close to the surface (see figure 4-a). We simulated the brain shift under two different scenarios: the first one mimicked the brain sliding towards the back of the skull under gravity as illustrated in figure 4-b and the second simulates a protrusion deformation that can occur due to brain swelling after the craniotomy as illustrated in figure 4-c. From the simulated data we extracted the 3D vessel centerlines to perform the matching where we added Gaussian noise with a standard deviation of 5 mm and a 5 clustering decimation on a subset of points composing the centerlines. We calculated the error as the Euclidean distance between the center of mass of the simulated tumors and the ones measured using our method. The results reported in figure 5 show the registration error for each tumor. Depending on the location of the tumors, we measured a distance error ranging from 0.63 mm to 2 mm with a vessel correspondence of 50% which corresponds to a 81% and 68% brain-shift correction respectively (The brain-shift correction is measured as the difference between the initial error and the corrected error divided by the initial error).

In addition to tumor locations, we measured the registration error on the brain surface using Hausdorff distance. As shown in Figure 6, the maximum measured Hausdorff distance was 3.2 mm with a vessel correspondence of 50%. We can also observe that gravity-induced deformation generates less errors than deformations caused by brain swelling at the site of the craniotomy.

We furthermore downgraded the matching quality by decreasing the number of vessel correspondences, from 90% to 10% to demonstrate the capability of the model to interpolate deformation where no image information is available. As shown in figures 5 and 6 the registration error on the tumors and the brain surface decreases monotonically as the number of correspondences increases. In the worst case where only 10% of the correspondences are used, we measured a maximum Euclidean distance between tumors’ center of masses of 3.61 mm for tumor 1, 2.20 mm for tumor 2 and 4.46 mm for tumor 3, and a maximum Hausdorff distance of 10.19 mm for the brain surface.

We tested our method on a subset of centerlines nodes to get closer to a craniotomy case where only a small part of the brain is visible. This subset is composed of 8 nodes regularly distributed over the centerlines. We plotted in Figure 7 the actual registration error between the estimated tumors’ positions and the ground truth positions for both type of deformation. The errors are more important compared to the global registration but are still acceptable with a minimal distance of 1.08 mm and a maximal distance of 2.61 mm. Tumor 2 located near the visible surface vessel subset exhibits the lowest distance errors while the tumor 3 located deeper w.r.t the selected surface shows the highest distance errors.

3.2. Real Data

We tested our method on real data using images acquired from a Carl Zeiss surgical microscope. The scope is equipped with two cameras with a video frame rate of approximately 30 frames per second. We used only the left image. The pre-operative mesh is first aligned manually on the image to obtain the rigid transform, then, the non-rigid refinement is performed, which will drive the whole brain mesh. The Figure 8 shows the non-registration process and the resulting mesh-to-image overlay.