# Efficient Needle Insertion Simulation using Hybrid Constraint Solver and Isolated DOFs

Claire Martin, Ziqiu Zeng, Hadrien Courtecuisse

► To cite this version:

Claire Martin, Ziqiu Zeng, Hadrien Courtecuisse. Efficient Needle Insertion Simulation using Hybrid Constraint Solver and Isolated DOFs. Eurographics 2023 - Short Papers, May 2023, Saarbrücken, Germany. 10.2312/egs.20231003. hal-04181906

HAL Id: hal-04181906

https://hal.science/hal-04181906

Submitted on 16 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient Needle Insertion Simulation using Hybrid Constraint Solver and Isolated DOFs

C. Martin<sup>1,2 ID</sup> , Z. Zeng<sup>1,2 ID</sup> , H. Courtecuisse<sup>1,2 ID</sup>

## Abstract

This paper introduces a real-time compatible method to improve the location of constraints between a needle and tissues in the context of needle insertion simulation. This method is based on intersections between the Finite Element (FE) meshes of the needle and the tissues. It is coupled with the method of isolating mechanical DOFs and a hybrid solver (implying both direct and iterative resolutions) to respectively generate and solve the constraint problem while reducing the computation time.

## CCS Concepts

* Computing methodologies → Physical Simulation; Real-Time Simulation;

## 1 Introduction

Needle-based procedures are beneficial for patients but are technically challenging to perform. Consequently, solutions for training are needed. Unlike traditional training methods, medical simulators allow training in many cases, from the most common to the rarest [TPRT00]. Moreover, trainees do not require the direct supervision of trained clinicians [SHZS*13], allowing saving time for the medical staff without compromising patients' safety. For this, simulators must be highly realistic and efficient regarding general modeling, interactions, and real-time performances.

Simulations of needle/tissue interactions have received extensive consideration (a recent survey can be found in [RM17]). Recent works [CNR*19] exposed the emerging use of haptics, along with virtual environments, in needle insertion simulation for training. The entire interaction process is highly nonlinear and usually includes several factors, such as the puncture force, the Stribeck effect (force drop after puncture), and the total force during both insertion and withdrawal as a non-linear function of insertion depth and velocity. An example can be found in [AKP10] which proposed a modified LuGre model. Okamura et al. [OSO04] consider the total force as the sum of stiffness (occurring before puncture), cutting (considered constant), and friction forces. The friction is described by a modified Karnopp model that covers both static and dynamic friction.

Considering interactive modeling, Roesthuis et al. [RVJM11] use springs to model needle-tissue interactions. Chentanez et al. [CAR*09] rely on Timoshenko beams' theory to predict the needle deflection inside soft tissues. Duriez et al. [DGM*09] proposed a constraint-based formalism to prescribe the relative displacements of the heterogeneous models of the needle (composed of beam ele-

ments with 6 DOFs) and a volume model (composed of tetrahedra or hexahedra), avoiding this way expensive remeshing operations.

The constraint solving for different contacts (unilateral, bilateral, frictional, etc.) is a general issue in Computer Graphics. The penalty methods [KKB18] are simple and fast but only result in approximate solutions at the end of time steps. It is also difficult to handle stable contacts with such methods. Incremental potential contact has proved to compute contact and friction very efficiently [LFS*20]. The constraint-based technique using Lagrange multipliers [DGM*09] gives accurate and stable solutions. After transferring the problems to the constraint space, they can be formulated as linear/non-linear complementarity problems according to the constraint type, which can be solved with different numerical methods [AE21]. Iterative methods such as projected Gauss-Seidel (PGS) [DDKA06, ZCC22] and Non-smooth Newton methods [MEM*19] are preferred in real-time simulations owing to the computing efficiency and the flexibility to control the time cost.

This work relies on the needle-tissue interaction model introduced in [DGM*09]. The location of constraints is improved by prescribing the relative displacements of all the DOFs involved during the insertion. It results from the intersections between the needle and volume meshes. Compared to previous solutions, our approach improves the accuracy and generates a realistic haptic rendering, but a significantly larger number of constraints are generated. To limit the computational cost, the method is combined with the IsoDOFs approach [ZCC22]. In addition, the large number of constraints raises convergence issues in solving constraint problems mixing bilateral and friction constraints. For this purpose, we introduce a hybrid solver, combining a direct resolution of bilateral constraints and the Gauss-Seidel method for friction.

<footnote n="1">¹ MIMESIS research team, INRIA, France</footnote>

<footnote n="2">² RDH/MLMS - ICube, University of Strasbourg, France</footnote>
## 2 Background

Timoshenko beams and a linear tetrahedral mesh are considered to model the needle and the volume, respectively. The co-rotational formulation of linear elasticity, described in [Fel00], treats large displacements of elements as small linear deformations coupled with rotations. It is used to model non-linearity while keeping real-time computational performances.

For both the needle and the tissue, the dynamic equation in the continuum mechanics can be formulated as a linear system, with the backward Euler method and a first-order approximation:

$$Ax = b + c \quad (1)$$

where the system matrix **A** gathers mass, damping, and stiffness, the vector **b** accumulates internal and external forces at the beginning of the time step, **x** is the unknown variation of velocities for the implicit integration, and **c** is the constraint forces (contact, friction...) contribution.

This work focuses on needle-tissue interactions, after puncture. In [DGM*09], as the needle progresses in the tissues, constraints are progressively generated at a distance *d* from the previous one along the needle shaft, as soon as this distance to the last constraint is reached by the needle tip. Each generated constraint group involves two bilateral constraints ($\lambda_b$) orthogonal to the needle shaft to laterally prescribe the needle-tissue relative motions, and one friction constraint ($\lambda_f$) along the needle shaft, gathered in a vector $\lambda$. The resulting system can be written as a Karush-Kuhn-Tucker:

$$\left\{ \begin{array}{l} \mathbf{A}_1 \mathbf{x}_1 - \mathbf{H}_1^T \lambda = \mathbf{b}_1 \\ \mathbf{A}_2 \mathbf{x}_2 + \mathbf{H}_2^T \lambda = \mathbf{b}_2 \\ \mathbf{H}_1 \mathbf{x}_1 - \mathbf{H}_2 \mathbf{x}_2 = \Delta \delta \end{array} \right. \quad (2)$$

where indices (1) and (2) denote respectively the volume and needle's models, $\Delta\delta = \delta^{t+h} - \delta^{explicit}$ is the gap between constraint violations after an explicit integration and at the end of the time step, and the constraint Jacobian **H** maps forces from the world to constraint space ($c_1 = \mathbf{H}_1^T \lambda$ and $c_2 = -\mathbf{H}_2^T \lambda$). One looks for $\lambda_{\text{bilateral}}$ such that $\delta_{\text{bilateral}}^{t+h} = 0$ for bilateral constraints, and $\lambda_{\text{friction}}$ follows a friction model derived from [DGM*09]. By eliminating the unknowns in Equation (2) we have:

$$\delta^{t+h} = \delta^{\text{free}} + \underbrace{\sum \mathbf{H} \mathbf{A}^{-1} \mathbf{H}^T \lambda}_{\mathbf{W}} \quad (3)$$

where **W** is the compliance matrix, $\delta^{\text{free}}$ is the constraint violation after "free motion" that describes physics dynamics without considering the constraints. When friction is added, the above equation describes a Non-Linear Complementary Problem (NLCP) where the state can change between stick and slip behaviors during the resolution depending on the force applied to the other constraints. Such NLCP can be solved by a PGS method [DDKA06], where the constraint function is verified in successive iterations. Finally, by including $\lambda$ into Equation (2), the constrained system is resolved.

An important limitation of this approach concerns the fact that the model relies on a distance parameter between constraints in order to save computation time and maintain stability. However, as indicated in [DGM*09], the accuracy of the interaction model decreases when the constraints are widely distributed. Conversely, too

tight constraints lead to over-constrained problems and raise stability and computation time issues. Indeed, building the compliance matrix **W** requires an extremely costly process that mainly depends on the dimensions of the constraints **c** and the mechanical DOFs **n**.

## 3 Methods

**Constraint definition:** To fully prescribe the needle-tissue relative displacements, all surrounding DOFs must be constrained (orange nodes in Fig. 1). To do so, we generate a new constraint at each intersection point between the edges of the needle and the triangles of the tetrahedral volume mesh. Constraints are then transferred to the mechanical DOFs of both models using barycentric coordinates.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Constraints are located at intersection points between the edges of the needle and the triangles of the tetrahedral mesh. The resulting constraints are irregularly spaced along the needle shaft.</figcaption></figure>

Consequently, no distance parameter between constraints is required, and the location of constraints is automatically adapted to the mesh resolution. However, fine meshes are expected to generate a large number of constraints, as the needle will intersect many elements. Regarding the considered interaction model, this will increase the size of the problem in constraint resolution. As exposed in Equation 3, building **W** is extremely costly.

**Compliance matrix assembling:** The constraint Jacobian **H** is sparse as each constraint impacts a limited number of DOFs. From this sparsity, Zeng et al. reformulated **W** in [ZCC22] to efficiently reduce its computational cost, also benefiting from a parallelized resolution on the GPU. Moreover, during needle insertion, the DOFs impacted by the constraints (also called isolated DOFs, or IsoDOFs) in the previously inserted part remain in the current time step. Therefore, the improved IsoDOFs method, which reuses common results in consecutive time steps, is exploited to further accelerate the resolution (see [ZCC22] for more details).

**Constraint solving:** Intersections occurring near a triangle vertex (for any mesh resolution) can result in close constraints. This will increase the condition number of the compliance matrix **W**, or even worse, lead to an over-constrained problem. Consequently, the PGS algorithm used to solve Equation (3) will hardly converge. On the other hand, the non-linear friction forces prevent the use of a direct solver since the friction status of constraints may change during iterations depending on the forces applied on other constraints.

To improve the constraint resolution step, a hybrid solver is used (see Fig. 2). The directions of all the bilateral forces are firstly included into **H**. The friction forces of the same constraint points are then considered. The resulting **W** is thus composed of two main blocks, $W_{\text{bilateral}}$ and $W_{\text{friction}}$. The coupling blocks (yellow blocks in Fig 2) describe the influence of forces on each other.
<figure><img src="image_3.png" /><figcaption>**Figure 2:** The bilateral constraints of the *n* constraint points are gathered in one single block of the compliance matrix **W**. The *n* friction constraints are then individually added to **W**. At each iteration of the Gauss-Seidel algorithm, the solver first inverts the bilateral block and then computes the friction forces depending on the result of the bilateral constraints resolution.</figcaption></figure>

<table><caption>**Table 1:** Maximum node displacement and amount of generated constraints depending on the considered insertion.</caption><thead><tr><th>Distance between constraints (cm)</th><th>Max node displacement (cm)</th><th>Generated constraints</th></tr></thead><tbody><tr><td>1</td><td>0.04167</td><td>32</td></tr><tr><td>0.5</td><td>0.02557</td><td>62</td></tr><tr><td>0.2</td><td>0.01991</td><td>152</td></tr><tr><td>0.1</td><td>0.01815</td><td>302</td></tr><tr><td>0.07</td><td>0.01815</td><td>430</td></tr><tr><td>0.04</td><td>0.01814</td><td>752</td></tr><tr><td>Intersection method</td><td>0.0175</td><td>206</td></tr></tbody></table>

An iteration *k* of the Gauss-Seidel algorithm starts with a direct resolution of all the bilateral forces. Yet, since a DoF may share several constraints, $W_{bilateral}$ is ill-conditioned. The problem is solved by computing a pseudo-inverse $W_{bilateral}^{+}$ (Equation (4)), which can be costly for large matrices. However, the constraints remain in much smaller amounts than the DOFs in the models: compared to the rest of the simulation, this cost remains negligible. Friction forces must then be computed individually due to their non-linearity, as usual Gauss-Seidel steps. Note that although the bilateral constraints are solved at once using the direct solver, their value may change depending on the frictional constraints at each Gauss-Seidel iteration due to the coupling between the constraints.

$$ \lambda_{\parallel \text{bilateral}}^{k+1} = -W_{\text{bilateral}}^{+} \delta_{\parallel \text{bilateral}}^{k} \quad (4) $$

where $\lambda_{\parallel \text{bilateral}}^{k+1}$ and $\delta_{\parallel \text{bilateral}}^{k}$ are respectively the force and violation vectors restricted to the bilateral constraints in iteration *k*.

# 4 Results and discussions

**Location of constraints:** To assess the improvement of constraint placement by the intersection process, a rigid needle was numerically inserted into a very soft beam-like gel (500 Pa) (see Fig. 3). A mesh of 50 nodes over 15 cm in the direction of the needle, against 6 nodes over 2 cm in the two other directions was used.

Insertions were performed with the intersection method, as well

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Bilateral constraints (green and red arrows) were generated between the rigid needle (purple) and the volume. Gravity was applied to generate a motion of the unconstrained volume nodes, producing unrealistic needle-tissue relative motions as circled.</figcaption></figure>

as the initial method with constraint spacing between 1 cm and 0.04 cm. After constraints were generated between the needle and the beam, gravity was applied, generating a motion of the soft beam while the rigid needle remained fixed. The displacements of nodes belonging to the beam and located around the needle shaft were measured. Table 1 shows that the distance model converges towards a maximal displacement of the volume nodes of 0.01815 cm with 302 constraints: all the nodes surrounding the needle are impacted by at least one constraint. Interestingly, our new constraint placement required 96 fewer constraints to impact the nodes around the needle, with a maximum node displacement of 0.0175 cm.

**Computational performances:** A needle was numerically inserted into a gel, following a defined curved trajectory. During the insertion, bilateral and friction constraints were generated from the intersection process. The time required to build **W** and the number of impacted DOFs were evaluated (see Fig. 4). Comparisons were made between the *IsoDOFs* and *reuse IsoDOFs* methods, exposed in [ZCC22]. Results show that only a minor portion of DOFs was newly constrained between consecutive time steps, thus decreasing the computation time required to build **W**: although our approach gives rise to a large number of constraints, the reuse of *IsoDOFs* allows the computation of **W** to be minimally impacted.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Average computation time to build **W** (in ms) over 50 time steps from the reuse and standard IsoDOFs methods, during a needle insertion. Amounts of reused and total IsoDOFs are compared. Most IsoDOFs are reused between consecutive time steps.</figcaption></figure>

Performances of the hybrid solver and the initial resolution process were then compared. The node spacing of the volumetric mesh is of the order of a centimeter. Constraints were regularly generated from the initial method. Simulations were run for constraint spacing decreasing from 0.012 m to 0.004 m, and the amount of Gauss-Seidel iterations was evaluated at the end of the insertion. The tolerance was set to $10^{-8}$; a maximum of 1000 iterations was allowed. With spaced constraints, the two solvers showed similar convergences, while the initial method hardly converged when constraints were closer than the order of the elements size (see table 2). The hybrid solver remained efficient, allowing the use of close constraints potentially generated by the intersection method.

**Friction force:** 44 mm deep needle insertions were performed at constant velocities. The force resulting from the friction constraints
<table><caption>**Table 2:** Average number of Gauss-Seidel iterations over 100 time steps, based on the constraint spacing and the considered solver</caption><thead><tr><th colspan="2">Constraint spacing (mm)</th><th>12</th><th>10</th><th>8</th><th>6</th><th>4</th></tr></thead><tbody><tr><th rowspan="2">Initial solver</th><th>Iterations</th><td>130</td><td>129</td><td>327</td><td>910</td><td>1001</td></tr><tr><th>Time (ms)</th><td>2.65</td><td>2.77</td><td>8.18</td><td>27.93</td><td>43.06</td></tr><tr><th rowspan="3">Hybrid solver</th><th>Iterations</th><td>134</td><td>130</td><td>134</td><td>142</td><td>133</td></tr><tr><th>Time (ms)</th><td>2.48</td><td>2.51</td><td>2.96</td><td>3.80</td><td>4.97</td></tr><tr><th>p-inv (ms)</th><td>0.02</td><td>0.02</td><td>0.03</td><td>0.04</td><td>0.10</td></tr></tbody></table>

is shown to be related to insertion velocity in Fig. 5. Although constraints are irregularly spaced along the needle shaft due to intersections with the tetrahedral mesh, the force profile remains smooth for constant velocities, creating a realistic haptic response for the user.

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Force profiles resulting from the friction model and constraint resolution methods for different insertion velocities. Intersection constraints were generated at the dotted lines.</figcaption></figure>

**Applications:** We have applied our method to some examples in medical simulations, and computed average measurements over 50 time steps. In 6a, the needle tip reaches a liver tumor, while 205 constraints (including collisions between the organs) impact 1238.7 isolated DOFs in average. With our approach, the simulator takes 5.65 ms to build W with the IsoDOFs method and 24.74 ms to solve the constraint with the hybrid solver. In a similar simulation with a kidney 6b, 136.80 constraints are generated, while 1169.70 isolated DOFs are impacted at the end of the insertion. The simulator takes 4.12 ms to build the compliance matrix and 6.56 ms to solve the constraints.

## 5 Conclusion and future work

In this work, a needle/tissue interaction method was developed to efficiently compute constraints between a needle and a tissue in the context of needle insertion simulation. This method was coupled with IsoDOFs and a hybrid solver, respectively accelerating the compliance matrix and the forces computations, in order to decrease the computation time required for the resolution step. In future work, we plan to integrate our model into an immersive virtual reality simulator for percutaneous procedure training.

**Acknowledgement:** This work was supported by French National Research Agency (ANR) within the project SPECULAR ANR-21-CE45.

## References

[AE21] ANDREWS S., ERLEBEN K.: Contact and friction simulation for computer graphics. *ACM SIGGRAPH 2021 Courses, SIGGRAPH 2021* (8 2021). 1

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Application of our approach in medical simulations.</figcaption></figure>

[AKP10] ASADIAN A., KERMANI M. R., PATEL R. V.: A compact dynamic force model for needle-tissue interaction. *2010 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBC'10* (2010), 2292–2295. 1

[CAR*09] CHENTANEZ N., ALTEROVITZ R., RITCHIE D., CHO L., HAUSER K. K., GOLDBERG K., SHEWCHUK J. R., O'BRIEN J. F.: Interactive simulation of surgical needle insertion and steering. *ACM Transactions on Graphics* 28, 3 (jul 2009), 1–10. 1

[CNR*19] CORRÈA C. G., NUNES F. L. S., RANZINI E., NAKAMURA R., TORI R.: Haptic interaction for needle insertion training in medical applications: The state-of-the-art. *Medical Engineering and Physics* 63 (2019), 6–25. 1

[DDKA06] DURIEZ C., DUBOIS F., KHEDDAR A., ANDRIOT C.: Realistic haptic rendering of interacting deformable objects in virtual environments. *IEEE Transactions on Visualization and Computer Graphics* 12, 1 (apr 2006), 36–47. 1, 2

[DGM*09] DURIEZ C., GUÉBERT C., MARCHAL M., COTIN S., GRISONI L.: Interactive simulation of flexible needle insertions based on constraint models. *Lecture Notes in Computer Science 5762 LNCS, PART 2* (2009), 291–299. 1, 2

[Fel00] FELIPPA C.: A systematic approach to the element-independent corotational dynamics of finite elements. Tech. Rep. January, College Of Engineeringuniversity Of Colorado, 2000. 2

[KKB18] KUGELSTADT T., KOSCHIER D., BENDER J.: Fast Corotated FEM using Operator Splitting. *Computer Graphics Forum* 37, 8 (2018), 149–160. 1

[LFS*20] LI M., FERGUSON Z. A. H., SCHNEIDER T., LANGLOIS T., ZORIN D., PANOZZO D., JIANG C., KAUFMAN D. M.: Incremental potential contact: Intersection- and inversion-free, large-deformation dynamics. *ACM Transactions on Graphics* 39 (7 2020). 1

[MEM*19] MACKLIN M., ERLEBEN K., MÜLLER M., CHENTANEZ N., JESCHKE S., MAKOVIYCHUK V.: Non-smooth Newton methods for deformable multi-body dynamics. *ACM Transactions on Graphics* 38, 5 (oct 2019). 1

[OSO04] OKAMURA A. M., SIMONE C., O'LEARY M. D.: Force modeling for needle insertion into soft tissue. *IEEE Transactions on Biomedical Engineering* 51, 10 (2004), 1707–1716. 1

[RM17] RAVALI G., MANIVANNAN M.: Haptic Feedback in Needle Insertion Modeling and Simulation. *IEEE Reviews in Biomedical Engineering* 10, c (2017), 63–77. 1

[RVJM11] ROESTHUIS R. J., VEEN Y. R. V., JAHYA A., MISRA S.: Mechanics of needle-tissue interaction. *IEEE International Conference on Intelligent Robots and Systems* (2011), 2557–2563. 1

[SHZS*13] SUTHERLAND C., HASHTRUDI-ZAAD K., SELLENS R., ABOLMAESUMI P., MOUSAVI P.: An augmented reality haptic training simulator for spinal needle procedures. *IEEE Transactions on Biomedical Engineering* 60, 11 (2013), 3009–3018. 1

[TPRT00] TASCHEREAU R., POULIOT J., ROY J., TREMBLAY D.: Seed misplacement and stabilizing needles in transperineal permanent prostate implants. *Radiotherapy and Oncology* 55, 1 (apr 2000), 59–63. 1

[ZCC22] ZENG Z., COTIN S., COURTEUISSE H.: Real-time fe simulation for large-scale problems using precondition-based contact resolution and isolated dofs constraints. *Computer Graphics Forum* (2022). 1, 2, 3