Original Article

Deformation invariant bounding spheres for dynamic active constraints in surgery

Proc IMechE Part H: J Engineering in Medicine 2014, Vol. 228(4) 350–361 Ó IMechE 2014 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/0954411914527440 pih.sagepub.com

Stuart A Bowyer and Ferdinando Rodriguez y Baena

Abstract Active constraints are collaborative robot control strategies, which can be used to guide a surgeon or protect delicate tissue structures during robot-assisted surgery. Tissue structures of interest often move and deform throughout a surgical intervention, and therefore, dynamic active constraints, which adapt and conform to these changes, are required. A fundamental element of an active constraint controller is the computation of the geometric relationship between the constraint geometry and the surgical instrument. For a static active constraint, there are a variety of computationally efficient methods for computing this relative configuration; however, for a dynamic active constraint, it becomes significantly more challenging. Deformation invariant bounding spheres are a novel bounding volume formulation, which can be used within a hierarchy to allow efficient proximity queries within dynamic active constraints. These bounding spheres are constructed in such a way that as the surface deforms, they do not require time-consuming rebuilds or updates, rather they are implicitly updated and continue to represent the underlying geometry as it changes. Experimental results show that performing proximity queries with deformation invariant bounding sphere hierarchies is faster than common methods from the literature when the deformation rate is within the range expected from conventional imaging systems.

Keywords Virtual fixtures, haptic rendering, proximity querying, space-partitioning trees, medical robotics

Date received: 30 October 2013; accepted: 7 February 2014

Introduction ‘Active constraints’ (commonly known as ‘virtual fixtures’) are control strategies used to provide assistance in human–robot collaborative manipulation tasks.1 Although they have been used in other areas, the primary focus of active constraint research is for robotassisted surgical procedures. In surgery, ‘regional’ active constraints can prevent a surgeon from accidentally penetrating and damaging delicate tissue regions,2 and ‘guidance’ active constraints can encourage a surgeon to accurately follow a plan defined preoperatively.3 The overwhelming majority of active constraint methods discussed within the literature assume that the environment, and hence the constraint geometry, is static throughout the surgical intervention. For procedures such as orthopaedic surgery, this is a valid assumption; however, for many others, it is not, and the tissue will move and deform throughout. Tissue movement and deformation can be due to physiological motion, such as respiration, tool–tissue interactions and other indirect effects of surgery. In these cases, ‘dynamic active constraints’ are required, such that the changes in the

patient’s anatomy are reflected in the geometries used to enforce the constraints.4–7 To apply active constraints, the constrained anatomy must be represented as a discrete computational object. Simple structures such as planes8 and threedimensional (3D) geometric primitives (i.e. cones, spheres and cylinders)9 can be used for this purpose; however, the coarseness of their approximation may make simple structures unsuitable when the required surgical margins are narrow and the surgeon is expecting to work in close proximity to the tissue structures which are protected by the constraint. When tight and accurate conformance to the patient’s anatomy is required, polygonal meshes are desirable for defining constraint geometries.2

Department of Mechanical Engineering, Imperial College London, London, UK Corresponding author: Ferdinando Rodriguez y Baena, Department of Mechanical Engineering, Imperial College London, London SW7 2AZ, UK. Email: [email protected]

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

Bowyer and Rodriguez y Baena

351

In order to compute robot commands which will enforce a constraint, the relative configuration between the constraint geometry and the surgical tool, typically the pair of closest points, is first required. When meshes are used to represent the constraint, closed-form, linear algebraic expressions cannot efficiently be used to establish this relationship, and the problem becomes one of proximity querying or collision detection. Li et al.,2 use a variant of the ‘oriented bounding box’ (OBB), bounding volume hierarchy to establish the relative tool-constraint configuration; however, a range of similar methods from within the computer graphics literature could also be used. A limitation of conventional bounding volume hierarchies is that they are rigidly constructed around the underlying mesh such that subsequent deformation of the mesh will likely lead to the invalidation of the hierarchy and incorrect results. In a dynamic active constraint, it can be expected that the mesh will undergo deformation, and therefore, this problem must be overcome. Some methods for rebuilding or refitting bounding volume hierarchies after deformations are proposed within the literature (as surveyed by Teschner et al.10); however, each rebuild or refit in these methods is computationally expensive. If the deformation rate (i.e. the discrete rate at which the surface geometry changes) approaches that of the traversal rate (i.e. the rate at which the tree is queried for a collision or proximity), then they become inefficient. To ensure accuracy in constraint placement, smoothness of surgeon interaction and overall system stability, the constraint geometry should update to track the tissue at as high a rate as is possible. The conventional deformable bounding volume hierarchies are subsequently unsuitable, and an alternative is introduced here. The proposed method is based on a conventional bounding sphere hierarchy; however, it is assembled in such a way that it is invariant to deformation in the underlying mesh. These ‘deformation invariant bounding sphere (DIBS) hierarchies do not require rebuilding or refitting, so that, even though some traversal efficiency is lost, when the deformation rate is high, the overall query time is up to an order of magnitude less than conventional methods. In the following section, relevant previous work is introduced. Subsequently, the surface ‘pinning’ of DIBS is discussed, as is the application of the method to surfaces which stretch during deformation. The methods proposed for constructing and traversing DIBS hierarchies are then detailed in the subsequent sections. The results of two comparative experiments with methods taken from the literature are next given along with a discussion of these results. Conclusions are then drawn in the final section.

Related work Dynamic active constraints Research to date into dynamic active constraints has focused on assisting in beating heart surgery. In Ren

et al.,11 a spatio-temporal distance field is produced from preoperative magnetic resonance images. This distance field is encoded into an artificial neural network, which can be used for proximity queries during constraint enforcement. Although this method is dynamic, it is based entirely on preoperative information and cannot adapt to intraoperative variation. In Navkar et al.,12 an adaptive dynamic active constraint is produced. In this work, four features within the ventricles are identified from intraoperative medical images, and a guidance pathway is constructed between them. Although this approach allows continuous updates, the reliance on relatively few control points makes the constraint inherently coarse. Dynamic active constraints are also proposed by Ryde´n and Chizeck,13 where they prevent contact with the surface of a beating heart. In this method, spherical constraints are attached to each point within a point cloud acquired in real time using a commercial range camera. By preventing a haptic proxy from penetrating into these spherical constraints, they are able to compute dynamic constraint forces. In Kwok et al.,4 dynamic active constraints for the body of a snake-like robot are presented. In this work, tubular constraint pathways are defined within the patient’s anatomy, through which the robot can move. For force computation, the position of the robot relative to a tubular constraint is computed using a graphics processing unit (GPU)-accelerated exhaustive search for the maximum separation between the robot body and the constraint.

Bounding volume hierarchies Bounding volume hierarchies are a well-established method for reducing the number of geometric tests required when assessing the spatial relationship between complex geometric objects such as in collision detection and proximity queries.14,15 The simple approach to computing these relationships involves exhaustively checking every possible pair of polygons from the object’s meshes. As complex geometries often contain many thousands of polygons, the inefficiency of this approach tends to be prohibitive. By constructing a bounding volume hierarchy for a complex geometric object, many pairs of polygons can be excluded from searches, causing significant speed improvements for traversal. Each bounding volume in the hierarchy encloses a set of the object’s polygons so that a single intersection/separation test can be performed against the bounding volume, potentially eliminating the need to check any of the polygons within. Various bounding volumes have been proposed for use in this way, and there is often a trade-off required between the properties of each. Bounding spheres16,17 are the simplest to represent computationally, and they are very efficient to examine during traversal; however, they do not generally fit the polygons which they contain particularly tightly and are therefore considered less efficient. Axis-aligned bounding boxes (AABBs)18

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

352

Proc IMechE Part H: J Engineering in Medicine 228(4)

are rectangular cuboids, where each edge is parallel to the base coordinate frame axes. AABBs tend to fit the underlying geometry tighter; however, each individual volume is less efficient to query. Oriented bounding boxes (OBBs)19 are similar to AABBs except they do not require that cuboid edges are parallel to the base axes. This freedom allows for even tighter bounding but individually they are slower to query. Discrete oriented polytopes (k-DOP)20,21 are polytopes with a different number of faces allowing them to be more flexible still. A common feature of all of these bounding volumes is that they are rigidly constructed around a model, and therefore, if the underlying model deforms, they require updating or reconstruction to ensure that they continue to represent the entire model accurately. A comprehensive survey of methods used for implementing bounding volume hierarchies for deformable models is given in Teschner et al.10 Van den Bergen,18 described an efficient AABB hierarchy update method in which the leaves are individually refitted and then a bottom-up algorithm refits each parent based on the boxes of their children. A bottomup method for bounding spheres is described by Brown et al.,22 where the locality of some deformations is exploited to minimise node updates. It was stated by Larsson and Akenine-Mo¨ller23 that if only a few of the deep nodes are reached during traversal, then the exhaustive update approach of Van den Bergen18 can be inefficient. As an alternative, they proposed a ‘hybrid’ approach whereby the top of the tree is updated bottom-up before traversal, and the bottom of the tree is updated when the information is requested during traversal. The majority of other approaches used for applying bounding volume hierarchies to deforming objects involve some kind of restriction upon the form of the deformation which is permitted. Bounded deformation (BD) trees24 can be efficiently updated when the object’s deformation can be described as an aggregation of displacement fields. In this case, the displacement fields themselves are used to update the bounding volumes rather than explicitly considering the surface points. Similarly, Spillmann et al.25 showed how linear and quadratic geometric deformations can be applied directly to update bounding sphere hierarchies. Unlike the methods described above, a DIBS hierarchy is a bounding volume hierarchy for deformable geometries, which does not require an update or reconstruction when deformation occurs. By simply updating the vertices within a deformed mesh, the DIBS hierarchy is implicitly updated and continues to allow accurate proximity queries for use in active constraints.

Surface ‘pinned’ bounding volumes A key aspect of the proposed method for constructing a deformation invariant bounding volume hierarchy is the way that the individual bounding volumes are

‘pinned’ onto the surface of the deforming object. This means that the bounding volumes follow the surface of the model as it deforms, as though part of the bounding volume was physically pinned onto the surface. The idea of having a topologically fixed hierarchy, where each bounding volume moves over time so that it tracks a certain portion of the model’s surface, is not a novel concept. This is the aim of many other deformable techniques such as that of Larsson and AkenineMo¨ller;23 however, these methods rely on manual repositioning and resizing of each bounding volume. Using ‘pinned’ bounding volumes, the repositioning is carried out automatically as the mesh changes, without any special consideration. Pinned bounding volumes can be implemented when the surface of the deforming model is represented using a set of Cartesian points such as the vertices of a polygonal mesh. Each one of these vertices provides a possible pinning location for one or more of the bounding volumes within the hierarchy. To pin a bounding volume to a specific vertex, the bounding volume is linked to an entry in a dynamic vertex table, rather than storing an explicit Cartesian coordinate, for its position. During traversal of the hierarchy, when the position of the bounding volume is required, the vertex is simply retrieved from the table and its current position is used. There are several limitations on the geometries and properties of bounding volumes, which can be pinned to a deforming surface in this way. For the purposes of clarity, only deformations where the external surface area of a model is constant will be considered here, and cases where this simplification does not hold will be discussed in the following section. An important consideration is that the surfaces of a model are only guaranteed to be flat at the lowest level of a hierarchy, where they contain a single polygon. The surface sections bounded by higher branch nodes can be curved, and therefore, deformation can cause flattening and, as a result, gaps in the surface coverage when pinning is used. An example of surface deformation causing this type of failure is shown in Figure 1(a)–(c). A possible solution to this problem is to use oversized spheres, based on the surface topology of the contained surface, for the bounding volumes. The radius of a conventional bounding sphere is the maximum Euclidean separation between the centre and any of the vertices within the mesh section. Instead of this, the maximum geodesic length between the centre and each vertex is used for the sphere’s radius (see Figure 1(d)). The geodesic length is the shortest distance between two surface points while following the surface profile. For a vertex v in the set of surface polygons P, the geodesic bounding sphere radius rg is a function of the sphere centre vertex vc and a geodesic separation function dg (  ) as follows rg = max dg (vc , v) v2P

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

ð1Þ

Bowyer and Rodriguez y Baena

353

Figure 1. An example of how surface deformation can cause tightly fitting pinned bounding volumes to fail (a–c) and how a conservative bounding sphere can overcome this (d). (a) Initial, undeformed mesh surface section with edges a, b, c and d. (b) Fitting and pinning of an AABB bounding volume to a mesh vertex. (c) Failure in coverage after surface deformation. (d) Bounding sphere with a geodesic radius computed from the maximum undeformed edge lengths in (a) rg = maxðkak + kbk, kck + kdkÞ. AABB: axis-aligned bounding boxes.

In a practical implementation of this method, a shortest path tree can be constructed using the edges and vertices of a mesh to approximate the geodesic separations to each vertex. A crucial property of spheres is that they are rotationally invariant with respect to the centre, and therefore, if the geodesic separation radius is used, then it will be guaranteed to cover the surface under all configurations. As this method uses spheres, and as the sphere centres have to be located on the surface of a model, each bounding volume will not fit the surface section as tightly as using other methods. It was shown by Gottschalk et al.19 that tightness of bounding was important and that bounding spheres do not perform as well as tighter fitting bounding volumes such as OBBs. However, this method is not solely concerned with tightly bounding the current geometry. Instead, it aims to bound a volume in space, which the model’s surface section will always remain inside as it deforms. An additional advantage of using bounding spheres is that they are considerably more efficient to query during hierarchy traversal than others as they are primarily considered as a single Cartesian point.

Assumption of limited surface stretch In the previous section, pinned bounding spheres were described where the surface area of a model remained constant. There may be some practical applications where this holds, but for the majority of cases, the deformation will involve some stretching of the model’s surface. To allow the extension of pinned bounding volumes to these cases, an assumption is made regarding the amount of stretch a surface will sustain, and this expected limit is used to inflate the pinned bounding spheres. It is these pinned and inflated bounding spheres that are referred to as DIBS. The limit on a model’s surface deformation might appear to be restrictive; however, soft tissue structures typically only stretch by a relatively small amount before being damaged (Fung26 gives the approximate distensibility of certain soft tissue structures as follows: mesentery 100–200%, ureter 60%, heart muscle 15%, arteries and veins 60%, skin 40% and tendons 2–5%. Additionally, Deck and Willinger27 state that 40% first principal strain in brain tissue causes a 50% probability of severe diffuse axonal injury). Therefore, DIBS can

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

354

Proc IMechE Part H: J Engineering in Medicine 228(4)

be effectively used in a range of clinical applications. In this work, the limit on a model’s surface deformation is described by a ‘limiting surface stretch factor’ s. The value of s states that after deformation, no pair of points on the surface can be further apart than s times their initial geodesic separation. Formally, for a deformation D, a surface is said to obey a limiting surface stretch if and only if sdg (p, q)5dg (D(p), D(q))

8p8q 2 P

ð2Þ

where p and q are vertices on the surface of the model and D(p) and D(q) are corresponding vertices on the surface of the deformed model. Once a suitable value for s has been established for a model, the final radius of the bounding sphere rs can be computed based on equation (1) rs = s max dg (vc , v) v2P

ð3Þ

If a conventional, static, bounding volume hierarchy was simply inflated in this way, it would not achieve large deformation invariance because each individual bounding volume’s position would still be fixed as it was set during initial construction. Oversized bounding volumes have been used to allow for small amounts of deformation between slow hierarchy refits,28 but these methods rely on temporary limits to the absolute deformation of a model such that each point in the model only moves by a limited amount from its initial predeformation position. This is a much more restricting concept than using them in combination with surface pinning, and many cases would require extremely large s values, making traversal inefficient. To summarise the previous two sections, the specific properties of DIBS are as follows:  



The centre of each bounding sphere is centred on the surface of the model at a pinning point. The unstretched radius of each bounding sphere is computed from the geodesic separation of the contained vertices using equation (1). The stretched radius of each bounding sphere is computed from the limiting surface stretch factor using equation (3).

The following section describes how a DIBS hierarchy with these properties can be constructed.

Hierarchical construction A top-down method was selected for the construction of DIBS hierarchies as this has previously been used for bounding spheres by O’Sullivan and Dingliana.17 In this top-down approach, the entire model is initially used as the root node. The model’s surface is then recursively subdivided into smaller and smaller surface sections, with new descendent nodes created for each, until single polygons remain at the leaf nodes. The choices of top-down construction and a binary

structure were only based upon trends within the literature, and further analysis is required to validate them. A bounding sphere is fitted to every node in the hierarchy during construction. There are two parameters required for each bounding sphere, the pinning vertex vc (i.e. a vertex on the surface section to which the sphere will be pinned) and a stretched geodesic radius rs . The computation of these parameters is described below for branch and leaf nodes.

Fitting bounding spheres to branch nodes To ascertain the pinning vertex (i.e. centre) for a branch node’s bounding sphere, the Cartesian mid-range u of the surface polygons P is first computed in the x, y and z Cartesian axes. The centre vertex is then taken as the vertex in P, which is closest to u. Numerically, vc is computed from (maxX(P) + minX(P)) 2 (maxY(P) + minY(P)) uy = 2 (maxZ(P) + minZ(P)) uz = 2 vc = arg min k v  uk

ux =

ð4Þ ð5Þ ð6Þ ð7Þ

v2P

where maxX and minX are functions returning the maximum and minimum vertex x values from the supplied polygon set. This is similarly defined for y and z. Alternative definitions for u, such as the mean and median vertex positions, were investigated; however, the mid-range was found to allow for the fastest traversal. To compute the stretched radius of the bounding sphere, the geodesic radius is required. The geodesic radius of a surface section is approximately equivalent to the distal node separation in a shortest path tree formed from its vertices and edges. Dijkstra’s algorithm29 is used to compute the shortest path tree for each surface section. The algorithm is initialised so that the vertices and edges in the surface section are used as vertices and edges in the graph, edge costs are the Euclidean separations between vertices and the centre vertex vc is used as the starting point. Once Dijkstra’s algorithm has returned the distal node separation from the polygon set, this can be used with equation (3) to produce the bounding sphere’s radius rs .

Fitting bounding spheres to leaf nodes Due to the comparatively high computational expense of polygon operations, it was decided that the leaf bounding spheres (i.e. those which contain a single polygon and have no child nodes) should be tight fitting, to prevent as many unnecessary polygon operations as possible. Tightened leaf bounding spheres were achieved using so-called elastic bounding spheres.

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

Bowyer and Rodriguez y Baena

355

Figure 2. Example of orienting a dividing plane so that it is orthogonal to the eigenvector associated with the largest eigenvalue of surface vertices’ covariance matrix.

The pinning vertex of an elastic bounding sphere is found in the same way as for branch nodes; however, no radius is stored during construction. Instead, if a leaf node is reached during hierarchy traversal, the algorithm computes the maximum Cartesian separation between the pinning vertex and the other vertices in the polygon and uses this value (without any stretching) as the radius of the bounding sphere. Establishing the maximum vertex separation for complex polygons during traversal could be costly; however, this is compensated for by the reduction in number of polygon operations. This elastic method of refitting leaf bounding spheres when required is similar to top-down hierarchy updates discussed in Larsson and Akenine-Mo¨ller.23

Binary subdivision of surface sections If the node is not a leaf, then once a bounding sphere has been fitted, the recursive construction algorithm subdivides the surface section into two parts for the creation of child nodes. The subdivision is performed using a single, oriented, dividing plane. Taking the concept from that used in OBB,19 the dividing plane is oriented so that it is orthogonal to the eigenvector associated with the largest eigenvalue of the surface vertices’ covariance matrix (Figure 2). The dividing plane is positioned along the maximum eigenvector so that it is exactly midway between the maximum and minimum vertex positions when projected onto the maximum eigenvector. The implementation of this division process iterates through the list of surface polygons and allocates each to one of the two child lists based on their scalar projection onto the maximum eigenvector. Formally, the polygon classification is performed as follows   (~ pmax + p~min ) Pchilda = p 2 Pparent : p~5 ð8Þ 2   (~ pmax + p~min ) Pchildb = p 2 Pparent : p~ \ ð9Þ 2

where p is a polygon on the parent surface section Pparent , p~ is the scalar projection of the centre of polygon p onto the maximum eigenvector and p~max and p~min are the maximum and minimum scalar projections, respectively. A result of using geodesic surface separations is that the surfaces considered must be contiguous. The geodesic separation between two disconnected points in a non-contiguous surface is undefined, and therefore, the above description for fitting bounding spheres would fail. To overcome this problem, a check of contiguity after each subdivision is performed, and orphaned sections are swapped between siblings. Additionally, the necessity of geodesic surface separations means that all raw meshes to be searched with a DIBS hierarchy must be fully contiguous to start with.

Real-time proximity queries To compute active constraint forces, the proximity between the tool and the constraint geometry, and the corresponding pair of closest points, must be computed from the DIBS hierarchy. The method used, Algorithm 1, is a recursive depth-first search and is similar to that for conventional bounding sphere hierarchies as proposed by Quinlan.16 The initial parameters for this algorithm are the DIBS hierarchy D and the current tool geometry T, and it returns the minimum separation sep and the closest point pair hcpD , cpT i. It should be noted that the framework presented here is specifically for cases where the representation of the tool (potentially as a proxy) and the constraint surface should not collide, and their relative configuration can be described based solely on minimum proximities between them. If the tool and mesh do collide, the algorithm will identify that the proximity is zero and it will return one of the, potentially infinite, locations at which their surfaces intersect. If full contact information is required, then the DIBS hierarchy can also be used for collision detection, similarly to the method described by Van den Bergen;18 however, this is not the focus of the research presented here.

Experiments and results To validate the effectiveness of using DIBS hierarchies for proximity queries, two experiments were carried out. Standard refittable AABB hierarchies,18 hybrid refittable AABB hierarchies23 and exhaustive searches were implemented for comparison, as these are widely used and do not rely on hardware acceleration. The first experiment was entirely simulated and involved randomly deforming a surface mesh while randomly positioning a tool geometry relative to it and then performing a proximity query between them. In the second experiment, the proximity between a timedependent beating heart model and a human-controlled tool was continuously computed.

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

356

Proc IMechE Part H: J Engineering in Medicine 228(4)

Algorithm 1 Recursive proximity query function procedure DIBS-Proximity-Query D, T m ‘ n root(D) return Recursive-PQ(m, n, T) end procedure procedure Recursive-PQ(m, n, T) sep ‘ Null hcpD , cpT i if leaf(n) then poly getAssociatedPolygon(n) findClosestPts(poly, T) hcpD , cpT i sep separation(cpD , cpT ) else i 0 for all c 2 children(n) do pp getPinningPt(c) rad getRadius(c) dibsseps(i) separation(pp, T)  rad i i+1 end for dibsseps sortAscending(dibsseps) while front(dibsseps) \ m do c associatedChild(front(dibsseps)) hcs, c1, c2i Recursive-PQ (m, c, T) if cs \ m then m cs hcs, c1, c2i hsep, cpD , cpT i end if dibsseps popFront(dibsseps) end while end if return hsep, cpD , cpT i end procedure These experiments were performed using single 3D geometric primitives, namely, spheres and capsules (cylinders with rounded ends), for the tool geometries and manifold, triangular, surface meshes for the deforming objects. The proximity queries were implemented in a single C++ thread and were run on a 3.0 GHz Core i7 CPU (Intel Corporation). An abstract base class was used to implement the hierarchy traversal algorithms. This class was inherited by each of the different bounding volume hierarchies, thus allowing for an unbiased comparison of their effectiveness.

Randomised positioning and deformation experiment For this experiment, four surface meshes of a human brain were used to represent constraint geometries. The four meshes were constructed from the same initial surface, but they were each sampled at different resolutions to create mesh sizes across four orders of magnitude. These meshes ranged in size from 500 to 500,000 polygons and are shown in Figure 3. For each proximity query, the mesh in use was modified so as to induce both local and global random surface

8 Initialise minimum separation 8 Initialise search node structure

8 Initialise separation 8 Initialise closest points 8 Leaf node

8 Branch node 8 Compute child DIBS separations

8 Recursively search ordered children

deformation. The deformation was produced using a combination of random, independent, sinusoidal deformations and linear scalings in the x, y and z axes. The rigid tool was positioned and oriented randomly within close proximity of the deforming mesh such that there was a mean separation of approximately 40 6 30 mm. After each random deformation and tool positioning, a proximity query was made between the tool and surface using DIBS, AABB, hybrid AABB and exhaustive searches. Based on the brain tissue strain limit given by Deck and Willinger,27 the DIBS hierarchy was constructed with a conservative limiting stretch factor of 2 (i.e. allowing for a tissue stretch of up to 100%). Experiments were used to ascertain the optimum level for splitting bottom-up and top-down updates of the hybrid hierarchy for each mesh size. In this experiment, the deformation rate was equal to the traversal rate, such that the total search time was calculated from the sum of the refit time and the traversal time. Each of the four meshes was randomly deformed and searched 10,000 times with the four methods. The results of these tests are shown in

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

Bowyer and Rodriguez y Baena

357

Figure 3. Multi-resolution brain mesh used within the randomised validation of the DIBS method: (a) 500 polygons, (b) 5000 polygons, (c) 50,000 polygons and (d) 500,000 polygons.

Table 1. Comparative performance for proximity queries between randomly positioned spheres and deforming meshes of differing sizes. Mesh size

Method

Refit time (ms)

Traversal time (ms)

Polygon separation tests

Total search time (ms)

500 polygons

DIBS AABB Hybrid AABB Exhaustive DIBS AABB Hybrid AABB Exhaustive DIBS AABB Hybrid AABB Exhaustive DIBS AABB Hybrid AABB Exhaustive

0:00 6 0:00 0:04 6 0:01 0:01 6 0:00 0:00 6 0:00 0:00 6 0:00 0:52 6 0:04 0:11 6 0:01 0:00 6 0:00 0:00 6 0:00 5:55 6 0:39 1:44 6 0:16 0:00 6 0:00 0:00 6 0:00 75:09 6 1:05 19:93 6 0:43 0:00 6 0:00

0:02 6 0:01 0:00 6 0:00 0:02 6 0:01 0:02 6 0:01 0:06 6 0:02 0:01 6 0:00 0:06 6 0:03 0:20 6 0:02 0:24 6 0:09 0:02 6 0:01 0:16 6 0:08 1:92 6 0:10 0:63 6 0:25 0:09 6 0:05 0:60 6 0:35 18:23 6 0:47

52 6 20 14 6 7 14 6 7 500 6 0 125 6 46 26 6 15 26 6 15 5000 6 0 253 6 121 67 6 45 67 6 45 50,000 6 0 639 6 325 187 6 134 187 6 134 500,000 6 0

0:02 6 0:01 0:05 6 0:01 0:03 6 0:01 0:02 6 0:01 0:06 6 0:02 0:52 6 0:04 0:17 6 0:03 0:20 6 0:02 0:24 6 0:09 5:57 6 0:39 1:60 6 0:18 1:92 6 0:10 0:63 6 0:25 75:18 6 1:05 20:53 6 0:56 18:23 6 0:47

5000 polygons

50,000 polygons

500,000 polygons

DIBS: deformation invariant bounding spheres; AABB: axis-aligned bounding boxes. The total search time for each test was calculated from the sum of the refit time and the traversal time. The results shown are the means and standard deviations of 10,000 tests.

Table 1 for the sphere primitive and Table 2 for the capsule primitive. The graphs in Figure 4 show the trends for total search times across the four meshes. A one-way analysis of variance test was applied to the total search times for each mesh size, and all were found to be significantly affected by the search method (p \ 106 ). Tukey’s honest significant difference test was then used to perform pairwise comparisons between each search method, and all were found to be significantly different from all others at the 99.9% confidence interval.

Beating heart experiment To assess the DIBS hierarchies in a scenario closer to conventional active constraints, an experiment was performed using a simulation of a beating heart (taken from Kwok et al.30). The heart simulation was a surface mesh of 10,272 triangles, with each vertex’s position defined based on a sinusoidal interpolation between three control frames. In practice, a periodic motion, such as that of a beating heart, presents efficient

opportunities for precomputed bounding volume hierarchies or distance fields for each time step (as in Ren et al.11); however, for the purposes of this experiment, it was posited that the deformation could not be predicted and therefore real-time updates were required. In this experiment, the simplified tool geometry was positioned relative to the beating heart by a human operator with the use of a Phantom Omni (SensAble Technologies, Inc.) haptic device, without force feedback. The human operator moved the capsule both around and through the beating heart, while the bounding volume hierarchies were used to compute the closest point pairs. Based on the estimated 15% distensibility of heart muscle,26 a conservative limiting stretch factor of 1.5 was sufficient for this experiment. The supplementary video (See Supplementary material online) demonstrates the procedure during these experiments, an image from which is shown in Figure 5. To ensure the integrity of the experiment, all four search methods were evaluated at each time step, and therefore, they all operated upon the same deformation and position configuration. The mean search frequency

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

358

Proc IMechE Part H: J Engineering in Medicine 228(4)

Table 2. Comparative performance for proximity queries between randomly positioned capsules and deforming meshes of differing sizes. Mesh size

Method

Refit time (ms)

Traversal time (ms)

Polygon separation tests

Total search time (ms)

500 polygons

DIBS AABB Hybrid AABB Exhaustive DIBS AABB Hybrid AABB Exhaustive DIBS AABB Hybrid AABB Exhaustive DIBS AABB Hybrid AABB Exhaustive

0:00 6 0:00 0:04 6 0:01 0:01 6 0:00 0:00 6 0:00 0:00 6 0:00 0:52 6 0:04 0:11 6 0:01 0:00 6 0:00 0:00 6 0:00 5:55 6 0:39 1:44 6 0:16 0:00 6 0:00 0:00 6 0:00 75:09 6 1:05 19:93 6 0:43 0:00 6 0:00

0:04 6 0:02 0:04 6 0:02 0:05 6 0:03 0:11 6 0:03 0:11 6 0:04 0:08 6 0:04 0:12 6 0:07 1:01 6 0:06 0:32 6 0:16 0:17 6 0:11 0:30 6 0:19 9:48 6 0:34 0:90 6 0:42 0:45 6 0:31 0:92 6 0:60 86:79 6 1:97

55 6 22 12 6 7 12 6 7 500 6 0 131 6 54 23 6 15 23 6 15 5000 6 0 268 6 171 59 6 49 59 6 49 50,000 6 0 664 6 407 163 6 132 163 6 132 500,000 6 0

0:04 6 0:02 0:08 6 0:02 0:06 6 0:03 0:11 6 0:03 0:11 6 0:04 0:59 6 0:06 0:24 6 0:07 1:01 6 0:06 0:32 6 0:16 5:72 6 0:41 1:74 6 0:25 9:48 6 0:34 0:90 6 0:42 75:54 6 1:09 20:85 6 0:73 86:79 6 1:97

5000 polygons

50,000 polygons

500,000 polygons

DIBS: deformation invariant bounding spheres; AABB: axis-aligned bounding boxes. The total search time for each test was calculated from the sum of the refit time and the traversal time. The results shown are the means and standard deviations of 10,000 tests.

2

2

10 DIBS AABB Hybrid AABB Exhaustive

1

10

Total search time (ms)

Total search time (ms)

10

0

10

−1

10

−2

10

DIBS AABB Hybrid AABB Exhaustive

1

10

0

10

−1

10

−2

2

10

3

4

5

10 10 10 Mesh size (number of polygons)

6

10

10

2

10

(a)

3

4

5

10 10 10 Mesh size (number of polygons)

6

10

(b)

Figure 4. Log–log plots of total search time against mesh size for the four proximity query methods: (a) sphere tool primitive and (b) capsule tool primitive. DIBS: deformation invariant bounding spheres; AABB: axis-aligned bounding boxes.

for proximity queries across the entire experiment is given in Table 3.

Discussion The results shown in Tables 1–3 illustrate that using DIBS hierarchies in proximity queries, between continuously deforming and rigid geometries, can be more efficient than two widely used bounding volume hierarchy methods and exhaustive searches. In Figure 4, the total search times for DIBS hierarchies are shown to be lower than those for the other methods across all four mesh sizes. Additionally, the trend for the DIBS hierarchies has a lower gradient than for the other methods, such that increasing mesh sizes gives DIBS an increasing benefit when compared to the others. For the two larger meshes, the DIBS

hierarchy was the only method of the four which returned from a proximity query, on average, faster than 1000 Hz. For haptic applications, such as active constraints, 1000 Hz is generally considered to be the control frequency required for smooth and stable interactions,31 and therefore, this improvement is significant. To create bounding spheres which are invariant to deformation, it is necessary to make them oversized and therefore looser in bounding the objects which they represent. The consequence of this is shown in Tables 1 and 2, where both the traversal times and the number of polygon separation tests for the DIBS hierarchy are often greater than those for the AABB methods. The reduced efficiency of traversing a DIBS hierarchy, when compared to an already refitted AABB hierarchy, means that its intended application must be considered before it is used.

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

Bowyer and Rodriguez y Baena

359 Table 3. Comparative performance for proximity queries between a capsule and a deforming heart mesh. Proximity query method

Frequency (Hz)

DIBS AABB Hybrid AABB Exhaustive

7222 873 1819 503

DIBS: deformation invariant bounding spheres; AABB: axis-aligned bounding boxes. Results shown are the mean produced from a total of 23,081 queries.

Figure 5. Image from the proximity query experiment between a beating heart mesh and a human positioned capsule. The red line identifies the pair of closest points.

speed of the AABB methods will increase, while the DIBS speed will remain the same. The result of this is that as the refit rate reduces, AABB hierarchies will eventually become faster than DIBS hierarchies. The point at which the refits are rare enough such that the AABB hierarchy becomes more efficient than DIBS is defined by a crossover refit ratio hc , which is a function of the DIBS traversal time tdt , the AABB traversal time tat and the AABB refit time tar

0.35 Sphere Capsule

Crossover refit ratio (ηc)

0.3

hc =

0.25 0.2 0.15 0.1 0.05 0 102

104 Mesh size (number of polygons)

106

Figure 6. Plot of the crossover refit ratios hc , above which DIBS is more efficient. The computed values are based on the mean traversal and refit times taken from Tables 1 and 2. No value is shown for the capsule with the 500 polygon mesh because the DIBS traversal time was actually shorter than the AABB traversal time such that DIBS is faster for all refit ratios. For this small mesh, the search times are so short that this is thought to be an anomalous result.

In the tests shown here, the object deforms and the hierarchy is updated at every search. That is to say that the ratio between refits and traversals is equal to one. This ratio is referred to as the refit ratio h h=

hierarchy refits per second hierarchy traversals per second

ð10Þ

If the object deforms only intermittently, and the refits are required less frequently, or if multiple tools are considered at each step, then the refit ratio will be less than one. When this happens, the cost of each refit can be spread across several traversals, and the average

tdt  tat tar

ð11Þ

The experimental values for the crossover refit ratios, based on the mean times in Tables 1 and 2, are shown in Figure 6. For the smallest mesh, the crossover refit ratio is relatively high (approximately one refit in three traversals) limiting the number of applications where DIBS hierarchies would be advantageous. However, exhaustive searches are almost efficient for these small meshes, and, in practice, a bounding volume hierarchy would be unnecessary. It is for larger meshes where the crossover refit ratio drops and DIBS hierarchies are of most benefit. For the largest mesh, the crossover refit ratio is less than 0.01, meaning that if the mesh deforms more often than once in 100 traversals, then DIBS hierarchies are more useful than AABB hierarchies. For a conventional traversal rate of 1000 Hz, this translates to a deformation frequency of 10 Hz, which is well within the range of most conventional tracking and imaging systems. Although Figure 6 shows that DIBS hierarchies can be slower to search than conventional methods when the deformation rate is low, DIBS hierarchies are not necessarily unsuitable in these cases. When the refit ratio is less than one, the time for each AABB refit is spread across several searches to guarantee a consistent control rate. The effect of this is that there is a latency between an input deformation and the enforcement of new constraints based on that deformation. If the motion of a constrained region is high, then this delay can cause inaccuracies in the constraint positioning, instabilities or poor surface rendering for the user. All of these problems can be mitigated with the immediate updates achieved using DIBS hierarchies.

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

360

Proc IMechE Part H: J Engineering in Medicine 228(4)

Conclusion and future work

Supplementary Material

To enforce a dynamic active constraint, it is necessary to compute the relative configuration between a continuously deforming tissue region and a rigid surgical tool. DIBS hierarchies have been presented as a new form of bounding volume hierarchy which, once constructed, continues to cover the underlying geometry, while it undergoes deformation. DIBS hierarchies allow tissue–tool relationships to be computed more efficiently than using a conventional, refittable, bounding volume hierarchy. The proposed method was used to perform proximity queries against a polygonal mesh and was compared to two widely used bounding volume hierarchies, which can be efficiently refitted after deformation. Both a fully randomised experiment and a human-controlled experiment showed that the DIBS hierarchy can be searched faster than these alternatives when the deformation rate is sufficiently high. The results show that for a mesh with 500,000 vertices, which is typical in medical imaging applications, the DIBS hierarchy is faster than the investigated alternatives if the model deforms faster than 10 Hz (as could be acquired from a conventional imaging system). Additionally, as the rate of imaging and real-time tissue modelling techniques improve, the benefit from using DIBS hierarchies will increase further. Results show that deformation invariance, which is achieved by pinning and inflating conventional bounding spheres, causes hierarchy traversals to be somewhat less efficient, when compared to AABB-based methods. However, it was shown that, for the cases considered, the decreased traversal efficiency is more than compensated for by the increased overall efficiency from not requiring hierarchy refits. There are various modifications which could be made to the method presented here to improve its efficiency, usefulness and applicability. Using Dijkstra’s algorithm to compute geodesic separations from mesh edges will likely mean that the path is longer than that of the physical surface, reducing the traversal efficiency. Using more advanced geodesic computation methods, such as those of Novotni and Klein32 and Surazhsky et al.,33 could improve this. Additionally, if the deformation of a model is described by a combination of vertex displacement functions, then ‘lazy evaluation’ of the displaced vertices could be incorporated into the DIBS hierarchy traversal, significantly reducing the number of vertex positions which must be updated during each search.

The online video is available at http://pih.sagepub.com/ supplemental-data

Declaration of conflicting interests The authors declare that there is no conflict of interest. Funding This work was supported by EU Grant FP7-ICT-20096-270460 on Active Constraints Technologies for Illdefined or Volatile Environments (ACTIVE).

References 1. Bowyer SA, Davies BL and Rodriguez y Baena F. Active constraints/virtual fixtures: a survey. IEEE T Robot 2014; 30(1): 138–157. 2. Li M, Ishii M and Taylor RH. Spatial motion constraints using fixtures generated by anatomy. IEEE T Robot 2007; 23(1): 4–19. 3. Bettini A, Marayong P, Lang S, et al. Vision-assisted control for manipulation using virtual fixtures. IEEE T Robot 2004; 20(6): 953–966. 4. Kwok K-W, Tsoi KH, Vitiello V, et al. Dimensionality reduction in controlling articulated snake robot for endoscopy under dynamic active constraints. IEEE T Robot 2013; 29(1): 15–31. 5. Bowyer SA and Rodriguez y Baena F. Dynamic frictional constraints for robot assisted surgery. In: Proceedings of the IEEE World Haptics Conference, Daejeon, Korea, 14–17 April 2013, pp.319–324. New York: IEEE. 6. Petersen JG and Rodriguez y Baena F. A dynamic active constraints approach for hands-on robotic surgery. In: Proceedings of the 2013 IEEE/RSJ international conference on intelligent robots and systems, Tokyo, Japan, 3–7 November 2013, pp.1966–1971. New York: IEEE. 7. Bowyer SA and Rodriguez y Baena F. Dynamic frictional constraints in translation and rotation. In: Proceedings of the IEEE international conference on robotics and automation, Hong Kong, China, 31 May–7 June 2014. New York: IEEE. 8. Park S, Howe RD and Torchiana DF. Virtual fixtures for robotic cardiac surgery. In: Proceedings of the 4th international conference on medical image computing and computerassisted intervention, Utrecht, The Netherlands, 14–17 October 2001, pp.1419–1420. London: Springer-Verlag. 9. Prada R and Payandeh S. On study of design and implementation of virtual fixtures. Virtual Real 2009; 13(2): 117–129. 10. Teschner M, Kimmerle S, Heidelberger B, et al. Collision detection for deformable objects. Comput Graph Forum 2005; 24(1): 61–81. 11. Ren J, Patel RV, McIsaac KA, et al. Dynamic 3-D virtual fixtures for minimally invasive beating heart procedures. IEEE T Med Imaging 2008; 27(8): 1061–1070. 12. Navkar NV, Deng Z, Shah DJ, et al. Visual and forcefeedback guidance for robot-assisted interventions in the beating heart with real-time MRI. In: Proceedings of the IEEE international conference on robotics and automation, Saint Paul, MN, 14–18 May 2012, pp.689–694. New York: IEEE. 13. Ryde´n F and Chizeck HJ. Forbidden-region virtual fixtures from streaming point clouds: remotely touching and protecting a beating heart. In: Proceedings of the IEEE/ RSJ international conference on intelligent robots and systems, Vilamoura, 7–12 October 2012, pp.3308–3313. New York: IEEE. 14. Ericson C. Real-time collision detection (ed Cox T). San Francisco, CA: Morgan Kaufmann Publishers Inc., 2004. 15. Lin MC and Gottschalk S. Collision detection between geometric models: a survey. In: IMA conference on

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015

Bowyer and Rodriguez y Baena

16.

17.

18.

19.

20.

21.

22.

23.

24.

361

mathematics of surfaces, Birmingham, UK, 31 August–2 September 1998, pp.37–56. Southend-on-sea: The Institute of Mathematics and its Applications. Quinlan S. Efficient distance computation between non-convex objects. In: Proceedings of the IEEE international conference on robotics and automation, San Diego, CA, 8–13 May 1994, vol. 4, pp.3324–3329. New York: IEEE. O’Sullivan C and Dingliana J. Real-time collision detection and response using sphere-trees. In: Proceedings of the spring conference on computer graphics, Budmerice, Slovakia, 28 April–1 May, 1999, pp.83–92. Bratislava: Comenius University. Van den Bergen G. Efficient collision detection of complex deformable models using AABB trees. J Gr Tool 1998; 2(4): 1–13. Gottschalk S, Lin MC and Manocha D. OBBTree: a hierarchical structure for rapid interference detection. In: Proceedings of the 23rd annual conference on computer graphics and interactive techniques, SIGGRAPH ‘96, New Orleans 4–9 August 1996, pp.171–180. New York: ACM. Klosowski JT, Held M, Mitchell JSB, et al. Efficient collision detection using bounding volume hierarchies of kDOPs. IEEE T Vis Comput Gr 1998; 4(1): 21–36. Zachmann G. Rapid collision detection by dynamically aligned DOP-trees. In: Proceedings of the IEEE 1998 virtual reality annual international symposium, Atlanta, GA, 14–18 March 1998, pp.90–97. New York: IEEE. Brown J, Sorkin S, Bruyns C, et al. Real-time simulation of deformable objects: tools and application. In: Proceedings of the 14th international conference on computer animation, Seoul, Korea, 7–8 November 2001, pp.228–258. New York: IEEE. Larsson T and Akenine-Mo¨ller T. Collision detection for continuously deforming bodies. In: Eurographics, 5–7 September 2001, pp.325–333. Oxford: Blackwell Publishers Ltd. James DL and Pai DK. BD-Tree: output-sensitive collision detection for reduced deformable models. In:

25.

26. 27.

28.

29. 30.

31.

32.

33.

Proceedings of the annual conference on computer graphics and interactive techniques, SIGGRAPH ‘04. Los Angeles, USA, 8–12 August, 2004, pp.393–398. New York: ACM. Spillmann J, Becker M and Teschner M. Efficient updates of bounding sphere hierarchies for geometrically deformable models. J Vis Commun Image R 2007; 18(2): 101– 108. Fung YC. Biomechanics: mechanical properties of living tissues. 2nd ed.New York: Springer-Verlag, 1993. Deck C and Willinger R. Improved head injury criteria based on head FE model. Int J Crashworthiness 2008; 13(6): 667–678. Mezger J, Kimmerle S and Etzmub O. Hierarchical techniques in collision detection for cloth animation. Journal WSCG 2003; 11(1): 322–329. Dijkstra EW. A note on two problems in connexion with graphs. Numer Math 1959; 1: 269–271. Kwok K-W, Vitiello V and Yang GZ. Control of an articulated snake robot under dynamic active constraints. In: Proceedings of the international conference on medical image computing and computer-assisted intervention (Lecture Notes in Computer Science), Beijing, China, 20–24 September 2010, vol. 6363, pp.229–236. Berlin/Heidelberg: Springer. Zhuang Y and Canny J. Haptic interaction with global deformations. In: Proceedings of the IEEE international conference on robotics and automation, San Francisco, CA, 14–28 April 2000, vol. 3, pp.2428–2433. New York: IEEE. Novotni M and Klein R. Computing geodesic distances on triangular meshes. In: Proceedings of the 10th international conference in Central Europe on Computer Graphics, Visualization and Computer Vision, Plzen-Bory, Czech Republic, 4–8 February 2002, pp.341–347. Vaclav Skala: Union Agency. Surazhsky V, Surazhsky T, Kirsanov D, et al. Fast exact and approximate geodesics on meshes. SIGGRAPH ’05. Los Angeles, USA, 2–4 August, 2005, pp.553–560. New York: ACM.

Downloaded from pih.sagepub.com at MICHIGAN STATE UNIV LIBRARIES on June 12, 2015