PHYSICS-BASED
M
ODELING


AND
SIM
ULATIONS
IN
M
EDICINE
PV177
Sem
inary Group
Igor Peterlik
October 2015
OUTLINE
§Motivation for physics-based simulations

‣ training, planning, navigation
§Examples of models and simulations

‣cataract surgery: training ophthalmologists
‣cryoablation planning: pre-operative tool for interventional radiologists
‣augmented reality framework for hepatic laparoscopy
§Basic concepts of modeling

‣Images vs. Models: reconstruction of models from images
‣conceptual, mathematical, physics-based models
2
M
OTIVATION
From
Training to
Navigation
COMPUTER-BASED MEDICAL SIMULATION
§Main areas of interest 

‣procedural training: practical and ethical considerations
‣pre-operative planning and rehearsal
‣per-operative guidance
§Different requirements on each level

‣Increasing levels of complexity as we get closer to the operation room
4
2010 2014 2018
Procedural Training Pre-operative planning Intra-operative guidance
MAIN CHARACTERISTICS
§procedural training

‣interventions in eye surgery, catheter
‣realistic, interactive (visual and haptic rendering), generic models
§pre-operative planning

‣liver, kidney resection, deep-brain surgery
‣realistic, not necessarily interactive, patient-specific models
§intra-operative navigation

‣catheter, needle insertion navigation, laparoscopic augmented reality
‣realistic, interactive, robust, patient-specific
5
NOTE: HAPTIC DEVICE
6
‣“3D mouse” with force feedback
‣allows for touching (haptein)
virtual objects
‣often necessary for training as
visual perception is not sufficient
(e.g. cutting of tissue)
‣main issue: high refresh rate needed to guarantee the fidelity of rendering
‣usually 1000 Hz is reported (although the required minimal frequency
rather depends on the mechanical properties of objects being rendered
‣other issues: stability, passivity (might depend on the quality of device)
FROM TRAINING...
7
Input /

Haptic deviceGeneric Model
Visual Feedback
8
Patient Specific Data
Tracking
device
Visual Feedback
Augmented view
...TO INTRA-OPERATIVE ASSISTANCE
Per-operative imaging
Robotic

device
EXAM
PLES
From
Training to
Navigation
TRAINING: CATARACT SURGERY
§metabolic changes crystalline lens fibers (loss of sight)

§several types of surgery

‣phacoemulsification: standard in developed countries (ultrasonic)
‣extra-/intra-capsular cataract extraction (ECCE, ICCE)
‣manual small incision cataract surgery (MSICS)
‣lens is extracted through a tunnel which is watertight (if created properly)
‣lens capsule is intact
‣outcomes comparable to phacoemulsification
‣much lower cost ($50 vs. $2500) and time (5 to 15 minutes)
‣requires very high dexterity of the pharmacologist
10
HELP ME SEE PROJECT
§large impact of cataract in the third world

‣estimated 20 million children, 100 million adults blind
§mission:

‣train a large number of ophthalmologists (30,000 in 2030)
‣use a virtual simulator for training
‣shifting paradigm
§large call for project

‣first prototypes tested by skilled ophthalmologists
‣consortium (InSimo, SenseGraphics, MOOG)
11
MSICS SIMULATOR PROTOTYPE
12
MSICS SIMULATOR PROTOTYPE
13
PLANNING: NEEDLE INSERTION FOR CRYOABLATION
§interventional radiology

‣destruction of a tumor using a (steerable) hollow needle
‣argon is used to freeze the tumor by forming an ice-ball (ice-rod)
‣insertion/placement of the needle plays a crucial role
‣avoid important objects during insertion (vessels)
‣insert the needle so that the ice-ball covers the tumor (+safety margin)
‣multiple needles inserted (synergy effect)
§actually, two stages are studied

‣needle insertion planning
‣prediction of ice-ball formation
14
INSERTION PLANNING
15
INSERTION PLANNING
16
ITERATIVE PROCESS
17
ICE-BALL FORMATION PREDICTION
18 H.Talbot et al. Interactive Planning of Cryotherapy Using Physically-Based
Simulation Proc. Medicine Meets Virtual Reality, 2014

19
NAVIGATION: AUGMENTED REALITY IN
LAPAROSCOPY
§laparoscopy: minimally invasive approach (keyhole surgery)

‣operation through small incisions
‣surgeon follows the intervention through camera (mono/stereo)
§pre-operative data available 

‣e.g. pre-operative abdominal CT
‣however, the actual position during surgery is often different (e.g.
supine vs. flank vs. prone position)
‣huge deformation occurs in abdominal cavity (mainly in patient with
higher body mass index)
‣surgeon has to create a mental image
LAPAROSCOPIC HEPATECTOMY
20
21
AUGMENTED REALITY I
Best 

paper
§track the image acquired by the camera (Computer Vision)

§drive the model (from CT data) during the interaction
N. Haouchine, J. Dequidt, I.P., E. Kerrien, M.-O. Berger, S. Cotin.
Image-guided Simulation of Heterogeneous Tissue Deformation For
Augmented Reality during Hepatic Surgery. In ISMAR proc. 2013 

the least reliable 3D point in y point set. We denote q the 3 m
tor formed by all q values, assuming the quality is isotropic at
h point. Note that other measures of reliability could be considd,
such as euclidean distances of the matched descriptors or the
envalues of the covariance matrix on the reconstructed points.
We used two sorts of stereo endoscope, the first one consists of
o mounted endoscope from Karl Storz Endoscopy and the secd
one is the stereo endoscope from the Da Vinci robot illustrated
4. Both camera lens generate radial distortion and have relae
small baselines (respectively 16mm and 6mm). The cameras
calibrated following Zhang [40] approach and lens distortions
rectified before performing the tracking. We take the assumpn
that the stereo endoscope is fixed, which is the case when the
geon manipulate the surgical instruments.
ure 4: Stereo endoscope used: on (left) a stereo endoscope of
Da Vinci robot and (right) different views of the two mounted enscope
from Karl Storz Endoscopy.
BIOMECHANICAL MODEL
his section we provide a description of the biomechanical model
d to compute the deformations of the liver. Before giving details
etrahedral model employed for parenchyma, we focus on model
vascularization and mechanical coupling between this two. Since
final composite model is heterogeneous and anisotropic due
he vascular structures, we finally describe the solution process
ed on a direct solver, still allowing for real-time performance.
correspondences between stereo images are obtained with a nearest
neighbour criterion on feature descriptors coupled with the epipolar
constraint to filter out outliers. After triangulation, we obtain
a sparse set of m 3D points, denoted by 3 m coordinate vector
y. Examples of point correspondences and reconstructed 3D points
from laparoscopic images are shown in Fig. 3.
as fP = JT fvi where J is a 3NP ⇤6NV Jacobian matrix of the mapping
between the nodes of parenchyma and vessels [7]. The global
stiffness K matrix is then computed as K = KP +J⌅KVJ.
Concerning the tumor, since this work focus only on small
tumors, we assume that its influence on the overall mechanical
behaviour is negligeable and therefore the coupling with the
parenchyma is only geometric. However the tumor can easily be
modelled as a real mechanical model with different properties from
those of the parenchyma, which will not affect the performance of
our method.
4.4 Numerical Solution
In this paper, a quasi-static scenario is considered, i. e. the actual
shape of deformable object under applied forces is computed using
the finite element formulation without dealing with the dynamic
properties of the tissue. Therefore, in each step of the simulation,
a linear system given by Eq. 1 is resolved. A wide range of direct
and iterative solvers has been proposed in the past to solve such
a system of equations emerging in the physics-based modeling of
deformable bodies. In case of homogeneous systems in which the
finite element formulation results in well-conditioned matrices, iterative
solvers such as conjugate gradients have proven to be efficient
techniques converging rapidly to the optimal solution. However,
in our case, the final matrix K gathers mechanical contributions of
both the parenchyma and vessel walls. As the experiments reports
a significant difference in stiffness of these two components (e. g.
see [38]), the composite system results in poorly conditioned matrix.
In this case the convergence of iterative solvers becomes an
important issue.
For this reason, we rely on direct LDL solver which requires
an explicit factorization of the system matrix. Although the solver
imposes more strict limitations on the size of the system being resolved,
it still provides a stable real-time solution applicable to the
problems considered in the scope of this paper.
5 NON-RIGID REGISTRATION
5.1 Initial registration
A correct initial registration is an important step in the framework
since it can significantly impact the estimated position of the tumor.
For that purpose, special care must be taken during the initialization
phase. Similar to related works [36], our initialisation is done manually
through a Graphical User Interface following these steps:
• the real camera parameters acquired from camera calibration
are loaded on the virtual camera.
• the 3D model of the liver is manually aligned on the first
pair of laparoscopic images based on salient geometrical landthe
model and the 3D features extracted from the stereoscopic images.
Only features that intersect the liver surface after ray-casting
are kept, the features that do not belong to the liver are filtered out
from laparoscopic images.
5.2 Non-rigid Registration
The non-rigid registration can be seen as an optimization problem
between the three dimensional features recovered from laparoscopic
images representing the liver and the biomechanical model
derived from preoperative CT data. In this paper we propose to
perform this registration as an error-minimization problem that accounts
for the internal energy of the biomechanical Ei and the tracking
energy Et. Deriving this energy shows that an extremum (minimum)
is reached when the internal forces equal the tracking forces.
Thus, the internal forces are expressed as:
fi(x) = ReK·(Re
⌅
x x0
) (2)
where K represents the stiffness matrix, Re the co-rotational matrix
and x,x0 are vectors of size 3n representing the position of the n
degrees of freedom of the mechanical model, respectively at any
time t and time t = 0.
Figure 6: Non-rigid registration: one frame of the non-rigid registration
where red spheres represent the projected features y0 from the
initialisation step, green spheres represent the tracked features y and
green lines represent the forces linking both point sets.
We propose to handle the non-rigid registration by adding external
stretching forces induced by the tracking step. The tracked
3D features y represent how visible parts of the liver are moving.
Due to noise or other reconstruction inaccuracies, point locations
in y are often regularized by generating a smooth and dense surface
model that is later used as boundary conditions for the biomechanical
model [34]. Instead, we propose here to incorporate such
or two order of magnitude higher than
to impose the boundary conditions disimulation
remains stable even at large
ULTS
dation remains very challenging. In our
omplex since neither qualitative results
idate the deformation of internal strucaparoscopic
images. In order to asses
oach, we compared our simulation rental
scenarios. First, we demonstrate
a whether and where an heterogeneous
geneous one for the prediction of tumor
n a realistic phantom liver to quantitaween
the simulation and a ground truth.
d on an actual laparoscopic procedure
, allowing us to qualitatively estimate
orm in a real surgical environment.
(b)
6.1 Experimentation with computer-generated data
We evaluate the impact, on the registration, of using a heterogeneous
model instead of a homogeneous model (which can be seen
as providing similar results as an advanced geometric approach).
This is done by calculating the euclidean distance between the estimated
tumor location in the cases of homogeneous and heterogeneous
deformations (see 7). We also measure this influence depending
on the location of the tumor in the liver, at three different
locations: 1) close to the the point of interaction in order to quantify
local deformation, 2) away from the point of interaction to quantify
global behaviour, and 3) in the middle of the vascular network to
assess its influence. The simulations are generated using the SOFA
framework [2].
Results illustrated in figure 8 show that taking into account the
vascular network impact the tumor location. We can notice a difference
in distance of about 15 mm in the case where the tumor is
located close to the deformation. We also notice that even if the
tumor is located far from the point of interaction, it remains influenced
by the vascular network with a distance of more than 3 mm.
However, when the tumor is very close to the boundary conditions,
the impact of the vascular network is considerably reduced, which
is an expected result.
sults against three experimental scenarios. First, we demonstrate
with computer-generated data whether and where an heterogeneous
model differs from an homogeneous one for the prediction of tumor
location. Second, we rely on a realistic phantom liver to quantitatively
measure the error between the simulation and a ground truth.
Third, our approach is tested on an actual laparoscopic procedure
performed on a human liver, allowing us to qualitatively estimate
how our approach could perform in a real surgical environment.
(a) (b)
(c) (d)
framework [2].
Results illustrated in figure 8 show that taking into account th
vascular network impact the tumor location. We can notice a di
ference in distance of about 15 mm in the case where the tumor
located close to the deformation. We also notice that even if th
tumor is located far from the point of interaction, it remains influ
enced by the vascular network with a distance of more than 3 mm
However, when the tumor is very close to the boundary condition
the impact of the vascular network is considerably reduced, whic
is an expected result.
Figure 8: Impact of the vascular network on the tumor deformatio
depending on its position in the liver: the distance between the tu
mor using homogeneous and heterogeneous biomechanical mod
is important locally (red) and globally (blue) and less important whe
the tumor is constrained by the vessels. The meshes illustrate th
distance between the position of the tumor in a homogeneous and
heterogeneous case for each location in the liver.
6.2 Experimentation with liver phantom data
We believe that performing a CT scan of the a phantom liver befor
and after a deformation is an ideal way of defining a ground trut
(a) (b)
(c) (d)
(e) (f)
Figure 7: Computer-generated data with (left) the liver at rest and
(right) the liver after deformation: in (a) and (b) the volumetric mesh
composed of tetrahedra, in (c) and (d) the beams generated along
the vessels described in Section 4 an in (e) and (f) the heterogenous
liver including the vascular network in wireframe.
Figure 8: Impact of the vascu
depending on its position in th
mor using homogeneous and
is important locally (red) and g
the tumor is constrained by th
distance between the position
heterogeneous case for each l
6.2 Experimentation wit
We believe that performing a C
and after a deformation is an
for the location of an interna
surgical instruments as well
large image artifacts in an in
conditions on the liver are diffi
fluence the organ motion and
a validation protocol using a
prior knowledge of its mecha
geometry of the liver and vasc
tual patient liver (but with a 1
is placed outside the scanned
using a wire attached to the ph
quisitions.
The scenario involved the
22
AUGMENTED REALITY II
Best 

paper
N. Haouchine, J. Dequidt, I.P., E. Kerrien, M.-O. Berger, S. Cotin.
Image-guided Simulation of Heterogeneous Tissue Deformation For
Augmented Reality during Hepatic Surgery. In ISMAR proc. 2013 

TRAINING, PLANNING, NAVIGATION: COMPARISON
§medical training

‣example: cataract surgery
‣realistic, interactive (visual and haptic rendering), generic models
§pre-operative planning

‣example: needle insertion planning in cryoablation
‣realistic, not necessarily interactive, patient-specific models
§intra-operative navigation

‣example: augmented reality for laparoscopic surgery
‣realistic, interactive, robust, patient-specific
23
Physics-based
Performance
Specificity
24
§Starting point: simulation for training 

§Ongoing transition towards

‣pre-operative planning of procedures
‣intra-operative navigation/guidance
2010 2014 2018
Procedural Training Pre-operative planning Per-operative guidance
Planning of needle insertion Augmented reality for liver surgery
complexity
Physics-based
Performance
Specificity
TRAINING, PLANNING, NAVIGATION: SUMMARY
SYNERGY
25
§Parameter identification of simulation models 

§Patient-specific modeling for real-time simulation

§Medical robotics

§Continue working on all objectives as they influence each other

§And transfer as many things as possible
planning
guidance
training
clinical
products
R&D
SOFASIM
26
COURSE OUTLINE
‣ 23/10: Motivation, context, examples. Images vs. Models.
‣ 30/10: Geometry: creating a mesh. Gmsh, CGAL, Paraview, SOFA.
‣ 06/11+?: Kinematics, kinetics, linear elasticity. Finite element method.
‣ 20/11: Modeling a simple quasi-static deformable object. First simulation in
SOFA: linear solvers, mappings, rendering.
‣ 27/11: Including non-linearities: co-rotational and hyper-elastic models. Nonlinear
solvers, convergence.
‣ 04/12: Dynamics: explicit vs. implicit time integration methods.
‣ 11/12: Advanced topics: contacts, interaction, visual and haptic real-time.
‣ 18/12: Discussion, perspectives, practicals…