Projection-Based Collision Detection Algorithm for Stereoelectroencephalography Electrode Risk Assessment and Re-Planning

Surgical planning is crucial to Stereoelectroencephalography (SEEG), a minimally invasive procedure that requires clinicians to operate with no direct view of the brain. Decisions making involves different clinical specialties and requires analysis of multiple multimodal datasets. We present a DepthMap tool designed to localize, measure, and visualize surgical risk, and an AlternativeFinder tool, designed to search for alternative trajectories maintaining adherence to the initial trajectory with three different re-planning strategies: similar entry, similar target, or parallel trajectory. The two tools transform the 3D problem into the 2D domain using projective geometry and distance mapping. Both use the graphics processing unit (GPU) to create a 2D depth image used by DepthMap for measurement and visualization, and by AlternativeFinder to find alternative trajectories. Tools were tested with 12 SEEG cases using digital subtraction angiography. DepthMap was used to measure vessel distance. AlternativeFinder was then used to search for alternatives. Computation time and displacements of the entry and target points for each trajectory and adherence strategy were recorded. The DepthMap tool found vessels in 118 initial trajectories (out of 145). Ninety alternative trajectories were found to meet the required avascular constraints (average 820K alternatives evaluated per initial trajectory). The average computation time was 449 ms per initial trajectory (77 ms when alternatives were found). The tools presented helped clinicians examine and re-plan SEEG trajectories to avoid vascular risks using three adherence strategies. Quantitative measurement of the adherence shows the potential of this tool for clinical use.


I. INTRODUCTION
Around one-third of epileptic patients do not respond to antiepileptic drugs [1] and could be treated with epilepsy surgery [2]. Stereoelectroencephalography (SEEG) is a stereotactic technique used in complex drug-resistant The associate editor coordinating the review of this manuscript and approving it for publication was Songwen Pei . epilepsy cases to evaluate a possible surgical intervention [3]- [6]. This is an invasive technique that consists of the placement of several rectilinear depth electrodes. It is used to measure electric fields inside the brain with the multiple contacts in a row placed in each trajectory [7].
SEEG is an invasive technique that is performed on selected patients to answer clinical questions when a clear chance for a surgical solution is feasible, and previous noninvasive information is not enough to elucidate proposed options for depicting the epileptogenic zone (EZ) and epileptic network. This noninvasive patient information includes seizure semiology, especially the initial symptoms that appear at seizure start, neuropsychology testing, MRI, video EEG recordings, and frequently PET scanning and ictal and interictal (and SISCOM) information. It may be difficult to understand the onset area from video-EEG surface electrodes recordings (mostly when it is deep as for instance if it originates from the insula, the orbital frontal, the mesial structures as the cingulum), or to elucidate the side of seizure onset in a rapidly generalizing seizure.
The implantation of SEEG requires careful planning of the trajectories. The SEEG trajectory plans evolve from an initial draft which evolves in successive degrees of refinement until a final surgical plan is obtained. SEEG has a strong multidisciplinary component [8], [9], with epileptology and neurosurgery being the most involved clinical specialties. In this document, we use the term initial plan to denote a concept similar to the 'epileptologist plan' [10], [11] or 'the strategy of the implantation' [12], [13] mentioned in related literature.
The initial steps of the SEEG workflow that is carried out at Hospital del Mar, Barcelona requires the epileptologists to suggest an initial trajectory plan based on an EZ localization hypothesis based on non-invasive clinical information. Following that initial plan from the epileptologist, neurosurgeons modify the initial plan looking for safe trajectories, taking into account the new information incorporated about cerebral vascularization (e.g., DSA or CT scans) and adapt that plan for surgical safety [10]. The work presented here starts right after the epileptologist has provided its initial trajectories, and it is aimed at aiding neurosurgeons with tools to early detect and discard high-risk trajectories and find alternative ones when required.
There is consensus that the most important risk to avoid is intracranial bleeding [3], [6], [14] which can cause significant damage or death. Special attention is taken to avoiding vessels at the beginning of the trajectory in the cortex, as they are associated with greater bleeding risk [11]. Surgical planning is patient-dependent and has been described as a time-consuming and inefficient process [15]. It may involve the visual inspection of multiple 2D slices from 3D imaging datasets so that a cylinder shape that surrounds the electrode, commonly referred to as the security zone (SZ) is free of vascular structures to be avoided. During the rest of the manuscript, the SZ is defined by the diameter of the cylinder that has the trajectory as its main axis. The SZ used depends on the precision of the implantation device (e.g., surgical robot) and institution/clinician criteria.
Previous work in stereotactic automatic planning has been presented for Deep Brain Stimulation (DBS) [16], needle biopsy [17], [18], SEEG [19], [20] and other straight-access neurosurgical interventions [21], [22]. For procedures where only the tip is active (e.g., DBS, needle biopsy, or laser ablation) a common approach is to compute the optimal entry point given the desired target defined by the user. However, this strategy is not entirely sufficient for SEEG electrodes, which measure electrical data along its entire length.
Specifically for SEEG planning, the Milan group presented a method that takes into account the distance from vessels, sulci avoidance, and penetration angle based on the user selection of one target point and an entry region [19]. The method was later modified so that the user selects a rough entry and a rough target [23]. This work was further developed to optimize trajectories [11] computing for each trajectory all possible entry-target combinations from within a user-defined entry and target permitted areas, where the exact number of candidate entry and target points was defined by image resolution. In its evaluation study clinicians were asked to rank qualitatively (either 'good', 'acceptable' or 'discarded') the adherence of the alternative trajectory to the initial one, both at the entry and the target regions. In a later work from the same group, a new method was presented with an increased computational speed of 160±102 s/trajectory [15]. To increase safety, a Maximum Intensity Projection (MIP) method was presented which takes into consideration the first centimeter of the trajectory [24].
Epinav is another software platform for automated SEEG planning, reaching interactive rates for finding alternative entry points for a given target point. Then a method is presented which is dependent on the skull mesh segmentation results to initialize the search [20]. More importantly, the method searches for the best entry point for a given target but cannot be used to search for the best target point for a desired entry. It makes fast computations using the GPU (5000 candidate entry points in 250ms). Epinav also has a visualization tool that provides a graph for each trajectory in which the horizontal axis represents the length of the electrode, and the vertical one represents the distance to the nearest critical structure in any direction. Later in [25], a multiple trajectory planning algorithm was presented which, takes into consideration the amount of gray matter traversed to increase electrode efficacy, apart from considering electrode collisions. The reported time was below a minute for 7-12 electrode plans. Regarding alternative path computations, results are displayed with a color-coding directly on the skull surface, similarly to other tools published for DBS [26], needle biopsy [17], [18], and general keyhole neurological interventions [21], [27]- [29].
Zelmann et al. [30] presented another method that also considers the location of individual contacts to maximize recording volume while constraining the trajectories to safe paths, considering a wide variety of constraints. Only 3 target structures per hemisphere were considered (1 in the amygdala and 2 in the hippocampus) and a distance map is computed from the surface inwards to represent the greater importance of measuring from the center of the structure. Gaussian sampling was used to sample possible target points from segmented deep target structures and possible entry points from the entry regions to favor the evaluation of trajectories VOLUME 9, 2021 traversing central areas. On average, 7 minutes per optimized trajectory was reported.
Most of the described approaches focus on adhering to the distal target (e.g., hippocampus). However, the proximal entry point of the brain is sometimes critical in SEEG. This entry point may be crucial to measure from an identified structural lesion (e.g., abnormal brain tissue) or a superficial functional cortical region such as Broca's or Wernicke's [9], [31]. Although adherence to the entry or the target has been used for qualitative alternative evaluation, but to the best of our knowledge it has not been used as an input to the optimization. Furthermore, adherence to the insertion angle is yet another form of adherence that might be useful for SEEG approaches (such as [32]) in which all trajectories are as parallel as possible to each other. Thus, there is a need for an automated trajectory planning tool that takes as input different approaches to maintain adherence to the initial plan, and which performs at interactive rates.

II. METHOD
In this section, the two tools used by our method, the DepthMap and the AlternativeFinder, will be described, as well as their integration with a GUI and its experimental evaluation. A combination of projective geometry and distance map computation is used to transform the 3D problem into the 2D domain to simplify the problem and speed up computation. The tools will be presented to avoid volumetric no-go zones (i.e., vessels) which we will refer to as the mask volume. This mask volume is obtained by registering and merging all DSA datasets into one and then performing a threshold segmentation. Those DSA datasets contain venous and arterial information of each patient (two for each implanted brain hemisphere). Nevertheless, the described functionality can be used with other datasets and can be adapted to consider go zones and/or mesh inputs. The description also assumes that the datasets are contained in a scene graph data structure and are registered via relative transformations [33], [34].

A. DepthMap TOOL
The DepthMap tool accepts two inputs: a mask volume containing no-go zones and a trajectory. The tool is designed to early detect and visualize critical structures within the SZ of the trajectory. It works by projecting the portion of the mask volume contained inside the trajectory SZ into 2D space to create a depth image, which contains depth information for each rendered pixel. Later, that depth image is evaluated to find 'cylindrical SZ vs. no-go zone' collisions. This procedure shares similarities with the MIP implementation presented in [24]. A preliminary version of this tool was presented in [35] and mentioned in [34].
To obtain the depth image, a rendering pipeline needs to be configured. First, a mask volume is positioned in object space with a transformation parametrized by several DICOM tags (e.g., pixel spacing, image orientation). Then, it is multiplied by its model transform T model , obtained by concatenating all its ancestor relative transformations contained in the scenegraph, which again modifies its pose (i.e., translation and orientation) placing it in world space (i.e., registered). The view transform T view is then applied, leaving data in camera space, which is centered at the entry point facing the target point. The next step is to apply an orthographic projection T proj defined by the smallest square prism which fits the cylindrical security zone of the trajectory (Fig. 1). After that perspective division takes place, and finally, the viewport transform T viewport maps the result to the final 2D image (resolution of 128 × 128 pixels).
After the rendering, the depth image (i.e., z-buffer) is retrieved. A non-linear look-up-table (LUT) is used to map each value of the depth image from black (far clipping plane) to a solid color (near clipping plane). The pixels inside the circumference that are not black (i.e., depth different than the far clipping plane) represent vessels inside the SZ. The algorithm takes advantage of the volume rendering capabilities provided by VTK [36].
An approximate formulation of the order in which each step is applied can be represented by the following equation: Radial distance to critical structures is preserved under the orthographic projection, allowing for the computation of the distance from closest vessel to electrode axis by measuring the closest non-null pixel to the 2D depth image center. The span of the trajectory can be subdivided (e.g., entry region and rest), allowing for the computation of vascular distances on each subpart of the trajectory. Please note that in this tool and the rest of the manuscript, distance to vessels is expressed in diameter, instead of radius, to comply with the SZ definition used in some robotic implantation systems.

B. AlternativeFinder TOOL
This tool also takes as input a mask volume containing no-go zones and a trajectory. Its purpose is to find alternative trajectories (i.e., trajectories without vessels inside the SZ) that adhere to a given initial trajectory. The metric to decide what is closer (or more adherent) is selected by the user from among three available options: trajectories with similar entry, similar target, or parallel (with the same angle relative to world coordinates).
In contrast with the DepthMap tool, where each depth image computation allowed for the computation of one trajectory, the Trajectory Finder computes for each render as many trajectories as pixels contained in the 2D depth image (currently 128×128 pixels, resulting in over 16K alternatives/depth image).
First, the Signed Mauer distance map algorithm [37] is used to create a distance map from the binary image. A transfer function is created and configured with a Boolean predicate function, mapping 1 (fully opaque) to values smaller than the SZ radius of the trajectory, and 0 (transparent) to those above it. When visualizing the distance map with that transfer function, a fattened version of the vessel tree appears. The reason to perform this operation is that finding vessels inside the security zone (this time defined by a capsule instead of a cylinder) is equivalent to finding the intersection between line segments and the threshold distance map described (Fig. 2).
The parametrization of the rendering is the same as in the DepthMap tool, apart from the location of the camera, and the use of orthographic (parallel) or perspective (similar entry/target) projection (Fig. 3), which depends on the adherence strategy selected: • Similar entry: a perspective camera located at the entry point facing towards the target point. Field-of-view is set to 15 • .
• Similar target: same as the previous one but located at the target point facing the entry point.
• Parallel: same parameters as in the previously described step. Resolution and lateral frustum bounds could be modified to affect trajectory density and search space. After the rendering, depth is recovered, where each pixel with a value of 0 represents an avascular trajectory with no vessels inside its SZ, (P2 in Fig. 4). Then, the closest 0-value-pixel from the image center (P1 in Fig. 4), if found, is used to compute its associated trajectory. For that, all pose transformations described for the DepthMap tool must be applied to P2 in reverse order.
In the similar entry strategy, when an alternative is not found in the first render, the search continues moving the  (2) were radial distance RD is 0.1 mm and ω is ≈ 135.5 • , an approximation of the golden ratio in the angular form. With: Where: A new set of trajectories are evaluated for each iteration, and the search stops as soon as one avascular trajectory is found, or a maximum number of iterations (n max = 32) is reached. The similar entry points lie on an Archimedean spiral (Fig. 5) and are visited in order (from the inside to the outside), with higher density near the initially selected entry (Fig. 6). The similar target strategy provides equal functionality but moving the target instead. When more than one mask volume or mesh is available (e.g., DSA volumes for arteries or veins considered separately), each dataset is projected individually, and the resulting depth images are aggregated before the closest risk-free trajectory selection. The same approach is used when considering different SZ requirements along the length of the trajectory. This situation is depicted in Fig. 5, where two frustums (an orange frustum for the entry region and a black frustum for the rest of the electrode) are used to search for vessels at different subsections of the trajectory.

C. GUI: INTEGRATION WITH A SURGICAL PLANNER
The method has been integrated into SYLVIUS, a SEEG surgical planning system designed specifically for epilepsy surgery [34]. The DepthMap tool is controlled by the user with the interface depicted in Fig. 7. Here, information obtained from a DSA acquisition is used to provide the vascular structures. For each trajectory, the interface lists the avascular SZ diameter (i.e., twice the radial distance to the axis). The interface also signals collision with other SZ, which is computed by projecting the meshes of every other trajectory (green screen when no collisions, greyscale otherwise).
When the user clicks on any vessel on the DSA image seen along the trajectory of the electrode, the 3D and tri-planar views get centered at that point for user inspection. For that, the transformations described in the previous section must be computed in reverse order, starting from the 2D image coordinates and its corresponding z-value, all the way back to world space. Please note that this requires the depth value, which is normally lost in normal rendering and MIP.
If a trajectory is found to be unsafe with the DepthMap tool, the AlternativeFinder tool can be configured with different avascularity requirements for the entry region (first 6 millimeters) and the rest of the trajectory. It can also be instructed to search for alternatives with a specific adherence strategy. When a trajectory is selected, only close by vessels are depicted in the 3D view to facilitate visual inspection.

D. EXPERIMENTAL DESIGN
To retrospectively evaluate both tools, a set of 12 patient cases requiring SEEG implantation was used. The study was approved by the local ethical committee on clinical  investigations. All patients or their legal representatives provided written informed consent for the use of the data for research and publication. The cases were from Hospital del Mar, Barcelona (part of EpiCARE, the European Reference Network for epilepsy treatment) and contained preliminary and unrevised plans imported from two robotic platforms: the ROSA R (Zimmer Biomed Inc.) and neuromate R (Renishaw plc.). Electrode coordinates were imported referred to T1-weighted double-contrast Gadolinium T12c images (Philips Medical Systems, 256×256×100 image resolution). The cases were then co-registered within SYLVIUS with their DSA datasets (Philips Medical Systems, 256 × 256 × 256 resolution) and stored in a scene-graph data structure [34]. The registered DSA datasets (up to 4 for each case) were fused to produce a single 3D volumetric image. A threshold to segment vessels was then manually selected by a trained neurosurgeon, obtaining the vessel mask volume.
For each initial trajectory, the DepthMap tool was used to locate vessels inside the SZ. When vessels were found, Trajectory Finder was used to find alternatives. DepthMap was used once again to measure the alternative avascular diameter for comparison. The experiment was conducted on a laptop equipped with an Intel R Core i7-4810MQ CPU at 2.80GHz, and an NVIDIA Quadro K3100M GPU. Collision with other electrodes was not considered in the experiment.
The avascularity constraint was set to 7 mm in diameter for the first 6 mm of the trajectory (entry-SZ), and 5 mm diameter along the rest of the electrode. The number of electrodes that were already safe under that criteria, the ones with no alternative found, and the ones with alternative paths along with its type were recorded. To measure the adherence to the initial plan for each re-planning strategy, the Euclidean distance from the original and alternative entry and target points was computed. The re-planning computation time was also recorded for each trajectory and adherence strategy.

III. RESULTS
The alternative trajectories obtained following the described experimental design are shown in Fig. 8, arranged by the re-planning strategy used. For each alternative found, the avascular diameter measured for the initial (black) and the alternative (red, green, and blue) is shown, with an arrow connecting them. Distance to vessels is expressed as an equivalent avascular SZ diameter (twice the radial distance) for compatibility with the SZ convention used in the ROSA R and neuromate R robotic planning software. A histogram with the Euclidean distance between original and alternative entry and target points is displayed for each re-planning strategy.
When an alternative was found, the average computation time was 77 ms per trajectory, with 75% of cases under 85 ms, and a maximum computation time of 765 ms. When no alternative was found, the time ranged between 8 and 889 ms. The average number of trajectories computed per millisecond was 1781 (with a central 50% of 1638-2028). Information regarding the total number of electrodes found in each case, its type, and computation times can be found in Table 1.
To evaluate a different use-case, an additional experiment was conducted aimed at maximizing distance-to-vessels. For each electrode, both the SZ and entry-SZ were measured with DepthMap, and later, the AlternativeFinder was run iteratively with increasing diameter of entry-SZ and SZ requirements (in steps of 1mm, favoring the entry-SZ over the SZ), until a maximum was found. Please note that when the diameter of the SZ is over 10 mm, it is presented as 10 mm in the graph, as that is the maximum diameter that can be measured with the DepthMap in the described configuration. In this additional experiment, every trajectory obtained enhanced avascular SZ values (Fig. 9).

IV. DISCUSSION
SEEG is a complex procedure requiring collaboration from several clinical specialties. The presented work fits in the workflow right after the epileptologist has provided its initial trajectories, and it is aimed at aiding neurosurgeons, with tools to early detect and discard problematic trajectories and find alternative ones when required. The AlternativeFinder tool is designed to adhere to different aspects of an initially planned trajectory. The reason we have not yet addressed the definition of that initial trajectory is that we have not found a practical way to that map seizure semiology (an important criterion for establishing the EZ localization hypothesis) into 3D space due to its non-image nature and complex symptom classification.
An experiment was presented to evaluate the search for the closest avascular trajectory given two fixed expected SZ values (usually the same for the whole plan). In the experiment, the user clicks and waits for the answer. Computation time is within acceptable values for an interactive tool, with response times well below one second. The presented method uses only hard constraints [11], [25], and it can early abort a search when no possible alternatives are found for one of the constraints (e.g., vessels at the entry), yielding response times as low as 8 ms for the parallel-alternatives computation.
As shown in Fig. 8, the similar entry approach produces smaller displacements on the entry than on the target. The opposite happens in the similar-target approach, demonstrating adherence to the initial trajectory in different ways. From all cited research, only De Momi et al. [11] addressed the issue of adherence to the entry and the target points and provides a qualitative evaluation of the results. Our study introduces a quantitative evaluation of the adherence adding yet another adherence option with the parallel strategy, obtaining alternative trajectories with the same insertion angles.
In Fig. 7 the histogram for parallel alternatives is identical for the entry and target point displacement, which was expected. This method may be used for parallel implantations [32], [38], and it offers similarities with the original SEEG carried out in the Sainte Anne Hospital, France [4]. The parallel search described in this work is too strict and does not consider trajectories whose angles deviate slightly that of the initial trajectory, and it will be the subject of future work. Visualization of surgical risk in a cortical mesh is valid for procedures where only the tip of the trajectory is relevant (e.g., DBS, needle biopsy) but it would not be useful for plotting risks associated with a fixed entry and multiple target points. Our visualization shares similarities with the MIP approach presented by the Milan group [24], but offers some extra advantages. The DepthMap tool records depth along the trajectory which is used by the user interface to let the user interactively inspect the tri-planar view of the original image simply by clicking on the desired part of the 3D vessel. This visualization is more natural to the user's eye than MIP, as objects closer to the user's eye occlude distant ones. Lastly, our tool scans the whole length of the trajectory, instead of the first centimeter. Our visualization provides radial positioning of the risks with respect to the trajectory axis, which cannot be visualized in the distance graph widget provided within Epinav [20]. This in turn provides more precise depth information that could be very useful during the implantation. Both visualization widgets are complementary and will be the subject of future work.
A recent publication examines 8 different heuristics used by different tools to optimize SEEG [39], where the first one was the maximizing distance from the vasculature. The second one was to avoid sulci (which we have tested using sulci information provided by Freesurfer [40]) but the main reason to do so is to avoid vessels, and we have opted to use DSA for that constraint. Furthermore, sulci avoidance conflicts with the next listed heuristic which is to maximize gray matter sampling. The drilling angle with the skull is another common heuristic, which our tool handles differently. Instead of providing the system with only a target point, we require an initial trajectory to be specified, and we guarantee that the trajectory does not deviate more than a certain angle.
The presented computations have the advantage of working directly with the original volumes, not requiring an intermediate mesh representation of the patient, although currently still require manual thresholding of the DSA. It can also be used with meshes, either with a distance map computation or directly projecting the mesh, in which case special care needs to be taken to prevent rendering of the inside of the mesh, which is empty and could return false negatives.
Another contribution of our method is that it is fully parametrizable, independently of the characteristics of the input datasets (e.g., cortical mesh resolution). This allows flexible configuration (search angle, number of trajectories to scan) through modification of specific parts of the rendering pipeline, as described in the DepthMap and AlternativeFinder subsections. When considering DSA imaging, finding avascular trajectories becomes more challenging than with a T1weighted MRI scan with gadolinium enhancement, due to the increased number of vessels. DSA is an invasive technique with its own associated risks, and its benefits are under active investigation [41]. As discussed, visualization and assisted planning are highly dependent on vessel segmentation quality [23], which is an area of active research [42].
A pilot study has been conducted with clinicians from Hospital del Mar to evaluate the use of the described methods in a prospective setup, showing that the tool might be useful for early risk assessment. The tool has proven useful in finding avascular trajectories both from preliminary trajectories and for more definitive ones, where the algorithm performs very subtle modifications to the trajectories. The most prevalent cause of false-positive alternatives has been incorrect vessel segmentation. Another undesired effect happened when the proposed alternative moved to a different brain gyrus. The described method has been integrated into the SYLVIUS stereotactic surgical planning platform. A patent has been filed for the Trajectory Finder method that searches for alternative trajectories. Future work will concentrate on exploring better segmentation of the DSA and other less invasive imaging modalities, computing electrode reports on tissue traversed, enhancing visualization and interaction, and providing preset avascular trajectory search within a certain brain region obtained from brain parcellations.

V. CONCLUSION
We have presented two automated tools that allow clinicians to examine and re-plan SEEG trajectories for avoiding risky structures while maximizing adherence with three possible adherence strategies, both performing at interactive speeds. A quantitative experiment has been conducted to measure the different re-planning strategies yield the expected adherence results. A graphical interface designed for trajectory evaluation in the clinical environment was implemented to allow the user to evaluate the tools. The presented tools work directly with volumetric images and can also be used with meshes. These tools provide early alarms and fast suggestions but do not replace the manual inspection of the trajectories, as they rely on the quality of the segmentation. From 1988 to 1999, he was a member and a Research Staff with the Institute of Systems Science, Singapore. There he developed the Dextroscope, a VR interface for the manipulation of 3D medical imaging data. In 2000, he co-founded Volume Interactions Pte Ltd., a Singapore spin-off company setup to commercialize the Dextroscope technology. Volume Interactions applied the Dextroscope technology to radiology and neurosurgery, creating products for pre-operative planning and image-guided surgery. The company was bought by the Bracco Group, in 2002, had operations in Europe, China, and USA. Back in Spain, he co-founded Galgo Medical SL, Barcelona, in 2013, a software company dedicated to providing 3D medical imaging surgical planning solutions. He leads the development of the arrhythmia care and neurosurgery lines. He is the author of more than 20 articles, and more than 25 inventions. His research interests include surgical planning based on medical imaging, and virtual reality. UPF, and has more than 300 publications in peer-reviewed scientific journals and conferences. His research interests include image processing and computer vision for medical image analysis, computer assisted surgery, and computational simulation for predictive medicine.
Prof. González Ballester was a Program Chair of the IEEE International Symposium on Biomedical Imaging 2019. VOLUME 9, 2021