A practical framework for generating volumetric meshes of subject-specific soft tissue [Peeters2013].pdf

(2649 KB) Pobierz
The Visual Computer manuscript No.
(will be inserted by the editor)
A practical framework for generating volumetric meshes of
subject-specific soft tissue
Pieter Peeters
·
Nicolas Pronost
Abstract
Studying human motion using musculoskele-
tal models is a common practice in the field of biome-
chanics. By using such models, recorded subject’s mo-
tions can be analyzed in successive steps from kinemat-
ics and dynamics to muscle control. However simulating
muscle deformation and interaction is not possible, but
other methods such as a
finite element
(FE) simulation
are very well suited to simulate deformation and inter-
action of objects. In this paper we present a practical
framework for the automatic generation of FE ready
meshes based on subject-specific segmented MRI data.
The proposed method resolves several types of data
inconsistencies: noise, an incomplete dataset and self-
intersections. This paper shows the different steps of
the method, such as solving overlaps in the segmented
surfaces, generating the volume mesh and the connec-
tion to a musculoskeletal simulation.
Keywords
Virtual human
·
Finite element simula-
tion
·
Musculoskeletal simulation
·
Subject-specific
modeling
·
Computational geometry
1 Introduction
The simulation of human motion using anatomically-
based musculoskeletal models is of interest in many
P.W.A.M. Peeters
Utrecht University (The Netherlands), Games and Virtual
Worlds
N. Pronost
Utrecht University (The Netherlands), Games and Virtual
Worlds
Tel.: +31-302537578
Fax: +31-302534619
E-mail: nicolas.pronost@uu.nl
fields including computer graphics, biomechanics, mo-
tion analysis, medical research and virtual character an-
imation. Much research has been done to understand
the musculoskeletal system, and so there is a large am-
ount of data [1] describing the mechanics of muscle,
the geometric relationships between muscles and bones,
and the motions of joints. In the medical field, the
neuromuscular system has been studied to get a bet-
ter understanding of movement disorders in patients
with cerebral palsy, stroke, osteoarthritis and Parkin-
son’s disease. Thousands of patients have been studied,
recording their neuromuscular excitation patterns both
before and after treatment. However, the detailed un-
derstanding of the function of each of the elements of
the musculoskeletal system remains a major challenge.
Researching these diseases in real-life experiments
has the following limitations. First, many important
variables are hard to measure in an experiment, such as
the force generated by a muscle. Second, it is difficult
to deduce cause effect relationships in complex dynamic
systems based on experimental data alone. These prob-
lems that arise when analyzing experimental data can
be solved for a large part by combining the experimen-
tal data with a neuromuscular simulation framework.
Neuromuscular simulation allows one to study the dif-
ferent facets of neuromuscular activity, specifically the
cause-and-effect relationships between muscular exci-
tation patterns, muscle forces and resulting motion of
the body. It can integrate theoretical models describing
the anatomy and physiology of the musculoskeletal sys-
tem and the mechanics of multi-joint movement. Simu-
lations also enrich experimental data by providing esti-
mates of important variables such as muscle and joint
forces, which are difficult to measure experimentally.
Another approach to simulating motion is the use
of
finite element analysis
(FEA). Finite element anal-
ysis is a method to simulate a complex environment
2
Peeters and Pronost
by representing it as a set of small elements which are
interconnected by a means of (differential) equations,
which define the properties of the environment. Typ-
ical data sources for creating FEA setups are volume
scans of subjects, using techniques such as
magnetic
resonance imaging
(MRI) or
X-ray computed tomogra-
phy
(CT). Using such a volumetric dataset of a subject,
a detailed FE simulation setup, consisting of a set of ele-
ments that describe bones, muscles and tendons, can be
created. The main problem of these datasets is the inac-
curacies which arise from and are inherent to the data
acquisition process. Even the best acquisition methods
such as done for the visible human dataset need some
steps to refine the data [2]. Since FEA offers such an ad-
vantage in simulating interaction between muscles and
bones over a more high-level type of simulation as mus-
culoskeletal simulation, the question arises whether we
could create a bridge between these types of simulation
environments.
The aim of this study is to propose a practical frame-
work capable of preparing the raw volumetric dataset
of a subject to an FE simulation of the deformation of
its soft tissue. The framework has also been designed
in a way that connection to musculoskeletal simulations
of the same subject be facilitated. We demonstrate our
method on rigid and soft bodies of the lower limb with
a particular attention to the knee area. We believe that
our framework is a step towards the automatic simu-
lation of muscle deformations in virtual humans in a
predictive manner through the interaction of the mus-
cle anatomy with applications in computer graphics.
Similarly, the method can simulate a complex muscle
structure so that muscle function can be investigated
in biomedical engineering.
2 Related work
The biomechanics community has ongoing efforts to
create detailed human musculoskeletal models. Although
recent models [3,4] provide accurate muscle parameters
for the whole body, they by definition do not provide
subject-specific geometrical data needed to simulate de-
tailed muscular models. These models use lumped pa-
rameter 1D muscle models that do not account for mus-
cle geometry. In musculoskeletal representation physi-
ological parameters such as muscle lengths and mus-
cle forces have been of primary interest, and the re-
alistic visualization has played a secondary role. Over
the years, neuromuscular simulation has evolved from a
fragmented community where each research group de-
veloped their own simulation software into a more col-
laborative environment. This has partly been the re-
sult of the efforts by the creators of the OpenSim plat-
form [5,3], by providing an open-source platform that
lets the user develop a wide variety of musculoskeletal
models that can be shared due to the open nature of the
software. We also observe a recent tendency to trans-
fer these musculoskeletal models and techniques to the
computer animation community [6].
Many attempts have been made at simulating mus-
cles in high detail using the method of finite element
simulation [7–9], for example, to investigate intramus-
cular pressures [10]. Teran et al. have designed a frame-
work for extracting and simulating musculoskeletal ge-
ometry from the visible human dataset [2,11]. In their
study they had to use a motion from a different sub-
ject, since the visual human dataset is obtained from a
deceased subject and no motion capture has been pre-
viously performed. The same dataset was used in ref.
[12] to create a 3D model of the human leg, specifically
for visualization of deformations, and incorporated also
the rendering of muscle fibers using textures.
In order to simulate interactions between rigid and
soft bodies within a common framework, two-way cou-
pling methods have been investigated. For example,
Shinar et al. [13] proposed a time integration scheme
for solving dynamic elastic deformations in soft bod-
ies interacting with rigid bodies. Despite heavy com-
putation times required for the finite element simula-
tion, their framework can handle two-way coupled con-
tact, collision, stacking, friction, articulation, and PD
control. Real-time performances have been achieved by
Kim and Pollard [14] by using reduced deformable mod-
els and linear-time algorithms for rigid body dynam-
ics. The primary focus of their work was to develop
a fast physically based simulation system for skeleton-
driven deformable body characters. In [15] Stavness et
al. focused on modeling and simulating subject-specific
anatomical structures composed of both hard and soft
tissue components. They demonstrated their framework
on the dynamics of the jaw tongue hyoid complex where
the deformation of the 3D FEM tongue model is driven
by muscle activations. Their work on a solver insists on
the need to integrate constraints directly into the ve-
locity solution.
In [2] the segmentation of anatomical entities was
performed by creating a level set representation of each
entity relevant for the simulation. The signed distance
function required for the level set procedures and to
generate the triangulated surface was the fast marching
method [16]. They use slice-by-slice contour sculpting
to repair problem regions. First they manually exam-
A practical framework for generating volumetric meshes of subject-specific soft tissue
3
Fig. 1
Our framework prepares subject-specific segmented MRI datasets for the FE simulation of soft tissues.
ine each slice visually to check for and eliminate errors.
Level-set smoothing techniques are applied afterwards,
such as motion by mean curvature [17] to eliminate
any further noise. In this work, we used MRI data of a
much lower resolution than the visual human dataset.
Also, the segmentation process we used is automatic
and therefore needs a more rigorous repair process.
Blemker and Delp [18] used an MRI dataset to cre-
ate simulations of several sets of muscles to study the
variation in moment arms across fibers within a mus-
cle and later to predict
rectus femoris
and
vastus in-
termedius
fiber excursions [19]. They used MRI of live
and cadaver specimen, and manually segmented the ar-
eas of interest. They also created a fiber map for each
muscle of interest, based on template fiber geometries
morphed to each muscles target fiber geometry. They
used a manual segmentation process instead of an au-
tomatic segmentation as is used in this work. While
this method of segmentation provides a surface mesh
of higher quality, the method is not automatic. The
segmentation method used in this project comes from
Schmid and Magnenat-Thalmann [20]. This method is
based on an earlier work of Gilles et al. [21], who pro-
posed a registration and segmentation method for clini-
cal MRI datasets based on discrete deformable models.
The main contribution of this study to the field of
musculoskeletal and FE simulations is the semiauto-
mated pipeline, that resolves data inaccuracies, such as
muscle overlaps and self-intersecting muscles. Further
developments of the framework will eventually allow for
the connection between musculoskeletal and FE simu-
lations, and thereby building a step towards a unified
platform for subject-specific motion analysis.
This paper is structured as follows. In Section 3 we
detail the main features of our framework manipulating
the data sources and preparing the final FE simulation.
In Section 4 and section 5 we respectively discuss the
limitations of the framework and summarize the find-
ings of the study and suggesting possible research di-
rections.
3 Practical framework
3.1 Overview
Our framework consists of several key steps (see Fig-
ure 1). First the lower limbs of a subject are scanned in
an MRI machine. The acquisition protocol used in our
experiments is presented in Appendix B. Then the lower
limbs are segmented using the algorithm of Schmid and
Magnenat-Thalmann [20]. The segmentation algorithm
produces closed surfaces of muscles and bones and the
attachment sites where the muscle connects to the bones
(Section 3.2).
The closed surfaces from the segmentation result
cannot be converted directly to volume meshes as the
segmented data may contain segmentation artifacts and
noise (Sections 3.3, 3.5 and 3.6) or be incomplete (Sec-
tion 3.4). These inconsistencies are first solved before
the surfaces can be given to the mesh generator (Sec-
tion 3.7). To complete the FE ready meshes, we finally
create attachment specifications that can be used by a
FE solver to constraint the movement of the soft tissues
(Section 3.8).
3.2 Segmenting MR images
To illustrate the features of our framework, we use a
segmented dataset resulting from a template method
that uses a minimal energy optimization to fit a tem-
plate muscle or bone to the acquired MRI data [20].
Since the method is dependent on forces, balancing
the weights of the non-penetration constraints and the
other constraints can be a difficult if not impossible
task. Therefore, the resulting surfaces may suffer from
intersections between surfaces and self-intersections, wh-
ich we resolve in this study. Since the density of the
vertices of the shapes is uniformly spread, the resulting
shapes have nearly the same property. The segmenta-
tion algorithm also includes tendon and attachment lo-
cations specified by vertex indices in the resulting mus-
cle meshes.
4
Peeters and Pronost
The segmentation algorithm produces surfaces of
bones and muscles. Surfaces are represented as a list
of vertices
V
=
{v
i
: 1
≤
i
≤
n
V
}
and a list of faces
F
=
{f
k
: 1
≤
k
≤
n
F
},
each face
f
k
= (i
k
, . . . , i
k
f
)
n
k
1
consisting of a sequence of indices in the vertex list.
A surface
S
=
{V,
F
}
is a pair of a vertex list com-
bined with a face list. All surfaces are consisting of
only triangles, so each face has always three vertices
f
k
= (i
k
, i
k
, i
k
).
1 2 3
An attachment site
A
is defined in the segmenta-
tion result as indices in a corresponding vertex list
V
in the order they are defined in the muscle files:
A
=
{i
1
, . . . , i
n
A
}.
Since some tasks in the pipeline change the amount
and ordering of the vertices in the data files, we process
the attachment sites to be index-invariant by defining
the attachment areas by geometry. This is explained in
Section 3.8. The tendons are defined in the same way
as the attachments, and for those we also need to con-
vert the indices to a geometric representation. This is
the protocol that we used to specify attachments and
tendons but any equivalent representation is possible
and can be used with minimal adjustments to our al-
gorithms.
Fig. 2
Smoothing results. Left is before and right is after the
smoothing operation.
3.4 Generating missing tendons
It is usual that after segmentation of the anatomical en-
tities from MRI some tendons are not present or that
tendons are present but their attachment to the mus-
cle is not defined. Due to MRI acquisition protocols,
location and shape of muscles, typical missing tendons
on the lower body are observed for the
gastrocnemius
muscle and the
soleus
muscle, which are attached to
the femur, tibia or fibula at the top of the muscle and
are missing their connection to the foot. Since these
muscles are valuable for a simulation of the lower body,
we developed an algorithm that automatically edits the
surface mesh to incorporate a tendon structure, includ-
ing an attachment and tendon specification expected
further along in the pipeline.
The tendon generator gets a list of muscles that do
not have tendons. For these muscles, the ‘bottom’ part
of the muscle is detected, and is extruded towards the
foot. In our example muscles, the bottom direction is
the negative direction along the superior-inferior axis,
which in our dataset is the
−z
direction. It should be
adapted to other scanning configurations, for example
by detecting the principal axis of the tibia bone.
z
As a first step of the algorithm, the
z-value v
f oot
max
of the top vertex of the foot
S
f oot
is determined so
z
z
that
v
f oot
max
> v
i
, v
i
∈
V
f oot
. For each muscle
S
muscle
,
the faces on the bottom part are recursively traversed,
z
starting with the faces connected to the vertex
v
muscle
min
with the lowest
z-value,
creating a set
F
bottom
of faces
belonging to the bottom part of the muscle. The set
F
bottom
is defined as follows.
z
F
bottom
0
=
{f
k
:
v
min
∈
f
k
}
af
F
bottom
t+1
=
F
bottom
t
∪ {f
k
:
f
k
∩
F
bottom
t
=
∅,
angle(f
k
)
< θ
muscle
}
−1
(f
n
angle(f
k
) =
cos 180π
k
·−z)
af
Here
f
k
is the set of faces adjacent to
f
k
,
angle(f
k
)
the function providing the angle between a face and the
xy-plane, f n
k
is the face-normal of face
f
k
and itera-
tions of
F
bottom
are denoted as
F
bottom
#
. Summarized,
we collect the set
F
bottom
of adjacent faces that have
3.3 Smoothing the surfaces
First, we have to apply a smoothing method to re-
move the high frequency artifacts from the segmenta-
tion algorithm. In this work we applied the method of
Taubin [22], since it is fitted to our type of biological
surfaces, which are closed surfaces with a nearly uni-
form vertex distribution. This method does not intro-
duce shrinking and is also the most suitable option with
regards to implementation complexity.
The algorithm is a low-pass filter, with three param-
eters:
λ, µ
and the iteration count
N
. The low-pass filter
properties are the pass-band frequency
k
P B
, the pass-
band ripple
δ
P B
, the stop-band frequency
k
SB
, and the
stop-band ripple
δ
SB
. The parameters are related to the
low-pass filter properties as follows.
λ<
1
k
SB
(1)
(2)
(3)
(4)
1
1
+ =
k
P B
λ µ
(λ
−
µ)
2
/(−4λµ)
N
N
<
1 +
δ
P B
< δ
SB
[(1
−
λk
SB
)(1
−
µk
SB
)]
Figure 2 shows the result of the operation on the
gastrocnemius
muscle using the parameters given in
Appendix B.
A practical framework for generating volumetric meshes of subject-specific soft tissue
5
Fig. 3
Left: The
soleus
muscle without tendons connecting
to the foot. Right: The
soleus
muscle with an attachment to
the foot generated automatically.
an angle with the
xy-plane
lower than
θ
muscle
. For our
missing tendons we used a value
θ
muscle
of 60 degrees.
After determining the faces that comprise the bot-
tom of the muscle, we start the extrusion. For each face
f
k
∈
F
bottom
, we set the
z-value
of each of its vertices to
z
z
z
v
iz
=
v
i
+ (v
f oot
max
−
v
muscle
min
). Then we determine
the edges
E
border
on the border of
F
bottom
. The vertices
in the border edge set
E
border
are duplicated and are set
to their original position. The gap between the original
border and the duplicated border is filled with regular
spaced vertices and corresponding faces
F
gap
. All the
faces from
F
bottom
and the newly added faces
F
gap
are
combined into a new set
F
tendon
=
F
bottom
∪
F
gap
, and
are together locally smoothed according to the method
described in Section 3.3. In Figure 3 it is shown how the
newly generated tendon for the
soleus
muscle compares
to the version without tendon.
The newly generated piece of muscle is also saved
as being tendon material by immediately generating a
tendon convex hull for the new vertices, as described in
Section 3.8.
This correction is not an ideal solution but the main
objective is to attach the muscle to the bone. Nonethe-
less, it is not advised to rely on simulations performed
on an area of interest where tendons have been gen-
erated this way. In our experiments on the knee joint,
only tendons near the foot needed to be generated, so
we are confident that our simulations are relevant in
regards to this issue.
3.5 Resolving self-intersections
The surface data provided by the segmentation algo-
rithm can contain self-intersections. The intersected faces
are inverted parts of the volume that the closed surface
describes. As these parts occur from bad segmentation,
our goal is not to untangle them as it would typically be
done in cloth simulation [23]. Indeed, that would add
information to the surface that should not be present
in the dataset, and it would be generated without pur-
pose. The only semantically correct way of resolving
them is by removing them in a way that the final sur-
face is still closed. We solved this issue by adapting
the algorithm proposed by Jung et al. [24] originally
designed to remove self-intersections from a raw offset
triangular mesh. The algorithm uses a region growing
approach, keeping a list of valid triangles. Starting with
an initial seed triangle, the algorithm grows the valid
region to neighboring triangles until it reaches triangles
with self-intersection. Then the region growing process
crosses over the self-intersection and moves to the ad-
jacent valid triangle. Therefore the region growing tra-
verses valid triangles and intersecting triangles adjacent
to valid triangles only. The region growing process is
complete when no more faces have to be added to the
valid region. We will describe here the adaptations we
made to the algorithm of Jung et al. to fit it to our
purpose.
As in the original algorithm, we remove faces that
have edges smaller that a predefined
ε
e
and reconfig-
ure faces that have angles smaller than
ε
α
. All faces in
each surface are checked for these two properties and
solved accordingly. For faces where an edge
e
p
= (v
i
, v
j
)
is smaller than
ε
e
, that edge is collapsed by setting
v
i
= (v
i
+
v
j
)/2 and
v
j
is removed. Then the faces
containing
v
j
are updated to refer to
v
i
. But after ap-
plying this operation a number of times, the geometry
might contain some new topological problems. The sur-
face can contain two triangles that are opposite to each
other. We have adapted the algorithm so that these
faces are detected and removed, which is a valid op-
eration, since it is a part of the surface that defines
no volume. For a triangle
f
k
that has an angle smaller
than
ε
α
, we do not remove any edges, but we change
its configuration. First we determine the longest edge
e
u
= (v
i
, v
j
), and find the triangle
f
m
sharing that edge.
We then have
f
k
=
{v
h
, v
i
, v
j
}
and
f
m
=
{v
g
, v
j
, v
i
}.
We now swap the edge
e
u
, such that
f
k
=
{v
h
, v
i
, v
g
}
and
f
m
=
{v
g
, v
j
, v
h
}.
Since the method originally was developed for raw
offset triangular meshes, selecting a seed triangle was
easily determined by calculating the convex hull
V C
⊂
V
, and picking any triangle
f
∈
F
which has at least
vertex
v
i
∈
V C.
In our case however, triangles taking
part in the convex hull can also be invalid, as is the
case in the
vastus intermedius
muscle. For this reason,
we added a constraint to the seed triangle selection.
First, we calculate the center of mass of the vertices of
the surface, by a simple averaging of the vertex posi-
tions.
v
i
p
com
=
v
i
∈V
(5)
|V |
Since the surface meshes in this work are of approx-
imately of uniform density (see Section 3.2), we have

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika:

Zgłoś jeśli naruszono regulamin