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 ﬁnite element modelingand allows for a fast update of a large number of tetrahedral elements.The speed-up is obtained by the use of a speciﬁc 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 ﬁrst 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
phacoemulsiﬁcation
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 ﬁrst steps toward asimulator that could help to solve this training bottleneck.
30
Unlike during phacoemulsiﬁcation, 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 stiﬀer 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 stiﬀer than at the periphery (blue).
Needs and contributions:
In this context of tradeoﬀ between accuracy andcomputation time, we need: (1) to account for precise contact response in highly
40
constrained cases and for strong stiﬀness inhomogeneities. We note that theseinhomogeneities create ill-conditioned systems that are more diﬃcult to solve.(2) to use a large number of nodes and elements in the ﬁnite 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 stiﬀnesses 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
stiﬀness
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
deﬁnite, a popular algorithm to eﬃciently solve this problem is the Conjugate
Preconditioner-based contact response and application to cataract surgery 3
Gradient (CG) iterative solver. However, it suﬀers 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: ﬁrst 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 stiﬀness fac-
75
tor. This stiﬀness is diﬃcult to determine, as it must be tuned accordingly to thestiﬀness of the object in contact. This becomes particularly an issue when deal-ing with inhomogeneous objects. Stiﬀ 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 stiﬀness, 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 conﬁguration 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 inﬂu-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 signiﬁcantly 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 ﬁeld 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 stiﬀ kernel and softer boundaries, whilethe meshes of these structures contain elements of potentially very diﬀerent 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 conﬁguration 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 diﬀerence would be reduced after several time-steps, using
125
a more accurate compliance produces better intermediate conﬁgurations 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. Stiﬀer 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 conﬁguration 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 ﬁrst 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 ﬁnal 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

Search Related

Positive psychology and its application to heICTs Application to Language Teaching and LeaPositive psychology and its application to heapplication of GIS and RS to town planning isPilgrimage and travel to the Holy LandCombating desertification and adaptation to CTask-based language teaching and learningenvironmental geology and application of eartApplication to Data AnalysisITproject and application estimationand meas

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