P H S Perrusi

P Baksic

H Courtecuisse

Interactive Finite Element model of needle insertion and laceration

Keywords:

This paper introduces an interactive model of needle insertion, including the possibility to simulate lacerations of tissue around the needle. The method relies on complementary constraints to couple the Finite Element models of the needle and tissue. The cutting path is generated from mechanical criteria (i.e. cutting force) at arbitrary resolution, avoiding expensive remeshing of Finite Element meshes. Complex behavior can be simulated in real time such as friction along the shaft of the needle, puncture and cutting force resulting from interactions of the needle with the tissue. The method is illustrated both in an interactive simulation of a needle insertion/cutting and in a robotic needle insertion in liver tissue during the breathing motion.

CCS Concepts

•

Introduction

Needle-based approaches are more and more chosen as a viable option for various medical applications ranging from biopsies to surgical treatment (such as Radio Frequency Ablation of liver tumors). Yet, despite many benefits for patients it also increases the technical difficulty for practitioners, pushing the need for training solutions. Simulation has become a powerful tool in learning [START_REF] Courtecuisse | Simulation-based medical teaching and learning[END_REF], especially in the medical field where the practice of surgical procedures is essential to reduce the post-operative damages. Among the difficulty to accurately place needles, additional constraints lies in the fact that during the insertion needles may lacerate and damage surrounding tissues due to physiological motions of organs. Although the lacerations are limited when the needles are held by a practitioner (avoiding applying too much stress on the needle inside tissues), they are particularly observed when needles are held by a stiff robotic system [BJMB * 14], which is why some robotic devices were conceived to avoid this behavior [PBB * 09]. For percutaneous applications, the extents of the cuts are usually not sufficient to separate the tissue. Indeed, the compression of surrounding organs and the limited extension of the cut surface prevent the organs to mechanically separate in two parts in the abdominal cavity. The present method relies on this observation and introduces a new model for needle insertion and lacerations. The proposed method relies on an intermediate plan defined independently from the needle and tissue mesh resolutions. The rupture and the extension of the cut plane are triggered by the amount of interaction forces between the needle and the tissue. The cut plane is used to simulate the additional degrees of freedom of the needle inside the tissue, resulting in an interactive model of the cut, without any need for expensive remeshing operations.

Related Works

One of the first models for needle-deformable tissue interactions was proposed by DiMaio et al. [START_REF] Dimaio | Needle insertion modelling for the interactive simulation of percutaneous procedures[END_REF]. Relying on a linear model, the method allows simulating stick-slip behaviors. Significant amount of work was also performed by Misra et al. [APM * 16], including experimental validations with robotic systems, in order to steer bevel-tip needles. From a numerical point of view, large scale simulations involving multiple organs and complex tissue interactions were proposed for training. Chentanez et al. [CAR * 09] presented a FE approach based on the beam's theory to predict needle deflection. Yet, the method relies on expensive remeshing operations which raises additional difficulties to handle topological changes and deformations as the needle penetrates the tissue, especially in an interactive environment. Duriez et al. [DGM * 09] introduced complementary constraints for the simulation of a needle insertion, avoiding this way remeshing operations. The method allows simulating complex phenomena such as tissue deformations, needle-tissue friction, and puncture force in real time. Fast GPU-based implementation has also been proposed [Pat12] for both nonlinear and linear elasticity. From a general point of view, the simulation of cut and topological modification has been widely studied. Methods can be classified into three parts, namely the mesh-based methods, the meshless methods and hybrid methods [WM18]. Mesh-based methods simulate the deformations using FE and the cutting process with topological changes. Among those changes, the easiest one is the element deletion, which consists in deleting an element that crosses the cutting path [FDA05]. However, element deletion induces volume loss in the mesh and an approximative cutting surface. Other methods use element duplication for each element crossing the cutting path [MBF04] but the volume of the mesh is still modified. To alleviate these issues, some methods proposed to re-mesh locally in order to fit the cutting path [PUC * 15]. However, although having a better cutting surface, it may create ill-shaped elements that lead to instabilities during the simulation. Meshless methods use particles and their interactions with neighbors to compute the deformation. In [CLLZ17], the cutting process is simulated when the forces increase above a threshold separating the two parts around the cutting plane. The main advantage is the reduced computation time but meshless methods lack of geometrical information that leads to a less efficient collision detection. Hybrid methods combine the advantages of both using a meshed surface and particles to compute the deformations. In [SLLW18], the authors reconstruct the surface mesh using an interactive triangulation algorithm during the cut. The cutting process can also be simulated using XFEM [KBT17], that will enrich the shape functions of the elements to introduce discontinuities. Although lacerations represent an important concern from the medical point of view, very few methods have been proposed for the simulation of the cut by the needle. Bui et al. recently proposed to rely on nonconforming meshes [BTB19], but the method still raises stability and computational issues in an interactive context. The contribution of this paper is a hybrid method where the cut is modeled through complementary constraints. The method is an extension of [DGM * 09] introduced for needle insertion. Based on mechanical criteria the method generates an intermediate plane, mapped into the volume mesh, describing the cutting path. Interactions between the needle and the cut surface are then handled by mechanical constraints, allowing this way to properly simulate the effects of the cut with a minimal amount of topological changes.

Constraint-based model of needle insertion

This section introduces the principles of the method inherited from [DGM * 09]. The needle is described as a set of serially linked beams composed of two nodes, each node having 6 DoFs (position and rotation). The model of the needle is based on the Timoshenko's formulation, which relies on the beam's theory. The volumic model is meshed with linear tetrahedral elements and follows the Hook's law. In addition, to handle large displacements both models are combined with the co-rotational approach. After the discretization, the internal forces can be written as a nonlinear function F(q) where q are the nodal positions of the FE meshes. Note that the proposed solution is generic and could be applied to other nonlinear models. An implicit integration of Newton's second law is chosen in order to increase the stability during interactions. Lagrangian Multipliers are used to enforce constraints between the needle and the tissue. The non-linear problem formulation is explained in [AGM * 19]. The problem is solved with a single Newton-Raphson iteration (see [CAK * 14] for justification). The linearized problem can be represented by the following equation:

Wλ = δ free (1)
where W is the schür complement of the system, λ are the constraint forces and δ free the free-motion interpenetration. Using complementary constraints both λ and δ free are unknown. Indeed, con-straints used in this work need to satisfy the Signorini's law:

0 ≤ δ ⊥ λ ≥ 0 (2)
The above equation means that either the objects are distant δ > 0 and no contact force is applied λ = 0, or instead constraints are violated and must be solved δ = 0 with a non-null constraint force λ > 0. In addition Coulomb friction can be added to model the friction of the needle with the tissue, producing this way a Non-Linear Complementary Problem (NLCP).

The NLCP is solved with a modified Gauss-Seidel (GS) algorithm described in [CAK * 14]. Although the convergence of the GS is generally slow, the size of W (being the number of constraints) remains small, offering this way the possibility to use the method in an interactive context. The interaction model is composed of 3 types of atomic constraints:

1. Bilateral Atomic Constraints φ are holonomic constraints φ(q) = 0, used to impose the relative positions of constrained points on the orthogonal plane defined by the constraint normal.

The combination of two Bilateral constraints enforce the points to follow a line at the intersection of the two planes. 2. Unilateral Atomic Constraints χ are used to enforce Signorini conditions χ(q) > 0. A force is generated if and only if the constrained point moves on the negative side of the plane defined by the constraint normal. 3. Friction Atomic Constraints ψ resist against the motion of the needle. The complementarity constraint defines two states: stick when there is no relative motion due to static friction and slip when the relative motion is not null (dynamic friction). After the linearization of equation (1), each atomic constraint is associated with a direction n c s (called normal) and a Lagrangian Multiplier λ c s corresponding to the force needed to satisfy the atomic constraint, where c ∈ {φ, χ, ψ}. The needle penetration and insertion phases (Figure 1a and1b) are similar to [DGM * 09]: when the needle is located outside, a penetration constraint (χp, ψp, ψp) is applied between the needle tip and the closest point located on the triangulated surface. Friction is also added in the tangential plane. Inside the tissue, insertion constraints (ψ i , φ i , φ i ) are added along the path taken by the needle. A new constraint is created as soon as the needle's tip is located at a distance d i of all the other constraints. Atomic constraints φ i prevent lateral motion of the needle along the inserted trajectory, and ψ i adds friction along the needle shaft.

Needle-tissue cutting and laceration

The cutting phase starts when the average lateral force reaches a cutting force threshold fc, denoted as tearing criteria:

t i > n i * fc, with t i = ∑ n φ i λ φ i (3)
where n i is the number of insertion constraints and t i the initial direction of the cutting surface. For each insertion constraint, a point cloud pt of the cut surface is computed by pt = p i + dt * t i where p i are the position of the insertion constraints and dt a parameter controlling the extent of the cutting surface. The Point-Cloud Library (PCL) [MRB09] is used to generate a triangular mesh S from positions of the point cloud. First, a polysubmitted to EUROGRAPHICS 2021. If several constraints are associated to the same needle's edge, only the closest one is kept in order to limit the computational expense.

A second mechanical criterion is introduced to extend the cutting path. If any contact force λ χ c associated to Unilateral Atomic Constraints χc between the needle and the border of S reaches a tearing force threshold ft , the cutting surface is extended in the direction of the contact normal -n χ c at distance dt . The cutting surface (Figure 1c) being extended at different times of the simulation, it must be defined independently from the current deformation in order to get consistent triangulation of the point cloud over the various extents of the cut. For this purpose the points of S are mapped into the tetrahedral mesh using barycentric coordinates. For each new point introduced in the cutting surface S, the position pr is computed as follows:

pr = 4 ∑ k=1 α k q(l) k ( 4
)
where l is the index of the tetrahedron containing the point and α k the associated barycentric coordinate in the deformed state. q(l) k is the nodal position of the tetrahedron l at the initial position. The set of points pr provides a static and undeformed description of the cutting surface. These positions are used to generate the triangulation with PCL, enforcing this way a constant triangulation for all the history of the cut while taking into account new points of the cut. Finally, equation ( 4) is also applied with the current nodal position of the tissue, providing the position of the cutting surface in the deformed state on which are applied the interaction constraints.

Results and performances

The parameters of the simulations are the distance between constraints d i during needle insertion, the distance dc between points of the cut surface, and the geometrical parameters for triangular mesh generation θ min , θmax and θs. In addition, mechanical parameters include the Young modulus E and the Poisson ν ratio of both the needle and the volume. The penetration force fp, the cutting force fc, the tearing force ft , the friction at the surface µp and during the insertion µ i . Although, the patient-specific identification of the parameter goes far beyond the scope of this paper, each of them is easily parametrizable in order to obtain the desired behavior for a training simulator. The method is first evaluated according to the computation time. Since no remeshing is necessary to perform the cut, the computation time is almost not impacted by the model (see Tab. 1). A first simulation is performed involving a volume mesh composed of 3072 tetrahedron elements with E = 500kPa and ν = 0.4 to simulate the tissue. The rigid needle's model is 32cm long composed of 16 beam elements with 0.8cm radius, E = 100GPa and ν = 0.4. In a first phase, a straight needle insertion 17mm places the needle inside the tissue with fp = 0.8N. Then lateral motion of 8mm is applied to the needle base and tissue tearing is observed when the cutting criteria fc = ft = 1N is reached. A second experiment is performed where the needle is controlled with a haptic device. During the insertion, the user can feel the various constraints (penetration, friction). When lateral motion is imposed to the needle's base, the cut is performed (red surface in Tab. 1) and the user can feel the cutting force necessary to extend the cut. The computation time to solve constraints remains constant for any step of the insertion. During cutting, the most expensive operation is the triangulation of the point cloud, representing on average 2 ms. Yet, this procedure is called only for new added points, i.e. every 19.8 time steps on average. At the end, the point cloud was composed of 110 points and a 184 element triangle mesh was generated. The GS requires on average 2.54 iterations during the insertion and 80 during the cuting phase. Since the matrices remains very small, the simulation maintains a frame rate of 65 FPS. 1: Performances for the various steps of the method. The size of the cube is 10 cm 3 . dc = d i = 0.5 mm.

Steps

An important advantage of the method lies in the fact that the cut surface is generated from mechanical information and not just geometric data. To enhance the importance of this contribution, the submitted to EUROGRAPHICS 2021.

model is employed during a robotic needle insertion inside the liver during the respiratory motion (see Fig. 2). A robotic motion has been precomputed in order to compensate for breathing motion (see [BCDB20]). The needle insertion is performed with both a flexible (E f lexible = 50GPa) and a rigid (E rigid = 200GPa) 12cm long needle following the same precomputed motion of the robot. As expected, the flexible needle deforms and does not reach the threshold to cut the tissue. Instead, the rigid needle introduces significantly higher stress in the volume, resulting in the cut of the tissue. The cutting path is therefore extended allowing the needle to slide along the cut direction.

6 Conclusion, discussion and future work This paper introduces an interactive model of needle/tissue insertion with the possibility to cut the tissue using mechanical criteria. The cutting path is obtained from the triangulation of a point-cloud whose positions are derived from interaction forces between the needle and the tissue. The method is introduced in the context of training with generic parameters. In terms of limitations, since the cutting surface is generated on the current deformed state of the tissue, the angle between triangle faces may vary significantly when remapped to the rest state. This is particularly observed when the tissue is significantly deformed and when rapid movements are imposed to the needle during the cutting phase. In this situation, the bilateral constraints can shift their directions rapidly when the needle is sliding through the cut plane and lead to stability issues in the simulation. However, the deformation of organs during percutaneous applications usually remains small, limiting the appearance of such instabilities. Sudden drops in the interaction forces were also observed when the contact force reaches the cutting threshold λ χ c ≥ ft to extend the cutting surface. However, this effect can be mitigated with the addition of friction constraints along the cutting surface. Finally, a future direction lies in the possibility to tune the parameters based on data provided by a force sensor installed on a robotic system. Drops in the force profile can then be detected to initiate the cut producing this way a patient specific model of the cut.

( a )Figure 1 :

 a1 Figure 1: Needle insertion and cutting: first the needle penetrates the tissue. Then needle/tissue constraints prevent the relative displacements of the models. Finally, in case of lateral motion, the cutting plane is generated from mechanical forces. 3 types of atomic constraints are used: Bilateral atomic constraints (blue), Unilateral atomic constraints (yellow) and Friction atomic constraints (red).nomial fitting, using Moving-Least Squares, is applied to estimate a normal direction for each point with respect to their neighbors. Then a greedy algorithm is used to triangulate the point cloud into a triangular mesh based on geometrical parameters. These parameters are the triangle's minimal angle θ min , their maximal angle θmax and maximum surface angle between triangles θs. During the cutting phase, all the insertion constraints along the needle are removed, except the first insertion point to avoid tearing the surface. Instead, a cutting constraint (φc, χc) is added for all the needle's points whose projection is located inside the cutting surface S. The cutting constraints are composed of one Bilateral Atomic Constraint φc between needle's points and their respective closest point on S. This constraint enforces the needle to slide along the cutting path. In addition, one Unilateral Atomic Constraint χc is applied between each point located on the border of the cutting surface of S and their respective closest point along the needle shaft. If several constraints are associated to the same needle's edge, only the closest one is kept in order to limit the computational expense. A second mechanical criterion is introduced to extend the cutting path. If any contact force λ

Figure 2 :

 2 Figure 2: External view of the needle insertion with robotic assistance simulation (left). Flexible needle close view (upper right). Rigid needle close view (lower right).

Table

		Atomic Constraint Total
		cst.	Solving (ms) (ms)
	Penetration	1	2 (13.2%)	15.1
	Insertion	36	3 (19.4%) 15.47
	Cutting	12	4 (24.3%)	12

submitted to EUROGRAPHICS 2021.

Acknowledgement: This work was supported by French National Research Agency (ANR) within the project SPERRY ANR-18-CE33-0007 and the Investissements d'Avenir program (ANR-11-LABX-0004, Labex CAMI).