1
Igor Peterlik
October 2015
Modelling of Living
Systems: Medical Images
and Physical Models
OUTLINE
• Motivation
– comparing images and models
• Biomechanical modeling
• Building a Model: Geometry
– segmentation and meshing
• Building a Model: Physics
– elastic formulation, parameters
– discretization, finite elements
• Example: Linear elasticity over P1 elements
• Numerical Solution
– direct and iterative solvers
– towards dynamics
2
FROM IMAGESTO MODELS
(AND BACK)
3
m
IMAGE ACQUISITION
4
MRI
US
CT
DIFFERENT MODALITIES
5
n the
this
rmae
3D
marks
Feao
the
en to
ature
The SURF detector provides a measure of reliability for each extracted
feature as the determinant of the Hessian matrix. We take
advantage of this measure to associate each reconstructed 3D point
with a quality q defined as the average reliability of its two corresponding
features in the left and right images. Quality values are
normalized so that q = 1 for the most reliable 3D point and q = 0
for the least reliable 3D point in y point set. We denote q the 3⇥m
vector formed by all q values, assuming the quality is isotropic at
each point. Note that other measures of reliability could be considered,
such as euclidean distances of the matched descriptors or the
eigenvalues of the covariance matrix on the reconstructed points.
We used two sorts of stereo endoscope, the first one consists of
two mounted endoscope from Karl Storz Endoscopy and the second
one is the stereo endoscope from the Da Vinci robot illustrated
in 4. Both camera lens generate radial distortion and have relative
small baselines (respectively 16mm and 6mm). The cameras
are calibrated following Zhang [40] approach and lens distortions
are rectified before performing the tracking. We take the assumption
that the stereo endoscope is fixed, which is the case when the
surgeon manipulate the surgical instruments.
Figure 4: Stereo endoscope used: on (left) a stereo endoscope of
the Da Vinci robot and (right) different views of the two mounted endoscope
from Karl Storz Endoscopy.
.
sentation of medical data
of scanners
MEDICAL IMAGING: PROS
•direct output of a scanning machine
• although already post-processed...
•doctors are used to look at images
• familiar representation
•huge source of information
• about geometry AND physics (elastography...)
•statistical evaluation
•(dis)-similarity
• visually-based
• stat/math-based
6
MEDICAL IMAGING: CONS
•noise
• loosing information due to modality
•noise
• loosing information due to motion
•noise
• loosing information due to post-processing
•a set of pixels/voxels
• no explicit physical meaning (although might
carry enough information about physical
parameters)
7
IMAGE PROCESSING
•filtering
• smoothing, denoising, edge-detection
• see Slicer3D (an open-source software)
•comparing: similarity metrics
• mean absolute differences
• summed squared differences
• normalized cross-correlation
• mutual information (different modalities)







8
MEDICAL IMAGE REGISTRATION
•minimizing dissimilarity
•many criteria
• inter/intra-patient
• single/multi-modal
• slice-slice, volume-volume,
slice-volume
• rigid, affine, deformable
• intensity/feature based
• model-based
• others...
•real-time in per-operational
scenario
9
target. One can observe that the algorithm captured the surface shift and the ventricular
thinning very accurately. The gray-level mean square difference between the target scan
and the deformed original scan on the image regions covered by the mesh went down
from 15 to 3. However, one can also notice that the left ventricle (lower one on the
Figure) was not able to fully capture the thinning. This is due to the approximate model
of the lateral ventricles we used in this experiment.
a) b) c) c)
Fig. 3. Slice 29 of a) initial scan b) target scan c)initial scan deformed using our algorithm
d)difference between target scan and deformed initial scan.
Figure 6 shows orthogonal cuts through the target intraoperative scan with transparently
overlayed color-coding of the intensity of the deformation field. The arrows show
the actual displacement of the nodes of the mesh. The extremely dense vector field in
the neighborhood of the lateral ventricles is due to the adaptive refinement of the mesh
at these locations.
Figure 5a shows the obtained deformation field overlayed on a slice of the initial
scan, and Figure 5b shows the same slice of the initial scan deformed with the obtained
deformation field. Several landmarks have also been placed on the initial scan
(green crosses) and deformed onto the target scan (red crosses), and these last landmarks
have also been overlayed on the target scan for comparison with the actual deformed
anatomy.
Similar landmarks as those shown on Figure 5 have been placed on 4 different slices
where the shift was most visible, and the distance between deformed landmarks and target
landmarks (not represented here for better visibility) have been measured. The surface
based landmarks on the deformed scan were within 1mm of the landmarks on the
target intraoperative scan. The errors between the landmarks placed in between the midsagittal
plane and the cortical surface were within 2-3mm from the actual landmarks.
The largest errors were observed at the level of the mid-sagittal plane and ventricles,
which can be explained by the fact that the surface matching of the ventricles was not
perfect. Nevertheless, the algorithm reduced the distance between landmarks in the iniFerrant
et al (2001): Registration of 3-D intraoperative MR images of the brain using a finite-element
biomechanical model.
(a) Source image (supine) (b) Warped image (c) Target image
(flank)
Fig. 5. Illustration of the accuracy of the registration for a cut in the source, warped
and target volume.
The supine and flank configurations are displayed on Fig. 6ab showing an
important deformation of the liver and surrounding tissues due to the impor-
MODELS
• A model is an abstract structure that uses mathematical language to describe
the behaviour of a system.
•typical examples of models of living systems:
•electrophysiological model: describes electrical properties of tissues
•e.g. electrophysiological model of heart
•model of fluid dynamics: describes behaviour of liquids
•e.g. cardiovascular fluid mechanics (blood circulation)
•biomechanical model of an organ: describes elastic (plastic) behaviour of tissues
•e.g. hyperelastic model of liver
•the mathematical language is usually based on differential equation
•since the behaviour usually means “a change of state”
10
BIOMECHANICAL MODEL: EXAMPLE
• medical image registration of volumes with important deformations
• two volumes taken at different configurations (pre-/intra- operational data)
• the goal is to align (register) the two images, i.e. match the voxels
11
• model based solution: an “energy-minimization” problem:
•an “error energy” given by difference between the two data (similarity
metric, difference in feature positions)
•an “elastic energy” given by a regularization term provided by an elastic model
12
BIOMECHANICAL MODEL: EXAMPLE
Vessel Tree in Body
• the Bezier tree can be immersed into the 3D tetrahe
IMAGE-MODEL COUPLING
13
ADVANCED MODELLING
14
ADVANCED MODELLING
15
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 

BUILDING A BIOMECHANICAL MODEL
FROMTHE IMAGE DATA
•two aspects: geometry (domain) and physics (formulation and parameters)
•the two aspects are closely interconnected
•geometry:
•type of the geometry structure is given by the nature of the problem and the
physical formulation (e.g. the basic “unit” is a tetrahedral element with 4 nodes)
• particular realisation is extracted from the image (e.g. the domain covered by the
elements is given by the shape of the organ)
• physics:
•formulation is given by a set of differential equations solved over the geometric
domain (e.g. finite element formulation of hyper-elasticity over linear tetrahedra)
• particular behaviour is determined by the physical parameters, usually obtained
by a measurement [invasive, non-invasive] (e.g. stiffness of the liver parenchyma)
16
BUILDING A BIOMECHANICAL MODEL:
GEOMETRY
17
GEOMETRY DISCRETIZATION: MESH
•usual geometric representation is given by a mesh (discretization of domain)
•a set of (connected) elements of given dimensionality and type
•1D: line mesh, beam mesh, spline mesh
•2D: triangular- and quad-mesh, shell mesh
•3D: tetrahedral mesh, hexahedral mesh
•mixed meshes
•used in computer-aided design (CAD) for decades
•many mesh generators from CADs
•commercial solutions (Ansys)
•open-source: GMsh,TetGen
18
MESHING OF MEDICAL IMAGES
• classical mesh generation from images consists of two steps
• segmentation: delimitation of the domain of interest in the image
• meshing: discretization of the segmented domain
19
IMAGE SEGMENTATION I
• manual segmentation (time consuming in 3D)
• semi-automatic methods:
– basic: histogram-based, edge detection, region-growing
– PDE-based: active contours (snakes, subject ofTP), level-sets methods
– graph-based segmentation: using graphs flows and cuts
20
IMAGE SEGMENTATION II
• atlas-based methods
– probabilistic methods (mean shape and possible variations)
• methods based on training
– neural-networks
• many open-source programs: ITKSnap, 3DSlicer,TurtleSeg
21
MESHING OF SEGMENTED DOMAIN I
• two-step approach:
– first step: generate surface representation (triangular mesh) of the segmented
domain (e.g. marching cubes)
– second step: generate 3D volume mesh from the surface mesh (e.g.TetGen
computing tetrahedral mesh from surface triangular mesh stored in STL)
– surface meshes can be very dense or with holes: reparation must be performed
before the second step (e.g. MeshLab)
• direct approach:
– direct generation of 3D volume mesh from the segmented domain: CGAL.org
– can be problematic for sharp features (usually not crucial in medical imaging)
and correct separation of boundaries (can be a problem, solution exists but is
not implemented…)
22
23
MESHING OF SEGMENTED DOMAIN II
BONUS:VARIATIONAL IMAGE MESHING
24
6 IEEE TRANSACTIONS ON MEDICAL IMAGING, TMI-2010-0299 PREPRINT
(a) (b) (c) (d) (e) (f)
Fig. 5. Mesh optimization on a 2D MR image slice (a) of brain ventricles. Initial (b) and optimized (c) discretizations with 59 nodes; initial (d) and
optimized (e) discretizations with 111 nodes. The finer optimized mesh is seen as overlaid on image (f).
(a) (b) (c) (d) (e)
Fig. 6. Mesh optimization on a 2D CT image slice (a) of the kidney. Initial (b) and optimized (c) discretizations with 61 nodes; and initial (d) and optimized (e)
discretizations with 338 nodes.
GOKSEL AND SALCUDEAN: IMAGE-BASED VARIATIONAL MESHING 9
modulus distribution, E(x), within the domain M is known a
priori. There exist several methods in the literature for the
acquisition and derivation of tissue elasticity, see [39] for
a review. Consequently, the material stiffness matrix in the
element can be formulated as C(x) = E(x) Cj. Substituting
this in (10) yields:
Ej
strain(uj
) =
1
2
ujT
BjT
CjBj
uj
j
E(x)dx . (14)
For the discretization within each element to be optimal, these
two energy formulations in (13) and (14) should be equal,
leading to:
˜Ejvj
=
j
E(x)dx (15)
which is satisfied when ˜Ej is the mean of the distribution
within element ⇥j
.
To demonstrate mesh optimization from mechanical tissue
features, the method was applied to prostate elastography
images acquired using the vibro-elastography technique of
Salcudean et al. [24]. Elastography is the technique in which
tracked localized displacements in response to a mechanical
excitation allow for the identification of mechanical tissue
properties [39], [40]. For the purpose of this paper, a 2D
sagittal transfer function image of the prostate is meshed. The
prostate, which is typically stiffer that its surrounding, is seen
in Fig. 11. The optimized meshes and their corresponding
discretizations are also shown in this figure.
D. Connectivity and Node Updates
For the case where the objective function is purely geometric,
i.e. = 0, EG has a simple algebraic (quadratic)
definition as in (3). This can also be observed for relatively
small values of as in Fig. 3(b). For EG alone, a geometrical
closed-form expression of its critical-point within each 1(a)
(b)
Fig. 11. Meshing of a transversal (a) and a para-sagittal (b) slice from
prostate vibro-elastography. The prostate is the darker oval structure in the
center.
is found using numerical integration over the voxels, during
which the enclosure of voxels by elements is determined using
• Orcun Goksel and Septimiu E.
Salcudean, "Image-Based Variational
Meshing", IEEE Trans Medical Imaging
30(1):11-21, Jan 2011.
• direct generation of meshes from the image
– no segmentation needed
– initial regular mesh is adapted to the image
– works for limited range of intensities
DISCRETISATION METHODS
•wide range of algorithms
• Tesselations (tiling of a plane)
• Delaunay triangulations [DT]
(no point of the triangulation
lies inside any circumcircle of
any triangle of the triangulation
• Voronoi diagrams (dual to DT)
25
DISCRETIZATION QUALITY
•quality of elements is a crucial in physics-based applications (vertex Jacobian)
• degenerated elements result in numerical instability (singularity of the
Jacobian)
•various measures of element quality:
• smallest angle/largest angle (2D)
• dihedral angle (3D)
• determinant of vertex Jacobian
• ratio of inscribed/circumscribed radii
• others (edge ratio, Frobenius aspect etc)
26
defined an e⇥ciency index ⇤ [13] as
⇤ = exp
1
ne
neX
e=1
⇤e
!
with ⇤e = le 1 if le < 1 and ⇤e =
1
le
1 if le ⇥ 1. The e⇥ciency inde
and should be as close as possible to ⇤ = 1.
For measuring the quality of elements, various element shape measure
literature [33, 27]. Here, we choose a measure based on the element rad
between the inscribed and the circumcircles.
If K is a triangle, we have the following formula
K = 4
sin ˆa sinˆb sin ˆc
sin ˆa + sinˆb + sin ˆc
,
ˆa, ˆb and ˆc being the three inner angles of the triangle. With this defin
triangle has a K = 1 and degenerated (zero surface) triangles have a K
For a tetrahedron, we have the following formula:
K =
6 6 Vk
4X
i=1
a(fi)
!
max
i=1,...,6
l(ei)
,
with VK the volume of K, a(fi) the area of the ith
face of K and l(ei) the
C. GEUZAINE, J.F. REMACLE
aluate the adequation between the mesh and the prescribed mesh size field, we
iency index ⇤ [13] as
⇤ = exp
1
ne
neX
e=1
⇤e
!
(2)
1 if le < 1 and ⇤e =
1
le
1 if le ⇥ 1. The e⇥ciency index ranges in ⇤ ⇤ [0, 1]
as close as possible to ⇤ = 1.
g the quality of elements, various element shape measures are available in the
27]. Here, we choose a measure based on the element radii ratio, i.e. the ratio
cribed and the circumcircles.
ngle, we have the following formula
K = 4
sin ˆa sinˆb sin ˆc
sin ˆa + sinˆb + sin ˆc
,
g the three inner angles of the triangle. With this definition, the equilateral
K = 1 and degenerated (zero surface) triangles have a K = 0.