Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Acquisition of human tissue elasticity properties using pressure sensors
(USC Thesis Other)
Acquisition of human tissue elasticity properties using pressure sensors
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Acquisition of Human Tissue Elasticity Properties Using
Pressure Sensors
by
Danyong Zhao
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(COMPUTER SCIENCE)
May 2022
Copyright 2022 Danyong Zhao
I would like to thank my advisor, Prof. Jernej Barbiˇ c a lot, for the inspiration and
support during my PhD study.
I would also like to thank my parents for encouragement and support,
the committee memebers for constructive comments and feedbacks,
my labmates Hongyi Xu, Yijing Li, Bohan Wang and Mianlun Zheng for discussions
and help,
Computer Science Department of USC for funding and help.
ii
Table of Contents
Dedication ii
List of Tables vi
List of Figures vii
Abstract ix
Chapter 1: Introduction 1
1.1 Elasticity Theory on Human Simulation . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Deformations and Elastic Forces . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Isotropic Elastic Material . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.3 Analytical Solution for Continuum Linear Elasticity . . . . . . . . . . . . 5
1.2 Human measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Chapter 2: Related Work 9
2.1 PDEs for the displacement field . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Measuring Young’s Modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Material Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Chapter 3: Human Measurements in Vivo 12
3.1 Overview of our approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2 Design of the Measuring Device . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.3 Human Mesh Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Chapter 4: Fast Simulation Using Symmetric Facts 20
4.1 Cylindrical Coordinate and Axial symmetry . . . . . . . . . . . . . . . . . . . . . 20
4.2 Gaussian Points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.3 2D Finite Element Method on 3D mesh . . . . . . . . . . . . . . . . . . . . . . . 22
4.4 Validation on 3D object . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Chapter 5: Contact Force Modeling 27
5.1 Collision Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.2 Constraint-based Normal Contact Force . . . . . . . . . . . . . . . . . . . . . . . 28
iii
5.3 Frictional Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.4 Line Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.5 Reduced Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Chapter 6: Material Optimization 32
6.1 Spline-based material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.1.1 Valanis-Landel Material . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.1.2 Natural Cubic Spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
6.2 Conjugate Gradient Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.3.1 Foam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.3.2 Silicone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
6.3.3 Convergence on mesh subdivision . . . . . . . . . . . . . . . . . . . . . . 40
6.3.4 Extension: Estimate the depth of the material sample . . . . . . . . . . . . 42
6.3.5 Measurement on Human . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
Chapter 7: Application: Designing Ergonomic Furniture 46
7.1 Introduction of Ergonomic Design . . . . . . . . . . . . . . . . . . . . . . . . . . 46
7.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
7.3 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
7.3.1 Modeling the human body . . . . . . . . . . . . . . . . . . . . . . . . . . 53
7.3.2 Modeling supporting surfaces . . . . . . . . . . . . . . . . . . . . . . . . 55
7.3.3 Nonlinear elastostatics coupled to the rigid core . . . . . . . . . . . . . . . 56
7.3.4 Computing the static equilibrium . . . . . . . . . . . . . . . . . . . . . . 58
7.3.5 Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.4 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7.4.1 Ergonomic Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7.4.2 Objective Function Gradient . . . . . . . . . . . . . . . . . . . . . . . . . 65
7.4.3 Deformable Supporting Surface . . . . . . . . . . . . . . . . . . . . . . . 66
7.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
7.5.1 Human sitting in a chair . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
7.5.2 Human standing on a shoe sole . . . . . . . . . . . . . . . . . . . . . . . . 69
7.5.3 Human grasping a mouse . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
7.5.4 Human sitting in a deformable chair . . . . . . . . . . . . . . . . . . . . . 71
7.5.5 Smoothing Strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
7.6 Physical validation of our method . . . . . . . . . . . . . . . . . . . . . . . . . . 74
7.7 Comparison to “Sit & Relax” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.8 Comparison to “Pose to Seat” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Chapter 8: Conclusion 84
References 87
iv
Appendices 93
A Contact Gradients and Hessian of Contact Energy . . . . . . . . . . . . . . . . . . 94
A.1 Vertex anchor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
A.2 Edge anchor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
A.3 Face anchor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
v
List of Tables
3.1 Samples on the human arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
6.1 Optimization results, numbers of optimizing iterations and timings of the four ma-
terial samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
7.1 Ergonomic scores and optimization timings . . . . . . . . . . . . . . . . . . . . . 68
7.2 Ergonomic scores using different smoothing values . . . . . . . . . . . . . . . . . 76
vi
List of Figures
1.1 Locality in poking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.1 Overview of the poking devices . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Relationship between the poking indentation and distance captured by laser . . . . 15
3.3 Contact force and human skin deformations captured using our experiment . . . . 16
3.4 Ultrasound scan images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5 Contact force and human skin deformation captured from the experiment . . . . . 18
3.6 Poking locations on the right arm . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.7 Tetrahedral mesh generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.1 Symmetry ofϕ component . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.2 Positions and weights of Gaussian points . . . . . . . . . . . . . . . . . . . . . . . 23
4.3 Validation using 3D simulation on the bunny . . . . . . . . . . . . . . . . . . . . . 26
4.4 Experiment set up on the bunny . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.1 Distribution of theλ values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
6.2 Man-made material samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.3 Fitting and validation results of the Soma Foama 25 material . . . . . . . . . . . . 39
6.4 Fitting and validation results of the Soma Foama 15 material . . . . . . . . . . . . 40
6.5 Fitting and validation results of the softer Dragon Skin FX Pro material . . . . . . 41
6.6 fitting and validation results of the stiffer Dragon Skin FX Pro material . . . . . . . 41
6.7 Validation of 2D simulation method . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.8 Fitting results on different height depths . . . . . . . . . . . . . . . . . . . . . . . 43
6.9 Fitting results of the best depth guess . . . . . . . . . . . . . . . . . . . . . . . . . 43
vii
6.10 Optimization results of human tissue material . . . . . . . . . . . . . . . . . . . . 45
7.1 Optimizing the shape of the supporting surface (chair) to maximize ergonomics of
contact against a human . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
7.2 Modeling and optimizing the contact between a human and a man-made object . . 54
7.3 Contact between the human body and a supporting surface . . . . . . . . . . . . . 55
7.4 Bidirectional contact model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.5 Sensitivity of our ergonomic cost to surface bumpiness . . . . . . . . . . . . . . . 64
7.6 Optimized pressure maps and ergonomic scores (chair). . . . . . . . . . . . . . . . 70
7.7 The rigid core model of the foot . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
7.8 Optimized pressure maps and ergonomic scores (shoe sole) . . . . . . . . . . . . . 72
7.9 Optimized pressure map and ergonomic scores (mouse) . . . . . . . . . . . . . . . 73
7.10 Optimized pressure maps and ergonomic scores (deformable chair) . . . . . . . . . 75
7.11 Optimized shapes and ergonomic scores over the iterations of the optimization,
using different smoothing values . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
7.12 Manufactured models and the force sensor . . . . . . . . . . . . . . . . . . . . . . 78
7.13 Output voltages of the reference sensor across the 78 measurements . . . . . . . . 78
7.14 Comparison between simulation and real measurements . . . . . . . . . . . . . . . 80
7.15 Comparison to Sit & Relax [68] . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.16 Comparison to Pose to Seat [69] . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
viii
Abstract
Physically based simulation of the human body in three dimensions is important in many appli-
cations in computer graphics, animation, virtual reality, virtual commerce, ergonomics and virtual
medicine. Finite Element Method (FEM) is a robust and reliable approach to simulate deformable
dynamics of three-dimensional elastic structures. However, for quality simulation that matches the
behavior of real human tissues, FEM needs accurate material properties that correctly model real
relationships between the strains and stresses in the human tissue. This thesis presents methods to
capture such nonlinear materials for the human musculo-skeletal tissue (skin, fat, muscles) in vivo,
through carefully designed poking experiments, force meters, lasers and ultrasound measuring de-
vices. From our experiments, we obtain ground-truth relationships between the contact force and
skin deformation. We then fit material models that best approximate the acquired real-world data.
First, we design a measuring device that can simultaneously capture the skin contact force and
the skin deformation, consisting of a force meter and a laser distance measuring device. We design
sliding and pivoting joints to rigidly attach the force meter to the laser device, so that diverse human
body locations (arm, hand, belly, etc.) can be measured ergonomically and reliably. We also use
an ultrasound device to capture the depth of the human subcutaneous fat at different locations;
enabling us to generate a 3d model of the fat layer, and optionally also the muscle layer, of the
human body.
Second, we propose a novel approach to equivalently convert 3D FEM simulations into 2D
simulation, suitable for our material capture. This method permits us to greatly speed up our
material optimizations, without losing any accuracy. We validated this approach by comparing the
simulation result from 2D equations and the 3D traditional equations. We propose to use natural
cubic splines to parameterize the three separable scalar elastic energy density functions based on
ix
the Valanis-Landel material model. Based on our novel 2D simulation method and the spline-
based non-linear isotropic material, we present an efficient method to compute the gradient of the
objective function for optimization and use the conjugate gradient optimization method to optimize
the material.
Lastly, we use our acquired materials to design the geometric shape of rigid supporting surfaces
to maximize the ergonomics of physically based contact between the surface and a deformable
human. We model the soft deformable human using a layer of FEM deformable tissue surrounding
a rigid core, with measured realistic elastic material properties, and large-deformation nonlinear
analysis using our material capturing and optimizing method. We define a novel cost function to
measure the ergonomics of contact between the human and the supporting surface. We give a stable
and computationally efficient contact model that is differentiable with respect to the supporting
surface shape. This makes it possible to optimize our ergonomic cost function using gradient-based
optimizers. We 3D-print the optimized shoe sole, measure contact pressure using pressure sensors,
and demonstrate that the real unoptimized and optimized pressure distributions qualitatively match
those predicted by our simulation.
x
Chapter 1
Introduction
Human body simulation is an interesting research problem in many different industries such as film
and games, especially when one needs to compute human skin deformations when the character is
in contact with other objects. To better model the human body, many models for human anatomy
have been proposed. In most anatomy models, human body has skin, fat, fascia, muscle and the
bone layers. To improve skin dynamics, each layer may have different material properties and dif-
ferent simulation methods. Some methods, such as Finite Element Method and Projective Dynam-
ics, are based on physics laws, while some are non-physical such as Position-base Dynamics. This
thesis uses Finite Element Method (FEM) to simulate the human body. It builds the relationship
between elastic forces and object deformations. We also aim at designing a good material model
for the various human tissue layers so that the FEM simulation matches the force-deformation re-
lationships captured from real experiments when the human body is poked by different objects.
We also propose new devices to measure the properties of human tissue.
1.1 Elasticity Theory on Human Simulation
Three-dimensional FEM is commonly used in solid deformable objects including human charac-
ters because FEM obeys physical laws. When the character is in collision, penalty-based method
and constraint-based method are two popular ways to model the contact force. The constraint-
based method builds a set of equations and inequalities so that the entire mesh is collision-free.
1
The penalty-based method applies contact force on vertices that penetrate into another surface, and
the magnitude of the contact force depends on the penetration depth. Usually, the penalty contact
force is linear in the penetration depth for simplicity. Then FEM is used to compute the human de-
formation given the contact force and the elastic force, as well as the material properties of human
tissue. There are multiple material types in common simulation usage, such as St.Venant-Kirchhoff
(StVK), Neo-Hookean and Ogden. As human tissue is largely incompressible, the volume preser-
vation is another important consideration.
1.1.1 Deformations and Elastic Forces
The deformation of a deformable solid object can be described by a function⃗ x=φ(
⃗
X), where
⃗
X
is the rest position and⃗ x is the deformed position of the object. As is standard, we further define
deformation gradient F, and first Piola-Kirchhoff stress P for hyperelastic materials as
F =
∂⃗ x
∂
⃗
X
, (1.1)
P=
∂Φ(F)
∂F
, (1.2)
where the function Φ is determined from elastic material properties. We will discuss it in more
detail in Section 1.1.2. In our work, the simulation mesh is always a tetrahedral mesh and the
deformation gradient of one tetrahedron can be computed by
F =
x
1
− x
0
x
2
− x
0
x
3
− x
0
X
1
− X
0
X
2
− X
0
X
3
− X
0
− 1
. (1.3)
Irving’s work [1] proposed a method to compute elastic forces on a vertex i of a tetrahedron
as g
i
=− P(A
1
N
1
+ A
2
N
2
+ A
3
N
3
)/3, where A
j
N
j
are the area-weighted normals of the three faces
that are incident to vertex i at the rest shape. Then for the entire simulation, the elastic force f
e,i
on
vertex i is computed by adding g
i
from all incident tetrahedra.
2
1.1.2 Isotropic Elastic Material
Due to symmetry of isotropic materials, the strain energy density function of isotropic elastic
materials can be defined as Φ(λ
1
,λ
2
,λ
3
) :R
3
7→R, whereλ
i
are singular values of the deformation
gradient F. The first Piola-Kirchhoff stress P can also be computed from the singular values of F
by
P= P(F)= UP(
ˆ
F)V
T
, (1.4)
where F = U
ˆ
FV
T
and diag(
ˆ
F)=(λ
1
,λ
2
,λ
3
) [1].
For hyperelastic isotropic materials, invariant-based strain energy will be more convenient to
compute the elastic force and stiffness matrix. The three invariants I
C
,II
C
,III
C
can be computed
from the three singular valuesλ
1
,λ
2
,λ
3
I
C
= trace(C)=λ
2
1
+λ
2
2
+λ
2
3
, (1.5)
II
C
= C : C=λ
4
1
+λ
4
2
+λ
4
3
, (1.6)
III
C
= det(C)=λ
2
1
λ
2
2
λ
2
3
, (1.7)
where C= F
T
F. The strain energy density function of many commonly used isotropic materials
can be written using the invariant-based function. For StVK, the function is
Φ(I
C
,II
C
,III
C
)=
1
8
λ
Lam´ e
(I
C
− 3)
2
+
1
4
µ
Lam´ e
(II
C
− 2I
C
+ 3), (1.8)
where λ
Lam´ e
and µ
Lam´ e
are the Lam´ e constants. For Neo-Hookean material, the strain energy
density is
Φ(I
C
,II
C
,III
C
)=
1
2
λ
Lam´ e
(I
C
− 3)− µ
Lam´ e
(I
C
− 3)logJ+
1
2
µ
Lam´ e
(logJ)
2
, (1.9)
3
where J=
√
III
C
. For Mooney-Rivlin model, the strain energy density is
Φ(I
C
,II
C
,III
C
)= C
01
(
I
C
− II
C
2
J
− 4
3
− 3)+C
10
(I
C
J
− 2
3
− 3)+ D
1
(J− 1)
2
, (1.10)
where C
01
, C
10
and D
1
are material constants. For Yeoh model, the strain energy density is
Φ(I
C
,II
C
,III
C
)=
3
∑
i=1
C
i
(I
C
− 3)
i
, (1.11)
where C
i
are material constants. Not all isotropic materials have an invariant-based strain energy
density function. Due to the symmetry in all the principle streches λ
1
, λ
2
, λ
3
, the strain energy
density function of many materials can be simplified into three seperated functions (Valanis-Landel
hypothesis) f(x),g(x),h(x) as
Φ(λ
1
,λ
2
,λ
3
)= f(λ
1
)+ f(λ
2
)+ f(λ
3
)+ g(λ
1
λ
2
)+ g(λ
2
λ
3
)+ g(λ
3
λ
1
)+ h(λ
1
λ
2
λ3). (1.12)
For linear co-rotational material, the seperated functions are
f(x)=
1
2
λ
Lam´ e
(x
2
− 6x+ 5)+µ
Lam´ e
(2− 2x),
g(x)=λ
Lam´ e
(x− 1)+µ
Lam´ e
(x− 1),
h(x)= 0. (1.13)
For Ogden material, the seperated functions are
f(x)=
N
∑
i=1
µ
i
α
i
(x
α
i
− 1),
g(x)= h(x)= 0, (1.14)
whereµ
i
andα
i
are material constants.
4
1.1.3 Analytical Solution for Continuum Linear Elasticity
If we only consider linear elasticity, when the system reaches its static equilibrium, Navier pre-
sented equations for three dimensional displacement field u in a homogeneous isotropic solid as
1
1− 2ν
∂(▽· u)
∂x
+▽
2
u
x
= 0, (1.15)
1
1− 2ν
∂(▽· u)
∂y
+▽
2
u
y
= 0, (1.16)
1
1− 2ν
∂(▽· u)
∂z
+▽
2
u
z
= 0, (1.17)
where ν is the Poisson’s ratio. In a special case, if the displacement field is axis-symmetric, a
cylindrical coordinate system can be used to reduce the 3D displacement field into a 2D partial
differential equation [2, 3] as
1
1− 2ν
∂(▽· u)
∂r
− u
r
r
2
+▽
2
u
r
= 0, (1.18)
1
1− 2ν
∂(▽· u)
∂z
+▽
2
u
z
= 0. (1.19)
If we assume that all contacts are frictionless and the poking depth is uniform, the boundary con-
ditions are
u
z
(r,D)= 0, (1.20)
u
r
(r,z)= u
z
(r,z)= 0,r→∞, (1.21)
u
z
(r,0)= x,r≤ a, (1.22)
(
∂u
r
∂z
+
∂u
z
∂r
)|
z=0
= 0,r> a, (1.23)
(
1
1− 2ν
(▽· u)+
∂u
z
∂z
)|
z=0
= 0,r> a, (1.24)
where a is the radius of the contact circle, x is the poking depth and D is the depth of the tissue.
5
Ragione [4] gives solutions to the effective stiffness, using the Hankel transform theory [5] as
F
x
=
aEG(
D
a
)
2(1− ν
2
)
, (1.25)
where F is the total poking force and G is a function that only depends on D/a. Ragione did not
give the analytical solution for function G but gives a numeric method to compute G. When D/a
is very large, G(∞)= 4 and Equation 1.25 converges to the standard Boussinesq solution
F
x
=
2aE
(1− ν
2
)
. (1.26)
When the depth of the solid is very shallow, G(D/a→ 0)→ a/D and we have
F
x
=
a
2
E
2(1− ν
2
)D
. (1.27)
1.2 Human measurements
Accurate geometric and material properties of human tissue can significantly improve the qual-
ity of soft-body human simulations. There are two different types of human measurement, one
is measuring tissue outside of the body (in vitro), and the other is measuring humans in vivo. In
medical research, measuring human tissues in vitro from cadavers is commonly conducted because
different human tissues can be detached and measured more accurately. However, in vivo measure-
ments have the promise to be much more accurate, and they make it possible to personalize the
measurements to each specific human. They are also easier to perform than in vitro measurements.
This makes in vivo measurement attractive in computer graphics. In this thesis, we design devices
and procedures to measure humans in vivo. Our devices are able to capture the human shape,
tissue depths and force-deformation curves for the human body, which can be used in subsequent
physically based simulation.
6
From our experiments, we observe that the skin deforms very locally when poked ( 1.1). This
implies that other locations that are far from the poking center are almost not affected by the
poking at all. Because of this, we can poke the human sample by sample, and measure the force-
deformation curves for different locations independently.
Figure 1.1: Locality in poking: The figure shows that when the arm is poked by a cylinder, only
the skin near the cylinder deforms, whereas most parts still remain at their rest shape. The left
two images are the result of a poking simulation whereas the right two images are at the rest shape
when the cylinder just touches the skin surface. In the left-bottom picture, we can clearly see that
when poked by a cylinder, the indentation is very local and most of the limb does not deform;
compare to the rest shape on the right.
1.3 Thesis Overview
This thesis first proposes a new device to measure human tissue stiffness. It works by combining
a force meter and laser measurements to capture contact forces and human skin deformation data
stimuteneously. We also use ultrasound to capture depths and deformations of interior human
tissues such as fat and muscles.
7
To improve the human skin simulations, we also use a 3D scanner to obtain the rest shape of
the human body. With the depth data from the acquisition experiment, we generate the geometry
of the human body and tetrahedralize the mesh. Then we propose a new material model and
optimize the material properties for each poking location so that the force-deformation curve from
the simulation is close to the one acquired from the real poking experiment.
To adapt to different materials, we use the conjugate gradient optimization method. We try to fit
the material properties by fitting the three scalar funtions in Valanis-Landel material model using
natural cubic splines so that the simulation result could be close to the ground-truth deformation-
force curves obtained through the poking experiment. Then we try to validate this set of material
properties and the whole simulation model by performing another poking experiment at the sample
locations; but with different poking shapes. Our goal is to show that the simulation results match
the poking experiment even under contact shape not used during training.
Finally, we introduce an application of acquiring the human tissue elasticity propertis, design-
ing the geometric shape of rigid supporting surfaces to maximize the ergonomics of physically
based contact between the surface and a deformable human. We divided the human body into the
elastic part that uses the acquired elasticity properties by our method, and a rigid core. We mod-
eled the contact using uni-directional sliding springs and define a novel cost function to estimate
the ergonomics. We also validate our results by scanning a real human subject’s foot, 3D-printing
the optimized shoe sole, measuring contact pressure using pressure sensors, and demonstrate that
the pressure distributions qualitatively match those predicted by our simulation.
8
Chapter 2
Related Work
2.1 PDEs for the displacement field
In the scope of 3D linear elastic theory, Navier presented a displacement formulation for homo-
geneous, linear, isotropic materials. For the common case where indentations occur on a circular
contact surface, the system is axisymmetric, and cylindrical coordinates are commonly used to
simplify the 3D PDEs into a 2D PDE in terms of radius and depth. Hankel transformation is then
used to solve the PDE [6, 7]. The Boussinesq problem [8] simplifies the human tissue into an
infinite half space made of a linear homogeneous material. Analytic solutions of under several
boundary conditions, such as poking depth, contact pressure and contact shape, are given in [9–
11].
In addition to the global homogeneous problem, the analytical solution for layers of homo-
geneous isotropic materials has also been explored. Under a uniform vertical circular load, [12]
gives the solution for three-layer materials. For an infinite layer of finite depth under a uniform
vertical circular deformation, [13] gives the analytical solution for effective stiffness for shallow
depths, and approximate solutions for arbitrary layer depths.
9
2.2 Measuring Young’s Modulus
Material properties of human tissues are important in many fields. Young’s modulus is a key tis-
sue property, establishing the relationship between stress an strain. To measure Young’s modulus
on humans or animals, the tissue is assumed to be an isotropic incompressible homogeneous hy-
perelastic material [14]. The suitability of this assumption has been confirmed experimentally
for many types of tissues such as human breast [15, 16], human spleen [17, 18] and human
liver [19–21]. We note that there are some types of human tissues that are compressible such
as lungs [Sobin:1988:CEF, 22]. Many medical experiments perform a very shallow poking [23–
25] so that the formulas derived from the Boussinesq problem can be used to compute the Young’s
modulus. Besides computing Young’s modulus from the Boussinesq solution, using optimization
methods through physically based simulation [26, 27] is another approach to fit the Young’s mod-
ulus to the experiment data.
2.3 Material Design
For humans and other biological tissues, we can assume isotropic behavior when poked. For a
small deformation, the simplest material is linear co-rotational material [28, 29]. This material
behaves linearly, but extracts the local rigid rotation using a polar decomposition [30]. Under large
deformations, commonly used non-linear material models include Neo-Hookean [31], St. Venant-
Kirchhoof(StVk) [32], Mooney-Rivlin [33, 34], Ogden [35] and Yeoh [36, 37] materials. The
strain energy produced by these well-known material models are all polynomial in the deforma-
tion. In [38], they generalized these polynomial materials, and proposed a spline-based polynomial
material that is C2-continuous so that the material non-linearity can be controlled, and implicit in-
tegrations used for simulation. For non-polynomial strain energy materials, Fung-Demiray and
Veronda-Westmann materials use exponential strain energy, whereas Arruda-Boyce and Gent ma-
terials are based on a molecular chain network. [14] compares many types of isotropic material
10
model and gives analysis of their suitability for different types of tissues such as spleen, liver,
kidney, bran breast and lung, for different species including pigs, cattle, rats and humans.
Anisotropic material properties have also been researched in several references [39, 40]. Non-
linear anisotropic material are also used in medical simulation for soft tissues under a large de-
formation and transversely isotropic behavior [41–43]. These anisotropic materials require two
directions of equal stiffness. When all three direction give different stiffnesses, one obtains or-
thotropic materials, suitable for wood, plants and muscles [44, 45].
11
Chapter 3
Human Measurements in Vivo
3.1 Overview of our approach
In this chapter, we give the design of our material measuring device, and the method to generate the
simulation mesh geometry. The goal of our measurements is to obtain the deformation-force curve
that can be used as a ground truth for human simulations. In our poking experiment, we select
several sample positions on the human skin and poke them, in the normal direction. We use a force
gauge and a (class 2) laser to capture the contact force and skin deformation simultaneously. Such
a setup is able to minimize errors caused by human movements. We also use a force gauge and an
ultra-sound device to measure the human tissue deformation under the skin.
To improve synchronization, we always stay at the “just-touching” position for two seconds af-
ter the experiment started. Then, we slowly perform the poke until reaching some maximum force
(typically≈ 8N). Then we slowly pull the device out of contact and pause at the “just-touching”
position for two seconds again, before breaking contact. In this way, we can easily align the contact
force with deformation and plot the deformation-force curves.
After fitting all deformation-force curves and scanning the rest shape of the human body, we
generate a simulation tetrahedral mesh and ground-truth deformation-force relationships that the
simulation should achieve at the sample positions.
12
3.2 Design of the Measuring Device
The goal of our measuring procedure is to capture deformations and contact forces simultaneously.
The simplest part is to measure the contact force using a force gauge. Our contact forces do not
exceed 10N, and therefore, we choose PCE-LFG 100 [46] with an accuracy of 0.01N, which is
sufficient for our experiments. This force meter also has a replaceable poking tip so that we can
control the shape of contact. Since the force gauge can only display the magnitude of the contact
force, we are always carefully poking the human skin along the normal direction, so that we can
know the three components of the total contact force within the poking area.
Figure 3.1: Overview of the poking devices: We use a force gauge to measure the contact force
along the surface normal. There is a replaceable poking tip and an optional extender mounted
on the force gauge (left part of the figure) so that we can control the contact shape. Here, we
3D-printed several poking cylinders with different radii. The force gauge also supports mounting
additional devices; and we mounted wooden sliders and joints (top part) to hold the laser setup (the
small square box at the center). When adjusting the slider, we can control how far the laser beam
is from the poking location. The joint is to make the shooting direction adjustable so that we can
ensure the laser beam direction is parallel to the poking direction.
Measuring the skin deformation is more difficult than measuring contact forces. Because it is
almost impossible to keep the human body perfectly still, a single positional sensor is typically
not enough to capture the relative skin deformation. If a set of position sensors are placed on the
human body, the synchronization would be another challenge. Here, we propose a new method
to directly capture the skin deformation using a laser device. Laser is commonly used to measure
13
the distance and we mount a FASTUS laser measurement sensor [47] on the force gauge. In our
design, we aim the laser at a fixed human tissue position, located close to the poking location but
sufficiently far so that it does not deform by the poking. Since the poking deformation is very local
(Fig. 1.1), such a position is not difficult to find. There are two benefits of using a laser sensor. One
is that the accuracy can be as high as 20µm. The other is that as long as the direction of the laser
is parallel to the direction of poking, the captured deformation is always accurate regardless of the
normal direction. However, to reduce the error caused by the deformations of the laser contact
point, we always try to aim at a location whose normal is similar to the poking direction.
The laser housing needs to be position at a distance of 5cm to 15cm to the laser contact point.
We designed two sliding wooden frames so that the distance between the force gauge and the laser
is adjustable in both horizontal and vertical directions. Because of this design, our device can be
used to poke multiple types no matter how the local human geometry is shaped. Also, due to the
requirement that the poking direction must be parallel to the laser shooting direction, we add a
joint between the two wooden frames, so that we can adjust the laser optical axis to be parallel to
the force meter.
Since we are measuring human tissues in vivo, the deformation of the inner human tissue layers
is more challenging to capture as it is invisible to lasers. We propose an approach to measure the
inner layer using ultrasound. We rigidly mount the ultra-sound device onto the force gauge so
that we can simultanesously capture the contact force and the depth of each layer from ultra-sound
images. Since the depth is always measured from the skin layer, the human movement does not
affect the results, which is a good property to have. From the skin deformation measurements and
depth data, we can compute the deformations of each layer interface.
The synchronization of the contact force and distance data captured by the laser is facilitated
by both devices having wired connections to a computer. This makes it possible to record the
data in real time. We need to synchronize the start time and the capture rate for both devices.
The ultra-sound device outputs the data in the form of an image. The X-axis of the ultra-sound
image represents time. We start the force gauge and ultra-sound scan at approximately the same
14
Figure 3.2: Relationship between the poking indentation and distance captured by laser: The
black lines represent the human surface we are poking (not necessarily flat) and the blue rectangle
represents the poking force meter. The green box represents the laser housing; it aims the red laser
onto the human surface. The poking force meter and the laser housing are rigidly connected. The
top two images illustrate a successful measuring case where the laser optical axis is parallel to the
poking direction. Comparing the top-left image and the top-right image, the force meter indents
the human surface by d, and the laser also moves d in the shooting direction. In this successful
case, the laser is aiming at the same location on the human surface during poking, and the normal
at the laser contact does not need to be the same as the poking location. The bottom images show
failure cases where the poking direction is not parallel to the laser optical axis.
time, and wait for two seconds before starting to poke. After the poking is done, we also maintain
the touching position for about two seconds before stopping the two devices. In this way, we
can see the start and end of the poking clearly on both curves and temporally align the data. We
post-process the curves using low-pass filtering to remove vibrations and noise.
3.3 Human Mesh Generation
To generate the simulation mesh, we first use Artec Eva 3D scanner to scan the rest shape of
human body. Generally, the mesh captured from a scanner is polygon soup, and the Artec Eva
software will always try to reconstruct the surface into a high-quality 2-manifold mesh. However,
the resolution is too high, causing slow simulation times. Thus, we need to resample the surface
15
Figure 3.3: Contact force and human skin deformations captured using our experiment: The
left image shows the contact force captured over time, and the right image shows the voltage data
output by the laser. We pre-tabulate the relationship between the laser output distance and laser
output voltage so that we can convert the voltage curve into the distance curve. The distance at the
start of poking represents the rest shape. We can compute the skin indentations at any time from
the distance curve.
mesh by first finding 20 vertices on the mesh that are as far among each other as possible. Then we
compute an isosurface mesh from the initial 20 sample vertices by the method described in [48].
If we scanned the whole human body, another method to generate a plausible human surface is
to use a high quality template of human mesh [49]. We first adjust the pose of a template human
mesh so that it is close to our scanned result. Then, a non-rigid ICP for surface registration [50] is
performed so that the shape of the high quality template mesh can match the shape captured from
the 3D scanner. After generating the surface mesh, inner surface meshes are needed to approximate
the interface layers of different human tissues. One approach is to assume contact tissue depth.
Another is to generate inner layers with different tissue depths.
Inner layers with uniform-depth can be generated by displacing the surface layer with a
uniform depth. We first compute a signed distance field of the surface mesh [51]. We then use
marching cubes [52] for a constant non-zero depth to generate an inner layer mesh. In the second
step, the Delaunay resampling method mentioned in the previous paragraph can be used to generate
a good inner layer mesh.
Inner layers with spacially-varying depth: From the ultra-sound images, we can obtain the
depths of bones, muscles or fat at the sample locations in the rest shape. Then, we compute
Lapacian weights for all surface vertices over the sample locations so that the tissue depths can be
16
Figure 3.4: Ultrasound scan images: The horizontal axis represents time and the vertical axis
represents tissue depth. Black color means a weak ultrasound reflection signal, blue color means
intermediate reflection signal, and white color means denotes strong reflection signal. In the two
images, we can clearly see that there are several human tissue layers. The interface depths change
over time during the poking process. Here, we assume that the first interface is between the fat
and the muscle while the last interface is between the muscle and the bone. We use Catmull–Rom
splines to fit the depths over time and draw the splines (green) on the ultrasound images. The red
dots are the control points of Catmull–Rom splines we selected on the image. The yellow lines at
the endings of the splines are the tangent directions of the splines at the two endings.
interpolated to the entire mesh. We displace all vertices along the normal direction by the depth
values, and use the ElTopo library to generate the surface of the interface layers.
Finally, we use TetGen [53] to generate tetrahedral meshes for simulation from good 2-manifold
surface meshes. Since we assume that human bones are rigid, we do not need to mesh the bone in-
terior. We can remove those tetrahedrons by computing winding numbers [54] of each tetrahedron
gravity center over the bone interface layer mesh.
3.4 Results
We first use the poking device to poke the arm and the human back. From the the indention curve
and the contact curve, we manually determine the the poking starts and poking maximums. By
selecting the data points between this time frame, we obtain the indentation-contact force profile.
Some profiles are quite linear and some show that the human tissue becomes stiffer under a larger
deformation (Figure 3.5). To help generate the simulation mesh and the outer surface that has
been captured by the 3D scanner, we need to generate an approximate inner surface layer, using
17
Figure 3.5: Contact force and human skin deformation captured from the experiment: The
blue dots are the sample data read directly from the poking device. The red line is a fitted curve of
all sampled data points, processed using a low-pass filter.
the distances between the bone and skin surfaces. We used the ultra-sound device to measure
32 locations (Figure 3.6) around the arm and processed the ultra-sound scan images to get tissue
depths at sample locations (Figure 3.4).
Position
Shoulder ———————¿ Wrist
1 2 3 4 5 6 7 8
A
Fat depth / mm 8.033 6.621 2.908 3.9 4.843 6.255 4.79 3.536
Bone distance / mm 24.3 33.71 36.22 23.2 18.703 13.211 9.81 6.9875
B
Fat depth / mm 6.57 9.456 8.557 3.954 4.895 4.791 4.53 4.32
Bone distance / mm 23.46 26.443 27.594 19.75 26.44 21.684 18.86 29.37
C
Fat depth / mm 6.203 10.387 6.9875 0.973 2.6465 3.64 5.05 4.163
Bone distance / mm 31.516 28.12 29.32 9.34 11.85 36.22 20.74 10.91
D
Fat depth / mm 5.67 4.37 7.772 3.22 5.47 5.68 6.15 3.483
Bone distance / mm 10.23 29.53 28.17 20.9 34.184 16.977 16.7677 28.43
Table 3.1: We sampled four sides of the human arm. For each side, we sample 8 points with
similar distances to each other. From the ultrasound image, we measure the fat depth and the
distance between the bone and the skin surface. The location for each label is shown in Figure 3.6
.
We use the distances between the bone and the skin surface at the sample points, and then use
BBW weights to generate a smooth depth field on the arm. An inner surface layer is generated
using the ElTopo library. The tetrahedral mesh is generated using TetGen. We also compare the
shape between the uniform depth mesh and the space-depth mesh (Figure 3.7).
18
Figure 3.6: Poking locations on the right arm: The green dots represent the poking locations.
We sampled 4 sides of the arm and several (8) locations 1-8 at each side, from the shoulder to the
wrist.
Figure 3.7: Tetrahedral mesh generation: The left image shows a uniform-depth tetrahedral
mesh generated using the distance field and the mesh resample method. The right image shows the
spatially-varying tetrahedral mesh generated using ElTopo.
19
Chapter 4
Fast Simulation Using Symmetric Facts
4.1 Cylindrical Coordinate and Axial symmetry
To accelerate the simulation and the optimization process, we want to reduce the 3-D simulation
into a 2-D case. First, we assume that the material we are measuring is isotropic, or at least
locally isotropic. Second, we assume that the object will deform locally in response to the poking,
therefore, we can approximate the object as an infinite plane with a certain uniform depth. Third,
we assume that we always use a cylinder to poke the object on the normal direction. With these
assumptions, the deformation and the contact force will become axis-symmetric. If we use a
cylindrical coordinate system (r,ϕ,z), the deformation and the contact force will be independent
on ϕ. Figure 4.1 illustrates how we simplify a 3D cylinder mesh into a 3D rectangle grid by the
symmetry of the circle. Then the total elastic energy can be written as
ε =
Z
Ω
Ψ(λ
1
(F),λ
2
(F),λ
3
(F))dV
=
Z
2π
0
dϕ
Z
h
0
Z
R
0
Ψ(λ
1
(F),λ
2
(F),λ
3
(F))· r· dr· dz
= 2π
Z
h
0
Z
R
0
Ψ(λ
1
(F),λ
2
(F),λ
3
(F))· r· dr· dz, (4.1)
20
x
z
y
r
φ
(r, φ, 0)
Figure 4.1: Symmetry ofϕ component: The top figure shows the cylindrical coordinate system
we are using and the 3D cylinder mesh we want to simulate. The deformations of the vertices are
the same when the ϕ components are the same. The bottom figure is the 2D rectangle mesh we
simpified using the symmetry of the ϕ-axis.
where F is the deformation gradient, λ
1
, λ
2
and λ
3
are the principle stretches of F and Ψ is the
elastic energy density function for the isotropic elastic material.
21
4.2 Gaussian Points
To effectively compute the elastic energy and the elastic force, we used Gaussian point method [55]
to approximate the integral of the elastic energy, we divide the planeϕ= 0 into a grid mesh. Inside
each rectangle element, we sampled 9 Gauss quadrature points as Figure 4.2 shows. Therefore, the
integral of the elastic energy can be approximated as
ε = 2πw
k∑
k
Ψ(λ(F
k
))· r
k
, (4.2)
where w
k
is the weight of the Gaussian quadrature point and r
k
is the r coordinate of that point.
4.3 2D Finite Element Method on 3D mesh
To derive the deformation gradient F
k
for each Gaussian quadrature point in the cylindrical coor-
dinate, we need to first consider a general case. For any point in the cylinder, the deformation can
be interpolated as
u(x,y,z)=
∑
i
(u
r
i
· e
r
+ u
z
i
· e
z
)φ
i
(r,z), (4.3)
where u
r
i
and u
z
i
are the deformation of the grid vertex on the ϕ = 0 plane, e
r
=(cosϕ,sinϕ,0),
e
z
= (0,0,1) and φ
i
(r,z) and the weight of the sample vertex at the position. Because of the
symmetry, the shape function φ
i
is also independent on ϕ. Here we use the bilinear interpolation
method to define the shape function φ
i
.
Then we can have the deformation gradient as
F(x,y,z)= I
3× 3
+
∂u(x,y,z)
∂(x,y,z)
= I
3× 3
+
∑
i
· ∂[(u
r
i
· e
r
+ u
z
i
· e
z
)· φ
i
]
∂(r,ϕ,z)
· [
∂(x,y,z)
∂(r,ϕ,z)
]
− 1
(4.4)
22
Point X Y Weight
1 − p
3/5
p
3/5 25/81
2 0
p
3/5 40/81
3
p
3/5
p
3/5 25/81
4 − p
3/5 0 40/81
5 0 0 64/81
6 − p
3/5 0 40/81
7 − p
3/5 − p
3/5 25/81
8 0 − p
3/5 40/81
9
p
3/5 − p
3/5 25/81
Figure 4.2: Positions and weights of Gaussian points: The upper figure shows how the 9 Gausian
points are positioned in one rectangle of the 2D rectangular grid mesh. The four circles at the four
corners denote the vertices on the 2D rectangular gird mesh and the 9 black dots denote the 9
Gaussian points inside one rectangle. Suppose the bottom-left of the rectangle is positioned at
(-1,-1) and the top-right of the rectangle is positioned at (1,1), the positions and weights of the 9
Gaussian points are listed in the bottom table. In simulation, the positions of the Gaussian points
will be transformed according the position of the rectangle while the weights always keep the same.
Since all Gaussian quadrature points are on the ϕ = 0 plane, we can simplify the deformation
gradient of Gaussian quadrature points as
F
k
(r
k
,z
k
)= I
3× 3
+
∑
i
u
r
i
∂φ
i
∂r
0 u
r
i
∂φ
i
∂z
0 u
r
i
φ
i
r
k
0
u
z
i
∂φ
i
∂r
0 u
z
i
∂φ
i
∂z
(4.5)
23
In the formula, we can see that the deformation gradient can always be divided into a 2× 2 and a
1× 1 submatrices. Thus when performing the SVD operation on F
k
, λ
3
is the center element and
we only need to actually perform SVD on a 2× 2 matrix to getλ
1
andλ
2
. Whenλ
1
andλ
2
are very
close, there would be a singularity problem in the computation of the stiffness matrix K. To solve
this situation, we perturb the values ofλ
1
andλ
2
to make sure the two values have a difference of
at least 10
− 6
.
The contact force on the grid vertex is
f =
∂ε
∂u
= 2π
∑
k
ω
k
· r
k
· ∂Ψ
∂u
= 2π
∑
k
ω
k
· r
k
· ∂Ψ
∂λ
· ∂λ
∂F
k
· ∂F
k
∂u
, (4.6)
whereλ =[λ
1
,λ
2
,λ
3
]
T
. From Equation 4.5, we know that
∂F
k
∂u
i
T
=
∂F
k
∂u
i
(r,φ,z)
· ∂(r,φ,z)
∂(x,y,z)
=
∂φ
i
∂r
0
∂φ
i
∂z
0
φ
i
r
k
0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0
∂φ
i
∂r
0
∂φ
i
∂z
· (
∂(x,y,z)
∂(r,φ,z)
)
− 1
. (4.7)
24
The term
∂Ψ
∂λ
depends on the property of the material and we will discuss it in the following
Chapter 6. The term
∂λ
∂F
k
is the standard SVD gradient. The stiffness matrix is
K=
∂
2
ε
∂u
2
= 2π
∑
k
ω
k
· r
k
· ∂
2
Ψ
∂u
2
= 2π
∑
k
ω
k
· r
k
· (
∂F
k
∂u
)
T
· ∂
2
Ψ
∂F
2
k
· ∂F
k
∂u
= 2π
∑
k
ω
k
· r
k
(
∂F
k
∂u
)
T
· [(
∂λ
∂F
k
)
T
· ∂
2
Ψ
∂λ
2
· ∂λ
∂F
k
+
3
∑
m=1
∂
2
λ
m
∂F
2
k
· ∂Ψ
∂λ
m
]· ∂F
k
∂u
. (4.8)
Similiar as in Equation 4.6, the term
∂
2
Ψ
∂λ
2
depends on the property of the material and the term
∂
2
λ
m
∂F
2
k
is the standard SVD hessian.
4.4 Validation on 3D object
To test our fitted material property on some non-axis symmetric simulations, we also built a 3D
simulator using Finite Element Method. In this simulation, we used the material property fitted for
the silicone. We also manufactured a Stanford Bunny using the same silicone we manufactured the
cube. We used a square object to poke the bunny. To make the simulation stable, we used penalty-
based method to model the contact and implicit backwards Euler to simulate the dynamics. The
contact force is recorded when the mesh is in the force equilibrium. Figure 4.3 shows the simulated
force is closed to the captured force.
25
Figure 4.3: Validation using 3D simulation on the bunny: The figure show that the simulated
contact forces closely match real contact forces for different poking depths.
Figure 4.4: Experiment set up on the bunny: These two pictures show how we poke the bunny
using a square head in the capturing experiment from two views.
26
Chapter 5
Contact Force Modeling
5.1 Collision Detection
To apply contact force on the human simulation mesh, collision detection is needed to determine
where the contact happens. In this thesis, we do not detect human body self-collision but only
detect inter-collision between the human body and the poking object. Since both objects are rep-
resented by triangle meshes, the naive approach is test intersections for every triangle pair [56,
57]. Due to the large number of triangles in the human mesh, we can speedup the technique using
axis-aligned bounding volume hierarchies [58].
We also use another analytical method to do collision detection that is optimized for our specific
poking application. Because our contact shape is a circle, we know the exact poking center. So we
just need to detect vertices that are within the contact radius by calculating the Euclidean distance
between each vertex and the poking center.
27
5.2 Constraint-based Normal Contact Force
We use constraint-based methods to compute the normal contact force. In our poking experiment,
our contact shape is always a circle. Therefore, we set linear constraints that force all contact
vertices to lie on the same plane as the contact plane
Cq= rhs, (5.1)
where C is a constant matrix, q is the deformation, and rhs is a constant that is related to the poking
depth. Since the deformation of the human skin is quite local due to poking, we rarely observe
contact vertices going in and out of the contact circle under different poking depths. Therefore, to
speedup the simulation, we only perform collision detection once and never change the colliding
vertex set. To combine the linear constraints with the force static equilibrium equations, we use a
Langrange multiplierλ, and the linear system becomes
K C
T
C 0
∆q
λ
=
F
e
rhs− C(q
0
)
, (5.2)
where F
e
is the current elastic force, and K is the gradient of the elastic force (stiffness matrix) at
the current deformation. After solving this equation, we add∆q to the deformation of the object
u. According to force equilibrium, we can compute the normal contact force from the contact
constraints by
f
n
=− C
T
λ. (5.3)
Since we use the same contact vertex set in the poking process, and the normal of the contact
plane never changes, the topology of the system matrix does not need to be modified, which helps
improve the simulation efficiency.
28
5.3 Frictional Contact
We assumed that the friction force between the contact plane is large enough to prevent any tan-
gential sliding from the rest position. In this case, if we know the poking depth, we can prescribe
the positions of all contact vertices. Thus we divide all vertices into three groups, the free vertices
u
F
, the contact vertices u
C
and the attached vertices u
A
. The contact vertices contains the vertices
that are in contact on the top, and the attached vertices that are attached to the ground at the bottom.
Since the deformation of u
C
is known, we need to solve u
F
from a constrained system that
minimizes the total elastic energy
min
u
F
ε(u
F
,u
C
,u
A
). (5.4)
The gradient and hessian of the elastic energy have been derived above. We can use Newton-
Raphson and line search technique to solve the optimization problem. After u
F
is solved, the total
simulated contact force for the poking is
f
S
=
∑
i
f
C
i
. (5.5)
5.4 Line Search
Implicit methods are usually more stable than explicit methods. In a full implicit solver, the stiff-
ness matrix K depends on the deformation ∆x, and the relationship is non-linear in most cases.
To solve such a non-linear system, line search [59] is used so that the non-linear system can be
simplified by iteratively solving a series of linear systems.
In static equilibrium, elastic energy will be minimized under the contact constraints. If we
use the linear force model, the stiffness matrix K is constant, and the deformation change∆q can
be solved using Equation 5.2. However, since our system is non-linear, ∆q may not minimize
the elastic energy. Sometimes, if we apply the deformation change, the elastic energy may even
29
increase. A common solution is to accept the direction of∆q, but find a proper step size so that the
new deformation can minimize the elastic energy
min
s
E(q
0
+ s∆q), (5.6)
where E(q) is the elasticity energy, s is the step size and q
0
is the current deformation. Since q
0
and∆q are known, we can treat this as a one-dimensional optimization problem,
F(s)= E(q
0
+ s∆q). (5.7)
There are two line-search methods we are using. The first one, and the simplest one, is back-
tracking. In the backtracking algorithm, we first try s= 1 and evaluate the elasticity energy. Then
we keep reducing s by a factorτ< 1, and detect the minimum elasticity energy. We keep reducing
s until Armijo–Goldstein condition [60] (Equation 5.8) is fulfilled, or after a maximal number of
iterations
F(0)− F(s)> c(∆q·▽ E(q)),, (5.8)
where c is a scalar parameter. However, this simple backtracking method may fail. When it fails,
for all s> 0 searched, the elastic energy F(s) is always larger than the original elastic energy
F(0). This happens because we only discretely sampled the step size s. Brent’s method is better at
searching a proper step size. The main idea of Brent’s method is to assume that the energy function
can be fitted by a quadratic function locally.
To integrate the static equilibrium equations using a linear search method, we need to enforce
the constraints once the poking depth changes. We modify the positions of vertices on the contact
area by moving them along the normal direction of the contact plane, until their positions are
30
exactly on the plane (Cq
0
= rhs). Once the constraints are satisfied, the following iterations of line
search will always guarantee good constraints,
C· (q
0
+ s∆q)= Cq
0
+ s· (C∆q)= rhs. (5.9)
5.5 Reduced Simulation
Since the deformation of the human skin is very local, most of the simulation does not contribute
to the physical simulation, but remains in its rest shape. Thus, to accelerate the simulation compu-
tation in the haptic cycle, we can virtually remove the geometry that is far from the poking center,
and only simulate a local simulation mesh. In our experiments, we observe that an active radius of
1cm is enough to keep the simulation realistic.
31
Chapter 6
Material Optimization
6.1 Spline-based material
6.1.1 Valanis-Landel Material
Many traditional works only capture the material property in the linear region. However, we ob-
serve that with even with a shallow poking (20% of the total sample depth), the deformation is
large enough that the material is likely to behave non-linearly. Figure 6.1 shows the distribution
of the deformation gradient principle stretchλ when the poking depth is 20% of the total sample
depth.
To capture the non-linearity of the material, we assume that the isotropic material can be mod-
elled as Valanis-Landel’s material [valanis:1967:SFH]
Ψ= f(λ
1
)+ f(λ
2
)+ f(λ
3
)+ g(λ
1
λ
2
)
+g(λ
2
λ
3
)+ g(λ
3
λ
1
)+ h(λ
1
λ
2
λ
3
), (6.1)
32
Figure 6.1: The figure shows the distribution of the λ values: The upper figure is the histogram
of{λ
1
,λ
2
,λ
3
} values. The middle figure is the histogram of {λ
1
λ
2
,λ
2
λ
3
,λ
3
λ
1
} values. The bottom
figure is the histogram of {λ
1
λ
2
λ
3
} values.
33
where f , g and h are scalar functions. To get the complete formulas for Equation 4.6 and Equa-
tion 4.8 in Section 4.3, we have
∂Ψ
∂λ
i
= f
′
(λ
i
)+λ
j
g
′
(λ
i
λ
j
)+λ
k
g
′
(λ
i
λ
k
)+λ
j
λ
k
h
′
(λ
i
λ
j
λ
k
), (6.2)
∂
2
Ψ
∂λ
2
i
= f
′′
(λ
i
)+λ
2
j
g
′′
(λ
i
λ
j
)+λ
2
k
g
′′
(λ
i
λ
k
)+(λ
j
λ
k
)
2
h
′′
(λ
i
λ
j
λ
k
), (6.3)
∂
2
Ψ
∂λ
i
∂λ
j
=λ
i
λ
j
g
′′
(λ
i
λ
j
)+λ
i
λ
j
λ
2
k
h
′′
(λ
i
λ
j
λ
k
), (6.4)
where{i, j,k}={1,2,3}.
6.1.2 Natural Cubic Spline
We used natural cubic splines to model these three functions. For example, to fit the function f
using a spline, we have a set of control points (x
k
,y
k
) and{x
k
}
0≤ k≤ n
is sorted ascendingly. The
natural cubic spline function can be written as
f(x)= f
k
(x)= a
k
(x− x
k
)
3
+ b
k
(x
k+1
− x)
3
+ c
k
(x− x
k
)+ d
k
(x
k+1
− x),
x
k
≤ x≤ x
k+1
,i= 0...n− 1. (6.5)
In addition, natural cubic spline should belong to C-2 curves, that is, f(x), f
′
(x) and f
′′
(x)
should be continuous, which satisfies
f
k
(x
k+1
)= f
k+1
(x
k+1
), (6.6)
f
′
k
(x
k+1
)= f
′
k+1
(x
k+1
), (6.7)
f
′′
k
(x
k+1
)= f
′′
k+1
(x
k+1
). (6.8)
34
If we denote z
k
= f
′′
k
(x
k
),k= 0...n− 1 and z
n
= f
′′
n− 1
(x
n
), we always have z
0
= z
n
= 0 as the
boundary condition of the natural cubic spline and z
1
...z
n− 1
can be solved by a tridiagonal system
v
1
h
1
h
1
v
2
h
1
h
2
v
3
h
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. h
n− 2
h
n− 2
v
n− 1
·
z
1
z
2
z
3
.
.
.
z
n− 2
z
n− 1
=
u
1
u
2
u
3
.
.
.
u
n− 2
u
n− 1
, (6.9)
where h
k
= x
k+1
− x
k
, v
k
= 2(h
k− 1
+ h
k
) and u
k
= 6((y
k+1
− y
k
)/h
k
− (y
k
− y
k− 1
)/h
k− 1
). After
solving for{z
k
}, we can compute all coefficients of the natural cubic spline by
a
k
=
z
k+1
6h
k
, (6.10)
b
k
=
z
k
6h
k+1
, (6.11)
c
k
=
y
k+1
h
k
− z
k+1
h
k
6
, (6.12)
d
k
=
y
k
h
k
− z
k
h
k
6
. (6.13)
We sampled{x
k
} log-uniformly in a range [λ
min
,λ
max
]. A typical λ
min
can be 2
− 1
and a typical
λ
max
can be 2
1
as Figure 6.1 shows. We can adjust the values of y
k
to control the spline and
the material. Similarly, for the spline functions g and h the x
k
of the sample points are sampled
log-uniformly in the range shown in Figure 6.1.
6.2 Conjugate Gradient Optimization
The goal of the optimization is to find the correct values for the three natural cubic spline functions
so that the difference between the simulated contact forces and the captured contact forces can be
minimized. We define the objective function ε as
35
min
µ
ε =
∑
i
|| f
C
i
(µ,u
F
i
(µ))− f
E
i
||
2
, (6.14)
whereµ is the control point positions of the splines in the material energy density function and
u
F
is the deformation of the simulated mesh. The gradient is
dε
dµ
=
∑
i
2( f
C
i
− f
E
i
)
T
(
∂ f
C
i
∂µ
+ K
CF
i
· ∂u
F
i
∂µ
) (6.15)
When the simulation is in the static equilibrium, we have
f
F
(µ,u
F
(µ))= 0. (6.16)
If we compute the derivative ofµ and apply the chain rule, we can have
d f
F
dµ
=
∂ f
F
∂µ
+ K
FF
· ∂u
F
∂µ
= 0. (6.17)
To calculate∂ f
C
/∂µ in Equation 6.15 and∂ f
F
/∂µ in Equation 6.17, we can compute the derivate
of Equation 4.6 that
∂ f
∂µ
= 2π
∑
k
ω
k
· r
k
· (
∂F
k
∂u
)
T
· (
∂λ
∂F
k
)
T
· ∂
2
Ψ
∂λ∂µ
. (6.18)
The term(
∂F
k
∂u
)
T
and(
∂λ
∂F
k
)
T
has been defined in Section 4.3 and from Equation 6.2, we can derive
∂
2
Ψ
∂λ∂µ
=
∂ f
′
(λ
i
)
∂µ
+λ
j
∂g
′
(λ
i
λ
j
)
∂µ
+λ
k
∂g
′
(λ
i
λ
k
)
∂µ
+λ
j
λ
k
∂h
′
(λ
i
λ
j
λ
k
)
∂µ
(6.19)
36
Since∂ f
′
/∂µ,∂g
′
/∂µ and∂h
′
/∂µ are very similar, we take∂ f
′
/∂µ as an example that
∂ f
′
∂µ
=
∂ f
′
k
∂µ
= 3
∂a
k
∂µ
(x− x
k
)
2
− 3
∂b
k
∂µ
(x
k+1
− x)
2
+
∂c
k
∂µ
− ∂d
k
∂µ
,
x
k
≤ x≤ x
k+1
,i= 0...n− 1. (6.20)
From Equations 6.10 to 6.13, we can derive that
∂a
k
∂µ
=
∂z
k+1
∂µ
1
6h
k
, (6.21)
∂b
k
∂µ
=
∂z
k
∂µ
1
6h
k+1
, (6.22)
∂c
k
∂µ
=
∂y
k+1
∂µ
1
h
k
− ∂z
k+1
∂µ
h
k
6
, (6.23)
∂d
k
∂µ
=
∂y
k
∂µ
1
h
k
− ∂z
k
∂µ
h
k
6
, (6.24)
∂y
∂µ
= I
(n+1)× (n+1)
. (6.25)
From Equation 6.9, we can get
∂z
∂µ
=
v
1
h
1
h
1
v
2
h
1
h
2
v
3
h
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. h
n− 2
h
n− 2
v
n− 1
− 1
· 6
1
h
0
− 1
h
0
− 1
h
1
1
h
1
.
.
.
.
.
.
.
.
.
1
h
n− 2
− 1
h
n− 2
− 1
h
n− 1
)
1
h
n− 1
(6.26)
With the gradient of the objective function, we can use conjugate gradient optimizing method to fit
the spline material. To make the optimizer converge faster, we will first fit the Young’s modulus.
37
6.3 Results
We used two types of man-made material, foam and rubber to demonstrate our capturing and
optimizing method. For each type, we also manufactured a softer one and a stiffer one. Our
optimization is able to fit the material with one force-indentation profile and we can validate the
material property by another poking shape that is not used in the optimization. We also validate the
material property on a 3D mesh and a non-circular contact shape. We also show that our method
can be adapted to estimate the height of the material sample for some applications. All examples
were computed on a 3.00 GHz INtel Xeon i7 CPU E5-2687W v4 processor. We used all the 48
threads on our CPU to compute the static equilibrium for each simulation iteration. Table 6.1
shows the change ofε in Equation 6.14 after our optimization method and the time it spent on the
optimization process. In Table 6.1, the initialε is the value ofε after fitting the Young’s modulus
and the final ε is the value ofε after the non-linearity of material is fitting using our optimization
method. The #iterations also refer to the number of iterations for fitting the non-linearity of the
material property using natural cubic splines. It always only takes one iteration to fit the Young’s
modulus as fitting the Young’s modulus is an unconstrained quadratic optimization
Figure 6.2: Man-made material samples: The left figure is the foam material sample we manu-
factured and the right figure is the rubber material sample.
38
Material name initialε final ε #iterations time
Soma Foama 25 0.026 0.006 44 127 mins
Soma Foama 15 0.011 0.002 54 23 mins
Softer Dragon Skin 0.181 0.038 14 34 mins
Stiffer Dragon Skin 4.805 0.299 10 58 mins
Table 6.1: Optimization results, numbers of optimizing iterations and timings of the four
material samples.
6.3.1 Foam
We manufactured two samples that are made by Soma Foama 25 and Soma Foama 15. For each
sample, we made a 10cm x 10cm x 1cm cube. The appearance of this material sample is shown
in the left of Figure 6.2. We used some cylinders whose radii are 1.03mm, 2.06mm, 3.01mm,
4.05mm and 5.06mm to poke the material sample and measured the contact forces. We used
the one cylinder to fit the material property and then use other cylinders to validate the material
property. Figure 6.3 and Figure 6.4 indicate our method can fit the material propery of these foam
samples with one cylinder poking and the validations on other cylinder shapes are also good. For
Soma Foama 25, the value of
√
ε is reduced by 53%, from 0.160 to 0.076, and for Soma Foama
15, the value of
√
ε is reduced by 56%, from 0.104 to 0.046.
Figure 6.3: Fitting and validation results of the Soma Foama 25 material: These two figures
are the fitting and validation results of the Soma Foama 25 material. We used the cylinder with
5.06mm radius to fit the material property of Soma Foama 25. The left figure is the captured and
simulated force-indentation curves for the initial guessed material property. The right figure is
the fitting and validation results of the final converged material property after optimization. The
dashed lines are simulated relations between the contact forces and the poking depths.
39
Figure 6.4: Fitting and validation results of the Soma Foama 15 material: These two figures
are the fitting and validation results of the Soma Foama 15 material. We used the cylinder with
3.01mm radius to fit the material property of Soma Foama 15. The left figure is the captured and
simulated force-indentation curves for the initial guessed material property. The right figure is the
fitting and validation results of the optimized material property.
6.3.2 Silicone
In the example, we manufactured two 10cm x 10cm x 0.5cm samples using two versions of Dragon
Skin FX Pro, one is softer and one is stiffer. Generally, the rubber behaves stiffner than the foam
mentioned in the previous section. The appearance of this material sample is shown in the right
of Figure 6.2. We did the same measurement experiment by using one or two cylinders to fit the
material and the data from other radii to validate the captured material property. Figure 6.5 and
Figure 6.6 indicate our simulated forces can approximate the captured forces on different sizes of
contact areas. For the softer Dragon Skin FX Pro, the value of
√
ε is reduced by 54%, from 0.426
to 0.194, and for the stiffer Dragon Skin FX Pro, the value of
√
εn is reduced by 75%, from 2.192
to 0.547.
6.3.3 Convergence on mesh subdivision
To validate our simulation method on the 2D mesh, we want to compare the simulated force com-
puted from our method to a classic finite element method on the 3D tetrahedral mesh. Since we ran
our 2D simulation on a grid mesh, we made a tetrahedral mesh for a cuboid and poke the cuboid
at the center of the top surface. We observed that the subdivision level of the mesh will affect the
40
Figure 6.5: Fitting and validation results of the softer Dragon Skin FX Pro material: These
two figures are the fitting and validation results of the softer Dragon Skin FX Pro material. We used
the cylinder with 4.05mm radius to fit the material property of softer Dragon Skin FX Pro. The
left figure is the captured and simulated force-indentation curves for the initial guessed material
property. The right figure is the fitting and validation results of the optimized material property.
Figure 6.6: Fitting and validation results of the stiffer Dragon Skin FX Pro material: These
two figures are the fitting and validation results of the stiffer Dragon Skin FX Pro material. We
used the cylinder with 3.01mm and 5.06mm radii to fit the material property of stiffer Dragon
Skin FX Pro. The left figure is the captured and simulated force-indentation curves for the initial
guessed material property. The right figure is the fitting and validation results of the optimized
material property.
41
simulation force so we only need to compare the converged contact forces from these two methods.
Figure 6.7 shows as the subdivision levels of those two meshes go up, the simulation forces will
converge to a similar value.
Figure 6.7: Validation of 2D simulation method: The figure show the simulated contact forces
on a 2D mesh and a 3D mesh converge to a similar value as the subdivision levels go up.
6.3.4 Extension: Estimate the depth of the material sample
For some applications, the total depth of the layer to be measured may also unknown, such as
capturing the human tissue property in vivo. In the example, we used the captured contact force
from all poking radii to fit the material property and we demonstrate that this method can infer the
depth of the material sample. We built several optimizers with different guessing of the sample
depths from 3mm to 9mm. From Figure 6.8, we can see that the 5mm-depth got the lowest dif-
ference between the simulation and reality. 5mm is the ground truth depth of our manufactured
material sample. Figure 6.9 shows the simulated force-indentation curves are close to captured
force-indentation curves for all poking radii.
42
Figure 6.8: Fitting results on different height depths: This figure displays how the value of the
objective function decreases over the conjugate gradient optimizing iterations for different depth
guessing. In the figure, the 5mm depth got the lowest objective function value and this is the same
as the groundtruth.
Figure 6.9: Fitting results of the best depth guess: The left figure is the initial guess of the
material property and the right figure is the final converged result. The solid lines are captured re-
lations between the contact forces and the poking depths. The dashed lines are simulated relations
between the contact forces and the poking depths.
43
6.3.5 Measurement on Human
We got an IRB approval from the university so we are also permitted to measure the human body.
We did the poking experiment by the customized hand-held poking device and optimized the ma-
terial properties of the human tissue using the optimization method described above. The fitting
results shown in Figure 6.10 demonstrate our method can work on the human body and the simu-
lated force-indentation curve can match the reality.
44
Figure 6.10: Optimization results of human tissue material: These two figures show our simula-
tion result using the fitting material property can qualitatively match the experimental data captured
by our poking device. The blue circles denote the experimental data captured by our poking device
and the red curve gives our simulation result. The depth of tissue at Position-009 (dorsum back) is
20mm and the depth of tissue at Position-041 (buttock) is 37mm.
45
Chapter 7
Application: Designing Ergonomic Furniture
7.1 Introduction of Ergonomic Design
Humans routinely interact with supporting surfaces – we sit on them (chairs), lie on them (beds
and sofas), stand on them (footwear), and work all day with our hands on them (mice and key-
boards). It is important for our health and well-being to shape these surfaces to be comfortable
and ergonomic [61]. Traditionally, such contoured surfaces are shaped by craftspersons who draw
upon generations of design knowledge to form shapes that are likely to be comfortable for a wide
range of people, often with human testers in a feedback loop to refine the design. However, the
traditional design process does not scale easily to customized support for individual people: it is
expensive, time-consuming, requires extensive interaction and iteration with the target user, and is
usually not based on fine-grained biomechanical modeling of the human body. At the same time,
customized fabrication methods such as 3D printing have become more and more accessible to the
general public [62]. Yet as long as customized design remains difficult, such fabrication methods
remain grossly under-utilized.
A promising path to closing this gap is to develop computational methods that can automat-
ically synthesize supporting surfaces whose contact ergonomics are personalized to a specific
person. In the geometric shape design literature, shapes have been optimized for various objec-
tives, typically aesthetics, visual plausibility and crude mechanical function [63]. A few methods
46
Initial guess Iteration 10 Iteration 20 Iteration 26 (nal)
vertex
displacements
contact
pressures
Human in chair
Figure 7.1: Optimizing the shape of the supporting surface (chair) to maximize ergonomics
of contact against a human. We model the deformations of the human soft tissue using the Finite
Element Method. Given any shape of the supporting surface, we calculate the contact pressures
using physically based simulation, and score the ergonomics using an ergonomic score function
that penalizes localized contact forces. Our gradient-based optimizer then optimizes the supporting
surface chair shape to minimize the score.
have tried to explore human interaction with the object, using non-physics-based data-driven mod-
els [64–67], or a precomputed approximate pressure map for an inaccurate contact distribution (a
negative image mold) that is used only as a fixed importance map for geometric support synthe-
sis [68, 69].
In this paper, we propose an automatic, physics-based technique to optimize supporting sur-
faces for specific people, with known geometric (bodyshape) and elastic (tissue) properties. Our
method is based on the computational contact mechanics between the human body and the sup-
porting surface. While humans are obviously biomechanically complex with intricate muscle and
fat geometry, we make our optimization tractable by modeling human body parts as consisting of
a deformable layer (with measured material properties) surrounding a rigid core. The supporting
surface is modeled as rigid. Based on ergonomics literature, we define a meaningful and tractable
ergonomic cost function that scores the pressure distribution exerted by the support on the hu-
man. We estimate this contact pressure using large-deformation FEM deformable object elasticity.
Importantly, and specifically for our optimization use-case, we design our cost function to be dif-
ferentiable with respect to the parameters defining the geometry of the supporting surface. Hence,
the latter can be efficiently optimized using gradient descent in order to minimize the cost.
47
We apply our pipeline to design optimal contact surfaces of chairs, shoe insoles, and computer
mice. We use a custom acquisition pipeline to gather elasticity maps of the human body parts
in contact with these surfaces, and use these material properties in our simulation. We validate
our results against baselines, as well as in an experiment with an actual fabricated surface, to
show that our pipeline is accurate and robust. In practice today, ergonomic surfaces are rarely
designed completely using computer models; but instead real prototypes are manufactured, real
people evaluate them by interacting with them, and then the design is iterated. Our work can
shorten the design loop because it can provide a good initial design for this practical subsequent
design process, and/or makes it possible to virtually prototype the design iterations.
Our technical contributions include a stable, easily computable, and differentiable ergonomic
function to measure the ergonomics between a rigid object and a nonlinear elastic deformable
object. We also give a phenomenological simulation model of the human subskin soft tissue that
is sufficiently complex to be able to express realistic tissue deformations, yet not overly complex
so that it can still be simulated efficiently and iterated in a design loop. The model consists of
a rigid core surrounded by and coupled to a layer of soft nonlinear elastic tissue. We give a
stable method to simulate the resulting coupled rigid-deformable human mechanical system in
contact with an external object (a supporting surface), with spatially varying contact regions and
pressures. Finally, we give an efficient method to optimize the shape of the object to maximize the
ergonomic comfort of the contact (Figure 7.1). To the best of our knowledge, our work is the first
to successfully incorporate an ergonomic objective, nonlinear large deformation FEM elasticity,
and distributed unilateral contact into ergonomic shape design of man-made objects. Distributed
contact is defined as contact with (potentially) several simultaneous contact sites, each distributed
over a non-zero surface area [70]. Contact is unilateral if it can only push but not pull, i.e., it
has same characteristics as real-world contact. Note that distributed unilateral contact is the most
general and most realistic contact between two 3D objects, and is significantly harder to process
and optimize than single-point contacts, bilateral attachments, or permanent conforming contact.
Furthermore, we also contribute a new experimental method to measure real distributed contact
48
pressures between two curved surfaces. We demonstrate how this can be reliably achieved using
standard, inexpensive resistance-based pressure sensors. We use our method to experimentally
validate that our optimized pressure distributions qualitatively match real-world pressures.
7.2 Related Work
Shape design: There are various functional and aesthetic goals which shape design has attempted
to address, and many methods explored to reach these goals. Topology optimization aims to add
or remove material from a structure to optimize its properties. Often this is done by optimizing
a dense grid of voxels to maximize the stiffness or robustness of a structure while minimizing
its weight, in the presence of either static [71, 72] or stochastic [73, 74] external forces. Other
methods optimize a level-set representation of the geometry [75]. While capable of producing
intricate structures, topology optimization is computationally expensive, and is not generally used
with nonlinear material models.
Other shape optimization methods keep the user in the loop and rely on precomputation to
generate suitable exploration spaces. Shugrina et al. [76] adaptively sample the design space of
a parametric object, mapping out which parts of the space satisfy constraints (such as manufac-
turability). Then, during runtime, users can quickly explore the design space, while ensuring they
remain in the valid constraint region. Yumer et al. [77] allow high-level geometrical editing by
relying on crowd-sourced comparisons of semantic attributes. This data is used to learn a mapping
between semantic attributes and geometry. While capable of producing a wide variety of useful
shapes, the amount of pre-computation that would be required to capture the fine variation between
different humans is immense. Instead we rely on optimization based on human models.
Traditionally, shape optimization methods operate with static constraints (e.g. Umetani et
al. [78], Wang et al. [79]), although some methods consider dynamic properties of the object
being optimized. There exist ways to design spinning [80], standing [81], flying [82] and float-
ing objects [83]. Dynamics-Aware Coarsening [84] can speed up the design of compliant (large
49
deformation) objects with dynamic motion, such as jumping. In our work, we focus on optimizing
the human ergonomics of contact, as opposed to the object’s dynamic or static properties. Several
publications in computational fabrication design geometry or material properties under equilib-
rium constraints [85, 86], typically using the implicit function theorem, as in our work. However,
previous work has not incorporated unilateral contact, or proposed a methodology to validate such
design by measuring real contact pressures on curved surfaces. Recently, Montes et al. [87] opti-
mized the patterns of tight-fitting clothing so that the fit to the human is as ergonomic as possible.
Because they modeled tight-fitting clothing, they could assume that the contact is permanent and
conforming (essentially a bilateral attachment). In contrast, we model distributed unilateral con-
tact between the human and the external object. As such, our method can optimize the ergonomics
of contact between the human and arbitrarily shaped, not necessarily tight-fitting objects. Like in
the real world, the specific contact sites and their sizes are not known a priori, but are determined
by our ergonomic score optimizer together with the human tissue deformations and the elasto-
static equilibrium of the human against contact forces and gravity. In the medical community,
researchers have evaluated the human body contact ergonomics for various body parts, such as feet
and buttocks [88, 89]. They modeled the individual musculo-skeletal anatomy [88, 89], which is
difficult to generate due to the required medical imaging, organ templates, and anatomy transfer
to the subject. Furthermore, unlike us, they did not aim to optimize the shape of external objects
in contact with the human, but rather merely calculated the contact strains for given fixed external
object shapes.
Furniture comfort evaluation: Considering how much time humans spend interacting with
furniture, evaluating furniture comfort is important for designers. There are general ergonomic
guidelines on how to choose and adjust chairs, e.g., from the US General Services Administra-
tion [90], but designing individual chairs to maximize comfort is still an active research area.
Different computational measures of comfort have been used, such as joint angles [91] and inter-
face pressure between the human and chair [92]. We choose to use interface pressure to quantify
model comfort, as it correlates well with user comfort ratings [92, 93]. Savonnet et al. [94] provide
50
an extensive review of various deformable body models (from twenty-seven different papers) that
have been created to predict the comfort of chair models or the development of pressure sores; they
conclude that “These models share a lack of validation and generally make little allowance for an-
thropometric diversity.” Generally, models aiming to predict pressure sores differentiate between
internal soft tissue (i.e., fat, muscle, and skin were modeled separately), while models to predict
seating comfort generally group the internal tissues together. Multiple material models have been
used to model the soft tissue deformation, including Mooney-Rivlin, neo-Hookean, Ogden, and
linear elasticity. We build our human model as a rigid-core (bones) surrounded by a deformable
layer (all soft tissue), using the neo-Hookean material model with compression resistance. The
depth of the soft tissue, and its material properties, are calibrated based on force measurements
from a user, in a process similar to [95].
Layered models are useful for animation of characters [96], where inner layers can control
large-scale behavior, and outer layers provide localized deformation due to contacts or other forces.
Pauly and Pai [97] investigated techniques for modeling contact between quasi-rigid objects. Ter-
zopoulos and Witkin [98] decompose deformation into reference (rigid) components and displace-
ment components. The reference component evolves according to rigid-body dynamics, allowing
the model to handle large deformations with linear elasticity. This was extended to deformable
reference shapes [99]. Implicit surfaces have been used for the outer layers [100]. A series of
papers [101, 102] present a fully coupled elastodynamic system of a rigid core surrounded by
a single layer of deformable tetrahedra. Recently, layered models were extended to articulated
“cores” [103, 104]. However, prior work has not used layered models for shape optimization
involving contact. In our problem, we only require the final converged state of the contact simu-
lation, and its gradient. We present a layered elastostatic model, which we found to be faster and
differentially more stable than elastodynamic models.
Ergonomic furniture: Previous studies have investigated the design of furniture based on
ergonomics. Saul et al. [105] developed a system which extruded 2D sketches into 3D furniture,
then physically simulated a rag doll (composed of rigid links) sitting on the shapes to check gross
51
characteristics such as sizing and balance of the furniture. The results from this quick simulation
can be used to adjust the furniture. This method is incapable of capturing fine differences between
humans, or pressure differences along the body. Zheng et al. [106] use a rigid skeleton model to
get rough body part lengths of users, and use ergonomics to infer contact configurations between
parts (e.g., the model’s back should be in contact with the seat back). The method then adjusts
the various rigid parts of the chair to support the model, and allows users to interactively edit
the model. This method can quickly show various chairs adapted for the target pose, but cannot
capture important interface pressure distributions. Perhaps most similar to our work is that of
Leimer et al. [68]. This method uses a rough estimate of contact pressure, akin to what would
happen if a body model was immersed in a negative image mold (so that every downward-facing
vertex was supported). This pressure distribution is used to fit a Catmull-Clark subdivision surface,
ensuring areas of high pressure are supported. While capable of automatically producing shapes
supporting various body poses, pressure induced by a negative image mold can be quite different
from that of real supporting surface and this is important for ergonomic design validation. In
recent work, Leimer et al. [69] gave an optimization approach whereby desirable contact pressure
distributions are computed by maximizing their smoothness and keeping the net forces and torques
on the human skeleton links as close to equilibrium as possible. Although plausible, their approach
relies heavily on geometric skinning-based deformation, which may be inaccurate. Unlike their
work, we model human tissue elastic properties and simulate the human in contact against the
supporting surface using computational contact mechanics of deformable objects. An important
key difference to Leimer’s work is that they require knowing a priori what human vertices are in
contact. However, real contact between a human and an external object is complex and spatially
varying (distributed). Simply postulating known contact areas leads to substantial innaccuracies,
as we experimentally show in Section 7.5. Our method determines the entire distribution of the
contact pressures automatically using contact mechanics in each iteration of shape optimization,
which leads to significantly more realistic contact pressure distributions. We compare to both
publications [68, 69] in Section 7.5.
52
7.3 Simulation Model
Our goal is to optimize the surface geometry of an external object (a supporting surface) to max-
imize the ergonomics of the contact against the human. The supporting surface is modeled as a
triangle mesh. The human body is modeled as a FEM deformable layer of varying thickness, sur-
rounding a rigid core that models the inner hard components such as bones. The local thickness
and stiffness of the deformable layer, as well as its surface geometry, is measured via actual human
studies, approved by the medical ethics committee (Institutional Review Board) at our institution.
Mass is distributed across the body according to realistic estimates. The deformable layer is mod-
eled as a tetrahedral mesh between the outer and inner surface of the body, and the volume inside
the inner surface (the “core”) is modeled as a rigid body. While this simplifies human anatomy,
such phenomenological models have been recently demonstrated to offer a good tradeoff between
the modeling complexity and invested computational costs, for human body contact mechanics
simulations [95].
Given a human body and a supporting surface, modeled as above, we define an ergonomic
cost function (§ 7.4.1) that estimates how uncomfortable the contact between the human and the
supporting surface is. Specifically, we follow the standard observation [92, 93] that high contact
pressures on the body are uncomfortable, and define our cost function to penalize such high contact
pressures. For each specific shape of the supporting surface, the spatially varying contact pressures
are computed using an FEM simulation (§ 7.3.3) that couples the rigid and deformable components
of the body model, and resolves the human contact against the supporting surface. In Section 7.4,
we then optimize the shape of the supporting surface, by calculating the gradients of the ergonomic
function with respect to the shape of the supporting surface.
7.3.1 Modeling the human body
Figure 7.2 gives an overview of how we model the human body in contact with a man-made
rigid object (the supporting surface). Our IRB-approved study acquired the human body shape
53
Figure 7.2: Modeling and optimizing the contact between a human and a man-made object.
The man-made object is fixed in space during each contact simulation, whereas the human body
is subjected to gravity. Our optimizer explores the space of man-made object shapes to maximize
contact ergonomics.
by scanning real persons using the Artec Eva scanner [107]. For our chair example, we scanned
the lower back, buttocks, torso and thigh. We approximated the sitting pose by asking the subject
to position themselves on all fours, with pillows supporting the torso. For our foot example, we
scanned the human foot. For our mouse example, we used a generic hand mesh purchased on-
line [Turbosquid]. Note that it is sufficient to scan and model the larger region of interest of the
human body; it is not necessary to scan the entire human body.
We measured human body elastic stiffnesses (Young’s modulus) and deformable layer depths
on contact regions using the technique described in the above chapters. We captured 60 samples
across the lower back, buttocks and thighs on the subject, in a pose and muscle activations approx-
imately similar to sitting. Subjects were fully clothed wearing tight-fitting pants. This exposed as
little privacy of the subject as possible. For each example, we measured the human in a pose simi-
lar to the one in the example, and with similar muscle activations. The measured Young’s moduli
ranged from 5,000 Pa to 15,000 Pa. After that, we interpolated captured soft layer stiffnesses and
depths to the whole mesh. We then used ElTopo [108] to push the surface mesh inward with the
54
interpolated depths, generating the inner surface of our human model. Finally we tetrahedralized
the volume between the outer and inner surface using TetGen [53] to create the tetrahedral mesh
for the FEM simulation.
7.3.2 Modeling supporting surfaces
Figure 7.3: Contact between the human body and a supporting surface. Left: human sitting
in a chair. Middle: foot contacting a shoe sole under human body gravity. Right: hand grasping a
mouse.
We model the chair by starting from a 0.5m× 0.2m× 0.5m cube, which we beveled into a
round arc in the knee region to avoid a sharp edge contact below the knee. Although this “chair”
shape is simplistic, it reasonably captures the essence of the contact between the human and a chair.
Our neutral shoe sole is flat, and was obtained by resizing a commercial shoe sole model purchased
online [Turbosquid]. The neutral mouse model was created by modeling two non-concentric
spheres, followed by a CSG union operation and a few iterations of Taubin smoothing [109]. This
created a reasonable initial “mouse” shape, reminiscent of ergonomic vertical mice such as the
Logitech MX Vertical, to serve as our baseline neutral model. All triangle meshes were sufficiently
subdivided to provide enough degrees of freedom for shape deformation. Figure 7.3 gives our
examples.
55
7.3.3 Nonlinear elastostatics coupled to the rigid core
The human is simulated via two-way coupling between a 6-DOF rigid object (the “rigid core”)
and nonlinear large-deformation elasticity of the deformable layer, modeled as a tetrahedral mesh
using FEM. The coupled rigid core + deformable system is subjected to gravity, whereas the rigid
supporting surface is not simulated and is pinned in space. This causes the FEM mesh to be in
distributed contact against the rigid supporting surface (Figure 7.2). Our model is designed so that
the static contact equilibrium under gravity is differentiable with respect to changing the shape of
the supporting surface.
The deformable layer is modeled using three-dimensional FEM nonlinear elasticity; we use
invertible isotropic hyperelastic neo-Hookean energy [1] with additional energy terms added for
volume preservation [110]. Let M∈R
3n× 3n
be the mass matrix, u∈R
3n
the vector of mesh vertex
displacements, and f
int
= f
int
(u)∈R
3n
and K= K(u)= d f
int
(u)/du∈R
3n× 3n
the internal elastic
forces and the tangent stiffness matrix, respectively. Assemble contact forces acting on mesh
vertices into f
contact
(u)∈R
3n
. Contact forces only apply to vertices on the exterior surface of the
tet mesh. We furthermore divide the tet mesh vertices into 2 groups: vertices that are attached
to the rigid core (by fixed constraints; the “Attached vertices”, denoted by A) and the rest of the
vertices (the “Free vertices”, denoted by F). We split the forces and matrices correspondingly into
the A and F components.
Consider an infinitesimally small displacement δu and infinitesimally small rigid body motion
of the coreδq=[(δx)
T
(δθ)
T
]
T
, where x is the position of the center of mass of the rigid core (in
world coordinates), and δθ ∈R
3
is an infinitesimally small rotation by angle ||δθ|| around axis
δθ/||δθ||. Because the attached vertices are constrained to the rigid object, we haveδu
A
= Bδq,
where
B=
I − e
h
1
.
.
.
.
.
.
I − e
h
a
. (7.1)
56
Here, h
i
is the vector from the rigid core’s center of mass to the current position of attached vertex
i, and integer a denotes the number of attached vertices. For a vector x∈R
3
, ˜ x denotes the skew-
symmetric matrix with the property ˜ xy= x× y for all y∈R
3
.
We can now state our elastostatic equilibrium equations. Both the elastic object and the rigid
object must be in equilibrium under gravity, contact to the supporting surface, and attachment
forces f
attach
between the rigid core and the deformable object:
K
FF
δu
F
+ K
FA
Bδq+ f
F
int
+ f
F
contact
= M
F
ˆ g, (7.2)
K
AF
δu
F
+ K
AA
Bδq+ f
A
int
+ f
A
contact
+ f
attach
= M
A
ˆ g, (7.3)
B
T
f
attach
+ w= 0, for (7.4)
ˆ g=[0 − g 0 0 − g 0 ... 0 − g 0]
T
∈R
3n
, and (7.5)
w=[F
T
τ
T
]
T
+[0 − mg 0 0 0 0]
T
. (7.6)
Here, ˆ g is the gravity acceleration vector (we use g=9.81), m is the mass of the rigid core, w∈R
6
contains the external forces and torques acting on the rigid body (the “wrench”), consisting of the
external force F∈R
3
and external torque (around the center of mass) τ∈R
3
acting on the rigid
core, plus gravity. Note that the torque of gravity is zero. In our chair and mouse examples, there
are no external forces except gravity. In the foot example, we model the force and torque of the rest
of the human body onto the foot, as explained in Section 7.5.2. Equations 7.2, 7.3 and 7.4 establish
the static equilibrium of the free vertices, attached vertices and the rigid core, respectively. By
57
eliminating f
attach
using Equation 7.3, we obtain the following linear system for an incremental
change in the elastostatic equilibrium:
I
B
T
K
I
B
δu
F
δq
=
I
B
T
M ˆ g− f
int
− f
contact
+
0
w
, (7.7)
δu=
I
B
δu
F
. (7.8)
7.3.4 Computing the static equilibrium
We perform a Newton-Raphson iteration where we repeatedly solve for δu and δq, followed by
a line search along the search direction(δu,δq). We are solving a difficult nonlinear contact me-
chanics problem, and as such the presence of the line search is critical. Furthermore, it is critical
to perform collision detection inside each iteration of the line search (below). We emphasize that
doing so is not commonly done in computer graphics applications, and has taken substantial exper-
imentation to discover. Without these algorithmic features, the shape design approach in this paper
converges poorly; but with the approach as presented here, we were able to consistently produce
stable and convergent results.
At each Newton-Raphson iteration, we evaluate the nonlinear internal elastic forces f
int
=
f
int
(u) and matrix K = K(u), and perform collision detection and contact resolution to calculate
the collision forces f
contact
(§7.3.5), and then solve Equation 7.7 for(δu,δq). The line search finds
η
min
≥ 0 that minimizes the change of the energy E
LS
away from the current configuration,
∆E
LS
(η)= E
elastic
(u+ηδu)− E
elastic
(u)+
E
contact
(u+ηδu)− E
contact
(u)− η(δu)
T
M ˆ g− η(δq)
T
w. (7.9)
58
Here, E
elastic
is the nonlinear elastic energy of the FEM mesh, and E
contact
is the collision energy
(§7.3.5). For each η, we perform collision detection again and calculate E
contact
(u+ηδu) with
respect to the new contact configuration. The last two terms are the potential energy of the external
forces and torques, and gravity, respectively. We use the Golden ratio algorithm [59] in our line
search. Note that the E
elastic
(u) and E
contact
(u) terms do not depend on η and can be omitted.
Finally, we update u→ u+η
min
δu and q→ q+η
min
δq.
In the initial configuration, our human is typically positioned slightly above the supporting
surface and therefore out of contact. To avoid singularities in the static equilibrium, we replace K
with K+αM whenever there are no contacts;α is a tunable parameter; we useα= 10
3
,10
6
,10
4
in
the chair, shoe sole and mouse examples, respectively. In all our examples, this only occurs during
the first few iterations while the human is “falling” under gravity, and is never needed again once
the first contact appears.
7.3.5 Contact
Figure 7.4: Bidirectional contact model is needed to correctly detect and penalize sharp features
(which are very un-ergonomic). A uni-directional contact model of the top object misses the spike
because vertex v is close the bottom object’s surface.
We resolve contacts using a custom-designed penalty method, by checking every vertex of the
surface of the human tet mesh against the volume enclosed by the supporting surface’s triangle
59
mesh. We use the penalty method because it permits us to define contact forces that are differ-
entiable (with respect to the collided point position) everywhere in the three-dimensional space
except on a set of measure zero, and whose derivatives we can readily analytically compute. This
enables us to define a differentiable ergonomic cost function (§7.4.1), calculate its gradient with
respect to the shape of the supporting surface (§7.4.2) and then perform supporting surface shape
optimization to maximize ergonomic comfort. Before we settled on our method, we tried model-
ing the volume enclosed by the supporting surface’s triangle mesh using a signed distance field,
and using it to define penalty forces [111]. We abandoned this approach because it was difficult to
define the gradient of the contact force with respect to changing the vertex positions of the support-
ing surface and human meshes. We note that recent approaches demonstrated how to differentiate
cloth in constraint-based contact [112]; however, these approaches were not demonstrated for a
contact between solid objects whereby the control parameters consist of the entire shape (vertex
positions).
We define our collision energy as follows. If human vertex a with world-coordinate position
p
a
collides with the volume enclosed by the supporting surface, we find the closest point on the
supporting surface; denote its position p
anchor
. Then, a contact spring is placed between p
a
and
p
anchor
with a collision energy of
E
contact
= k
contact
||p
a
− p
anchor
||
2
, (7.10)
where k
contact
is the penalty stiffness coefficient. We adjusted k
contact
by running a few static solves
against the neutral shape of the supporting surface, until we obtained one whereby the maximum
penetration depth in the static equilibrium is approximately 0.3mm. All our optimization results
have a maximum penetration depth of approximately 0.3mm. This value is small and on the order
of accuracy of optical surface scanners.
The penalty contact forces and their gradient can then be computed by the gradient and the hes-
sian of E
contact
. To accelerate the convergence of the static solver, it is very important to incorporate
into the gradient the fact that p
anchor
slides on the supporting surface as p
a
changes location. To
60
this end, we derive and use three different gradients, depending on whether p
anchor
is a vertex, on
the interior of an edge, or on the interior of a face of a triangle of the surface mesh of the supporting
surface (Appendix A).
The non-differentiable set of measure zero consist of the boundaries of V oronoi regions of
mesh vertices, edges and faces, and of the mesh itself. Although non-differentiable at those points,
differentiability elsewhere greatly accelerated our optimization. The non-differentiability is dealt
with by employing a line search, as explained in Section 7.3.4. We initially tried using a simpler
contact formulation where we omitted the sliding terms from the gradient. Although this still
gave us convergent optimization results, the sliding anchor method and contact gradients presented
herein gave us an order of magnitude speedup in our optimizations.
7.4 Optimization
We now describe how we optimize the shape of the supporting surface to maximize the ergonomics
of contact with the human.
7.4.1 Ergonomic Cost Function
As suggested by prior work [92, 93], contact pressures and their spatial distribution is an objec-
tive way to measure comfort. How can a given pressure distribution be converted into a single,
well-defined scalar cost function that both effectively measures ergonomics and can be efficiently
computed and differentiated for shape optimization? The first idea may be to integrate the contact
pressures,
R
Γ
pdS, over the contact areaΓ. However, this does not work because, for flat (or nearly
flat) contacts, this integral approximately equals the magnitude of the gravity force of the human
61
body, therefore having no differentiation power over different pressure distribution. After exper-
imenting with different formulations, we arrived at our definition of the ergonomic cost function
ζ,
ζ =
Z
Γ
human
p
3
dS+
Z
Γ
supporting surface
p
3
dS. (7.11)
Observe that pressures p are non-negative scalars by definition, hence no norm or absolute value is
needed on p. Note that we integrate both over the contact area on the human and that on the sup-
porting surface. In a continuous setting these contact areas are the same, and the pressure distribu-
tions are also the same due to Newton’s 3
rd
law. However, when Equation 7.11 is discretized over
the contact areas on the human and supporting surface, the two integrals in the equation discretize
differently. It is important to include both summation terms (“bi-directional contact”), otherwise
one of the contact surfaces is permitted to become spiky without any penalization (Figure 7.4).
Note that we use a uni-directional contact model (§7.3.5) in the static equilibrium calculations, but
a bi-directional contact model in our ergonomic cost function (Equation 7.11). This is because the
Hessian of E
contact
(Equation 7.10) is symmetric positive-definite for uni-directional contact, but
becomes indefinite under bi-directional contact, which decreases simulation stability (§7.3.4). In
the remainder of this paper, we will use f
uni
contact
and f
bi
contact
when referring to the uni-directional
contact forces on the human tet mesh vertices ( f
contact
from §7.3.5), and bi-directional contact
forces on the union of human tet mesh vertices and vertices of the triangle mesh of the supporting
surface, respectively.
Observe that our ergonomic cost function measures discomfort: meaning, higher values are less
favorable ergonomically. Although it would be in principle “cleaner” to define a score function that
is lower and higher when ergonomics is less and more favorable, respectively (say, by composing
ζ with x7→− x or x7→ 1/x), doing so unnecessarily complicates the formulation without a clear
gain. Hence, we adopted our “discomfort” measure ζ of Equation 7.11. Finally, the exponent
3 in Equation 7.11 could be replaced with some other suitable exponent such as, say 2 or 4 (see
62
Figure 7.5), but, as argued, not 1. Higher exponents penalize localized pressures more aggressively,
but are less stable due to discretization of contact.
Assemble vertex positions of the supporting surface into a vector x∈R
3m
. Our static equi-
librium contact pressures are an implicit function of x, and therefore so is ζ. After discretizing
Equation 7.11, we obtain
ζ(x)=
∑
i
s
i
¯ p
3
i
,
¯
P=(I− L)P, (7.12)
P= P
y(x),x
y(x)=[(u
F
)
T
(x), q
T
(x)]
T
, (7.13)
where we have joined the two integrals of Equation 7.11 into one sum: i runs over each contact
vertex on the human body mesh, and on the supporting surface (Figure 7.4). Quantity y contains the
converged displacements and rigid body transformation of the human body in static equilibrium, s
i
is the surface area associated with vertex i (either on human or supporting surface), and L is the 2× 2 block-diagonal matrix whose upper-left and lower-right blocks are the umbrella Laplacian matrix
of the human body and the supporting surface’s triangle mesh, respectively. Vector P contains
contact pressures; its i-th component p
i
= f
coll,i
/s
i
is the pressure computed at vertex i. Vector
P depends on both y and x because the contact pressures depend on positions of the vertices of
both the human and the supporting surface. In order to decrease mesh discretization artifacts, we
include a Laplacian smoothing factor I− L in the definition of ζ ; namely, vector
¯
P=(I− L)P
contains the smoothened pressures; and ¯ p
i
is its component at vertex i.
To promote smoother optimized supporting surfaces, we also add a smoothing score ψ(x) to
the cost function. We initially formulated our optimization without this term, but we observed a
small amount of noise, mostly in our mouse and foot examples when running from the negative
initial guess. We define the smoothing score as the difference between the original shape x and the
shape smoothed using Taubin smoothing [109],
ψ(x)=||x−
(I− µL
tri;3
)(I− λL
tri;3
)
m
x||
2
, (7.14)
63
Figure 7.5: Sensitivity of our ergonomic cost to surface bumpiness. We sit a human into a chair
that has a bump. We then vary the spatial frequency (“spikiness”) of the bump and compute our
ergonomic cost, under various cost exponents. For easier viewing, the scores are linearly re-scaled
so that they are equal at the lowest frequency. It can be seen that all exponents are discriminative.
where L
tri;3
is the three-dimensional umbrella Laplacian matrix of the supporting surface’s triangle
mesh, m is the number of Taubin smoothing iterations, andλ andµ are Taubin smoothing param-
eters. The smoothing score only depends on the shape of the supporting surface; contacts do not
affect this score. Finally, we define our combined objective function,
φ(x)=ζ(x)+γ· ψ(x), (7.15)
whereγ is the strength of smoothing. Parameterγ was determined experimentally, by comparing
ζ(x) andψ(x) on the initial guess, and settingγ so thatζ(x) andγψ(x) have approximately equal
magnitudes. We then explored several more values γ, computed by multiplication with 10
j
, for
64
a few integers j∈{− J,...,0,...J}. We used J = 4 in our examples. Further, we found that a
small amount of additional post-process smoothing improved both the visual appearance and the
ergonomic score in some scenarios (chair example in Section 7.5.1).
We minimize φ using nonlinear conjugate gradients combined with a line search. At each
iteration, our optimization first computes the gradient of φ with respect to x, and projects it against
previously explored search directions (as standard in nonlinear conjugate gradients [59]). This
produces the conjugate gradient search direction d
i
. We then perform a 1-dimensional line search
min
κ
φ(x
i
+κd
i
), (7.16)
where x
i
is the shape of the supporting surface at iteration i. The scalarκ≥ 0 is the step length in the
line search. Note thatζ must be re-evaluated for each x
i
+κd
i
, by computing the static equilibrium
against the supporting surface with vertex positions x
i
+κd
i
(§7.3). To simplify the optimization,
we use a pre-defined set of eight κ values, and test each one independently, which can be done by
launching eight independent static equilibrium calculations in parallel. We set the eight κ values
dynamically based on d
i
, so that the maximum positional change of the supporting surface vertices
equals 1µm, 10µm, 0.1mm, 0.2mm, 0.5mm, 0.7mm, 1 mm and 2 mm. Typically, the supporting
surface shape changes little as the optimizer “probes” different values ofκ. Therefore, in order to
speed up convergence, we re-compute the static equilibrium starting from the values of u and q
discovered while computing the previous iteration x
i
.
7.4.2 Objective Function Gradient
The gradient of ∂ψ/∂x is easy to compute because ψ is a quadratic function of x. The gradient
∂ζ/∂x can be computed using the chain rule,
∂ζ
∂x
=
∂ζ
∂ f
bi
contact
∂ f
bi
contact
∂x
+
∂ f
bi
contact
∂y
∂y
∂x
. (7.17)
65
The difficult term to compute in the equation above is
∂y
∂x
. We compute it using the implicit function
theorem [85–87], but adapted here for our rigid core elastostatic formulation. When the system is in
static equilibrium, the right-hand side of Equation 7.7 is zero. Letδx be an infinitesimal change of
the supporting surface shape. Then, the contact force changes infinitesimally by (∂ f
contact
/∂x)δx,
and the updated static equilibrium satisfies
ˆ
K
δu
F
δq
=−
I
B
T
∂ f
uni
contact
∂x
δx, (7.18)
where
ˆ
K is the system matrix of Equation 7.7. Finally, we combine Equations 7.17 and 7.18,
∂ζ
∂x
=
∂ζ
∂ f
bi
contact
∂ f
bi
contact
∂x
− ∂ f
bi
contact
∂y
ˆ
K
− 1
I
B
T
∂ f
uni
contact
∂x
. (7.19)
In order to compute the gradient, we therefore only need to solve a single linear system
ˆ
Kz=
∂ζ
∂ f
bi
contact
∂ f
bi
contact
∂y
T
. (7.20)
7.4.3 Deformable Supporting Surface
Our method can also be adapted in optimizing a deformable supporting surface. To do this, we first
add all DOFs of vertices on the deformable supporting surface to the system variable u
F
. Since the
contact force depends on the deformed positions of the vertices rather than their rest positions, we
need to modify Equation 7.13 to
P= P
y(x)) y(x)=[(u
F
)
T
(x), q
T
(x)]
T
. (7.21)
66
The gradient of the ergonomic score becomes
∂ζ
∂x
=
∂ζ
∂ f
bi
contact
∂ f
bi
contact
∂y
∂y
∂x
(7.22)
Similarly to rigid supporting surfaces, when the system is in static equilibrium, the right-hand side
of Equation 7.7 is zero, and we get
ˆ
K
δu
F
δq
=−
I
B
T
∂ f
int
∂x
δx. (7.23)
We therefore need to solve a linear system of equations equal in form to that of Equation 7.20.
7.5 Results
We give three examples demonstrating ergonomic contact between a human and a rigid supporting
surface (Figure 7.3). Our optimization is able to reduce the ergonomic score by large percentages
(Table 7.1), resulting in pressure distributions whereby the “hot-spot” pressures are reduced 2-3x,
and their carrying load evenly distributed across the entire contact surface. All examples were
computed on a 3.00 GHz Intel Xeon i7 CPU E5-2687W v4 processor. We used 8 threads for
optimization as we attempt 8 step lengths (κ) in one line search (Equation 7.16).
7.5.1 Human sitting in a chair
We designed an ergonomic chair for a specific person. We scanned and simulated the subject’s
lower back, buttocks and upper thighs region. The other parts of the body are modeled as rigid be-
cause they do not deform in our scenario. We acquired the depths and Young’s moduli at different
body locations, and created the rigid core by offsetting the external shape using the spatially vary-
ing depths [108]. The rigid core can rotate arbitrarily. Our chair is generated from a smoothed flat
block as initial guess. The human sits in the central region of the chair, and it is not our intention
67
Supporting surface Chair Shoe sole Mouse
Initial guess shape Flat plane Flat plane Negative image
Initial guess score 4.99× 10
11
2.99× 10
13
5.48× 10
7
Optimized score 1.39× 10
11
5.36× 10
12
1.37× 10
7
Percentage reduced 72.07% 82.10% 74.96%
Smoothing strength 1× 10
13
1× 10
19
1× 10
13
Optimization time 4 hours 19 hours 13 hours
#Human vtxs 10,231 14,060 77,749
#Supp. surface vtxs 26,891 5,894 13,614
#Contact vtxs 1,285 575 2,171
Table 7.1: Ergonomic scores and optimization timings. Human (deformable), supporting surface
(rigid) and contact vertex counts refer to the tet mesh, triangle mesh and the static equilibrium,
respectively. The human exterior surface is the tet mesh’s exterior surface.
to permit horizontal sliding. Therefore, we lock the two directions of horizontal translation, and
only permit the rigid core to translate along the gravity direction. Figure 7.6, top, shows that our
optimization can greatly reduce the ergonomic score (by 72.1%). The output pressure distribution
is also clearly better on the optimized chair. Note that we perform similar locking of translations
in all of our examples. Doing so prevents sliding and therefore we did not need to employ friction
forces in our contact model or the shape optimization, which makes the optimization simpler, faster
and more tractable. In this particular example, we found that adding a small amount of additional
Taubin smoothing (µ=− 0.53,λ = 0.5) to the final optimized chair smoothened out minor numer-
ical irregularities during optimization and further reduced the ergonomic score (from 3.62× 10
11
to 1.39× 10
11
). All of our chair examples have this smoothing applied before computing pressures
and ergonomic score. This postprocessing was not significant in the other examples (mouse, shoe)
and was omitted.
The above result was obtained starting from a flat chair as an initial guess. Another potential
initial guess is the “negative image” chair obtained by subtracting the human body from a flat chair
using CSG boolean subtraction. We create the negative image by penetrating the human into the
chair by 5cm. To remove the sharp edges caused by CSG, we use Laplacian smoothing. Figure 7.6,
bottom, shows that the negative image chair, when used directly as the chair shape, has a bad
68
ergonomic score because of the high pressures along the edges. Although we smoothed the edges,
the human body will still generate high pressures along the edges because of the deformation.
When starting from the “negative image” initial guess, our optimization method can reduce the
pressure substantially (92.3%), however the end result is inferior to optimization that started from
a flat chair (Figures 7.6). We also tried other penetrations: 1cm, 2cm, 3cm, 4cm, 6cm; they all
produced inferior results. The 6cm was unstable because the negative image was too deep, creating
a concave chair.
7.5.2 Human standing on a shoe sole
We purchased the triangle meshes of a human foot and a flat shoe sole from [Turbosquid]. We
tetrahedralized the human foot using our method described in Section 7.3.1. To simulate a real
human standing on the shoe sole, we model the presence of the upper human body using an external
gravity force equaling to half of the body mass. Note that the human body’s center of mass is
typically not on the vertical line through the rigid core’s center of mass; we model this effect by
adding an external torque applied to the rigid core (Figure 7.7). By tweaking this torque, we can
simulate the effect of the human body leaning forwards or backwards. Although our system does
not incorporate motor control, setting a specific value of this torque (either positive or negative)
will shift the contact pressure distribution on the shoe sole to compensate for the leaning body.
This effectively simulates a stably standing human, and is possible as long as the external torque is
consistent with a body center of mass being within the shoe sole surface area. Similar to the chair
example, we locked the two directions of horizontal foot translation. In addition, because we are
not interested in yaw and roll motion of the foot, we only allow the rigid core of the foot to rotate
around the pitch axis.
Similar to the chair model, we also use a flat shoe sole and a negative image as our initial
guesses (Figure 7.8). Our results demonstrate that the optimizer works well on the flat shoe sole
(82.1% ergonomic score reduction). When starting from the negative image, the score reduction is
smaller (29.0%). Overall, the best score was obtained starting from a flat shoe sole.
69
Initial guess: at chair Initial guess: negative image
Pressures (initial guess) Pressures (converged result)
Displacements
Displacements
Score vs optimization iteration
Score vs optimization iteration
Pressures (initial guess) Pressures (converged result)
Figure 7.6: Optimized pressure maps and ergonomic scores (chair). “Displacements” gives
the magnitude of vertex displacements between the initial guess and the converged results. When
using the flat chair initial guess (top), the ergonomic score was reduced from 4. 99× 10
11
to 1.39× 10
11
(72.1%); the maximum penetration depth of the optimized chair was 0.15mm. When using
the negative image initial guess (bottom), the ergonomic score was reduced from 2.05× 10
12
to
1.58× 10
11
(92.3%); the maximum penetration depth of the optimized chair was 0.16mm.
70
External gravity
External
torque
Figure 7.7: The rigid core model of the foot. Left: The outer layer and the inner layer of the foot.
The inner layer forms a rigid core of the foot. Right: The deformable tetrahedra form the outer
layer.
7.5.3 Human grasping a mouse
In this example, we explored how to create an ergonomic shape of the mouse as it contacts the
human hand. The mouse neutral shape was modeled using CSG operations on spheres, followed
by Taubin smoothing. We rigged the human hand using linear blend skinning, and posed it so that
it sits on the mouse in a natural configuration. Although we do not articulate the human hand in
our system, or model the biomechanics of the upper arm, our model does permit calculating the
contact pressures experienced by the hand due to gravity, as it rests on top of the mouse. We locked
all hand rotation axes and only permit the hand to translate up and down. Our optimizer can greatly
reduce the ergonomic score for both the neutral mouse (93.5%) and negative image (75.0%) initial
guesses. The negative image method produced a better ergonomic score than the neutral mouse
(Figure 7.9).
7.5.4 Human sitting in a deformable chair
We also demonstrate that our method can be adapted to designing the rest shape of a deformable
supporting surface. We used the same scanned human model as in the rigid chair example, and we
71
Initial guess: at shoe sole Initial guess: negative image
Pressures
(initial guess)
Pressures
(converged result)
Displacements
Pressures
(initial guess)
Pressures
(converged result)
Displacements
Score vs optimization iteration Pressure histogram
Score vs optimization iteration Pressure histogram
Figure 7.8: Optimized pressure maps and ergonomic scores (shoe sole). “Displacements” gives
the magnitude of vertex displacements between the initial guess and the converged results. When
using the flat shoe sole initial guess (top), the ergonomic score was reduced from 2. 99× 10
13
to
5.36× 10
12
(82.1%); maximum penetration depth of the optimized shoe sole was 0.30mm. When
using the negative image initial guess (bottom), the ergonomic score was reduced from 8.56× 10
12
to 6.08× 10
12
(29.0%); maximum penetration depth of the optimized shoe sole was 0.52mm. In
both results, it can be seen that the pressure histogram shifts to the left (i.e., to smaller pressures).
The bins at smaller pressures also became taller as the optimized shoe sole has a larger number
of contacts, each carrying a smaller contact pressure. The grey color in the pressure histogram
denotes the overlap of the blue color (negative image) and the pink color (optimized shoe).
72
Pressures
(initial guess)
Pressures
(converged result)
Displacements
Initial guess: neutral mouse Initial guess: negative image
Score vs optimization iteration Pressure histogram
Pressures
(initial guess)
Pressures
(converged result)
Displacements
Score vs optimization iteration Pressure histogram
Figure 7.9: Optimized pressure map and ergonomic scores (mouse). “Displacements” gives the
magnitude of vertex displacements between the initial guess and the converged results. When using
the neutral mouse shape (top) as the initial guess, the ergonomic score was reduced from 5.25× 10
8
to 3.34× 10
7
(93.5%); maximum penetration depth of the optimized mouse was 0.25mm. When
using the negative image shape (bottom) as the initial guess, the ergonomic score was reduced
from 5.48× 10
7
to 1.37× 10
7
(75.0%); maximum penetration depth of the optimized mouse was
0.13mm. Similarly to the shoe sole example, the pressure histograms shifted to the left (i.e., to
smaller pressures) and became taller, due to larger number of contacts each carrying a smaller
contact pressure; very high pressures disappeared (note: logarithmic pressure scale).
73
also used a smoothed flat box as the initial guess. Like in the rigid case, we locked the two direc-
tions of horizontal translations. The top part of Figure 7.10 shows that our optimization decreases
the ergonomic score of a deformable chair (50kPa Young’s modulus) by 63.1%. The output pres-
sure distribution is also clearly improved. The vertex displacements between the initial guess and
the converged result are smaller than in the rigid chair example because the deformable material of
the chair already accommodates some ergonomic “molding” to the human body. Because the rest
shape of the chair becomes less important for ergonomic design (as it exhibits smaller variability),
the number of iterations in our optimization becomes smaller.
We also modified the Young’s modulus of the deformable chair to 20kPa (2.5x softer) and
125kPa (2.5x stiffer), to investigate how the chair material properties affect the results (middle and
bottom parts of Figure 7.10). We observed that (as expected) the softer chair results in a lower (i.e.,
better) ergonomic score. Also, the rest shape of the soft chair changed less from its initial shape,
as the chair’s greater deformability partially enables a good ergonomic outcome on its own. More
rigid chairs, however, need a greater change in their rest shape to decrease the ergonomic score.
7.5.5 Smoothing Strength
In order to evaluate the effect of smoothing (γ) values to our optimization, we optimized the mouse
with different γ values, 0, 10
11
, 10
13
, 10
15
and 10
17
. Each γ gives a different optimized shape
(Figure 7.11) and ergonomic score (Table 7.2). The smoothing strength of 10
13
gives the best
optimized ergonomic score.
7.6 Physical validation of our method
We validated our simulation results using a physical experiment where we compared pressure dis-
tributions measured in the real-world to those predicted by our simulation method, both for the
initial guess supporting surface and the optimized supporting surface. In doing so, we designed a
new inexpensive experimental method to measure real distributed contact pressures between two
74
2.5x softer chair
Pressures (initial guess)
Pressures (converged result)
Displacements
Score vs optimization iteration
2.5x stiner chair
Pressures (initial guess)
Pressures (converged result)
Displacements
Score vs optimization iteration
Initial guess: at chair
Pressures (initial guess)
Pressures (converged result)
Displacements
Score vs optimization iteration
Young’s Modulus: 50kPa
Figure 7.10: Optimized pressure maps and ergonomic scores (deformable chair). “Pressures”
gives the pressure field on the deformed chair when the simulation reaches the static equilibrium.
“Displacements” gives the magnitude of vertex displacements between the rest shape of the initial
guess and the rest shape of the converged result. All experiments use the flat chair initial guess.
Left: under a Young’s modulus of 50kPa, the ergonomic score was reduced from 5.59× 10
11
to
2.06× 10
11
(63.1%); maximum optimized penetration depth was 0.20mm. Middle: softer chair
(25kPa). The ergonomic score was reduced from 2.46× 10
11
to 1.63× 10
11
(33.6%). Right: stiffer
chair (125kPa). The ergonomic score was reduced from 6.83× 10
11
to 2.62× 10
11
(61.7%).
75
Figure 7.11: Optimized shapes and ergonomic scores over the iterations of the optimization,
using different smoothing values. Top row: optimized shape with color-coded displacements.
Bottom row: ergonomic scores at each optimization iteration. The smoothing values used from the
leftmost column and the rightmost column are 0, 10
11
, 10
13
, 10
15
and 10
17
, respectively.
curved surfaces. Doing so is challenging because one has to somehow insert measuring equipment
between the two curved surfaces, but the presence of the sensors itself potentially changes the
distributed contacts, and thus interferes with the contact measurements. Initially, we attempted to
model the pressures using thin pressure measurement film (paper) from FujiFilm [113]. However,
this did not work well because the optimized shoe is a general surface with non-zero Gaussian cur-
vature (i.e., is not a developable surface). Consequently, when we placed the film in between the
foot and the optimized shoe, the film warped, which interfered with the pressure measurements.
Instead, we then used electrical resistance-based pressure sensors from Tekscan [114]. These sen-
sors have a shallow profile and minimally interfere with pressures. By designing proper electrical
circuits as described below, we were able to reliably measure the contact pressure distribution on
the entire shoe for less than $500.
Smoothing strength Ergonomic score Percentage reduced
0 5.04× 10
7
90.4%
10
11
4.70× 10
7
91.1%
10
13
3.40× 10
7
93.5%
10
15
1.41× 10
8
73.1%
10
17
7.94× 10
8
-51.1%
Table 7.2: Ergonomic scores using different smoothing values. The initial scores for the five op-
timizations are all 5.25× 10
8
. The entry in the last row is negative because the smoothing strength
is too large, causing the optimizer to overly smoothen the shape and compromise ergonomics.
76
First, we used the Alja-Safe material [115] to make a plaster model (Figure 7.12) of a real
human foot, and scanned it with a 3D scanner (Artec Spider [116]). We then ran ergonomic
shoe sole optimization against the real foot shape, as described in our paper. The ergonomic
score was reduced by 82.8%, and the pressure map was visibly improved (Figure 7.14). After
the optimization, we 3D-printed the optimized rigid shoe sole. We then asked the same subject to
stand on a flat rigid surface, and on the optimized shoe sole. We used the Tekscan [114] A101 force
sensor to measure real pressures at 78 sample locations on the shoe sole (Figure 7.12). This sensor
can measure force in the range from 0N to 44N, and its diameter is 3.81 mm, which is suitable
for our experiments. The Tekscan force sensor changes its resistance when a force is applied
based on a known linear profile provided in manufacturer’s datasheet. We used the manufacturer-
provided dedicated circuit to convert the force into a voltage signal and acquired the voltages into
the computer using a digital acquisition card. We used a thin electrical strap provided by Tekscan
to electrically connect the two pins on the force sensor to the Tekscan measuring electric circuit so
as to cause minimal interference with the mechanical contact. The connector provided by Tekscan
was too thick (∼ 2mm profile), and, as it is in immediate proximity of the sensor, found itself
“jammed” in between the foot and the shoe, interfering with the contact pressures. We removed
the connector and manually resoldered the wires, and then gently sanded the soldered electrical
connections down to a∼ 0.3mm profile using sand-paper. We measured each of the 78 locations
one by one. This kept our design inexpensive, and also reduced contact interference. The Tekscan
measuring circuit samples the contact force at the sensor at 10Hz. For each measuring location,
the subject stood in the shoe for 1 minute, which gave us 600 pressure points for each sample
location. The subject was standing on both legs and was told to keep a stable neutral standing
posture. Out of the 600 pressures, we discarded the smallest and largest 50 values to reduce the
noise. We then computed the average pressure and standard deviation at each sample location. We
then interpolated the pressures onto the vertices of the shoe sole using linearly precise biharmonic
weights [117], whereby each of the 78 sensor locations serves as a “handle”.
77
Figure 7.12: Manufactured models and the force sensor. Left: the plastic foot (photo not to
scale) manufactured using the Alja-Safe cloning process. This plastic foot is used only to stably
scan the foot 3D shape; it is otherwise unused in our method. Middle: the 3D-printed plastic
optimized shoe sole with sample locations marked. The reference force sensor is seen fixed at the
bottom of the shoe sole. Right: the relationship between the contact force and output voltage.
Voltages at reference force sensor
flat shoe optimized shoe
Voltage (V)
Figure 7.13: Output voltages of the reference sensor across the 78 measurements for the un-
optimized (flat) and optimized shoe. Variations are relatively small.
78
The contact pressures can be somewhat affected by variations in the human standing pose.
To address this, we placed a second reference force sensor with an independent circuit at a fixed
location on the shoe sole (Figure 7.12). We recorded a one-minute pressure signal for the refer-
ence sensor and computed its average pressure at the fixed location. For each of the 78 measuring
points, this gives us a primary pressure p
i
at the measuring location, and a pressure r
i
at the ref-
erence location. We then scale pressures p
i
with a scalar α
i
, where α
i
is chosen so that α
i
r
i
is
the same constant C for all i; this compensates for variations in the leaning postures across the 78
measurements. We choose constant C so that, after the pressuresα
i
p
i
are interpolated to the entire
shoe sole (using linearly precise biharmonic weights [117]), the total net contact force on the foot
equals half of the gravity force of the human subject. This convenient choice of C also avoids the
need for a separate calibration step to calibrate the output voltages versus sensor resistance. Gener-
ally, we want the reference sensor pressures to be close to each other when the primary measuring
sensor is placed at different locations, i.e., posture variations should be small; and this is indeed
the case (Figure 7.13). We plot the resulting pressures in Figure 7.14. It can be seen that (1) the
real and simulated pressures on the un-optimized shoe sole are qualitatively similar to each other,
and (2) the real and simulated pressures on the optimized shoe sole are qualitatively similar to each
other, and (3) the un-optimized shoe sole has substantial pressure hot-spots whereas the optimized
pressures does not have hot-spots and the pressures are spread over a large surface area.
7.7 Comparison to “Sit & Relax”
We compared our work with the “Sit & Relax” method proposed by Leimer et al. [68]. This method
conceptually immerses the supported body part into its negative image volume, in order to estimate
a pressure distribution across the surface. During optimization, the supporting surface is raised to
make contact with the surface of the body part, with each point’s affinity weighted by the corre-
sponding negative-image pressure. The smoothness term and an underlying Catmull-Clark spline
model act as counter-balancing regularizers. Since the pressure map is computed for a support
79
Figure 7.14: Comparison between simulation and real measurements. Top row: un-optimized
(flat) shoe sole. Bottom row: optimized shoe sole. Left column: pressure computed using sim-
ulation. Middle column: pressure obtained using real measurements. Right column: standard
deviation of the pressure obtained using real measurements during the 1-minute acquisition.
surface that is never actually realized (the negative image), and is frozen at the beginning of the
optimization instead of varying with the evolving support surface, the Sit & Relax approach leads
to un-ergonomic high-pressure regions, visualized in red in Figure 7.15. In contrast, our method
recomputes the pressure map in each iteration, using a physically realistic, differentiable model
of deformable-rigid contact, and achieves a much more evenly distributed pressure map without
sharp peaks. The ergonomic score (Equation 7.11) for the Sit & Relax chair is 8.40× 10
11
. In
contrast, the score for our optimized chair (starting from a flat surface) is 1. 39× 10
11
, or 83.4%
less. We also experimented with initializing our method with the optimized Sit & Relax surface.
This produced a score of 2.36× 10
11
, still not quite as good as with the flat initialization. We
conjecture that starting with an unbiased flat surface allows the gradient descent to reach a better
local minimum in this case. Note that these visualizations and scores assume that our simulation
model closely approximates real-world pressure distributions and hence can be used as a basis for
80
comparison of the two methods. We validated this in Section 7.6. Since the Sit & Relax optimiza-
tion involves no physical simulation in the inner loop, it is certainly orders of magnitude faster
than ours. Speeding up our optimization is an interesting avenue for future work; for example, by
parameterizing our supporting surface with Catmull-Clark splines.
Sit & Relax our method
displacements pressures
Figure 7.15: Comparison to Sit & Relax [68]. Left and right columns give the results optimized
with Sit & Relax and our method, respectively. “Displacements” gives the magnitude of vertex
displacements from the initial guess, and “pressures” are the contact pressures computed using
physically based simulation against the output of either method. Because our method incorporates
a physically based simulation, it gives a distribution with far fewer high-pressure points.
7.8 Comparison to “Pose to Seat”
A more recent method by Leimer et al. – “Pose to Seat” [69] – also optimizes supporting surfaces
via inferred pressure maps. The main contribution of this paper, compared to “Sit & Relax”, is
handling an articulated human model. For their contact model, they assume all downward-facing
vertices are in contact, effectively modeling the pressure distribution on a negative image mold
81
and keeping it unchanged throughout the geometric optimization of the support surface, just like
Sit & Relax. In determining the pressure distribution, they do not consider the other contact object
(the supporting surface), deformable elasticity, or contact simulation. In contrast, our method does
not handle articulation, but does model a (differentiable) physically-based contact pressure map
specific to a supporting surface, and the distributed pressure map is updated repeatedly in the inner
loop of our shape optimization.
Since these aspects (articulation vs distributed contact) are orthogonal aspects of the problem,
it is difficult to do an apples-to-apples comparison with this paper. Further, the support synthesis
step of Leimer et al. [69] (which uses the fixed pressure distribution as a precomputed importance
map but is otherwise completely non-physical) is designed to support the entire human body, using
a handcrafted template with pre-assigned patches per body part. This algorithm is highly specific
and it is not clear how it could be adapted to our scenarios.
Hence, we chose to compare to Leimer et al. [69] by studying the pressure maps produced by
following their assumptions, without considering articulations. To aid their method, we explicitly
computed more realistic sets of contact vertices for their method, by bringing the (non-elastic) foot
in contact with candidate supporting surfaces and marking the foot vertices within a small thresh-
old distance of the supports. In this fashion, we ran their method with two different sets of body
vertices in contact (Figure 7.16). We transferred their contact pressures to the shoe sole for visual-
ization and comparison to our method (Figure 7.16). Our experiment demonstrates that our pres-
sure maps are much closer to the real measured pressures than those obtained using [69], despite
our contact modeling improvements to their method (Figure 7.16). Considering all downward-
facing vertices to be in contact, as the original Leimer et al. “Pose to Seat” approach does, yields a
nearly uniform, even less realistic pressure distribution that is independent of the actual supporting
surface.
82
[Leimer 2020] our method
at sole optimzed sole
measured
Figure 7.16: Comparison to Pose to Seat [69]. The left, middle and right columns give the
pressure maps on the shoe sole computed by Leimer et al. [69], real-world experiments and our
optimization method, respectively. While Leimer et al. assume all downwards-facing vertices are
in contact, we aided their method by supplying more realistic contact regions (without considering
elasticity). We investigated two contact sets of vertices: “flat sole” and “optimized sole” give the
contact pressures between the foot and the flat shoe sole, and between the foot and our optimized
shoe sole. Because our method incorporates physically based simulation and our human foot has
a deformable layer, our method produces contact pressures that are qualitatively and quantitatively
closer to the measured pressures.
83
Chapter 8
Conclusion
This thesis proposed a method to measure the human elastic property including Young’s modulus
and material non-linearity. We designed a hand-held device that uses a force gauge and a laser
to capture the contact force and indentation depth simultaneously. Because the device is portable,
it can be easily used on the human body. To model the non-linearity of the material flexibly, we
used natural cubic spline to model the three scalar functions in Valanis-Landel material model. To
accelerate the optimization process, we proposed a novel idea that simplifies the 3-dimensional
simulation into a 2-dimensional system as the contact force and deformation can be symmetric
with the respect ofϕ component in the cylindrical coordinate system.
To measure the depth of human tissue in vivo, we proposed two methods. One is to use an ultra-
sound device and read the tissue depth from the scanned images. Another method is to use the same
optimzation method as estimating the human elastic porperty by having a series of reasonable depth
guesses and optimizing based on these guesses. Among the final converged optimized results, we
pick up the depth guess that has the best matching between the captured force-indentation curve
and the simulated one.
We also presented an application of our work in designing the shape of supporting surfaces to
maximize the ergonomics of their contact against the human body. While previous ergonomic work
mostly focused on geometric design, we contribute a physically grounded approach that models
elastic contact using large-deformation FEM nonlinear elasticity and uses measured real human
84
tissue properties. We validated our approach using 3D printing and confirmed that the optimized
pressure distributions qualitatively match those observed in the real world.
We modeled our human as a passive deformable object in a representative pose. While this
proved to be a physically plausible and computationally tractable model, method accuracy could be
improved by modeling articulated human bodies and adding active human motor control. To make
the optimization tractable and differentiable, we model contact using a penalty-based scheme. Al-
though our penalty forces are not differentiable at the point where a vertex leaves contact, the
non-differentiable set is small (has measure zero) and we were in practice able to mitigate the lack
of differentiability by taking sufficiently small steps in our line search. Constraint-based contact
models are a good avenue for future work. Our optimizer only requires the ergonomic score and
its gradient; we did not pursue Hessians. In our work, we are interested in contact whereby the
supporting surface touches the human in a chosen region of interest of the human body. Tangential
sliding is not desirable and we prevented it by locking the two tangential translational degrees of
freedom of the human. This also had an important pragmatic benefit of not needing a frictional
model. Friction forces are notoriously difficult to differentiate. Our optimization spends most of
the time in computing the static equilibrium against the supporting surface. Exploring contact
models, including friction, that permit rapid optimization and are readily differentiable, yet they
accelerate the optimization, is an important avenue for future work. While we launch our opti-
mization from multiple initial guesses, much like most optimization methods, our optimization
scheme is local and does not necessarily find the globally optimal solution. Our optimizer treats all
surface vertices of the supporting surface as the optimized degrees of freedom; lower-dimensional
deformation models (such as, for example, linearly precise biharmonic weights [117]) could be
used instead. In the future, more accurate internal human models could be pursued, such as those
involving anatomical shapes of muscles and subcutaneous fat. Because the human body is not
equally sensitive to pressure in all parts of the body, our ergonomic function could be equipped
with a “pain map” [94]. As a common choice in ergonomics literature, our ergonomic function in-
corporates contact pressures; but ergonomics is also related to posture comfort, such as muscle and
85
joint torques. Furthermore, while most surveyed ergonomic paper advocate for contact pressures
as the ergonomic metric, an alternative discomfort metric to explore is three-dimensional elastic
strain in the human body.
86
References
1. Irving, G., Teran, J. & Fedkiw, R. Invertible Finite Elements for Robust Simulation of Large
Deformation in Symp. on Computer Animation (SCA) (2004), 131–140.
2. Yang, F. Thickness effect on the indentation of an elastic layer. Materials Science and En-
gineering: A 358, 226–232 (2003).
3. Wang, D., Roesler, J. R. & Guo, D.-Z. Innovative algorithm to solve axisymmetric dis-
placement and stress fields in multilayered pavement systems. Journal of Transportation
Engineering 137, 287–295 (2010).
4. La Ragione, L., Musceo, F. & Sollazzo, A. Axisymmetric indentation of a rigid cylinder on
a layered compressible and incompressible halfspace. Journal of Mechanics of Materials
and Structures 3, 1499–1520 (2008).
5. Sneddon, I. N. Fourier transforms (Courier Corporation, 1995).
6. Yang, F. Thickness effect on the indentation of an elastic layer. Materials Science and En-
gineering: A 358, 226–232 (2003).
7. Maina, J. & Matsui, K. Elastic multi-layered analysis using DE-integration. Publications of
the Research institute for Mathematical Sciences 41, 853–867 (2005).
8. Boussinesq, J. Application des potentiels ` a l’´ etude de l’´ equilibre et du mouvement des
solides ´ elastiques (Gauthier-Villars, Imprimeur-Libraire, 1885).
9. Sneddon, I. N. The relation between load and penetration in the axisymmetric Boussinesq
problem for a punch of arbitrary profile. International journal of engineering science 3, 47–
57 (1965).
10. Johnson, K. L. & Johnson, K. L. Contact mechanics (Cambridge university press, 1987).
11. Love, A. E. H. A treatise on the mathematical theory of elasticity (Cambridge university
press, 2013).
12. Wang, D., Roesler, J. R. & Guo, D.-Z. Innovative algorithm to solve axisymmetric dis-
placement and stress fields in multilayered pavement systems. Journal of Transportation
Engineering 137, 287–295 (2010).
13. La Ragione, L., Musceo, F. & Sollazzo, A. Axisymmetric indentation of a rigid cylinder on
a layered compressible and incompressible halfspace. Journal of Mechanics of Materials
and Structures 3, 1499–1520 (2008).
14. Wex, C., Arndt, S., Stoll, A., Bruns, C. & Kupriyanova, Y . Isotropic incompressible hy-
perelastic models for modelling the mechanical behaviour of biological tissues: a review.
Biomedical Engineering/Biomedizinische Technik 60, 577–592 (2015).
15. Samani, A. & Plewes, D. A method to measure the hyperelastic parameters of ex vivo breast
tissue samples. Physics in Medicine & Biology 49, 4395 (2004).
16. O’Hagan, J. J. & Samani, A. Measurement of the hyperelastic properties of 44 pathological
ex vivo breast tissue samples. Physics in Medicine & Biology 54, 2557 (2009).
87
17. Carter, F. J., Frank, T. G., Davies, P. J., McLean, D. & Cuschieri, A. Measurements and
modelling of the compliance of human and porcine organs. Medical Image Analysis 5, 231–
236 (2001).
18. Davies, P. J., Carter, F. J. & Cuschieri, A. Mathematical modelling for keyhole surgery
simulations: a biomechanical model for spleen tissue. IMA Journal of Applied Mathematics
67, 41–67 (2002).
19. Chatelin, S., Oudry, J., P´ erichon, N., Sandrin, L., Allemann, P., Soler, L., et al. In vivo liver
tissue mechanical properties by transient elastography: Comparison with dynamic mechan-
ical analysis. Biorheology 48, 75–88 (2011).
20. Chui, C., Kobayashi, E., Chen, X., Hisada, T. & Sakuma, I. Combined compression and
elongation experiments and non-linear modelling of liver tissue for surgical simulation.
Medical and Biological Engineering and Computing 42, 787–798 (2004).
21. Nava, A., Mazza, E., Furrer, M., Villiger, P. & Reinhart, W. In vivo mechanical characteri-
zation of human liver. Medical image analysis 12, 203–216 (2008).
22. Toshima, M., Ohtani, Y . & Ohtani, O. Three-dimensional architecture of elastin and colla-
gen fiber networks in the human and rat lung. Archives of histology and cytology 67, 31–40
(2004).
23. McKee, C. T., Last, J. A., Russell, P. & Murphy, C. J. Indentation versus tensile measure-
ments of Young’s modulus for soft biological tissues. Tissue Engineering Part B: Reviews
17, 155–164 (2011).
24. Pailler-Mattei, C., Bec, S. & Zahouani, H. In vivo measurements of the elastic mechanical
properties of human skin by indentation tests. Medical engineering & physics 30, 599–606
(2008).
25. Tran, Q. T., Gerratt, B. R., Berke, G. S. & Kreiman, J. Measurement of Young’s modulus in
the in vivo human vocal folds. Annals of Otology, Rhinology & Laryngology 102, 584–591
(1993).
26. Bickel, B., Baecher, M., Otaduy, M., Matusik, W., Pfister, H. & Gross, M. Capture and
Modeling of Non-Linear Heterogeneous Soft Tissue. ACM Trans. on Graphics (SIGGRAPH
2009) 28, 89:1–89:9 (2009).
27. Pai, D. K., Rothwell, A., Wyder-Hodge, P., Wick, A., Fan, Y ., Larionov, E., et al. The hu-
man touch: measuring contact with real human soft tissues. ACM Transactions on Graphics
(TOG) 37, 58 (2018).
28. M¨ uller, M. & Gross, M. Interactive Virtual Materials in Proc. of Graphics Interface 2004
(2004), 239–246.
29. Civit-Flores, O. & Sus´ ın, A. Robust Treatment of Degenerate Elements in Interactive Coro-
tational FEM Simulations. Computer Graphics Forum 33, 298–309 (2014).
30. Chao, I., Pinkall, U., Sanan, P. & Schr¨ oder, P. A Simple Geometric Model for Elastic De-
formations. ACM Trans. on Graphics (SIGGRAPH 2010) 29, 38:1–38:6 (2010).
31. Bonet, J. & Wood, R. D. Nonlinear Continuum Mechanics for Finite Element Analysis, 2nd
Ed. (Cambridge University Press, 2008).
32. Barbiˇ c, J. Real-time Reduced Large-Deformation Models and Distributed Contact for Com-
puter Graphics and Haptics PhD thesis (Carnegie Mellon University, 2007).
33. Mooney, M. A theory of large elastic deformation. Journal of applied physics 11, 582–592
(1940).
88
34. Rivlin, R. S. & Saunders, D. Large elastic deformations of isotropic materials VII. Exper-
iments on the deformation of rubber. Philosophical Transactions of the Royal Society of
London. Series A, Mathematical and Physical Sciences 243, 251–288 (1951).
35. Ogden, R. W. Large deformation isotropic elasticity–on the correlation of theory and exper-
iment for incompressible rubberlike solids. Proceedings of the Royal Society of London. A.
Mathematical and Physical Sciences 326, 565–584 (1972).
36. Yeoh, O. H. Characterization of elastic properties of carbon-black-filled rubber vulcanizates.
Rubber chemistry and technology 63, 792–805 (1990).
37. Yeoh, O. H. Some forms of the strain energy function for rubber. Rubber Chemistry and
technology 66, 754–771 (1993).
38. Xu, H., Sin, F., Zhu, Y . & Barbiˇ c, J. Nonlinear material design using principal stretches.
ACM Transactions on Graphics (TOG) 34, 75 (2015).
39. Bower, A. Applied Mechanics of Solids (Taylor & Francis, 2011).
40. Bonet, J. & Burton, A. A simple orthotropic, transversely isotropic hyperelastic constitu-
tive equation for large strain computations. Computer Methods in Applied Mechanics and
Engineering 162, 151–164 (1998).
41. Picinbono, G., Delingette, H. & Ayache, N. Non-linear and anisotropic elastic soft tissue
models for medical simulation in IEEE Int. Conf. on Robotics and Automation 2001 (2001).
42. Ten Thije, R., Akkerman, R. & Hu´ etink, J. Large deformation simulation of anisotropic
material using an updated Lagrangian finite element method. Computer methods in applied
mechanics and engineering 196, 3141–3150 (2007).
43. Allard, J., Marchal, M., Cotin, S., et al. Fiber-based fracture model for simulating soft tissue
tearing. Studies in health technology and informatics 142, 13–18 (2009).
44. Li, Y . & Barbiˇ c, J. Stable Orthotropic Materials in Symp. on Computer Animation (SCA)
(2014), 41–46.
45. Li, Y . & Barbiˇ c, J. Stable anisotropic materials. IEEE transactions on visualization and
computer graphics 21, 1129–1137 (2015).
46. Instruments, P. User Manual PCE-DFG 500 Series Force Gauge https://www.pce-instruments.com/
().
47. Optex. Smallest displacement sensor in class https://www.optex-ramco.com ().
48. Boissonnat, J.-D. & Oudot, S. Provably good sampling and meshing of surfaces. Graphical
Models 67, 405–451 (2005).
49. Loper, M., Mahmood, N., Romero, J., Pons-Moll, G. & Black, M. J. SMPL: A skinned
multi-person linear model. ACM transactions on graphics (TOG) 34, 248 (2015).
50. Amberg, B., Romdhani, S. & Vetter, T. Optimal step nonrigid icp algorithms for surface
registration in 2007 IEEE Conference on Computer Vision and Pattern Recognition (2007),
1–8.
51. Xu, H. & Barbiˇ c, J. Signed Distance Fields for Polygon Soup Meshes. Graphics Interface
2014 (2014).
52. Lewiner, T., Lopes, H., Vieira, A. W. & Tavares, G. Efficient implementation of marching
cubes’ cases with topological guarantees. Journal of graphics tools 8, 1–15 (2003).
53. Hang Si. TetGen: A Quality Tetrahedral Mesh Generator and a 3D Delaunay Triangulator
2011.
89
54. Barill, G., Dickson, N. G., Schmidt, R., Levin, D. I. W. & Jacobson, A. Fast Winding Num-
bers for Soups and Clouds. ACM Trans. Graph. 37, 43:1–43:12. doi:10.1145/3197517.
3201337 (2018).
55. Wriggers, P. & Laursen, T. A. Computational contact mechanics (Springer, 2006).
56. Devillers, O. & Guigue, P. Faster triangle-triangle intersection tests PhD thesis (INRIA,
2002).
57. M¨ oller, T. A fast triangle-triangle intersection test. Journal of graphics tools 2, 25–30
(1997).
58. Van den Bergen, G. Efficient Collision Detection of Complex Deformable Models using
AABB Trees. J. of Graphics Tools 2, 1–14 (1997).
59. Press, W., Teukolsky, S., Vetterling, W. & Flannery, B. Numerical recipes: The art of scien-
tific computing third (Cambridge University Press, Cambridge, UK, 2007).
60. Armijo, L. Minimization of functions having Lipschitz continuous first partial derivatives.
Pacific Journal of mathematics 16, 1–3 (1966).
61. Putz-Anderson, V ., Bernard, B. P., Burt, S. E., Cole, L. L., Fairfield-Estill, C., Fine, L. J.,
et al. Musculoskeletal disorders and workplace factors. National Institute for Occupational
Safety and Health (NIOSH) 104 (1997).
62. Baudisch, P. & Mueller, S. Personal Fabrication. Foundations and Trends in Human-Computer
Interaction 10, 165–293 (2017).
63. Chaudhuri, S., Ritchie, D., Wu, J., Xu, K. & Zhang, H. Learning Generative Models of 3D
Structures. Computer Graphics Forum 39, 643–666 (2020).
64. Grabner, H., Gall, J. & Van Gool, L. What makes a chair a chair? in Proc. CVPR (2011),
1529–1536.
65. Kim, V . G., Chaudhuri, S., Guibas, L. & Funkhouser, T. Shape2Pose: Human-Centric Shape
Analysis. ACM Trans. on Graphics (Prof. of ACM SIGGRAPH 2014) 33 (2014).
66. Jiang, Y ., Koppula, H. & Saxena, A. Hallucinated Humans as the Hidden Context for La-
beling 3D Scenes in 2013 IEEE Conference on Computer Vision and Pattern Recognition
(2013), 2993–3000.
67. Hu, R., Yan, Z., Zhang, J., Van Kaick, O., Shamir, A., Zhang, H., et al. Predictive and
Generative Neural Networks for Object Functionality. ACM Trans. on Graphics (Proc. of
SIGGRAPH 2018) 37 (2018).
68. Leimer, K., Birsak, M., Rist, F. & Musialski, P. Sit & Relax: Interactive Design of Body-
Supporting Surfaces. Computer Graphics Forum 37, 349–359 (2018).
69. Leimer, K., Winkler, A., Ohrhallinger, S. & Musialski, P. Pose to Seat: Automated design
of body-supporting surfaces. Computer-Aided Geometric Design 79, 79:1–1 (2020).
70. James, D. & Pai, D. Modeling deformation of linear elastostatic objects. Haptic Rendering:
Foundations, Algorithms, and Applications, 395–420 (2008).
71. Bendsoe, M. & Sigmund, O. Topology Optimization (Springer-Verlag Berlin Heidelberg,
2004).
72. Liu, H., Hu, Y ., Zhu, B., Matusik, W. & Sifakis, E. Narrow-band Topology Optimization on
a Sparsely Populated Grid. ACM Trans. on Graphics (Prof. of ACM SIGGRAPH Asia 2018)
37, 251:1–251:14 (2018).
73. Langlois, T., Shamir, A., Dror, D., Matusik, W. & Levin, D. I. W. Stochastic Structural
Analysis for Context-aware Design and Fabrication. ACM Trans. on Graphics (Proc. of
ACM SIGGRAPH Asia 2016) 35, 226:1–226:13 (2016).
90
74. Ulu, E., Mccann, J. & Kara, L. B. Lightweight Structure Design Under Force Location
Uncertainty. ACM Trans. on Graphics (Proc. of ACM SIGGRAPH 2017) 36 (2017).
75. Van Dijk, N. P., Maute, K., Langelaar, M. & Van Keulen, F. Level-set methods for structural
topology optimization: a review. Structural and Multidisciplinary Optimization 48, 437–
472 (2013).
76. Shugrina, M., Shamir, A. & Matusik, W. Fab Forms: Customizable Objects for Fabrication
with Validity and Geometry Caching. ACM Trans. on Graphics (Proc. of ACM SIGGRAPH
2015) 34 (2015).
77. Yumer, M. E., Chaudhuri, S., Hodgins, J. K. & Kara, L. B. Semantic Shape Editing Using
Deformation Handles. ACM Trans. on Graphics (Proc. of ACM SIGGRAPH 2015) 34, 86:1–
86:12 (2015).
78. Umetani, N., Igarashi, T. & Mitra, N. J. Guided Exploration of Physically Valid Shapes for
Furniture Design. ACM Trans. on Graphics (Proc. of ACM SIGGRAPH 2012) 31 (2012).
79. Wang, Z., Song, P., Isvoranu, F. & Pauly, M. Design and Structural Optimization of Topo-
logical Interlocking Assemblies. ACM Trans. on Graphics (Prof. of ACM SIGGRAPH Asia
2019) 38 (2019).
80. B¨ acher, M., Whiting, E., Bickel, B. & Sorkine-Hornung, O. Spin-It: Optimizing Moment of
Inertia for Spinnable Objects. ACM Trans. on Graphics (Proc. of ACM SIGGRAPH 2014)
33, 96:1–96:10 (2014).
81. Pr´ evost, R., Whiting, E., Lefebvre, S. & Sorkine-Hornung, O. Make It Stand: Balancing
Shapes for 3D Fabrication. ACM Transactions on Graphics (SIGGRAPH 2013) 32, 81:1–
81:10 (2013).
82. Umetani, N., Koyama, Y ., Schmidt, R. & Igarashi, T. Pteromys: Interactive Design and
Optimization of Free-formed Free-flight Model Airplanes. ACM Trans. on Graphics (Prof.
of ACM SIGGRAPH 2014) 33 (2014).
83. Wang, L. & Whiting, E. Buoyancy Optimization for Computational Fabrication. Computer
Graphics Forum (Proc. of Eurographics) 35 (2016).
84. Chen, D., Levin, D. I., Matusik, W. & Kaufman, D. M. Dynamics-Aware Numerical Coars-
ening for Fabrication Design. ACM Trans. on Graphics (Proc. of ACM SIGGRAPH 2017)
34 (2017).
85. Miguel, E., Bradley, D., Thomaszewski, B., Bickel, B., Matusik, W., Otaduy, M. A., et al.
Data-Driven Estimation of Cloth Simulation Models. Computer Graphics Forum (Proc. of
Eurographics) 31 (2012).
86. P´ erez, J., Thomaszewski, B., Coros, S., Bickel, B., Canabal, J. A., Sumner, R., et al. Design
and Fabrication of Flexible Rod Meshes. ACM Trans. Graph. 34 (2015).
87. Montes, J., Thomaszewski, B., Mudur, S. & Popa, T. Computational Design of Skintight
Clothing. ACM Trans. on Graphics (Proc. of ACM SIGGRAPH 2020) 39 (2020).
88. Bucki, M., Luboz, V ., Perrier, A., Champion, E., Diot, B., Vuillerme, N., et al. Clinical
workflow for personalized foot pressure ulcer prevention. Medical engineering & physics
38, 845–853 (2016).
89. Luboz, V ., Bailet, M., Grivot, C. B., Rochette, M., Diot, B., Bucki, M., et al. Personalized
modeling for real-time pressure ulcer prevention in sitting posture. Journal of tissue viability
27, 54–58 (2018).
90. Administration, U. G. S. Ergonomic Seating Adjustment Guide https://www.gsa.gov/
cdnstatic/GSA\_Ergonomic_Seat_Adjustment_Guide.pdf. 2020.
91
91. Kyung, G. & Nussbaum, M. A. Specifying comfortable driving postures for ergonomic
design and evaluation of the driver workspace using digital human models. Ergonomics 52,
939–953 (2009).
92. Kyung, G. & Nussbaum, M. A. Driver sitting comfort and discomfort (part II): Relationships
with and prediction from interface pressure. International Journal of Industrial Ergonomics
38, 526–538 (2008).
93. De Looze, M. P., Kuijt-Evers, L. F. & Van Dieen, J. Sitting comfort and discomfort and the
relationships with objective measures. Ergonomics 46, 985–997 (2003).
94. Savonnet, L., Wang, X. & Duprey, S. Finite element models of the thigh-buttock complex
for assessing static sitting discomfort and pressure sore risk: a literature review. Computer
Methods in Biomechanics and Biomedical Engineering 21, 379–388 (2018).
95. Pai, D. K., Rothwell, A., Wyder-Hodge, P., Wick, A., Fan, Y ., Larionov, E., et al. The human
touch: measuring contact with real human soft tissues. ACM Trans. on Graphics (Proc. of
ACM SIGGRAPH 2018) 37, 58 (2018).
96. Chadwick, J. E., Haumann, D. R. & Parent, R. E. Layered Construction for Deformable
Animated Characters in Prof. of ACM SIGGRAPH 1989 (1989), 243–252.
97. Pauly, M., Pai, D. K. & Guibas, L. J. Quasi-rigid objects in contact in Proceedings of the
2004 ACM SIGGRAPH/Eurographics symposium on Computer animation (2004), 109–119.
98. Terzopoulos, D. & Witkin, A. Physically based models with rigid and deformable compo-
nents. IEEE Computer Graphics and Applications 8, 41–51 (1988).
99. Metaxas, D. & Terzopoulos, D. Dynamic Deformation of Solid Primitives with Constraints
in Proceedings of ACM SIGGRAPH 1992 (1992), 309–312.
100. Cani-Gascuel, M.-P. Layered Deformable Models with Implicit Surfaces in Proceedings of
Graphics Interface Conf. (1998), 201–208.
101. Galoppo, N., Otaduy, M. A., Mecklenburg, P., Gross, M. & Lin, M. C. Fast Simulation of
Deformable Models in Contact Using Dynamic Deformation Textures in 2006 ACM SIG-
GRAPH / Eurographics Symp. on Computer Animation (2006), 73–82.
102. Galoppo, N., Otaduy, M. A., Tekin, S., Gross, M. & Lin, M. C. Soft Articulated Characters
with Fast Contact Handling. Computer Graphics Forum 26, 243–253 (2007).
103. Kim, M., Pons-Moll, G., Pujades, S., Bang, S., Kim, J., Black, M. J., et al. Data-Driven
Physics for Human Soft Tissue Animation. ACM Trans. Graph. 36 (2017).
104. Romero, C., Otaduy, M. A., Casas, D. & Perez, J. Modeling and Estimation of Nonlinear
Skin Mechanics for Animated Avatars. Computer Graphics Forum (Proc. of Eurographics)
39 (2020).
105. Saul, G., Lau, M., Mitani, J. & Igarashi, T. SketchChair: An All-in-one Chair Design System
for End Users in Proceedings of the Fifth International Conference on Tangible, Embedded,
and Embodied Interaction (2011), 73–80.
106. Zheng, Y ., Liu, H., Dorsey, J. & Mitra, N. J. Ergonomics-Inspired Reshaping and Explo-
ration of Collections of Models. IEEE Trans. on Visualization and Computer Graphics 22,
1732–1744 (2016).
107. Artec3D. Eva Scanner, http://www.artec3d.com 2018.
108. Brochu, T. & Bridson, R. Robust topological operations for dynamic explicit surfaces. SIAM
Journal on Scientific Computing 31, 2472–2493 (2009).
109. Taubin, G. A Signal Processing Approach to Fair Surface Design in Proc. of ACM SIG-
GRAPH 1995 (1995), 351–358.
92
110. Kikuuwe, R., Tabuchi, H. & Yamamoto, M. An edge-based computationally efficient for-
mulation of Saint Venant-Kirchhoff tetrahedral finite elements. ACM Trans. on Graphics
28, 1–13 (2009).
111. Barbiˇ c, J. & James, D. L. Six-DoF haptic rendering of contact between geometrically com-
plex reduced deformable models. IEEE Transactions on Haptics 1, 39–52 (2008).
112. Liang, J., Lin, M. & Koltun, V . Differentiable Cloth Simulation for Inverse Problems in Ad-
vances in Neural Information Processing Systems (eds Wallach, H., Larochelle, H., Beygelz-
imer, A., d’Alch´ e-Buc, F., Fox, E. & Garnett, R.) 32 (2019), 772–781.
113. FujiFilm. Fujifilm Prescale Pressure Measurement Film, https://fujiprescalefilm.com 2020.
114. Tekscan. http://www.tekscan.com 2020.
115. AljaSafe. SmoothOn Inc. www.smooth-on.com. 2020.
116. Artec3D. Spider Scanner, http://www.artec3d.com 2018.
117. Wang, Y ., Jacobson, A., Barbiˇ c, J. & Kavan, L. Linear subspace design for real-time shape
deformation. ACM Trans. on Graphics (TOG) 34, 57 (2015).
93
Appendices
A Contact Gradients and Hessian of Contact Energy
Equation 7.10 defines the energy of collision and the collision force can be computed through its
gradient. To use an implicit integrator, the gradient of the force (the hessian of the energy) is
also needed. Suppose the three vertices of the triangle that contains the closest point to p
a
on the
supporting surface are p
1
, p
2
and p
3
. Then, the position of p
anchor
is a function of p
a
, p
1
, p
2
and
p
3
. Observe thatR
3
can be decomposed into 7 disjoint V oronoi regions, one for each vertex, edge
and the face interior. The query point p
a
is in one of the regions. As we perturb p
a
, we will remain
in the same V oronoi region (unless on the boundary). This means that the closest point will remain
the same vertex, slide on the edge, or slide on the face, depending on the specific case; and hence
the gradient and hessian formulas must be derived separately for each case. Although the gradient
is not continuous across V oronoi region boundaries, such a construction is vastly superior to not
making the gradient V oronoi-region aware; it made the difference of the optimizer converging
slowly to not at all, to consistently converging.
A.1 Vertex anchor
When p
anchor
is a vertex p
i
, we have
p
anchor
= p
i
, (A.1)
∂E
coll
∂p
a
= 2k
coll
(p
a
− p
i
)
T
, (A.2)
∂
2
E
coll
∂p
2
a
= 2k
coll
· I. (A.3)
94
A.2 Edge anchor
When p
anchor
is an interior point of an edge with endpoints p
i
and p
j
, we have
p
anchor
=
(p
a
− p
i
)· (p
i
− p
j
)p
i
+(p
a
− p
j
)· (p
j
− p
i
)p
j
||p
i
− p
j
||
2
, (A.4)
∂E
coll
∂p
a
= 2k
coll
(p
a
− p
anchor
)
T
I− (p
i
− p
j
)(p
i
− p
j
)
T
||p
i
− p
j
||
2
, (A.5)
∂
2
E
coll
∂p
2
a
= 2k
coll
·
I− (p
i
− p
j
)(p
i
− p
j
)
T
||p
i
− p
j
||
2
2
. (A.6)
A.3 Face anchor
When p
anchor
is an interior point of the face, we have
p
anchor
= p
a
− ((p
a
− p
1
)· n)n, (A.7)
∂E
coll
∂p
a
= 2k
coll
(p
a
− p
anchor
)
T
· nn
T
, (A.8)
∂
2
E
coll
∂p
2
a
= 2k
coll
· nn
T
, (A.9)
where n is the unit normal of the triangle.
95
Abstract (if available)
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Anatomically based human hand modeling and simulation
PDF
Plant substructuring and real-time simulation using model reduction
PDF
CUDA deformers for model reduction
PDF
Real-time simulation of hand anatomy using medical imaging
PDF
Perception and haptic interface design for rendering hardness and stiffness
PDF
Multi-scale dynamic capture for high quality digital humans
PDF
Closing the reality gap via simulation-based inference and control
PDF
Feature-preserving simplification and sketch-based creation of 3D models
PDF
Deformable geometry design with controlled mechanical property based on 3D printing
PDF
Traveling sea stars: hydrodynamic interactions and radially-symmetric motion strategies for biomimetic robot design
PDF
Scalable dynamic digital humans
PDF
Thermal analysis and multiobjective optimization for three dimensional integrated circuits
PDF
Inference of computational models of tendon networks via sparse experimentation
PDF
Compositing real and virtual objects with realistic, color-accurate illumination
PDF
Data-driven 3D hair digitization
PDF
Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures
PDF
Deformation control for mask image projection based stereolithography process
PDF
Design and use of a biomimetic tactile microvibration sensor with human-like sensitivity and its application in texture discrimination using Bayesian exploration
PDF
Computational fluid dynamic analysis of highway bridge superstructures exposed to hurricane waves
PDF
Hybrid methods for robust image matching and its application in augmented reality
Asset Metadata
Creator
Zhao, Danyong
(author)
Core Title
Acquisition of human tissue elasticity properties using pressure sensors
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Computer Science
Degree Conferral Date
2022-05
Publication Date
12/20/2021
Defense Date
10/12/2021
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
contact,elasticity,ergonomics,FEM,finite element method,force sensor,human skin,material measurement,OAI-PMH Harvest,optimization,shape design
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Barbic, Jernej (
committee chair
), Chen, Yong (
committee member
), Kuo, C.-C. Jay (
committee member
)
Creator Email
danyongz@usc.edu,zhdy1991@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC18807271
Unique identifier
UC18807271
Legacy Identifier
etd-ZhaoDanyon-10312
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Zhao, Danyong
Type
texts
Source
20211223-wayne-usctheses-batch-906-nissen
(batch),
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
elasticity
ergonomics
FEM
finite element method
force sensor
human skin
material measurement
optimization
shape design