Preconditioner-based contact response and application to cataract surgery

Please download to get full document.

View again

of 14
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
Document Description
In this paper we introduce a new method to compute, in real-time, the physical behavior of several colliding soft-tissues in a surgical simulation. The numerical approach is based on finite element modeling and allows for a fast update of a large
Document Share
Document Tags
Document Transcript
  Preconditioner-Based Contact Response andApplication to Cataract Surgery Hadrien Courtecuisse, J´er´emie Allard, Christian Duriez, St´ephane Cotin Shaman Group, INRIA Lille North Europe Abstract.  In this paper we introduce a new method to compute, in real- 5 time, the physical behavior of several colliding soft-tissues in a surgicalsimulation. The numerical approach is based on finite element modelingand allows for a fast update of a large number of tetrahedral elements.The speed-up is obtained by the use of a specific preconditioner thatis updated at low frequency. The preconditioning enables an optimized 10 computation of both large deformations and precise contact response.Moreover, homogeneous and inhomogeneous tissues are simulated withthe same accuracy. Finally, we illustrate our method in a subpart of asimulation of cataract surgery, which require to handle contacts with nonhomogeneous objects precisely. 15 1 Introduction Several methods have been presented for real-time bio-mechanical simulations of contacting soft-tissues. However, this work is motivated by a training simulatorfor cataract surgery, whose needs can not be addressed with existing methods.Thus, we first present the context and the motivations for this work before 20 highlighting the simulation needs and the contributions of our method. Context:  In developed countries, when the crystalline lens is being clouded overby a cataract, a therapy based on  phacoemulsification   is commonly practiced.Training simulation of this procedure has been studied [ ? , ? ], and a commercialsolution is provided by VRmagic. Another surgical procedure known as Manual 25 Small Incision Cataract Surgery (MSICS), requires low technology with almostequivalent results when performed by a well-trained specialist. A recent report[ ? ] shows that there is a great need of training for this surgery, to face the hugenumber of people that need care. This work provides the first steps toward asimulator that could help to solve this training bottleneck. 30 Unlike during phacoemulsification, MSICS consists in pulling out the crys-talline lens as a single piece through a small incision (see Fig. 1). This step gen-erates large deformations on both lens and eyeball. Moreover, for older patients,the center of the lens can be much stiffer than the periphery. This inhomogene-ity must be taken into account, as it impacts the successful completion of the 35 surgery. Indeed, it changes the mechanical behaviors and may require adaptingthe size of the incision. A training system for MSICS necessitates an accuratemodeling of these soft-tissues in a simulation executed in real-time.  2 Hadrien Courtecuisse, J´er´emie Allard, Christian Duriez, St´ephane Cotin IncisionIrrigating vectis LensIris   Scleral lip depression (a) Extraction of the nucleous (b) Simulated procedure Fig.1.  Nucleus extraction with an irrigating vectis. (a) motion of the tool, (b) simu- lated lens. Elements at the center (in red) are stiffer than at the periphery (blue). Needs and contributions:  In this context of tradeoff between accuracy andcomputation time, we need: (1) to account for precise contact response in highly 40 constrained cases and for strong stiffness inhomogeneities. We note that theseinhomogeneities create ill-conditioned systems that are more difficult to solve.(2) to use a large number of nodes and elements in the finite element formulation,in order to maintain the realism of the deformations.The contribution of this paper is to answer these needs (which probably exist 45 for other kind of surgery) in a real-time simulation. The method is based on animplicit integration scheme (to manage the large range of stiffnesses while keep-ing large time steps) and uses a preconditioner that is updated at low frequency.This preconditioner provides almost an exact solution to the deformation solver,allowing its fast convergence. Moreover, it provides a very good estimation of  50 the mechanical coupling between the contacting nodes. This allows for a precisesolving process based on accurate contact laws. 2 Related Works Many previous works involve the simulation of  soft body in contacts. We will 55 concentrate here on aspects related to the proposed strategies to solve the me-chanical and contact equations once they have been expressed.The core mechanical equations are based on Newton’s second law and werely on an linearized implicit integration scheme [ ? ], which is more stable whensudden contacts occur. Thus we update the next positions and velocities withthe following equations : Ma t + h  =  f  ( x t + h , v t + h )  v t + h  = v t  +  h a t + h  x t + h  =  x t  +  h v t + h  (1)  M +  h B +  h 2 K      A a t + h  = f  ( x t , v t ) +  h Kv t      b (2)where  M  is the  mass   matrix,  B  the  damping   matrix and  K  the  stiffness   matrix.Eq. (2) must be solved at each time-step, but continuously changes due tomaterial and geometrical non-linearities. As the system matrix is positive semi- 60 definite, a popular algorithm to efficiently solve this problem is the Conjugate  Preconditioner-based contact response and application to cataract surgery 3 Gradient (CG) iterative solver. However, it suffers from convergence issues for ill-conditioned matrices, which can appear for inhomogeneous materials or mesheswith varying element sizes. Preconditioning techniques consists in computing anapproximation P of matrix A which is easier to invert. Then, this preconditioner 65 can be used to solve (2) by relying on  P − 1 A , which is better conditioned andthus converge to an adequate solution in fewer iterations. However, computinga preconditioner adds two overheads: first to invert (or factorize) the precon-ditioner itself, then to apply it at each CG iteration. Several preconditionerscan be used, from simple diagonal matrices [ ? ] to costly but precise incomplete 70 Cholesky factorizations.When contacts occur, they must be taken into account in the above equationsystem. A common solution consists in using a penalty method, which handlescontacts by adding a contact force  f   =  kδ  n  at each contact point, where  δ   is ameasure of the interpenetration, n is the contact normal and  k  is a stiffness fac- 75 tor. This stiffness is difficult to determine, as it must be tuned accordingly to thestiffness of the object in contact. This becomes particularly an issue when deal-ing with inhomogeneous objects. Stiff penalty forces can either be considered asexplicitly-integrated external forces, which can introduce instabilities and wouldrequire lowering the time-steps. 80 Other methods rely on Lagrange multiplier [ ? , ? , ? ] to compute contact forcessuch that all intersections are removed at the end of each time step. The corecomputation of the algorithm involves solving a Linear Complementary Problem(LCP) deriving from Signorini’s law : δ  = HCH T  λ  +  δ 0  with 0 ≤ δ ⊥ λ ≥ 0 (3)where  H  is the jacobian of the contacts, and  C  is the  compliance   matrix,  λ  isthe contact force,  δ 0  and  δ  are the measure of interpenetrations before and aftercollision response. This equation means that, if   λ  is positive at the end of thetime step, then  δ   must be equal to zero, and vice versa.For explicits methods C is the inverse of the (often diagonalized) mass ma- 85 trix, whereas for implicit schemes it also involves damping and stiffness, i.e. C  = (  1 h 2 M +  1 h B + K ) − 1 for the scheme used in (2). This matrix changes atevery time steps, and its computation can be prohibitive for large deformablemeshes. It is possible to use an approximation of the compliance, because anapproximated local deformation in response to a collision is visually acceptable. 90 The most extreme approximation is to only consider the mass, as for explicitschemes. When it is used, contacts are corrected without any mechanical couplingbetween nodes. A more accurate method has been proposed in [ ? ], inspired bythe co-rotational formulation used in FEM to remove non-linearities introducedby rotations. C can be precomputed from the rest configuration and updated at 95 each time step based on a local estimation of the rotations. However, this ap-proximation can become inaccurate for large deformations. Moreover, it requiresstoring a dense matrix, with 9 n 2 values where  n  is the number of vertices of theobject, therefore preventing its application on detailed meshes.  4 Hadrien Courtecuisse, J´er´emie Allard, Christian Duriez, St´ephane Cotin 3 Inhomogeneous Deformable Model on GPU 100 Soft tissue deformations are modeled using a geometrically non-linear elasticformulation and a Finite Element Method (co-rotational model introduced in[ ? ]). The local rotation of each element is estimated in order remove the influ-ence of geometrical nonlinearities, allowing to account for large deformations. Tocompute the deformation in real-time with detailed meshes (more than 10,000 105 tetrahedral elements), we implemented a CG solver on the GPU [ ? ].Recently, we introduced a new approach [ ? ] to maintain a good approxima-tion of the inverse of the system matrix during the simulation. This methodexploits time coherency to update an exact factorization of the system matrixwithin a separated loop computed at a lower frequency. It only requires a few 110 simulation steps to provide a new factorization, thus constantly providing a goodapproximation usable as a preconditioner. It is further improved by estimatingrotations between the current position and the state used for the last factor-ization, similar to the co-rotational formulation. Even if the preconditioner isapplied on the CPU, the GPU can still be used for the remaining operations. 115 Finally, as it uses a sparse factorization, it supports simulations with a largenumber of elements.As a consequence the above method significantly reduces the number of iter-ations needed for the CG algorithm to converge. However, the paper only con-sidered simulations involving a single object without any contacts. In this paper,we extend this idea to the computation of contact responses, as is detailed in thenext section, and apply this new technique in the field of medical simulation. Inthe context of the cataract simulation, two deformable tissues need to be sim-ulated: the eye lens and sclera (Fig. 1) which both undergo large deformationscombined with multiple contacts during surgery. Furthermore, the crystallinelens is by nature inhomogeneous, with a stiff kernel and softer boundaries, whilethe meshes of these structures contain elements of potentially very different size(in particular near the incision located in the sclera). All of these characteristicslead to a very poor conditioning of the system matrix, and it is obvious that inthis case, preconditioning techniques could greatly improve the converge rate of the solver. For this application, we found the  LDL T factorization to be morestable than the LL T Cholesky factorization. Therefore, the complete expressionof the preconditioner we use is: P = R T t → t − ∆t  L t − ∆t D t − ∆t L Tt − ∆t  R t → t − ∆t  (4)where  R t → t − ∆t  is the current local rotations matrix, and  L t − ∆t D t − ∆t L Tt − ∆t the most recent factorization. 4 Preconditioner for Contact Response 120 Local mechanical coupling must be taken into account to compute accurate colli-sion response (see Fig 2). Indeed, while the constraint-based formulation will en-sure a contact-free configuration at the end of each time-step in all cases, the local  Preconditioner-based contact response and application to cataract surgery 5 deformation near contacts obtained without the correct compliance (Fig. 2(b))is invalid. While this difference would be reduced after several time-steps, using 125 a more accurate compliance produces better intermediate configurations and,more importantly for real-time simulations, allow larger time-steps. (a) homogeneous object (b) inhomogeneous objectwithout coupling(c) inhomogeneous objectwith coupling Fig.2.  Force repartition for homogeneous and inhomogeneous objects on which colli-sions are solved with and without mechanical coupling. Stiffer part are shown in red. To solve this issue, we propose to extend the work introduced in [ ? ], usingthe preconditioner as compliance matrix instead of a precomputed dense inversematrix. Indeed, the compliance matrix is only a scaled version of the inversesystem matrix: C =  h 2 A − 1 . Moreover, we showed in section 3 that we maintaina good approximation of  A − 1 by updating the preconditioner at low frequency.Thus, if  P remains a good approximation of the system, it can be re-used to buildthe compliance matrix. That way, we guarantee a contact-free configuration atthe end of the time step with almost the exact mechanical coupling of elementstaken into account. To build the LCP matrix, we combine eq. (3) and (4) to compute for each object: HCH T = h 2 HA − 1 H T ≈  H  RLDL T R T  − 1 H T ≈ h 2 J  LDL T  − 1 J T with  J = HR (5)This computation can be implemented in three steps. The first step consistsin applying the local rotations since the last update to  H . This operation isinexpensive because it is done by computing the product of a block-diagonal 130 matrix  R  with a sparse matrix into matrix  J  with the same sparsity structure.Next, we compute the product of  J T with the inverse of the factorization. This isachieved by computing rows independently, each requiring two triangular solveswith one row of   J . Finally, the resulting matrix is multiplied by matrix  J  toobtain the final contribution to the LCP matrix. 135 The number of solves in step 2 is proportional to the number of contacts.They can therefore become expensive in complex scenarios. However, each solveis independent. Thus, we propose to implement them on GPU, by computingseveral solves in parallel, reducing the level of parallelism that must be achievedwithin each solve. This maps nicely to the two-level SIMD architecture of to- 140 days GPU where synchronizations within a group of cores is fast, whereas global
Similar documents
View more...
Search Related
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks