![Fig. 1. Risk maps obtained using weights k1 = 0.1, k2 = 0.9. The computed risk values are normalized to the range of [0, 1]. Blue, white and red regions represent regions of low, medium and high risk, respectively. Black regions depict impermissible insertion areas. (a) Interactive planning by moving the target point within a region. Our high performance approach opens the possibility of interactive intra-operative visualization for computer assisted re-planning. (b) Moving the target point along a selected path. The shown risk maps are computed by placing the current target at the tip of the needle. The horizontal bar shows the normalized risk as a function of the percentage of covered pathway (i.e. 50% means that the needle is currently at half of the total length of the selected safe path). The insertion velocity can be planned according to the risk function along the selected path.](/figures/fig-1-risk-maps-obtained-using-weights-k1-0-1-k2-0-9-the-32cgnmr6.webp)
![Fig. 7. Comparison between voxel-based approach [3] versus our mesh-based approach in terms of computational time. The equivalence of input structures is as follows: R1 is a 32k Verts-65k Tris mesh (for our meshbased approach), or a 288x172x136 grid (for the voxelbased approach); R2 is a 62k Verts-126k Tris mesh, or a 576x172x136 grid; R3 is a 124k Verts-248k Tris mesh, or a 576x460x136 grid; R4 is a 185k Verts-370k Tris mesh, or a 576x768x136 grid; R5 is a 246k Verts-492k Tris mesh, or a 576x768x272 grid.](/figures/fig-7-comparison-between-voxel-based-approach-3-versus-our-2u60a9ma.webp)
IEEE COMPUTER GRAPHICS AND APPLICATIONS, VOL. XX, NO. XX, MARCH 2013 1
GPU-Accelerated Interactive Visualization and
Planning of Neurosurgical Interventions
Mario Rinc
´
on-Nigro, Nikhil V. Navkar, Nikolaos V. Tsekos, and Zhigang Deng, Senior Member, IEEE
F
Abstract—Advances in computational methods and hardware platforms
provide efficient processing of medical imaging data sets for surgical
planning. In the case of neurosurgical interventions that are performed
via a straight access path, planning entails selecting a pathway, from
the scalp surface to the targeted area, that is of minimal risk to the
patient. We propose a GPU-accelerated approach to enable quantitative
estimation of the risk associated with a particular access path at inter-
active rates. It heavily exploits spatially accelerated data structures and
efficient implementation of algorithms on GPUs. We evaluate the com-
putational efficiency and scalability of the proposed approach through
extensive performance comparisons, and show that interactive rates can
be achieved even for high-resolution meshes. Through a user study,
and feedback obtained from domain experts, we identify some of the
potential benefits that our high-speed approach offers for pre-operative
planning and intra-operative replanning of straight access neurosurgical
interventions.
Index Terms—GPU Acceleration, Neurosurgical Interventions, Risk
Maps, and Straight Access
1 INTRODUCTION
Planning of interventional procedures from pre-
operative imaging has led to the development of
new effective interventional paradigms. Among the
procedures that benefit from such pre-operative
planning are neurosurgical interventions such as
deep brain stimulation, biopsies, and shunt and
insertion of external ventricular drains (EVD). Planning,
often, entails processing of multi-contrast and multi-
modal imaging (MRI/CT) to map anatomical and
pathophysiologic features. A common practice is
manual selection and assessment of the suitability of
different entrance positions on the scalp surface of the
patient in order to select an appropriate insertion path.
To process, visualize and manipulate large volumes
of 3D imaging datasets, automated surgical planning
methods have been proposed in recent years [1], [2],
[3], [4], [5], [6]. The current trend is to move as much
of the planning tasks as possible into the operating
• M. Rinc´on, N. Navkar, and Z. Deng are with the Computer Graphics
and Interactive Media Lab and the Department of Computer Science,
University of Houston, Houston, TX 77204-3010.
• N. Tsekos is with the Medical Robotics Lab and the Department of
Computer Science, University of Houston, Houston, TX 77204-3010.
• E-mails: mrinconnigro2@uh.edu, nvnavkar|nteskos|zdeng@cs.uh.edu.
room to avoid errors that can appear because of
the discrepancies between pre-operative images and
intra-operative images. In this context, techniques for
integrating MRI, CT imaging, and surgical procedures
are being actively investigated [7], [8].
Recently, Herghelegiu et al. [9] proposed a multi-level
framework for planning biopsies of deep-seated tumors.
The system is meant to assist planning at all stages
from the selection of the target point within tumors, to
the selection of a safe entry point that minimizes the
risk of hit with vital structures. After selection of the
target point, the neurosurgeon defines a set of regions
on the surface of the skull, and the system displays a
color-coded map that defines the risk associated with
the paths within the region. The system also assists the
neurosurgeon in the validation of the path by displaying
information about regions in the path that require closer
inspection, specifically, the regions in which the path is
closer to any point in the brain vasculature. In terms of
visualization, the system displays a set of 2D slices and
volume rendered 3D views to assist the neurosurgeon in
the planning and validation of the selected pathway.
All the aforementioned works describe approaches to
find the optimal path for neurosurgical interventions
either automatically [4] or by assisting the decision-
making process of the neurosurgeon [5], [6], [9]. State-of-
the-art approaches use criteria defined by the operator
to generate color-coded projections of internal structures
of the brain on the scalp surface. The input includes a
geometrical representation of vital tissue (e.g., meshes
[2] or segmented volumetric data [3], [10]) that should
not be traversed by the path from a point on the scalp
surface through which the operator may insert a surgical
tool to the target point in the region of interest. With
those approaches, the operator needs to select the target
point carefully and precisely since the time required
by the algorithms for computing optimized paths is
generally long (e.g., [2]). Moreover, in most cases the
target (i.e., the vital structure) is not a single point but
rather a three dimensional (3D) structure or surface.
Therefore, accessing all different possible targeted points
within or on the region of interest requires orders of
magnitude longer processing time than the single target
point case [3], [2]. The possibility of performing this
IEEE COMPUTER GRAPHICS AND APPLICATIONS, VOL. XX, NO. XX, MARCH 2013 2
Fig. 1. Risk maps obtained using weights k
1
= 0.1,
k
2
= 0.9. The computed risk values are normalized to
the range of [0, 1]. Blue, white and red regions represent
regions of low, medium and high risk, respectively. Black
regions depict impermissible insertion areas. (a) Interac-
tive planning by moving the target point within a region.
Our high performance approach opens the possibility
of interactive intra-operative visualization for computer
assisted re-planning. (b) Moving the target point along
a selected path. The shown risk maps are computed by
placing the current target at the tip of the needle. The
horizontal bar shows the normalized risk as a function
of the percentage of covered pathway (i.e. 50% means
that the needle is currently at half of the total length of the
selected safe path). The insertion velocity can be planned
according to the risk function along the selected path.
process in real-time or near real-time allows interactive
visualization and planning, which would significantly
reduce the planning time and thus allow effective intra-
operative re-planning of neurosurgical interventions if
needed.
In this article, we propose a GPU-accelerated method
to process, visualize and plan neurosurgical interven-
tions at the interactive or near real-time speed as it
is required for intra-operative re-planning. It consists
of two main components: (i) Embedding the geometrical
structures that represent critical tissue areas, pertinent to
the procedure, into spatial data structures. This speeds up
computation of the geometric queries involved in the
estimation of the risk associated with paths. (ii) Im-
plementing those algorithms onto graphics processing units
(GPUs). This exploits the parallel nature of the problem,
while effectively handling the involved irregular work-
load. Through various comparisons (e.g., between our
approach and the CPU-based baseline implementation,
and between our approach and an existing voxel based
method), we demonstrate that besides its robustness
and interactive speed (e.g., two orders of magnitude
higher than the conventional method), our approach can
generate safer straight access paths for neurosurgical
interventions.
In practical terms, existing surgical planning ap-
proaches (e.g., the multilevel planning framework de-
scribed in [9]) could benefit from a high-performance
approach such as the one we propose, as it will open
the possibility of intra-operative planning. In addition,
the computational power might be beneficial in a pre-
operative scenario by allowing to explore and visualize
much larger planning spaces through fast processing,
and further avoid possible complications in the pro-
cedure. For instance, our high performance approach
would allow to include additional information and plan-
ning constraints (e.g. additional vital structures such
as functional brain regions, and consider the angle of
intersection of the tool with the targeted regions) with
no large increases in computation time.
To the best of our knowledge, this work is the first
to enable interactive estimation and visualization of
risk in straight access neurosurgical interventions with
mesh-based tissue representation. Our current clinical
workflow for using the herein proposed approach entails
interactively modifying (with a pointing device such as
a mouse) the location of the target point and visualizing
the effects on the risk maps. As a benefit of the high
speed of our approach, the operator can interactively
explore a much larger number of possible paths and
target points, thereby enhancing the effectiveness and
efficiency of the decision-making process involved in
neurosurgical interventions.
1.1 Applicability
Fig. 1(a) illustrates an example of how the proposed
approach can assist the operator to select both an optimal
insertion and a target point. In this particular paradigm,
the operator selects a target point within a pre-defined
region (e.g., into the core or onto the surface of a tumor)
and interactively visualizes the resulting risk map (color-
coded) to guide the selection of an optimal insertion
point. After a target point and an insertion point are se-
lected, the operator can also view the risk map associated
with the current position of the advancing needle (Fig.
1(b)). Such dynamic risk maps can significantly facilitate
IEEE COMPUTER GRAPHICS AND APPLICATIONS, VOL. XX, NO. XX, MARCH 2013 3
Fig. 2. Schematic view of the proposed approach. The
operator can interactively visualize the risk associated to
particular paths, and explore the effect of changing the
position of the target point.
decision-making for the operator. When the operator
selects a path, the path is usually further analyzed
and inspected by looking at imaging planes (MR slices)
orthogonal to the insertion path (also known as bird’s-
eye view). For stereotactic robot-assisted neurosurgical
interventions, such as in the NeuroArm and NeuroMate
systems, the velocity of the tooltip during insertion can
be controlled based on the risk function along the path.
The horizontal bar in Fig. 1(b) shows a visualization of
the current risk at each point in the selected pathway.
Each risk value along the pathway is computed by
setting the needle tip as current target. The risk along the
pathway increases non-linearly as the needle advances
and reaches the maximum at the final target position.
If the environment changes are known (when the tool
is near the vital tissue, etc.), the operator can adjust the
velocity of the insertion tool to advance more cautiously.
Indeed, as shown in Fig. 1(b), the high (interactive)
speed of this method allows the operator to see the exact
risk from the surrounding tissue when the tool is being
inserted, by simply placing the target point along the tip
of the tool and recomputing the risk map at each step.
2 METHODOLOGY
2.1 Definition of Planning Risk Maps
In the formulation of our problem, selection of safe
straight access paths entails minimizing a function that
quantifies their associated risks. We define such a func-
tion based on the following criteria: (1) the paths cannot
intersect any critical tissue, otherwise, this may result
in unacceptable patient injury; (2) the paths should be
as distant as possible from the critical tissue; and (3)
the path length (i.e. the distance between the entry
point and the target) should be as short as possible. In
general, additional selection criteria can be seamlessly
incorporated into this function, based on the particular
needs of the procedure.
Inspired by [2], we also represent structures seg-
mented from imaging as triangular meshes. Specifically,
we represent the scalp surface as a triangular mesh
m
s
= (V
S
, T
S
) and the critical tissue as another triangular
mesh m
c
= (V
C
, T
C
). We consider the set of candidate
paths as the line segments that extend from the vertices
of m
s
to a pre-defined target point p. To obtain the set
of permissible paths, we discard all the unsafe paths
that intersect with any triangle in m
c
. Then, the risk
associated with a particular path (belonging to the set
of the permissible paths) starting at a vertex v ∈ V
S
, is
computed as:
Risk(v, p) = k
1
||v − p|| − k
2
min
v
0
∈V
C
d(vp, v
0
) (1)
d(vp, v
0
) =
||v
0
− p|| if (v − p) · (v
0
− p) < 0
||v
0
− v|| if (p − v) · (v
0
− v) < 0
||(v −p)×(p−v
0
)||
||v −p||
otherwise
(2)
Here k
1
≥ 0, k
2
≥ 0 are user-specified weighting
parameters (see Fig. 3 for a visualization of the effect
of varying these parameters). All the risk map visualiza-
tions in this writing were created using weights k
1
= 0.1
and k
2
= 0.9, unless otherwise noted. The first term in
Eq. 1 accounts for the total length of the path, and the
second term accounts for the proximity of the path to
the critical tissue. As such, larger values of k
1
enforce
the selection of short paths, while larger values of k
2
enforce the selection of paths that are farther away from
critical structures. Brute-force computation of the set of
permissible paths, as well as the risks associated with
all the paths, turns out to be computationally-expensive
even for low-resolution meshes. In order to efficiently
perform the computation, we build acceleration spatial
data structures from the geometric primitives of the criti-
cal tissue mesh. This significantly reduces the amount of
computation required to obtain the risk map. An efficient
GPU implementation of our approach further enables the
computation of risk maps at an interactive rate even for
high-resolution meshes. Fig. 2 shows a schematic view
of our proposed approach. The operator can interactively
visualize the risks associated with particular paths, and
explore the effect of changing the position of the target
point.
Calculation of the risk maps involves two computa-
tionally expensive steps. The first is the identification
of all the paths that must be discarded, because they
make direct impact on the critical tissue. In our specific
formulation, it involves testing for intersection between
a line segment and all the triangles in m
c
. The second
step is the computation of the second term in the Eq. 1
above. This computation involves determining, among
all the vertices in m
c
, which one is the closest to a
specific line segment. Testing each path for intersection
against each triangle, and computing its closest distance
to each vertex are not affordable options for interactive
visualization and planning.
2.2 Acceleration Spatial Data Structures
To avoid a large bulk of the computation required to
obtain the risk map, without sacrificing precision, we
IEEE COMPUTER GRAPHICS AND APPLICATIONS, VOL. XX, NO. XX, MARCH 2013 4
Fig. 3. Visualizations of risk maps for different pairs of weights. (a) Position of the target point. (b) Enforcing a large
proximity: k
1
= 0.1 and k
2
= 0.9. (c) Enforcing the balance between length and proximity: k
1
= 0.5 and k
2
= 0.5. (d)
Enforcing the selection of short paths: k
1
= 0.9 and k
2
= 0.1.
Fig. 4. Bounding Volume Hierarchy. (a) is an example of the BVH of triangles, (b) is an example of the BVH of vertices,
and (c) shows bounding boxes for a 3D mesh (down to a depth of 3).
use acceleration spatial data structures. Specifically, we
opt for the use of Bounding Volume Hierarchies (BVHs)
[11], since they are inexpensive to compute, and they can
be quickly restructured to support meshes that may go
through dynamic deformations. It should be noted that
in this work we only consider static meshes, not general
dynamic meshes. A BVH is a partition of geometric
objects into a hierarchy that can be used to accelerate
geometric queries as the ones involved in our problem.
Formally, BVHs are binary trees, where each node rep-
resents a group of objects and the region in space that
they occupy, and children nodes are recursive partitions
of the group of objects within their parent nodes. Fig.
4(a)-(b) show 2D examples of BVHs, and Fig. 4(c) shows
bounding boxes for a 3D mesh (down to a depth of 3).
We construct two BVHs from the critical tissue mesh
m
c
: one holding its triangles and the other holding its
vertices. The BVH of triangles is used to compute the set
of permissible paths, and the BVH of vertices is used for
the computation of the proximity term (i.e., the second
term in Eq. 1).
A geometric query over the BVHs involves a traversal
of the data structure. To compute the set of permissible
paths we need to figure out whether each path intersects
a triangle in m
c
. Starting from the root, at each step of
the traversal, an intersection test against the children
nodes is performed. If a node is intersected, then the
traversal must be continued in a recursive manner over
the substructure for which the node is root. Otherwise,
the whole substructure can be safely ignored, because it
is guaranteed that none of the triangles that it holds can
be intersected by the line segment. Triangles are actually
contained within the leaf nodes of the tree structure, so
when a leaf node is reached, each of its triangles must
be tested for intersection against the line segment. The
traversal stops when a triangle hit by the segment is
found or when all the triangles within the intersected
nodes have been tested. Fig. 4(a) illustrates this idea for
the BVH of triangles. The BVH traversal starts with a test
of the box labeled AABB. As it is intersected, the traver-
sal continues with tests of its two children (i.e., AABB
1
and AABB
2
). AABB
1
is intersected, so both its children
must be tested as well. AABB
1,1
is not intersected by
the path (going further down that node can be safely
discarded), but AABB
1,2
is intersected. Since AABB
1,2
is
a leaf node, all the geometric primitives it contains must
be tested for intersection. In the case of AABB
2
, we can
see that only its child labeled AABB
2,1
must be tested
for intersection.
The computation of the proximity term (the second
IEEE COMPUTER GRAPHICS AND APPLICATIONS, VOL. XX, NO. XX, MARCH 2013 5
term in Eq. 1) takes place in a similar fashion over
the BVH of vertices. In this case, for each path we are
interested in determining what is the closest vertex (in
m
c
) to the line segment. Fig. 4(b) illustrates a query over
the BVH of vertices. For each path we determine which
of the children of a node is closer, and which one is
farther. In general the closest child is more likely to
contain the closest vertex than the other nodes. If the
closest vertex found within that node is closer to the
line segment than any point in the bounding box of the
farthest child node, then the farthest node can be safely
discarded. In our example, node AABB
2
is the closest
one, therefore the traversal is performed in a recursive
manner over this node. Within AABB
2
, child AABB
2,1
is
the closest, and it is also a leaf node (containing vertex
primitives), so the closest vertex within this node to our
path is found. Since the closest vertex is closer to our
path than the closest point in bounding box AABB
2,2
,
we can safely discard traversals through such node. In a
similar way we can discard traversals through AABB
1
.
2.3 GPU-based Parallel Implementation
In our proposed approach, the algorithm generates the
risk maps by performing the calculations discussed in
Section 2.2 for a large number of potential paths. Since
the calculations for each path are independent from
those of other paths, the computational task will be
ideal for implementation on massively parallel proces-
sors such as GPUs. However, in spite of the intrinsic
parallelizable nature of the problem, achieving a high-
performance on GPUs is not straightforward. The origin
of this is that the computations required for each path
are highly variable, and this irregular workload leads to
underutilization of GPU resources. To describe the pre-
sented GPU implementation, it is instructive to briefly
review its operational structure. A GPU consists of a
number of processing units, i.e. the streaming multipro-
cessors (SMs). In the streaming programming paradigm,
grids of threads are subdivided into groups known as
blocks. Different blocks are assigned for execution to
different SM in a static manner. i.e. after a block has been
assigned to an SM, it cannot be executed on another
SM. After being assigned, the threads within each block
are further subdivided into groups of 32 consecutive
threads, the warps. Threads within a warp are executed
in parallel within the SM in Single Instruction Multiple
Thread (SIMT) mode. In SIMT mode, each thread in a
warp executes the same instruction at the same time.
With respect to the problem of computing the risk
map, a naive scheme of assigning work to processing
units in CUDA or OpenCL would correspond to launch
the execution of a grid with as many threads as the
paths. Therefore, the computation for each thread would
be performed in a single thread. As a result, the high
variability of the required computations for different
paths will lead to resource underutilization. Practically,
the workload assigned to some blocks will need more
time than other blocks, eventually leading to a situation
where some of the SMs are idle for a portion of the
overall computation time. The difference between the
processing time required for different tasks can thus
lead to a high workload imbalance that translates into
performance degradation [12].
To obtain better performance, we balance the work-
load within the GPU using a centralized queue [12]. All
the paths to be computed are kept on a single queue
implemented as an array in the global memory. The head
of the queue is also kept in the global memory, and it is
visible to every thread in the grid. Instead of assigning
the computation of a path to a single thread, we make
long-running warps to retrieve batches of paths (tasks)
from the global queue. In our work, long-running warps
are those that only finish their execution when all the
paths in the global queue have been processed. This can
be appreciated when the warps process is considered:
when they are first launched, the first thread takes the
next task by updating the head of the queue through an
atomic addition operation (it needs to be atomic in order
to avoid race conditions). Then, each of the threads,
within the warp, processes one of the paths in the task.
When a warp completes this batch, it takes another task
of paths, or finish execution in case the queue is empty.
The processing of each path, first involves traversing
the BVH of triangles for determining whether it hits the
critical tissue, and then traversing the BVH of vertices to
compute the proximity term of the risk function. With
the centralized queue, none of the SMs is going to be
idle, as long as there are remaining tasks to be processed.
3 EVALUATION AND RESULT ANALYSIS
3.1 Computational Performance
To evaluate the efficiency of the proposed approach (i.e.,
to what degree it can be used in an interactive fash-
ion), we measured the time required for performing the
computation of the risk maps for meshes with different
resolutions. Specifically, the mesh sizes were varied as
follows: (i) for the scalp surface, the range of triangles
was 2k to 276k and that of vertices was 4k to 553k;
and (ii) for the critical tissue, the range of triangles was
2k to 276k and that of vertices was 46k to 300k. It is
noteworthy that the number of vertices in the scalp
surface mesh determines the number of paths to be
processed.
We implemented our proposed approach on a GPU
using CUDA, and timing was recorded using CUDA
time events. With this, the recorded time includes both
the execution time of the CUDA kernels, and the com-
munication overhead due to data transfer between host
and device. Specifically, it starts when the BVH data (as
well as the scalp surface mesh data) are transferred to the
GPU, and ends when the results from the computations
of the path length and proximity terms are transferred to
the main memory of the host. Our implementation was
tested on a mainstream commodity desktop computer