University of Southern California Dissertations and Theses

Artistic control combined with contact and elasticity modeling in computer animation pipelines

(USC Thesis Other)

Artistic control combined with contact and elasticity modeling in computer animation pipelines

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Artistic Control Combined with Contact
and Elasticity Modeling in Computer
Animation Pipelines
Yijing Li
Department of Computer Science
University of Southern California
Committee:
Dr. Jernej Barbiˇ c (Chair)
Dr. C.-C. Jay Kuo
Dr. Aiichiro Nakano
This dissertation is submitted to the
FACULTY OF THE USC GRADUATE SCHOOL
for the degree of
Doctor of Philosophy
Viterbi School of Engineering August 2019
Special thanks to my advisor, Prof. Jernej Barbiˇ c, for his enthusiasm, guidance and support.
I feel extremely grateful to have him as my PhD advisor.
I would also like to thank my loving parents for their encouragement,
the committee memebers for constructive comments and feedbacks,
and my lab mates Yili Zhao, Hongyi Xu, Danyong Zhao, Bohan Wang and Mianlun Zheng for
discussions and help.
Abstract
Although physically based deformations and contact are well-understood in science, their
application to practical computer animation pipelines poses several challenges. In computer
animation practice, characters, plants and inanimate objects are represented as triangle meshes,
and these meshes easily and often inter-penetrate, both during modeling and animation. It is
very difficult to design algorithms that resolve contact, yet simultaneously permit the artists to
retain direction and control over the animations.
This thesis introduces several methods to address intersections in various stages of computer
animation pipelines and improve their physical realism, all the while retaining artistic control
over the modeling and animation process. The first problem that we tackle is how to equip
standard shape modeling algorithms with robust collision response handling. Shape deformation
is a well-studied problem in computer graphics with numerous articles published on the topic
over the years. However, to date, there has been very little work on adding contact handling
to shape deformation modeling; and the existing methods do not work very well. To resolve
collisions when an artist is deforming a mesh, we give an interactive hierarchical reduced shape
editing algorithm. The algorithm can achieve interactive rates by adding constraints adaptively,
and by defining and computing a hierarchy of shape deformation bases. At runtime, a reduced
version of the as-rigid-as-possible (ARAP) energy is minimized to deform the shape interactively
subject to user handle constraints. Ours is the first algorithm to demonstrate interactive and
robust shape modeling in the presence of complex contacts and complex geometry.
Once the above shape modeling method has produced self-collision-free shapes, the next chal-
lenge to tackle is how to create controllable animations that properly resolve self-contact, all
the while following equations of motion as closely as possible. Our method can generate anima-
tions from scratch, or, more importantly in practice, correct existing non-physical animations
to preserve contact and add secondary physically based motion. The input motions may be
obtained from any method, such as keyframe animation, character rigging, 3D scanning, or
iv
geometric shape modeling. The input animation sequences can be non-physical, crude or even
incomplete. A tetrahedral mesh is first created to represent the embedded geometry and a fitting
process is then optimized to create a tetrahedral mesh animation that conforms to the input.
Then the system utilizes physically based Finite Element Method (FEM) simulation to add new
physical motions and resolve collisions. The artist provides weights for how much physically
based simulation should be allowed to modify the animation in any region of the model, and in
time. This permits smoothly turning physics on and off over space and time, making it possible
for the output to strictly follow the input, to evolve purely based on physically based simulation,
and anything in between.
Both the shape modeling method and animation modeling method discussed above require
self-intersection-free rest triangle meshes as inputs. If triangle meshes intersect already in the
input, then the intersecting parts will be “welded” together at the collision sites, precluding any
correct deformation shape modeling and animation. The third contribution of this thesis is to
overcome this limitation. We introduce an algorithm to analyze self-intersecting shapes and to
build non-welding self-intersecting tetrahedral meshes that conform to the input meshes. The
algorithm uses knowledge from algebraic and geometric topology to construct immersions that
permit meshing self-intersecting shapes without welding. We describe our algorithm in detail,
prove its correctness, and demonstrate its efficiency on several challenging examples.
The above three contributions permit the artist to create high-quality contact-aware controllable
animations that are compatible with realistic modeling and workflows in computer animation
practice. We demonstrate these algorithms using several realistic examples of controllable
contact between objects with complex geometry.
Table of contents
List of figures viii
List of tables xx
1 Introduction 1
1.1 Contact-Aware Editing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Enriching Animations With Physics . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Resolving Rest Shape Intersections . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Related Work 4
2.1 Collision Detection and Real-Time Response . . . . . . . . . . . . . . . . . . . . 4
2.2 Surface Mesh Modeling and Deformation . . . . . . . . . . . . . . . . . . . . . . 5
2.2.1 Modeling with Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2.2 Interactive Physically Based Modeling . . . . . . . . . . . . . . . . . . . 6
2.3 Multi-resolution hierarchies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Tetrahedral Mesh Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.5 Adding Physics on Animations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5.1 Simulation on Rigs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5.2 Simulation on Coarse Animations . . . . . . . . . . . . . . . . . . . . . 9
2.6 Animation Design and Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.7 V olumetric Mesh Tracking and Fitting . . . . . . . . . . . . . . . . . . . . . . . . 10
2.7.1 Inverse Caging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.8 Self-intersection Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.8.1 2D Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.8.2 3D Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Table of contents vi
3 Interactive Editing With Contact Handling 15
3.1 Shape Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Modeling in Full Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 Modeling in a Shape Subspace . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4 Multi-Resolution Shape Deformation Basis . . . . . . . . . . . . . . . . . . . . . 21
3.4.1 Single-Level Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4.2 Hierarchical Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.4.3 Assembling Bases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.4.4 Speed and Memory Considerations . . . . . . . . . . . . . . . . . . . . . 31
3.4.5 Basis Update Oracle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.4.6 Uniqueness of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.5 Collision Detection and Handling . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.1 Collision Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.2 Contact Constraint Formulation . . . . . . . . . . . . . . . . . . . . . . . 34
3.5.3 Updating the Contact Constraints . . . . . . . . . . . . . . . . . . . . . . 35
3.5.4 Static Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.6 Solving the Linear System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4 Adding Physics Using FEM Simulation 43
4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2 Tetrahedral Mesh Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.3 Tetrahedral Mesh Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.4 Constrained Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.5 Preserving High-Frequency Detail . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.6 Balancing Input Motion with Physics . . . . . . . . . . . . . . . . . . . . . . . . 53
4.6.1 Implicit Blending Force . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.6.2 Kinematic Blending . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.6.3 Generating Weights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.7 Collisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.7.1 Collision Mesh Creation . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.7.2 Collision Detection and Response During Simulation . . . . . . . . . . 57
4.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
63
Table of contents vii
5 Immersion of Self-Intersecting Solids and Surfaces 65
5.1 Topology of Self-Intersecting Meshes . . . . . . . . . . . . . . . . . . . . . . . . 66
5.2 Cell complex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2.1 Cells, Patches and Arcs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2.2 B-patches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.2.3 Computing the Cell Complex . . . . . . . . . . . . . . . . . . . . . . . . 73
5.3 Immersion graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.3.1 Simple immersions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.3.2 Rules for Graph Edges . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.3.3 Valid Graphs Give Immersions . . . . . . . . . . . . . . . . . . . . . . . 78
5.4 Algorithm to discover immersions . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.4.1 Running time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.4.2 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.5 Cell Tetrahedralization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.5.1 Virtual tets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6 Conclusion 97
6.1 Limitations and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
References 100
Appendix A Fast Right-Hand Side Evaluation for Reduced Shape Editing 110
Appendix B Immersion of Self-Intersecting Solids and Surfaces: Supplementary
Material 1: Topology 112
B.1 Topological definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
B.2 Topological proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
Appendix C Immersion of Self-Intersecting Solids and Surfaces: Supplementary
Material 2: Tetrahedralization 124
C.1 Inside-outside test without forming the region boundary surface . . . . . . . . . 124
C.2 Pseudo-normal on the piece boundary . . . . . . . . . . . . . . . . . . . . . . . . 125
C.3 Meshing self-touching cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
List of figures
2.1 A self-overlapping curve with two distinct immersions. Left: input curve.
Middle: decompositions of the two different immersions. Right: the immersion
method in this thesis is able to run on this 2D example and find both two
immersions. Green dots are handles used to pull out the mesh with ARAP
energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Input with multiple connected components in 2D. The input model consists
of two components, forming an egg-shaped hole. The immersion method in
this thesis is able to find the correct immersion and construct an “unglued”
simulation triangle mesh, which enables ARAP to pull the model apart. . . . . 12
3.1 Example of editing the chimpanzee with multiple contact sites. In the ex-
ample, the chimpanzee’s mouth was closed, forming a mouth self-contact. One
of its hands was also pulled to collide with its face, and the other to collide with
an ear. Upper-left: input triangle mesh, lower-left: tet mesh. Upper-middle:
deformed mesh with collision resolved, upper-right: activated bases: green
(L0), yellow (L1) and brown (L2). User handles are shown in blue. Lower-
middle and lower-right: zoom-in on the collision between the mouth and the
hand. Contact points are highlighted in red. The maximum computation time
at any frame in the method is 109× smaller to the maximum time in full-space
modeling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
List of figures ix
3.2 Comparison to Harmon’s method: The hierarchical method creates more
natural and smoother shapes. The example adds a handle to the armadillo and
pulls it to collide with the ground plane. The hierarchical method (left column)
then solves for smoothness and collisions together, whereas Harmon’s method
(right column) fixes collisions only as a post -process. With the hierarchical
method, the feet of the armadillo bend naturally in the initial iterations; slowly,
the entire body leans forward to minimize ARAP energy. Harmon’s method
projects the subspace control points to the closest collision-free states, moving
only control points at the feet up, creating unnatural shapes. . . . . . . . . . . . 17
3.3 Editing the Huangshan pine tree with collisions. Three handles (shown
in three yellow rectangles) are added on three branches to pull them against
other branches, creating self-collisions. The contact regions display translucent
branches and red colliding triangles. The hierarchical subspace algorithm
provides shorter and stabler computation times compared to the full method. . 18
3.4 Hierarchical bases provide a tradeoff between collision details and speed.
The top part (mustard color) of the mesh (top-left sub-figure) is pulled to cause
a self-contact against the bottom (scarlet color), as shown in the top-right sub-
figure. With level-1 local bases activated (top-right and middle-left in closer
view), collisions are resolved with visible penetrations and no deformation of
the bottom part, at 1× computation time. When level-2 local bases are added
(middle-right), the method is able to express the deformation of the bottom
part, at the cost of 1.19× of the computation time. When level-3 local bases
are added (bottom-left), the penetrations are reduced at the cost of 1.34× time.
When level-4 local bases are added (bottom-right), the penetrations are minimal.
With many modes activated, the method slows down to 3.10× time. . . . . . . . 22
3.5 Recomputing the Wang’s bases during collisions at runtime is extremely
slow, compared to the precomputed hierarchical bases. Octopus model with
20,165 vertices and 92,003 tetrahedra. . . . . . . . . . . . . . . . . . . . . . . . . 24
3.6 Examples of basis vectors at different levels of the hierarchy (second row).
Highest values are colored red and lowest are purple. Basis vectors in higher
levels (e.g., global level L0) have larger support while basis vectors in lower
levels (e.g., L2) have smaller support. Activated basis vectors due to self-
collision are displayed as colored points in the first row. They are green (L0),
yellow (L1) and brown (L2). The user handles are colored blue. . . . . . . . . . 25
List of figures x
3.7 Comparison between a different number of global manipulators: Fewer
global manipulators give un-natural motion when deformed only by the global
basis. In the left image, the chimpanzee’s left arm is not bent properly, compared
to the right image where a larger number of global manipulators smoothly
deform the arm at the elbow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.8 Building manipulators and Voronoi regions: This shows the process of form-
ing local regions to build local bases. From left to right, up to bottom: sample
level-0 (L0) manipulators (two in this case; orange); create V oronoi regions
(yellow and pink) for each manipulator; sample level-1 manipulators as children
of level-0 manipulators (green); extend each V oronoi region into overlapping
local regions (cross-hatched). The bottom-right image shows the approach
of using V oronoi regions directly as local regions. This approach produces
artifacts (Figure 3.9) and is therefore abandoned. . . . . . . . . . . . . . . . . . 28
3.9 Overlapping local region construction details, manipulator forest and ar-
tifacts of non-overlapping bases: Top-left: manipulators on part of the bunny
mesh (top-right). There are four level-0 (L0) manipulators (orange) and many
L1 ones (green). The local region (yellow) of the top-left L0 manipulator is
constructed by unioning balls centered on each vertex in the V oronoi region
of the L0 manipulator. V oronoi regions of different L0 manipulators are sepa-
rated by the blue lines. Dashed orange lines represent the actual local region
boundary on mesh vertices. Vertices inside the local region are thus covered by
the computed local basis of the L0 manipulator’s children manipulator group.
Note that this figure only shows a part of the bunny mesh. Bottom-left: the
manipulator forest corresponding to the manipulators in the top-left. Right:
artifacts occur if the local regions do not overlap. In this case, local regions are
just the V oronoi regions. The bunny’s ear was pulled to collide with its back.
Only L1 bases were activated to resolve collisions in order to demonstrate the
artifacts of non-overlapping L1 bases. Observe that the non-overlapping L1
bases cannot represent non-zero deformations across V oronoi boundaries, as
highlighted in the red rectangle. . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
List of figures xi
3.10 Static damping alleviates contact limit cycles: Here, the shape of the hand
of a character was edited. Left: Six frames performed without static damping
are superimposed with transparency to show the high-frequency vibration on
the thumb, highlighted in a red rectangle. The thumb looks blurred because it
is the superposition of several frames of a vibrating motion. Right: Six frames,
with static damping enabled (α = 0.5). Same (recorded) handle motion as in
the Left image. The thumb does not vibrate and the shape converges, therefore
the superposition is sharp. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.11 Comparison between using vs not using the additional thread to computeF
i
.
Klara cloth example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.12 Editing the bunny with self-contacts. Top row: the rest shape. Middle row:
edited shape with both ears in contact (left) and the same handle movement but
without contact handling (right). Bottom row: back view of the edited shape
with (left) and without (right) contact handling. . . . . . . . . . . . . . . . . . . 40
3.13 Editing the armadillo with self-contacts. Left: the rest shape. Upper-right:
the edited shape with one hand in contact with the head. Lower-right: same
handle movement but without contact handling. Note that a finger penetrates
into the head. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.14 Comparison between the ARAP energy and Projective Dynamics. Al-
though the energies are different, the output visual differences are minor in this
example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.15 Modeling cloth in contact with the body. The user adds two handles on the
cloth and adjusts the cloth shape while the subspace hierarchical algorithm
properly respects collisions with the body. Different levels of activated control
points are visualized in green (L0), yellow (L1) and brown (L2) in the right
figure. The user handles are shown blue and the collided manipulators are red.
11,518 cloth triangles, 32 contacts. Shape-editing runs at 43 FPS on average. . 42
List of figures xii
4.1 Preserving dynamic spatial high-frequency detail: Top: input animation of
the plant withering (aging), containing spatial high-frequency motion on the
leaves. The model and animation were downloaded from the Internet. Middle:
physically based simulation loses high-frequency detail. Bottom: the system
preserves high-frequency detail. Tetrahedral mesh can be coarse and is shown
overlaid. The fitting process can be seen as extracting the spatial low-frequency
part of the motion and modifying it using physically based simulation. The
high-frequency motion can be added back depending on needs. . . . . . . . . . 46
4.2 Self-intersections in the input mesh (human face). Left: the lips collide in
the highlighted region. Right: an interior view of the self-intersection. The
system accepts such ill-formed inputs using the method described in Chapter 5. 47
4.3 Completing a tracked face using physics: Due to occlusions, vision tracking
systems cannot capture all parts of a human head (mouth cavity, inner lips,
tongue, etc.). We utilize physically based simulation to produce the deforma-
tions of the untracked parts. Left (the input): tetrahedral mesh (18,391 tets),
the tracked triangle mesh (49,924 triangles) and the neutral complete triangle
mesh of the head. Middle, right (the output): the deformed tracked mesh and
the complete mesh computed using our method. 960 frames, 1.35 sec per frame. 47
4.4 Comparison of invertible StVK to corotational linear FEM. The input is a
severely stretched cube triangle mesh. With the same amount of total fitting
error, the shape produced by the corotational linear FEM (right) has self-
intersections along the cube edge, whereas the invertible StVK (middle) has no
such artifacts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.5 Comparison of different internal force evaluation strategies: A simple bar
vibrates with one end fixed. The yellow shape was obtained using u
0
k
= u
k
+∆tv
k
(Bergou’s initial guess). The gray shape was obtained using u
0
k
= u
k
(the other).
The initial displacement and velocity are the same in both methods. In the first
row,∆t= 0.003. Bergou’s initial guess exhibits a less-damped motion. In the
second row,∆t= 0.1. A position constraint is added to one vertex (shown as a
red dot). The other initial guess is stable under the large timestep, whereas the
Bergou’s is not. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
List of figures xiii
4.6 Preserving dynamic high-frequency detail: The difference between the in-
put vertex displacement ¯ q
i
and the interpolated displacement ¯ q
∗
i
of the fitted
triangle ¯ u
1
¯ u
2
¯ u
3
(left) is rotated by the rotation component R of the deforma-
tion gradient between the two triangles (2D example). It is then added to the
interpolated displacement q
∗
i
of the simulated triangle u
1
u
2
u
3
(right) to give the
final displacement q
i
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.7 High-frequency detail transfer failure produces spikes on the mesh (left)
when using tetrahedral deformation gradients to transfer detail from highly
compressed tetrahedra, as contrasted with the smooth shape (right) produced
by transfer with the rotation component (used in this thesis). An alternative
method to preserve high-frequency detail is to alter barycentric coordinates for
the triangle mesh at each frame. This method also suffers from such a “spike”
failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.8 Adding physics to an elephant walk cycle. Left-most column: the fixed Ω
′′
and freeΩ∖Ω
′′
region of the elephant tetrahedral mesh. Free regions include
the trunk, ears, belly and tail. Left four columns: Each column is a separate
animation with a different blending weight (Equation 4.11) between the artist
input and physics. Implicit blending method was used. From left to right, the
weights are decreased from following the input 100%, to a 100% physically
enriched result. InΓ∖Γ
′′
, the method produces secondary motion according
to the weights, whereas inΓ
′′
, the output is identical to the input. Right-most
column: kinematic blending with 50% physics. It produces similar results as
the implicit blending method, but suffers from less intuitive tuning of collision
handling, such as between the legs and belly. . . . . . . . . . . . . . . . . . . . 59
4.9 Adding physics to a keyframed sumo wrestler animation. Left-most col-
umn: the fixed Ω
′′
and freeΩ∖Ω
′′
region of the sumo tetrahedral mesh. Free
regions include the body, thighs, cheeks and derrière. The top row on the right
shows several frames of the input motion. The bottom row on the right shows
the result with secondary motion added. . . . . . . . . . . . . . . . . . . . . . . 61
List of figures xiv
4.10 Adding physics to a dragon animation. The fixed region consists of the
two back feet. The blending weights (third and fourth subfigures in the first
row) change through time to follow input at first, and add secondary motion
afterwards. Note that the weights on the back feet are ignored because they are
inΓ
′′
. The second row shows the resulting animation where the dragon rises
up. For cinematic effect, a particle-effect fire animation is added. Secondary
motion is present in the horns, whiskers, body and tail. . . . . . . . . . . . . . . 62
4.11 The method described produces more global motion than TRACKS with
40 automatically-generated regions. Two comparisons are made at each frame:
(1) the Euclidean distance between the displacement u
k
∈R
3n
and the input
frame ¯ u
k
∈R
3n
, and (2) the difference in the x-axis position of a selected dragon
vertex between the output and input. . . . . . . . . . . . . . . . . . . . . . . . . 62
4.12 Adding physics to a keyframed bird animation. The top row shows the fixed
Ω
′′
and freeΩ∖Ω
′′
bird regions. Free regions are the head, belly and wings,
except the wing tips which are fixed to drive the wings, similar to inverse
kinematics. The middle row shows several frames of the input motion. The
bottom row shows our result with subtle secondary motion added. This example
does not use blending. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.13 Comparison to as-rigid-as-possible energy: Left: the deformation produced
by our method. Right: the shape produced by as-rigid-as-possible energy [Sorkine
and Alexa 2007]. In order to improve the as-rigid-as-possible result, the neg-
ative weights in the as-rigid-as-possible energy are clamped; otherwise, the
as-rigid-as-possible method causes the obtuse mesh triangles to collapse. The
comparison were performed by setting the positions of the tracked meshΓ
′
as
constraints forΓ and minimized the as-rigid-as-possible energy. Lacking inte-
rior information, the as-rigid-as-possible energy resulted in a colliding shape at
extremely concave regions such as the mouth cavity and eyelids, whereas our
method produces a good result. . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.14 Adding physically based contact to artist animations: The soft ball is launched
at the elephant. The ears, trunk, belly and the tail are enriched with physics just
as in Figure 4.8. The system performs collision detection between the elephant
and the ball. Upon impact, contact forces are computed and added to both
objects, resulting in local contact deformations on both objects, all the while
the elephant still generally follows the original animation. . . . . . . . . . . . . 64
List of figures xv
5.1 Basic steps of the immersion method are illustrated in this 2D example. Given
the input shape (in this 2D example, a closed polygonal line / curve inR
2
), the
method extracts the cells, patches and arcs (Section 5.2), then duplicates and
connects cells to recover the immersion of a two-dimensional disk such that the
disk boundary immerses onto the input shape (Section 5.3 and 5.4). Note that
the largest cell self-touches, which is addressed in Section 5.2.2. To produce
the immersion of the disk, first triangulate the entire space covered by the input
shape, then create a submesh for each cell (Section 5.5). Some triangle vertices
are duplicated to avoid gluing parts of the cell together. For visualization
purposes here, a cyan line indicates that the two triangles joined by it are
actually the same triangle. Finally, submeshes for the cells are merged together
according to the cell connectivity recovered by the algorithm (Section 5.5). . . 68
5.2 Examples of self-intersecting curves inR
2
and surfaces inR
3
which are not
boundaries of immersions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.3 Cells, patches and arcs in 2D. Left: input 2D curve. Right: 4 cells (open
disks), 8 patches (open curves; Pi with different colors), four arcs (isolated
points). Cells and their B-patches are as follows: Cell 0: P0, P1, P2*, P3*;
Cell 1: P2, P3, P4, P6; Cell 2: P4*, P5; Cell 3: P6*, P7. Here, * denotes
that the B-patch orientation disagrees with the cell. Arcs and the topological
neighboring patches at those arcs are as follows: Arc 0: (P1,P2), (P1,P3); Arc
1: (P0,P4), (P2,P5); Arc 2: (P0,P6), (P3,P7); Arc 3: (P4,P7), (P5,P6). . . . . . 69
5.4 Cells, patches and arcs in 3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.5 Finding patch geometric neighbors on self-touching cells. Cell 0 in Fig-
ure 5.3 self-touches at Arc 0. First identify self-touching arcs (there is only
one in this example; bottom-left), then find the two components that should be
disconnected (bottom-right). Use them to identify the geometric neighbor pairs
at Arc 0 as (P2, P1) and (P1, P3). . . . . . . . . . . . . . . . . . . . . . . . . . . 71
List of figures xvi
5.6 2D examples of simple and non-simple immersions. Top: a disk-topology
abstract manifold
ˆ
S is immersed onto 2D. Its boundary
ˆ
M is immersed onto
M. For cell C
i
with winding number (WN) w, its pre-image, ˆ σ
−1
(C
i
) has
w connected components (PI-CC). The restriction of ˆ σ to each connected
component in ˆ σ
−1
(C
i
) is a surjective immersion onto C
i
, and the number of
sheets for the surjective immersion is 1. Therefore, this is a simple immersion.
Bottom: a ring-topology abstract manifold
ˆ
S is immersed onto 2D. Its boundary
ˆ
M is immersed onto a 2D multi-component curve, which creates three cells.
There are two possible immersions (each row, bottom-right). Both immersions
wind the 2D ring twice and connect it back to itself. Among the three cells,
the interesting one is the ring-shaped cell C
1
. In both immersions, C
1
has a
winding number of 2, but the number of connected components of its pre-image
is only 1 (highlighted in red in the PI-CC column). The restriction of ˆ σ to
ˆ σ
−1
(C
1
) (which has a sole connected component) is a surjective immersion
where the number of sheets is 2 (bottom-left). Intuitively, this means that ˆ σ
“wraps” ˆ σ
−1
(C
1
) onto C
1
twice. Therefore, ˆ σ is not a simple immersion. . . . 74
5.7 3D example of a non-simple immersion. Left: a torus mesh which winds
around the center axis twice, forming a 3D self-connecting cell. Intersecting
mesh faces are colored red. Right: the orange curve is the skeleton of the torus
mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.8 Immersion graph connectivity rules. Solid black curves denote “owned”
patches, dashed black curves denote “declined” patches, cyan lines are connec-
tions between nodes, and blue arrows specify patch orientations. Checkmarks
and red crosses denote that a rule is followed or violated, respectively. . . . . . 77
5.9 Examples of Rule 6: Top-left: one common case of Rule 6, where four patches
p,q,r,s and four nodes c,d,e, f from four cells join at the same arc. Bottom-left:
a rare case of Rule 6, where c and e are the same node and p and q are the same
patch. Right: examples of such situations. Red circles denote the four-patch
cases. Red rectangles denote the three-patch cases. Note that the rule does not
apply to arcs c
0
and c
1
highlighted with blue circles in the bottom-right figure,
because they both have a surrounding empty cell that does not belong to the
cell complex of M. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
List of figures xvii
5.10 An example of the execution of the immersion algorithm to construct an
immersion. The input is the 2D curve shown in the bottom-right of Figure 5.9.
Red, solid black and dashed black curves denote “undecided”, “owned” and
“declined” patches, respectively. At each iteration, a blue line represents the first
valid connection between two nodes, green lines are the subsequent connections
made by Rule 6 after the first connection, and cyan lines are previous connec-
tions. This is a complete example showing all 11 iterations of the algorithm to
produce the output immersion graph (bottom-right). . . . . . . . . . . . . . . . 82
5.11 Resolving the double-loop torus. A cube is added to cut the self-connecting
cell. Left to right column-wise: the model (collided faces in red) and the
tetrahedral mesh, the helper cube pulled away to show the torus mesh, volu-
metric ARAP deformation exposing the topology, the double-loop torus further
deformed into a single-loop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.12 Building the region graph. Top: input geometry, the steps to build the region
graph, and the final region graph with the assigned +,− labels. Blue arrows
are the triangle normals. Black empty circles are region graph nodes. Bottom:
nearly self-intersecting dragon. Naive meshing glues the mouth and back
and produces almost no motion in animation; whereas the immersion method
produces good dynamics with a coarse mesh. . . . . . . . . . . . . . . . . . . . 85
5.13 Examples of geometry inside a tetrahedron. Top row: 3 pieces partition the
tet into 4 regions. The embedded triangle mesh is shown dashed. Blue arrows
represent triangle normals. Regions labeled + are outside of the embedded
triangle mesh. Regions labeled − are inside of the triangle mesh, and are a
part ofΩ. Note that the left and right image have identical piece geometry,
but opposite +,− labels, due to the opposite piece orientation. Bottom row:
corresponding region graphs. The nodes of the+ and− regions are denoted by
empty blue and solid yellow circles, respectively. Each orange arrow represents
a directed edge corresponding to a piece. . . . . . . . . . . . . . . . . . . . . . . 85
5.14 Duplication of tet vertices to embed a comb. . . . . . . . . . . . . . . . . . . . 88
5.15 Yeahright-on-ground model. Top left: Yeahright dropped on the ground.
Collided faces are in red. Bottom left: generated tetrahedral mesh. Top right:
the Yeahright-on-ground is deformed using volumetric ARAP on the tetrahedral
mesh. Bottom right: closer look on how the immersion algorithm is able to
untangle the intersecting features. . . . . . . . . . . . . . . . . . . . . . . . . . . 90
List of figures xviii
5.16 Quintuple Torus model. Top-left: side view of the Quintuple Torus. Col-
lided faces are in red. Bottom-left: top view with holes visible. Top-middle:
Sacht’s implementation became stuck in the reverse flow due to the compli-
cated collisions, spending 8 hours without progress until manually stopped.
Bottom-middle: the immersion algorithm successfully generated a tetrahedral
mesh to embed the shape. Top-right: simulating the output mesh using FEM.
Bottom-right: the deformed tetrahedral mesh. . . . . . . . . . . . . . . . . . . . 91
5.17 Intermediate stages of building immersion graphs. The three rows corre-
spond to the examples of Sacht (small), yeahright-on-ground and quintuple
torus. Leftmost column: the first node added to the immersion graph. Rightmost
column: the final node added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.18 Unglued meshing of self-intersecting lips and mouth cavity: Column 1: side
view of the head and the self-intersecting mouth connected to an interior cavity.
Such self-intersecting modeling of the mouth cavity is very common in visual
effects and games, to add a plausible mouth appearance to characters. The
self-intersecting mouth is a common technology blocker for head simulation
algorithms. Column 2: head triangle mesh (top) and the resulting tetrahedral
mesh (bottom). Column 3: zoom to the self-intersecting lips. Column 4: the
tetrahedral mesh “unglued” the input self-intersection, enabling the simulation
to successfully open the mouth. . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.19 Self-intersection-aware tet meshing: The immersion method properly dupli-
cated and connected the tets to account for helix triangle mesh self-intersections.
Left: input triangle and tetrahedral mesh. Intersecting triangles are shown red.
Right: the output tet mesh has been successfully pulled apart using volumetric
ARAP [114]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.20 Tetrahedral meshing of self-intersecting procedurally generated geome-
try. The triangle mesh for this tree (328,152 triangles; 2,634 branches) was
generated using a procedural modeling method without regard for branch inter-
sections and self-intersections. The immersion method successfully generated
a simulation tet mesh that “unglues” all the branches. Because our tets can
be overlapping, our mesh can be relatively coarse (32,457) for a model of this
complexity, permitting faster FEM simulation. For comparison, naive meshing
glues branches and results in visibly suboptimal tree motion when simulated
using FEM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
List of figures xix
5.21 The immersion method produces topology-aware weights. Bounded bihar-
monic weights[57] are computed using the two handles “A” and “B” indicated;
the figure (left) shows the BBW weight function for handle “A”. BBWs com-
puted on our volumetric mesh (left) reflects the correct topology of the em-
bedded surface (we re-use the Sacht’s test model [97] from Figure 5.22). In
contrast, naive voxelization (right) glues the cylinders together, resulting in
unsatisfactory weights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5.22 Comparison to Sacht’s work: Part (A): Left to right: side and top view of
the model from [97] (intersections are in red), Sacht’s method result (failed in
the middle of the process), the tetrahedral mesh generated successfully by the
immersion method, and the tetrahedral mesh and its embedded mesh deformed
by volumetric ARAP, demonstrating the correct topology. Part (B): Left to
right: side and top view of the “small” input model; Sacht’s method produces
a highly deformed shape that results in a low-quality tetrahedral mesh for the
original surface; the immersion method produces a good volumetric mesh. . . 94
B.1 The generators ofπ
1
(C). In this case,C is a solid torus with 3 handles (k= 3),
viewed (on the left) under an orthographic projection from the top. . . . . . . . 119
C.1 Handling closest sites at piece boundaries. Left: a piece inside a tetrahedron
T
0
T
1
T
2
T
3
, where T
3
is occluded by the piece. Vector N
1
is the normal of the tet
face T
0
T
1
T
2
. Vertices v
i
are on the boundary of a piece. Vector n is the normal
of the triangle with the boundary edge v
1
v
0
. Top-right: degeneracy cases for
boundary edges. Bottom-right: tet face angle cases for boundary vertices. . . . 126
List of tables
3.1 Example Complexity for #mesh vertices and #tetrahedra, #average collision
constraints during contact, min/max basis size used during contact. . . . . . . 38
3.2 Runtime statistics for average computation time (at) broken down in terms of
constructing positionsp= Wv (at:pWv), constructing right-hand sides (at:rhs),
collision detection (at:col), and solving the linear system (at:solve). The av-
erage total computation time per iteration (at:total) and the average speedup
compared to the full method (at:sp) are also given. The table also reports
the maximum computation time (mt) during the manipulation (mt:total) and
its speedup (mt:sp). The implementation includes multi-threading with Intel
TBB; it parallelizes collision detection, rotation extraction, basis formation, and
computing mesh vertex positions from reduced coordinates. It also accelerates
the full method with the same amount of parallelism, for a fair comparison. No
GPU computation is employed. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3 Preprocessing Statistics for building the hierarchy (“hier”), clustering rota-
tions (“clust”), creating all basis vectors (“bases”), projectingL (“proj”). Col-
umn “total” gives the total preprocessing time for each model in Table 3.2. For
the Chimpanzee and Huangshan pine examples, rotation clustering is not used
because many clusters would be needed to accurately model deformations on
cylindrical-like shapes such as branches and limbs. Pre-processing such clusters
was relatively expensive and only gave marginal runtime speed improvements. 39
List of tables xxi
4.1 Simulation statistics. vtx, tri=#triangle mesh vertices and triangles; tet-vtx,
tet=#tet mesh vertices and tets; free tet-vtx, free tet=#tet mesh vertices and tets
in the free regionΩ∖Ω
′′
; frames=#graphical frames; sim-frames=#simulation
frames, including intermediate frames that were inserted for stable simulation;
fitting, sim = time for fitting the tetrahedral deformations and simulation for one
graphical frame, respectively;γ: parameter used in fitting. Intel Xeon 2.3GHz
2×6 cores, 32 GB memory. In the elephant example, 95% of the simulation
time is spent processing collisions. The head example only uses frame fitting;
hence, no simulation time is reported. . . . . . . . . . . . . . . . . . . . . . . . 58
5.1 Mesh complexity statistics. m.w.#: max winding number on the mesh. . . . . 89
5.2 Immersion algorithm timings on examples. t-cell, t-imm, t-mesh: times for
generating the cell complex, immersion algorithm and meshing. . . . . . . . . 89
5.3 Timing comparison to Sifakis’s algorithm. The total time spent in four key
CGAL exact arithmetic routines when using CSG operations as in [111], versus
the novel pseudo-normal and exact arithmetic Sutherland-Hodgman algorithm.
Yeah-right-on ground model. “Meshing total” is the time to generate the tet
mesh after the immersion algorithm has generated the graph, including CGAL
and other operations specific to each method. The “grand total” is the total time
from loading mesh M to producing the output tet mesh; it equals the immersion
algorithm time (16 sec; same for both methods), plus the “meshing total” time.
The libigl::mesh_boolean times are for tetrahedralization in Sifakis’s method;
the virtual tets method does not need it for tetrahedralization. it is used in the
immersion algorithm; where it occupies 12.5 sec of the running time. . . . . . 95
Chapter 1
Introduction
In computer graphics, there is an important problem that has to date not received substantial
analysis: how to efficently and robustly model physically based motion of three-dimensional
solids represented as polygon meshes (characters, tissue, plants, creatures, etc.) while simul-
taneously permitting artistic direction and control. This problem is especially difficult when
three-dimensional solids undergo contact or self-contact.
In a typical animation workflow, artists first generate surface representations of characters and
then bind skeletons or rigs onto the meshes to define how the mesh vertices move according to
the rig parameters. For facial performances, blendshapes are also used to create animations.
The tuning of rig parameters and blendshape weights often comes with large amount of tedious
manual work. With recent advances of mesh-based optimization, physically based simulation
and physically based contact handling, many virtual assets can be deformed or simulated with
ease. However, utilizing those techniques in combination with traditional computer animation
workflows remains an interesting problem. The goal of this thesis is to introduce several methods
to address collisions, self-collisions and complex nonlinear elasticity of three-dimensional
deformable models, while simultaneously retaining artistic control over the outputs.
1.1 Contact-Aware Editing
Interactive mesh editing techniques have been developed for designing new shapes by deforming
them from a neutral polygon mesh. Such techniques, however, often neglect potential collisions
1.2 Enriching Animations With Physics 2
during the manipulation. These collisions create visual artifacts and increase the difficulties
for physically based simulation methods in later stages. Complex physically based simulation
methods are able to resolve collisions accurately, but at the cost of long computation times. In
light of the requirement of intuitive control and fast computation during mesh editing tasks, this
thesis proposes a new reduced hierarchical method to address the challenge.
The method utilizes a reduced space of the original shape space to accelerate computation. At
each iteration, it reduces the well-known and widely used as-rigid-as-possible (ARAP) energy,
which describes how smooth the deformation is, subject to user handle constraints. To resolve
contacts at runtime, the method pre-computes a hierarchy of local bases, each of which is
smooth and gives additional degrees of freedom. When a contact site is detected, the method
adaptively adds new local bases to the system, increasing the degrees of freedom in the solver
to allow realistic local contact responses, while adding minimal computation overhead to the
application. With this method, users can control the tradeoff between contact accuracy and
interactivity.
1.2 Enriching Animations With Physics
Artists have many tools to create animations, ranging from motion capture devices using
computer vision to various deformers and rigs in computer animation software. To be compatible
with existing computer animation pipelines, this thesis proposes a system that corrects any
animation input (such as obtained using blendshapes or a skeleton rig) to more closely respect
physics and contact mechanics. The input of the system is either an animation of the entire
triangle mesh, or a subpart of it. Our system then generates physically based animations that
follow the input as closely as possible, while adding secondary motion and resolving collisions
and self-collisions.
The input animation is first converted into a simulation-compatible structure, such as a tetrahe-
dral mesh. The artist then chooses which part of the mesh to make physical, and assigns weights
to provide the balance between pure physics-driven motion and strict following of the input
animation. The system utilizes Finite Element Method (FEM) on deformable objects to simulate
the tetrahedral mesh with boundary constraints. Two methods, the implicit blending force and
kinematic blending are provided to balance physics with input according to artist-provided
weights. The physically based simulation includes collision detection and handling to remove
1.3 Resolving Rest Shape Intersections 3
penetration artifacts. Finally, the tetrahedral mesh animation is interpolated to the triangle mesh,
resulting in a triangle mesh animation that is collision-free and enriched with physically based
secondary dynamics.
1.3 Resolving Rest Shape Intersections
Both the interactive shape editing methods and physically based animation systems require
self-intersection-free rest triangle meshes as inputs. Common tetrahedral meshing methods
will glue the self-intersecting parts of the mesh, creating undesired tetrahedral mesh topology.
To address this limitation, this thesis introduces a topological approach to analyze a valid
self-intersecting mesh and build a compatible tetrahedral mesh for it. The resulting tetrahedral
mesh is self-intersecting, but shares the same topology as the triangle mesh, thereby allowing
the modeling and animation stages in the pipeline to function correctly. To rigorously solve
the problem, we study the mathematical concept of immersions, and give a deterministic and
constructive algorithm to determine if the input self-intersecting triangle mesh is the boundary of
an immersion. For near self-intersections, a robust algorithm is also given to properly duplicate
mesh elements and correctly embed the underlying geometry into the mesh element copies.
Both the self-intersections and near self-intersections are combined into one algorithm that
permits successful meshing at arbitrary resolution.
1.4 Thesis Overview
The rest of this thesis is organized as follows. Chapter 2 gives a general survey of prior work on
contact, modeling, shape and animation editing and simulation. Chapter 3 describes how the
reduced hierarchical basis is formed and how we perform contact-aware editing. In Chapter 4,
the details of tetrahedral mesh fitting and adding physics with simulation are given. Chapter 5
analyzes self-intersecting shapes and gives our algorithm to convert them into topologically-
consistent tetrahedral meshes. Finally Chapter 6 concludes the thesis and lists limitations and
directions for future work.
Chapter 2
Related Work
2.1 Collision Detection and Real-Time Response
To handle collisions in various graphics applications, one needs to first detect collisions. A
survey of real-time collision detection was given by Ericson [35]. In haptics and surgery
simulation, real-time contact feedback is widely studied [71]. Stable 6-DOF haptic manipulation
of rigid models can be achieved with implicit integration [93]. Rigid body and reduced
deformable object collisions can also be simulated at haptic rates [9, 64, 140]. Barbic and
James [9] employed a hierarchical point tree to provide an upper limit on collision detection
time. The interactive editing system described in this thesis exploits a similar level-by-level
collision detection, but demonstrates how to efficiently combine it with hierarchical deformation
bases for efficient shape editing. Allard et al. [ 2] introduced volume constraints that can control
contact response detail. Their method operates on the full model, whereas the hierarchical
method in this thesis uses a reduced model for faster computation. The method in this thesis can
achieve good computational speedups compared to the full method because it is designed such
that the reduced system matrix does not change unless there is a change of basis. Even then,
only the part corresponding to the changed basis changes. Talvas et al. [117] used nonuniform
pressure on volume constraints for fast hand simulation with contact. Their method can be used
on top of this hierarchical constrained system to simplify constraints.
For reduced simulation, several publications added more degrees of freedom into contact
regions [50, 120]. Different from their work, our method creates a hierarchy of bases to allow
2.2 Surface Mesh Modeling and Deformation 5
better control on tradeoff between contact resolution and speed. This iterative algorithm is
designed for shape editing and modeling, while theirs are used in reduced physically based
simulation. This method uses constraint-based techniques to handle collision, while they
use penalty forces that are known to be difficult in determining penalty stiffness. Recently,
Erleben [36] assessed various contact normal methods and concluded that no good local normals
are available. The method in this thesis shares the limitation but it is possible to alleviate the
problem by smoothing normals.
2.2 Surface Mesh Modeling and Deformation
Polygonal meshes are the most commonly used data structure for rendering and animating. The
most commonly used method to deform a mesh is the linear blend skinning, which requires a
low-dimensional space of bones inside a skeleton to control all the vertex positions on the mesh.
The linear blend skinning, and other deformers (Free-Form Deformation (FFD), blend shapes,
cage and nonlinear deformers in animation software such as Maya) define a direct nonlinear
mapping between the low-dimensional parameter space and the high-dimensional mesh vertex
positions. Another category of deformation techniques uses interactive optimization methods
to search for vertex positions that best minimize a given energy that usually describes the
smoothness of the mesh, subject to user constraints. These techniques can be seen as using
an indirect nonlinear mapping from the low-dimensional parameter space created by user
constraints to the high-dimensional mesh vertex positions, via energy minimization. A well-
known example is as-rigid-as-possible (ARAP) energy [114]. All those methods are designed
with efficiency in mind, enabling deformations of complex mesh with millions of vertices and
polygons. However, they fail to meet the needs of solving collisions for realisitic animations.
2.2.1 Modeling with Contact
To overcome the problem of contact during modeling, researchers proposed several solutions.
Vaillant et al. [127] used implicit distance fields to form a collision free mesh. Gain and
Dodgson [40] developed a self-intersection-free FFD scheme. Harmon et al. [49] resolved
contacts by minimizing space-time interference volumes, using the controls of the given
modeling tool. Their contact resolution method is essentially a post-processing step to modeling.
2.2 Surface Mesh Modeling and Deformation 6
Different from their approach, the interactive editing system described in this thesis integrates
contact handling directly into the modeling loop and provides the flexibility of adding localized
bases to address local contact handling. This method is compared to Harmon’s method in
Chapter 3, demonstrating that the method in this thesis is more stable and produces better
shapes.
2.2.2 Interactive Physically Based Modeling
The recent introduction of projective dynamics [20], its follow-up methods [22, 73] and
ADMM [88] make it possible to perform plausible real-time physically based simulation
on complex meshes. These methods are also directly applicable to interactive shape modeling,
and follow the same local/global optimization mechanism as the ARAP method. The method
in this thesis utilizes the ARAP method for interactivity, but adds one extra level of reduction
by optimizing in a hierarchical reduced space, allowing even more complex meshes to be
edited with contact. The basis hierarchy could be combined with projective dynamics / ADMM
methods to replace the ARAP energy. Note that Brandt et al. [22] used model reduction with
local/global optimization, too. Their work, however, cannot handle self-collisions, whereas the
method in this thesis works with both external and self-collisions. It requires a more complex
constraint handling mechanism (as presented in Chapter 3).
Generally, the method in this thesis differs from real-time deformable object simulation methods
with collision handling in that the specific application targeted is geometric deformation, which
does not require a per-frame accurate result, but focuses on interactivity and convergence. In
contrast, simulation methods focus on speed and physical plausibility. For shape modeling, one
needs speed, control, and gradual refinement of the shape. The system in this method satisfies
these requirements. It enables direct local control of shapes, gradually refines the shape, and
produces smooth and shape-aware deformations. As per collision constraints, the formulation
used is similar to other vertex-constraint-based simulation methods. The difference is that this
method can exploit the system structure to optimize computation (Section 3.5.2 and 3.6).
2.3 Multi-resolution hierarchies 7
2.3 Multi-resolution hierarchies
Multi-resolution hierarchies are commonly used for modeling, energy minimization and simula-
tion. Boier-Martin et al. [16] surveyed subdivision methods and multiscale modeling. Botsch et
al. [19] and Winkler et al. [135] performed hierarchical shape matching. Coarse meshes are
commonly used to accelerate nonlinear optimization [39, 76, 126]. Several researchers studied
multi-resolution FEM [26, 32, 44], and the multigrid method [92, 118, 143]. James and Pai [60]
modeled multi-resolution deformations using Green’s functions and the capacitance matrix
algorithm. Their method is limited to linear elasticity which produces artifacts under large
rotations. The interactive editing method in this thesis gives a multi-resolution basis hierarchy
suitable for geometric shape modeling with the ARAP energy in the presence of localized
contact.
Malgat et al. [75] gave a technique to combine two or more deformation models in separate
hierarchical layers. Their core innovation is a generally very useful framework whereby
deformations due to mechanical models at deeper layers are properly locally transformed
with time-varying affine transformations imposed by the ancestor layers in the hierarchy.
They describe how to create mass matrices and equations of motions consistent with these
assumptions; and give a dynamics filter to avoid redundant DOFs from different layers. The
method of this thesis formulates the total world-coordinate position of any mesh vertex as a
linear superposition of fixed shape vectors that do not rotate or affinely transform at runtime
(Equation 3.15). It therefore does not need the local affine transformability accommodated by
Malgat’s work. No filters are required and this basis hierarchy automatically avoids redundant
DOFs (Section 3.4.2). Malgat et al. [75] demonstrated applications to collision handling.
However, in their examples, they manually specified parts of the mesh whereby detailed
deformations are needed due to collisions. In this thesis, it activates and de-activates the local
detail automatically. It can be done efficiently when the collision sites are not known in advance.
2.4 Tetrahedral Mesh Modeling
Tetrahedral meshes are commonly used in shape editing and deformable object simulation.
Existing tools such as Tetgen [48] and NetGen [101] generate good-shaped tetrahedral meshes
given manifold surface triangle meshes. Other algorithms like [80] and [3] provide good
2.5 Adding Physics on Animations 8
tetrahedralization strategies that are suitable for simulation purposes. Most of prior methods
take as input intersection-free closed manifold surface triangle meshes.
For input geometry that is non-manifold, point clouds, or meshes with multiple intersected
components, one can convert the geometry into a distance field and apply surface meshing
techniques [17] to get a satisfying triangle mesh for tetrahedralization. Alternatively, one
can use Constructive Solid Geometry (CSG)-like operations to form a “union” of the mesh
and tetrahedralize it [53]. The disadvantage of such treatments is that the input may include
self-intersection artifacts. Self-intersections should be removed beforehand. However, both
the distance field and the union operation mistakenly “glue” the intersecting components. To
eliminate these artifacts before the tetrahedralization process, Sacht et al. [97] deformed a self-
intersecting mesh into an intersection-free mesh by leveraging conformalized mean-curvature
flow. However, we show that our method is faster, more robust and more general: it supports
high-genus shapes.
2.5 Adding Physics on Animations
To create realistic motions, physically based simulation is widely studied for decades. Many
researchers tried to integrate simulation with rigs to provide better simulation control for artists.
Another direction is to improve animation tools for direct manipulation like wiggly splines
used to create physical oscillatory motion [63]. The problem of simulating deformable objects
has been well-studied since the pioneering work of [123]. Several publications have improved
simulation versatility and efficiency, for example [8, 32, 77, 86], and robustness [55, 122].
2.5.1 Simulation on Rigs
Rigs are widely used in animation pipelines. A natural idea is to combine rigs with simulation
to create soft tissue motions on characters and animals. The simulated mesh is driven by the
motion of the embedded skeleton or character rig [25, 37, 41, 65, 66, 72, 77]. Shi et al. [104]
enhanced skeleton-driven animations with secondary motion extracted from sample sequences.
Hahn et al. [46, 47] formulated the equations of motion in the rig space so that some rig
parameters can evolve via physics. Different from these previous methods that generate soft
2.6 Animation Design and Control 9
tissue motion from the character rig, the simulation method in this thesis does not require a rig
input, and works with and corrects an existing triangle mesh animation.
2.5.2 Simulation on Coarse Animations
Several methods [14, 85, 96, 129] enriched coarse cloth animations with spatial high-frequency
effects such as wrinkles and folds. Among them, Bergou et al. [14] achieved this by running a
detailed mesh simulation constrained to the input coarse animation. The simulation method in
this thesis also uses constraint-based dynamics, but with a different design goal. It modifies the
coarse (spatially low-frequency) component of the animation, as representable by the chosen
simulation mesh, and preserve high-frequency motion, whereas Bergou’s method preserves
coarse motion and generates high-frequency detail. A comparison to Bergou’s method is made
in Chapter 4.
2.6 Animation Design and Control
Several work addressed the design and editing of elastic physically based simulations using
keyframes [10, 52, 54] or spacetime constraints [68, 69, 102]. Additionally, Kondo et al. [67]
allowed animation editing by both keyframes and trajectories. Barbiˇ c et al. [12] directed
physical simulations to given input trajectories leveraging LQR control in the reduced space.
Coros et al. [30] controlled deformable objects by changing their rest shapes. Barbiˇ c et al. [13]
introduced spacetime Green’s functions for interactive animation editing. These papers mainly
focus on direct user control. Although they can enrich animations with physics, they require
special models (e.g., model reduction or specific definitions of elastic forces), or operate under
specific assumptions. This thesis aims to achieve the goal of adding physical motions on
existing animations offline with a more general approach, thereby using common methods such
as Newton’s method and full-space FEM with boundary constraints.
2.7 V olumetric Mesh Tracking and Fitting 10
2.7 Volumetric Mesh Tracking and Fitting
When fitting a volumetric mesh animation to a triangle mesh animation, the surface vertices
of the volumetric mesh are usually over-determined by the given triangle mesh geometry,
whereas the interior vertices are in the opposite situation. Therefore, a smoothing term is usually
used to determine the positions of the interior vertices. Choi and Szymczak [29] proposed
fitting volumetric meshes to animated surfaces by minimizing a warped linear elasticity energy
and surface distance. Wuhrer et al. [138] used point clouds that are a subset of the surface
vertices of the volumetric mesh to infer material properties, forces and interior displacements.
Paillé et al. [94] compressed a volume grid to fit its boundary to the input surface. V olumetric
mesh tracking for 2D and 3D images is also used in biomechanics [91]. The volumetric mesh
fitting process in this thesis differs from prior work in that prior work tracked the surface of a
volumetric mesh. In contrast, the system of this thesis applies embedded simulation, a widely
useful technology in computer animation where the tracked triangle mesh is not necessarily on
the volumetric mesh surface, but may be embedded (deep) into the volumetric mesh, with the
goal of subsequent physically based animation.
2.7.1 Inverse Caging
Inverse caging is the process of finding a deformed cage (and sometimes its weights) that best
generates the given shape. The vertices of a volumetric mesh can be treated as the handles of
a cage that controls an embedded mesh. Lu et al. [74] and Savoye and Franco [99] adopted
Laplacian coordinates to deform or fit a cage. Chen et al. [28] used Green coordinates as cage
weights and transferred the deformation gradients of the input mesh to a cage. In order to
address over-determined systems occurring in the fitting problem, Thiery et al. [ 125] computed
the maximum volume submatrix of the interpolation matrix, to fit a cage animation to a triangle
mesh animation. Savoye et al. [100] converted performance mesh animation into cage-based
animation using a linear estimation framework. These methods focus on geometric caging
where the cage is typically a surface mesh, with applications in animation compression and
geometric editing. The goal of fitting in this thesis is to use a 3D solid physical energy to place
a 3D tetrahedral mesh with interior vertices, for physically based simulation. The advantage
of parameterizing the volume occupied by the input surface geometry is that it avoids interior
contact on triangle mesh deformation like collisions in concave regions (Figure 4.13). Barbiˇ c et
2.8 Self-intersection Analysis 11
al. [10] used a related method to generate volumetric mesh deformations from given triangle
meshes. However, they only used fitted deformations as animation keyframes. In contrast to
their model reduction approach which is limited to global low-frequency deformations, the
full-dimensional method applied in this thesis models local deformations, and preserves input
animation detail up to the volumetric mesh resolution.
2.8 Self-intersection Analysis
2.8.1 2D Geometry
Research on detecting and decomposing 2D self-intersecting shapes dates back to Shor and
Van Wyk [105], who gave a polynomial-time algorithm to find all possible immersions of a
2D shape which comes from stretching and bending a disk, by finding triangulations of the
input curve. Eppstein [33] analyzed the complexity of various immersion and embedding
problems and found out that it is NP-complete to determine several immersion and embedding
problems. Frisch [38] studied extending immersions of 2D circles into 2D disks and analyzed
existence and uniqueness of the problem. In computer graphics, Mukherjee [84] described
another method to solve the problem of 2D disk immersion. As acknowledged in the author’s
follow-up work [82], the original paper [84] did not provide details. Later, Mukherjee [83]
presented a way to interpolate boundary curves of disk immersions based on Shor and Van
Wyk’s triangulation. Mukherjee [82] gave a new method to solve the same disk immersion
problem by cutting the curve at “crest points”. These papers were limited to 2D. The common
self-intersections on 3D polygonal meshes in computer graphics was not studied by these
researchers.
As noted by Shor [105], a self-overlapping 2D curve can lead to more than one distinct
immersion (Figure 2.1). Both methods of Shor and Mukherjee [82] search for all possible
immersions of a given curve. Although the immersion method described in this thesis is
primarily for 3D applications, it is also directly applicable to 2D. When applied to immersions
of 2D disks in 2D, it is able to find all possible immersions, paralleling these previous results.
Different from previous approaches, it is applicable both to 2D and 3D immersions, and finds
all possible immersions on non-disks both in 2D and 3D (e.g., tori), and on multiple-component
inputs (e.g., with a hole in the shape, Figure 2.2).
2.8 Self-intersection Analysis 12
Figure 2.1 A self-overlapping curve with two distinct immersions. Left: input curve. Middle:
decompositions of the two different immersions. Right: the immersion method in this thesis is
able to run on this 2D example and find both two immersions. Green dots are handles used to
pull out the mesh with ARAP energy.
Figure 2.2 Input with multiple connected components in 2D. The input model consists of
two components, forming an egg-shaped hole. The immersion method in this thesis is able to
find the correct immersion and construct an “unglued” simulation triangle mesh, which enables
ARAP to pull the model apart.
2.8.2 3D Geometry
There has been further research on 3D immersions in computer graphics. Li [70] used Gauss
diagrams from knot theory to detect all distinct embeddings of a circle in 3D sharing the same
projection as the given self-overlapping curve. Weber and Zorin [132] introduced a method
to map a triangle mesh of disk topology to arbitrary domain with potentially self-intersecting
boundaries. Although operating in a 3D space, these methods have been limited to immersions
2.8 Self-intersection Analysis 13
onto 1D and 2D manifolds, whereas this thesis discusses volume immersions onto 3D manifolds.
Several mesh repair methods [5, 6, 24, 51] process 3D self-intersecting meshes, by gluing the
mesh using a union-like operation at the self-intersections, or forming the outer hull of the
model. In contrast, the immersion method in this thesis generates meshes that correctly separate
the self-intersecting parts. Sacht’s method [97] utilized conformalized mean-curvature flow
to resolve self-intersections of a 3D surface mesh. While their method can be used for both
2D and 3D immersions, the theory behind it is limited to disk or sphere topologies (no tori of
any number of handles), whereas the immersion method in this thesis is designed to process
meshes of any genus; we give a comparison to Sacht’s method in Section 5.6. Sacht’s work
aimed at “unwrapping” the input self-intersecting surface and did not directly aim to generate
a simulation tet mesh. Although an extension whereby a tet mesh can be generated in the
unwrapped space and then deformed forward to the world-space is mentioned, such a tet mesh
may suffer from poor quality in case the “unwrapped” mesh is flattened (see Section 5.6). In
comparison to Sacht, this thesis gives the method to generate the tet mesh in the world space
and as such the size and shape of tetrahedra is better controllable and suitable for subsequent
simulation. There is no aware of any work, in computer graphics or computational topology or
any other field, that has studied our 3D immersion problem onto volumes bounded by surfaces
of arbitrary genus. Recently, Mitchell [78] provided a method to create non-manifold level
sets for nearly self-intersecting and self-intersecting meshes. Their work addresses implicit
functions and focuses on self-collision resolution, whereas this thesis focuses on polygonal
meshes, generating an embedding into arbitrary user-provided tet meshes.
Note that several methods have utilized a self-intersecting volumetric mesh to embed an
intersection-free surface mesh with narrow features [79, 89, 109, 111, 121], or for mesh
cutting [111, 131]. The immersion method in this thesis can also operate on nearly self-
intersecting inputs, producing “un-glued” self-intersecting volumetric meshes for subsequent
embedded simulation. To do so, it follows the general spirit of Sifakis’s method [111]; but
greatly accelerate it by better managing the usage of exact arithmetic. In such applications, it
is critical to employ exact arithmetic for triangle vs tet intersections. Sifakis’s method uses
constructive solid geometry, whereas our method largely avoids exact queries of constructive
solid geometry, by demonstrating an exact arithmetic variant of the Sutherland-Hodgman tet vs
triangle clipping algorithm [116]. Tet vs triangle clipping has previously been used to re-embed
triangle meshes into dynamically subdivided tetrahedral meshes [136, 137]. The experiment
shows a 30× overall meshing and embedding speedup compared to the exact Constructive Solid
2.8 Self-intersection Analysis 14
Geometry (CSG) method of Sifakis in our examples. Teran [121] explained how to assign and
duplicate the tet vertices to connect the duplicated tets, which is re-used in this thesis.
Chapter 3
Interactive Editing With Contact
Handling
In this chapter, we present an efficient output-sensitive method for editing shapes undergoing
contact and self-contact. Geometric shape modeling is a fundamentally important task in
computer graphics and related fields. Previous methods have largely focused on how to model
shapes using various deformation energies, typically by prescribing a set of handles which the
user can employ to adjust the shape. In many applications, however, shapes undergo contact or
self-contact. There has been little work to address contact for geometric shape modeling. The
method described in this chapter will focus on this problem. The examples include corrective
shapes in pose-space deformation for character animation where limbs can penetrate bodies,
modeling shapes of trees, and editing the shape of clothing in contact with the human body. In
computer animation, this method could potentially also be used to improve blendshapes for
face animation, by resolving mouth and eye collisions.
Contact-awareness can in principle be added to any handle-based deformation method: simply
add a new handle to the system whenever a new contact is detected, constraining the handle
to some nearest collision-free position. The first issue with such a “naive” approach is that it
is suitable only for contact against external objects. For self-contact, it is not clear where the
displaced handles should be positioned. Naive placements leave gaps, cause penetrations or
produce non-smooth deformations at the collision sites. The second issue is that adding handles
is typically a very expensive operation, as it often requires re-calculating the shape functions, or
solving large systems of equations. For example, the widely-used bounded biharmonic weights
16
Figure 3.1 Example of editing the chimpanzee with multiple contact sites. In the example,
the chimpanzee’s mouth was closed, forming a mouth self-contact. One of its hands was also
pulled to collide with its face, and the other to collide with an ear. Upper-left: input triangle
mesh, lower-left: tet mesh. Upper-middle: deformed mesh with collision resolved, upper-right:
activated bases: green (L0), yellow (L1) and brown (L2). User handles are shown in blue.
Lower-middle and lower-right: zoom-in on the collision between the mouth and the hand.
Contact points are highlighted in red. The maximum computation time at any frame in the
method is 109× smaller to the maximum time in full-space modeling.
require re-solving an optimization problem on the mesh in order to add a new handle. This
is problematic in the presence of time-varying contacts where new contacts may appear or
disappear at every frame.
The method described in this chapter alleviates these issues by constructing a multi-resolution
hierarchy of shape deformation bases that is both rotationally invariant and linearly precise.
It enables smooth self-contact handling using a constraint formulation (Section 3.5.2). The
hierarchy starts with the constant and linear precision non-hierarchical subspace proposed by
Wang et al. [130], serving as the first level of the hierarchy. Then a novel multi-resolution basis
hierarchy is described that maintains constant and linear precision, and that can be used to
add more degrees of freedom locally to address contact. At runtime, smooth bases of various
localities can be efficiently added to the existing basis. Based on the current contact complexity,
17
Figure 3.2 Comparison to Harmon’s method: The hierarchical method creates more natural
and smoother shapes. The example adds a handle to the armadillo and pulls it to collide with the
ground plane. The hierarchical method (left column) then solves for smoothness and collisions
together, whereas Harmon’s method (right column) fixes collisions only as a post -process.
With the hierarchical method, the feet of the armadillo bend naturally in the initial iterations;
slowly, the entire body leans forward to minimize ARAP energy. Harmon’s method projects the
subspace control points to the closest collision-free states, moving only control points at the
feet up, creating unnatural shapes.
the runtime traversal can be interrupted at any level, making it possible to fit the calculation
into prescribed computational budgets.
The method requires a tetrahedral mesh to represent the volume inside the given polygon mesh.
Tetrahedral meshes are one kind of volumetric meshes that parameterize a given 3D volume.
The cell shapes of a tetrahedral mesh are tetrahedra, the simplest shape of all the ordinary
convex polyhedra. They are commonly used in modeling and simulation. In this thesis, this
type of volumetric meshes is chosen for its simplicity, although other volumetric meshes can be
used as well given modifications to corresponding algorithms are made. We put the details of
tetrahedral mesh creation into Section 4.2 because Chapter 4 contains related experiments.
This chapter also describes a new contact model suitable for geometric shape modeling. The
model is inspired by the active set method; but is designed for shape modeling, which is
18
Figure 3.3 Editing the Huangshan pine tree with collisions. Three handles (shown in three
yellow rectangles) are added on three branches to pull them against other branches, creating
self-collisions. The contact regions display translucent branches and red colliding triangles.
The hierarchical subspace algorithm provides shorter and stabler computation times compared
to the full method.
characterized by a lack of need for modeling dynamics, and a focus on smoothness, stability and
controllability. It utilizes the special system matrix format to accelerate computation, and also
uses multi-core computing, whereby the “updater” thread continually updates the internal data
structures to maintain fast update rates during modeling of transient contact. It applies hierarchy
and contact model to geometric modeling of static shapes with the widely used as-rigid-as-
possible (ARAP) energy [114]. This chapter presents several challenging examples of complex
geometry in self-contact (Figures 3.1,3.3) and contact with external objects (Figure 3.15). They
demonstrate that the hierarchical method outperforms full non-hierarchical shape modeling by
at least an order of magnitude. An experimental comparison is also made between contact model
to the state-of-the art for contact-aware shape modeling [49], and it demonstrates significant
improvements in shape quality (Figure 3.2).
3.1 Shape Modeling 19
3.1 Shape Modeling
The rest of the chapter first describes shape modeling in the full space (without reduction),
subject to arbitrary user and contact constraints. The widely-used as-rigid-as-possible (ARAP)
energy [114] is adopted to measure deformation smoothness. Although ARAP is the main
energy used in this work, it is possible to have other similar energies implemented. One of them
tested is the volumetric strain energy used in Projective Dynamics [20].
3.2 Modeling in Full Space
Let p∈R
n×3
contain the deformed mesh vertex positions. The user controls the shape by
prescribing the positions of one or more user handles placed at arbitrary mesh vertices,Sp=
s, where S∈R
h×n
is the selection matrix that selects h handle indices in p, and s∈R
h×3
are the handle positions controlled by the user. At each iteration, it performs the ARAP
local/global optimization [114]. The local step finds the best rotation around each vertex, via
polar decomposition. The global step solves a linear system
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
L S
T
C
T
S
C
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
p
λ
S
λ
C
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
b
s
c
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
, where (3.1)
L=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
L
L
L
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
, S=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
S
S
S
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
. (3.2)
Matrix L∈R
n×n
is the cotangent Laplacian matrix [114], and can be pre-factored. The right-
hand side matrix b∈R
3n
depends on local rotations R
i
(Appendix A). Vector p∈R
3n
stacks
the vertex positions at each dimension sequentially, p=[p
T
[∶0]
p
T
[∶1]
p
T
[∶2]
]
T
, wherep
[∶i]
represents
the i-th column of p. Similar conventions are followed for b and s. Throughout the chapter,
bold-face such as p∈R
3n
is used to denote a vector containing mesh vertex displacements (first
x-DOFs for all vertices, then y-DOFs, then z-DOFs), as well as matrices that operate on such
vectors. The sans-serif typeface such asp∈R
n×3
is used to denote arranging the displacements
or positions into a n×3 matrix. To resolve contact, planar contact constraints Cp= c are added,
where additional requirements are enforced so that the constraints are not pulling but only
3.3 Modeling in a Shape Subspace 20
pushing. This will be discussed in Section 3.5. The Lagrange multipliers are denoted byλ
S
and
λ
C
.
The local/global optimization is effectively a block-coordinate descent, which in practice
decreases the energy at each iteration. One benefit of the ARAP energy is that the x,y,z
dimensions are decoupled, enabling us to solve three pre-factored n× n systems in parallel.
However, since contact normals can point in arbitrary directions, and sliding constraints are
adopted for better contact resolution, the contact constraints couple all three dimensions. The
user can freely move the user handles while the algorithm deforms the rest of the mesh to
minimize the ARAP energy and resolve contact (Figure 3.3). Note that an alternative approach
would be to treat the positional constraints as Dirichlet boundary conditions, at the cost of
re-factoring the system each time a new constraint is added. Therefore, Lagrange multipliers are
chosen because they support non-positional (tangential sliding) constraints needed for contact,
and because they enable a faster system update method (Section 3.6).
3.3 Modeling in a Shape Subspace
The full optimization method given in Section 3.2 does not run at interactive rates for detailed
geometry. To make the modeling tractable, the hierarchical method in this chapter solves the
system in a hierarchical shape subspace. Other advantages of subspace methods are faster
convergence compared to a full method, and a decrease in the spiky artifacts around manipulated
handles. Suppose the mesh vertex positions can be expressed asp=Uv, whereU∈R
n×r
is the
shape basis, andv∈R
r×3
is the reduced coordinate. Equation 3.1 reduces to
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
˜
L
˜
S
T ˜
C
T
˜
S
˜
C
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
v
λ
S
λ
C
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
˜
b
s
c
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
, where (3.3)
˜
L=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
˜
L
˜
L
˜
L
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
,
˜
S=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
˜
S
˜
S
˜
S
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
, (3.4)
3.4 Multi-Resolution Shape Deformation Basis 21
for
˜
C= CU. Here, v=[v
T
[∶0]
v
T
[∶1]
v
T
[∶2]
]
T
, and similarly for
˜
b. The symbol ˜ x denotes the reduced
counterpart of x. The r×r subspace system matrix
˜
L=U
T
LU (3.5)
is much smaller to solve than L, assuming r≪ n. The subspace selection matrix is
˜
S=SU,
and
˜
b=U
T
b. Reduced models are known to be difficult for detailed collision handling. It
is addressed by dynamically inserting more degrees of freedom using the hierarchical bases
(Section 3.4).
Accelerating right-hand side computation: The time complexity of computing
˜
b is O(nr).
One can approximate
˜
b efficiently by exploiting the fact that the subspace vertex positions are
highly correlated. The hierarchical method in this chapter evaluated rotation clustering [56]
and cubature-like [4] methods, both properly adapted to the setting (details are in Appendix A).
Both methods provide fast right-hand side evaluation and are found to have similar runtime
performance and accuracy.
3.4 Multi-Resolution Shape Deformation Basis
This section describes how the hierarchical bases are created to allow a compromise between
collision detail and speed (Figure 3.4).
3.4.1 Single-Level Basis
A good shape deformation basis should be smooth and shape-aware. The deformation quality
and the user experience is greatly improved if the basis achieves both constant and linear
precision [130], i.e., if
U1
r
= 1
n
, (3.6)
¯ p=U¯ v, (3.7)
where 1
k
∈R
k
is a vector of ones, ¯ p∈R
n×3
contains the rest positions of mesh vertices, and ¯ v
contains the rest positions of manipulated vertices. Any affine transformation can be achieved
3.4 Multi-Resolution Shape Deformation Basis 22
Figure 3.4 Hierarchical bases provide a tradeoff between collision details and speed. The
top part (mustard color) of the mesh (top-left sub-figure) is pulled to cause a self-contact
against the bottom (scarlet color), as shown in the top-right sub-figure. With level-1 local
bases activated (top-right and middle-left in closer view), collisions are resolved with visible
penetrations and no deformation of the bottom part, at 1× computation time. When level-2 local
bases are added (middle-right), the method is able to express the deformation of the bottom part,
at the cost of 1.19× of the computation time. When level-3 local bases are added (bottom-left),
the penetrations are reduced at the cost of 1.34× time. When level-4 local bases are added
(bottom-right), the penetrations are minimal. With many modes activated, the method slows
down to 3.10× time.
once Equations 3.6 and 3.7 are satisfied,
pQ+1
n
t=(Uv)Q+(U1
r
)t=U(vQ+1
r
t). (3.8)
Given an affine transformation ( Q∈R
3×3
,t ∈R
1×3
), the transformation of the positionsp can
therefore be recovered by transforming the control points v. The basis generation method
3.4 Multi-Resolution Shape Deformation Basis 23
presented by Wang et al. [130] satisfies Equations 3.6 and 3.7, and therefore it is chosen as the
first level of the hierarchy.
For completeness, the computation of Wang’s basis is described below. Given a solid 3-manifold
¯
Ω∈R
3
made of 3D tetrahedra, Wang divides the tet vertices into two disjoint sets: “manipulators”
and “free” vertices. For each manipulator vertex, Wang’s method computes one basis vector
that describes how the free vertices move if one moves the manipulator vertex while keeping the
other manipulator vertices fixed. The word “manipulator” is used here to avoid a confusion with
the word “handles” which in this chapter denotes the user handles to control the shape. Given
r≥ 1 manipulator vertices with prescribed positions stacked in a matrixy∈R
r×3
, a smoothness
energy is minimized,
p= argmin
x∈R
n×3
1
2
trace(x
T
Ax) subject to Yx=y, (3.9)
wherep∈R
n×3
are the minimizing deformed mesh positions. The matrixA=K
T
M
−1
K∈R
n×n
gives a positive semi-definite quadratic form measuring the fairness of x, whereK∈R
n×n
is a
linearly precise cotangent Laplacian matrix, and M is the mass matrix. The selection matrix
Y∈R
r×n
selects the r manipulator positions. The solution of Equation 3.9 is
p=Wy= ¯ p+W(y− ¯ y), for (3.10)
W=Y
T
−
ˆ
Y
T
(
ˆ
YA
ˆ
Y
T
)
−1
ˆ
YAY
T
, (3.11)
where
ˆ
Y∈R
(n−r) ×n
is the selection matrix that selects the free degrees of freedom, and ¯ y are
the manipulator positions at rest. MatrixW∈R
n×r
gives the Wang’s basis. Each basis vector
is smooth and global, and assumes a value of 1 at exactly one manipulator, and 0 at all other
manipulators. The largest deformation magnitudes occur close to the manipulator that assumes
the value of 1. The basis vectors can be computed using Equation 3.11, by solving r linear
systems with the same sparse matrix
ˆ
YA
ˆ
Y
T
. Factoring the system matrix and solving for the
columns ofW can be done in parallel.
Wang’s method allows arbitrary sets of vertices to serve as manipulators. At first, the following
direct and “naive” approach was tried, and eventually abandoned in favor of the methods
presented in the rest of the chapter. The idea was to produce localized bases by adjusting the
number and distribution of Wang’s manipulators at runtime, based on the occurring collisions,
and recomputingW. Note that this “localization” is not exact, as all basis vectors are globally
3.4 Multi-Resolution Shape Deformation Basis 24
Figure 3.5 Recomputing the Wang’s bases during collisions at runtime is extremely slow,
compared to the precomputed hierarchical bases. Octopus model with 20,165 vertices and
92,003 tetrahedra.
supported across the entire mesh, albeit decaying to zero far away from the manipulators.
Therefore, thresholding is required to enforce any local support. Furthermore, the recomputation
of W is a major drawback because the basis functions need to be recomputed each time a
manipulator is added or removed. This is a global solve of size n, regardless of the intended
detail and basis content. It is slow and not practical to perform at runtime (Figure 3.5).
3.4.2 Hierarchical Basis
The hierarchical basis addresses these issues by pre-computing local bases and combining them
with the global basis at runtime, as visualized in Figure 3.1 and 3.6. At runtime, the global basis
is always activated. When collisions or other events require more detailed deformations in a
region, different levels of local bases of the region are activated, thereby adding more degrees
of freedom into the solver. The basis hierarchy was built using a hierarchy of manipulators.
The hierarchy structure is a forest, where a group of globally distributed manipulators serves
as the roots of the forest. In this hierarchy, higher manipulators correspond to more global
modes, covering deformations on larger areas, while lower manipulators correspond to more
local modes, covering deformations on smaller areas. The modes are grouped into global and
local bases, and so are the manipulators.
Global Basis
The global manipulator group is created using low-discrepancy sampling. The algorithm starts
at a random vertex, and then always picks the vertex that is the furthest from the already selected
3.4 Multi-Resolution Shape Deformation Basis 25
Figure 3.6 Examples of basis vectors at different levels of the hierarchy (second row).
Highest values are colored red and lowest are purple. Basis vectors in higher levels (e.g., global
level L0) have larger support while basis vectors in lower levels (e.g., L2) have smaller support.
Activated basis vectors due to self-collision are displayed as colored points in the first row. They
are green (L0), yellow (L1) and brown (L2). The user handles are colored blue.
manipulators and adds it into the selected manipulators, measured in semi-geodesic distances,
until the number of selected manipulators reaches a pre-defined value of B
0
. To avoid excessive
precomputation times on costly geodesic computations, The distances are measured by running
Dijkstra’s algorithm on a graph where each mesh vertex is a node, and two nodes are connected
by an edge if the two corresponding vertices share a common tetrahedra; hence “semi-geodesic”
distances as opposed to geodesic distances. The length of an edge is the Euclidean distance
between the two vertices. In this way, B
0
tree roots for the manipulator forest are computed first.
In Figures 3.8 and 3.9, red points are part of the root manipulators. This group of manipulators
are referred as the global group, or level-0 group. Its corresponding basis is computed using
Wang’s method on the entire mesh (Equation 3.11). The basis is called the global basis, or
level-0 basis. Then, the method divides the mesh into B
0
disjoint sets. These sets are the
V oronoi regions of the level-0 manipulators with respect to the semi-geodesic metric described
above. Note that the V oronoi regions are also used in meshless deformation models [37].
3.4 Multi-Resolution Shape Deformation Basis 26
Figure 3.7 Comparison between a different number of global manipulators: Fewer global
manipulators give un-natural motion when deformed only by the global basis. In the left image,
the chimpanzee’s left arm is not bent properly, compared to the right image where a larger
number of global manipulators smoothly deform the arm at the elbow.
Wang’s Bases With Exact Locality
Before the local basis are defined, it needs to first give a method to define and compute localized
Wang basis vectors, with exact localization and without any thresholding. In addition to free
and manipulated vertices, a new type of vertices is introduced, fixed vertices: these are vertices
where one wants Wang’s basis to be identically zero, such as, for example, outside of some
localized region in a hierarchy. The construction below works for arbitrary sets of manipulators,
free and fixed vertices. The idea is to make the fixed vertices into “manipulators” whose
positions are permanently set to their rest positions. Let Z be the selection matrix that selects
the fixed vertices. Then, the constraint in Equation 3.9 becomes
⎡
⎢
⎢
⎢
⎢
⎣
Y
Z
⎤
⎥
⎥
⎥
⎥
⎦
x=
⎡
⎢
⎢
⎢
⎢
⎣
y
Z¯ p
⎤
⎥
⎥
⎥
⎥
⎦
. (3.12)
Equation 3.11 now becomes
p= ¯ p+[W X]
⎡
⎢
⎢
⎢
⎢
⎣
y− ¯ y
0
⎤
⎥
⎥
⎥
⎥
⎦
= ¯ p+W(y− ¯ y), where (3.13)
[W X]=(I−
ˆ
Y
T
(
ˆ
YA
ˆ
Y
T
)
−1
ˆ
YA)[Y
T
Z
T
], (3.14)
3.4 Multi-Resolution Shape Deformation Basis 27
and therefore the formula for W is the same as the one given in Equation 3.11. Matrix W is
identically zero on rows corresponding to the fixed vertices, by construction. The submatrix of
the manipulator rows is the identity matrix.
Local Basis
The forest is grown by expanding its current leaf nodes, one at a time. This process converts
nodes that were previously leaves into interior nodes, and generates new leaf nodes. Suppose
nodeν
i
is currently a leaf, located at level i. Recall that given any set of tet mesh verticesV
and a subset ofV called sites, the Voronoi diagram is a decomposition ofV into disjoint sets
whereby all elements of a set share the same site as the closest site (in our pseudo-geodesic
metric). Now, define a V oronoi diagram whose set V itself is a V oronoi region determined
earlier in the recursion, namely the V oronoi region of the parent node of i. Let the sites be all
children of the parent of i (includingν
i
). The V oronoi region of siteν
i
is denoted byV
i
. The
method creates the children of ν
i
(level-(i+ 1) manipulators; the new leaves) by performing
low-discrepancy sampling inV
i
. The low-discrepancy sampling picks the vertex that is the
furthest from the already selected manipulators, and vertices outside ofV
i
. The method also
skips vertices already selected as manipulators to avoid manipulator redundancy. The result of
the sampling are B
i+1
new manipulators on level-(i+1) that are the children ofν
i
in the forest
hierarchy.
The nodes are expanded until all tet mesh vertices are assigned, or a maximum level is reached.
Deeper-level nodes have smaller support, and therefore more localized bases. In Figures 3.8
and 3.9, green points are manipulators at level-1, and are children of root-level manipulators. A
subset of the forest is shown in Figure 3.9. In practice, B
0
is larger than B
i
, for i≥ 1. On level
0, it is typically 30≤ B
0
≤ 50 for small meshes, and 100≤ B
0
≤ 200 for complex meshes. For
other levels, it is 4≤ B
i
≤ 10. A larger level-0 basis produces meaningful global deformation,
whereas smaller deeper-level bases give local deformation (Figure 3.7), without adding too
many degrees of freedom to the solver at runtime.
To create a level i+ 1 local basis for a local group which are the children of ν
i
, the initial
approach was to treat the vertices outside of the V oronoi region V
i
ofν
i
as fixed vertices and
build Wang’s basis on the children manipulators, as explained in Section 3.4.2. ν
i
and all its
ancestors are also included into the fixed vertex set, so that all manipulator bases on level i+1
vanish at the nodes on level i and higher. Because each Wang’s basis vector vanishes at the
3.4 Multi-Resolution Shape Deformation Basis 28
Figure 3.8 Building manipulators and Voronoi regions: This shows the process of forming
local regions to build local bases. From left to right, up to bottom: sample level-0 (L0)
manipulators (two in this case; orange); create V oronoi regions (yellow and pink) for each
manipulator; sample level-1 manipulators as children of level-0 manipulators (green); extend
each V oronoi region into overlapping local regions (cross-hatched). The bottom-right image
shows the approach of using V oronoi regions directly as local regions. This approach produces
artifacts (Figure 3.9) and is therefore abandoned.
fixed vertices (Equation 3.11), the local basis has zero values on vertices outside of V
i
. Note
that no thresholding is required; basis vectors are automatically zero, by construction, outside
ofV
i
. However, this basis is problematic because the basis vectors of children of ν
i
do not
overlap with their “cousins”, i.e., other bases from the groups at the same level. This means
that level-(i+1) bases cannot represent non-zero deformations that cross the boundary ofV
i
.
This issue is addressed by extending the V oronoi region into a larger local region that overlaps
with other nearby local regions at the same level (Figure 3.8), and then using the construction
from Section 3.4.2 with respect to this enlarged region. For each mesh vertex v inV
i
, we
compute the largest ball centered at v that does not include any manipulators from levels 0 to
i+ 1 located outside ofV
i
(Figure 3.9, top-left). The union of these balls is called the “local
region” (denoted in yellow in Figure 3.9, top-left). Then treat the vertices outside of the local
region as fixed vertices to compute the local basis for the children of ν
i
. The effect of this is
that the basis functions on the children ofν
i
extend to the nearby V oronoi regions (Figure 3.9,
top-left). In this way, the desired overlap is achieved while still keeping
˜
L sparse. Note that
3.4 Multi-Resolution Shape Deformation Basis 29
Figure 3.9 Overlapping local region construction details, manipulator forest and artifacts
of non-overlapping bases: Top-left: manipulators on part of the bunny mesh (top-right). There
are four level-0 (L0) manipulators (orange) and many L1 ones (green). The local region
(yellow) of the top-left L0 manipulator is constructed by unioning balls centered on each vertex
in the V oronoi region of the L0 manipulator. V oronoi regions of different L0 manipulators are
separated by the blue lines. Dashed orange lines represent the actual local region boundary on
mesh vertices. Vertices inside the local region are thus covered by the computed local basis
of the L0 manipulator’s children manipulator group. Note that this figure only shows a part of
the bunny mesh. Bottom-left: the manipulator forest corresponding to the manipulators in the
top-left. Right: artifacts occur if the local regions do not overlap. In this case, local regions are
just the V oronoi regions. The bunny’s ear was pulled to collide with its back. Only L1 bases
were activated to resolve collisions in order to demonstrate the artifacts of non-overlapping
L1 bases. Observe that the non-overlapping L1 bases cannot represent non-zero deformations
across V oronoi boundaries, as highlighted in the red rectangle.
the basis functions are still automatically zero, by construction, outside of the local region; no
thresholding is required.
3.4 Multi-Resolution Shape Deformation Basis 30
3.4.3 Assembling Bases
Finally, this section describes how the bases functions affect the final position of each tet mesh
vertex. For now briefly assume that only one local basis W
ℓ
is activated in addition to the global
basisW
g
; the hierarchical case can be easily generalized. The vertex positions are computed as
p=W
g
v
g
+W
ℓ
v
ℓ
=[W
g
W
ℓ
]
⎡
⎢
⎢
⎢
⎢
⎣
v
g
v
ℓ
⎤
⎥
⎥
⎥
⎥
⎦
. (3.15)
Here, [W
g
W
ℓ
] is the “working basis”. In Wang’s work, v represents both the manipulator
positions and the reduced coordinates of the system. However, in the basis hierarchy,v
ℓ
should
not be interpreted as the manipulator positions in the local region, because the global and local
bases overlap in each local region. Instead, a crucial insight that makes the entire construction
possible is to treat v
ℓ
as the displacements from positionsW
g
v
g
. This is demonstrated by the
following equation. The positions of manipulatorsp
ℓ
in each local region are
p
ℓ
=S
ℓ
(W
g
v
g
+W
ℓ
v
ℓ
)=S
ℓ
W
g
v
g
+v
ℓ
, (3.16)
where S
ℓ
is a selection matrix that selects verticesℓ from all positions. Note that S
ℓ
W
ℓ
is an
identity matrix because, by Wang’s basis construction, each basis function inW
ℓ
has a weight
of 1 on its own manipulator, and 0 on other manipulators in each local region.
Linear precision is still preserved in Equation 3.15. Any 3×3 linear transformation Q can be
produced in the subspace by transforming all the control points,
pQ=(W
g
v
g
+W
ℓ
v
ℓ
)Q=W
g
(v
g
Q)+W
ℓ
(v
ℓ
Q). (3.17)
Likewise, constant precision holds because any 1×3 translation vector t can be produced by
translating only the global control points,
p+1
n
t=(W
g
v
g
+1
n
t)+W
ℓ
v
ℓ
=W
g
(v
g
+1
r
t)+W
ℓ
v
ℓ
. (3.18)
Note that the above two equations hold becauseW
g
andW
ℓ
both satisfy Equations 3.6 and 3.7.
The hierarchical method used Wang’s basis due to its good shape modeling properties [130].
However, note that the multi-resolution basis construction could be employed with other shape
deformation basis construction methods. For example, one could use Bounded Biharmonic
3.4 Multi-Resolution Shape Deformation Basis 31
Weights [57] or Green’s functions computed using BEM [60]. In addition, if the basis weights
are linearly precise, so will be the multi-resolution bases.
3.4.4 Speed and Memory Considerations
Before modeling, the basis hierarchy can be pre-computed quickly using multi-threading.
The local bases are sparse and can be stored efficiently and loaded into memory at startup.
At runtime, when more degrees of freedom are needed, global and local bases are simply
concatenated into the working basis. In code, this can be done very efficiently simply by
assembling pointers to the modes. This provides good flexibility for contact handling. One can,
for example, trade collision resolution for speed by setting the maximum permitted activated
level at runtime. When changing the working basis, one needs to re-project the system matrix
using Equation 3.5. Since the entire basis hierarchy is pre-computed, it can also pre-compute
the projection ofL into any basis vector in the hierarchy, obtaining an n× n matrix
˜
L
all
. This
makes it possible to rapidly projectL into any working basis, simply by selecting proper rows
and columns of
˜
L
all
. Note that
˜
L
all
is block-sparse, since the basis vectors in different trees
overlap minimally. Alternatively, if memory is scarce and there are not many added basis
vectors, one can update
˜
L by exploiting
˜
L=
⎡
⎢
⎢
⎢
⎢
⎣
U
0
T
LU
0
U
0
T
LU
1
U
1
T
LU
0
U
1
T
LU
1
⎤
⎥
⎥
⎥
⎥
⎦
, (3.19)
where
˜
L is re-built using the previous basisU
0
and the newly added basis vectorsU
1
. One only
needs to computeU
1
T
LU
0
andU
1
T
LU
1
. Since
˜
L is symmetric, only its upper-triangular part is
computed.
3.4.5 Basis Update Oracle
At runtime, when new contacts appear, the hierarchical system adds local bases into the working
basis. Specifically, for a colliding vertex v on level-i, it adds the basis vectors of v and its sibling
nodes (the local Wang’s basis of the parent of v). it then repeats this process for the parent of v,
by adding the local Wang’s basis of the parent of parent of v. It continues this cascading process
until the top of the hierarchy is reached. The user can adjust the maximum permitted level of
3.4 Multi-Resolution Shape Deformation Basis 32
local bases, trading solver efficiency for collision resolution. For any node ν, it can remove the
basis vectors of its supporting region if all the descendant nodes have been contact-free for a
given number of frames.
3.4.6 Uniqueness of Solution
The local basis construction described in Section 3.4.2 ensures that no manipulator is assigned to
two or more basis vectors. It can be proven that any basis selected at runtime is always linearly
independent, and therefore the computed solution is unique. To prove linear independence, it is
sufficient to prove that the entire hierarchical basis is linearly independent, as any subset of a
linearly independent set is linearly independent.
Lemma: For each nodeν in the hierarchy, the only basis vectors that do not vanish atν are its
basis vector, and the basis vectors of nodes whose hierarchical level is closer to the root thanν.
Proof: Observe that the basis functions of siblings of ν and descendants of ν vanish at ν
by construction. For the other nodes at equal or larger hierarchical level than ν, revisit the
construction of overlapping regions in Section 3.4.2 (yellow sphere union in Figure 3.9): spheres
are grown until reaching nodes in sibling V oronoi diagram regions, and therefore basis functions
vanish at those nodes. ∎
Theorem: The entire hierarchical basis is linearly independent.
Proof: Suppose that a linear combination of the basis vectors equals zero at all mesh vertices.
For each node on level 0, the only basis vector that is non-zero is that corresponding to the node
itself; hence the linear combination coefficient of all nodes on level 0 must be zero. For a node
on level 1, by the Lemma, the only basis vector that is non-zero is that corresponding to the
node itself and the nodes on level 0, which have zero coefficients. Hence, the coefficient for
each node on level 1 must also be zero. By repeating this process, the coefficients of all basis
vectors in the entire hierarchy must be zero, i.e., basis vectors are linearly independent. ∎
3.5 Collision Detection and Handling 33
3.5 Collision Detection and Handling
The hierarchical system supports collisions with external static objects, and self-collisions. The
collision model uses signed distance fields to model external static objects and hashing tables
for self-collision detection. Collision constraints are then formed and added into the solver.
3.5.1 Collision Detection
The collision points can be any set of points sampling the surface of the object. In the examples
shown in this chapter, the set of mesh vertices is used; but one could employ a denser set by
sampling, say, more points on the surface triangles. For external collisions, the method queries
the signed distance value of each vertex to check whether it is in contact. For self-collisions,
the spatial hashing [124] is used, where the entire space is partitioned into uniform grid cells.
Each cell stores the deformed mesh vertices it contains, as well as the deformed tetrahedra
that overlap with it. This hash table is updated at every frame. For each cell, the method
checks the containing points against the containing tetrahedra. If it detects a point x inside a
tetrahedron that is not the tetrahedron containing this point in the rest configuration, then this
is a self-collision that will be resolved as follows. The method first computes the barycentric
coordinate of x in the tet, and then “pulls-back” x to its position ˆ x in the rest configuration, using
the tetrahedron’s inverse affine transformation. It then queries the rest shape signed distance
field to find a surface point ˆ x
s
that is closest to ˆ x [122]. Then transform ˆ x
s
forward using the
affine transformation of, and barycentric coordinates in, its containing tetrahedron, obtaining x
s
.
Improvement to Teran’s Method: In Teran’s method [122], the distance and direction from
x to x
s
are the penetration depth and normal. This pull-back method has a drawback where if the
collided area has a large shear, the x
s
may not be the closest point to x on the deformed surface.
To alleviate this issue, one trick can be used here: doing a local search on the surface triangles
around x
s
to find another point x
c
that is closer to x than x
s
. This local search is done by starting
at x
s
, it first computes the closest point x
t
to x on the surface triangle t where x
s
lies. If x
t
is the
same as x
s
, the local search ends. Otherwise, x
t
is closer, then the method starts the local search
around x
t
. The method also keeps track of visited triangles to avoid potential infinite looping.
The collision scheme is easy to implement and is embarrassingly parallel. Unless the hash table
degenerates (it did not in the examples), the query time complexity is O(n) when collision
3.5 Collision Detection and Handling 34
points are mesh vertices. Signed distance field queries, and the hash table traversal to read the
candidate (point, tetrahedron) pairs can be performed in parallel. Typically, collision detection
takes most of the time in the system. To further reduce collision detection costs, one can limit
the collision detection to only a subset of the tetrahedra at each frame, in a round-robin manner
to ensure all tetrahedra are checked in a fixed number of frames. In addition, the colliding
tetrahedra and nearly colliding tetrahedra (defined by a fixed maximum topological distance to
colliding tets) are always checked at each frame, regardless of which subsets they belong to.
3.5.2 Contact Constraint Formulation
The contacts in the hierarchical method are modeled as sliding constraints in the plane of
contact, as opposed to, say fixing a vertex to the closest surface point. The sliding flexibility
is tested to be essential for good shape modeling. Given an external collision on vertex i with
position x
i
= S
i
p∈R
3
and collision normal n
i
, it creates a constraint
n
T
i
(x
i
−x
s
i
)= 0, (3.20)
where x
s
i
is the closest point to x
i
on the surface of the external object. The selection matrix
S
i
∈R
1×n
selects the colliding vertex. Each contact contributes a single row to the collision
constraint matrix C and right-hand side c in Equation 3.1, as follows
C
single row
=[n
0i
S
i
n
1i
S
i
n
2i
S
i
], c
single row
= n
T
i
x
s
i
, (3.21)
where n
i
=[n
0i
,n
1i
,n
2i
]
T
. Given a self-collision on point x
i
, it creates a constraint
n
T
i
(x
i
−
2
∑
j=0
w
j
x
t
j
)= 0, (3.22)
where x
t
j
are vertex positions of the triangle that embeds x
s
i
, and w
j
are the barycentric coordi-
nates of x
s
i
. For self-collisions, the constraint rows of the system, and the right-hand side, are
C
single row
=[n
0i
S
t
n
1i
S
t
n
2i
S
t
], c
single row
= 0, (3.23)
whereS
t
=S
i
−∑
2
j=0
w
j
S
t
j
, andS
t
j
selects the triangle’s j-th vertex.
3.5 Collision Detection and Handling 35
After collision detection, all constraints are assembled and C and c are formed to add to
Equation 3.3. To avoid over-constraining the subspace, only constrained vertices are allowed
to be the currently active tree hierarchy nodes. This creates at most r constraints (r is current
basis size), while the entire system has 3r degrees of freedom. To resolve detailed collisions
that the current control points cannot represent, new local control basis vectors are added into
the system at runtime (Section 3.4). In addition, if the user moves the user handle into contact,
the system moves it out of contact to prevent further collisions.
3.5.3 Updating the Contact Constraints
This section describes how unilateral (i.e., non-sticking) contact is modeled. The approach
is inspired by the active set method [90], adapted here to shape collision editing, where there
is no dynamics, but shape quality is important. The approach ensures shape smoothness. At
the end of each solve, the method determines a set of contactsC that must persist to the next
iteration. This is done by checking the sign of the Lagrange multipliersλ
C
of all constraints. A
positive Lagrange multiplier means that there is a “sticking force” on the collision surface, in
which case the contact is not added toC. All other contacts are added toC. At the next iteration,
collision detection is performed again, forming new contacts. Contacts that are inC and that
also appear as new contacts, are simply treated as new contacts. Sometimes contact sites drift
on the surface after each iteration; this tends to happen in convex regions of the mesh under
contact. To encourage smooth sliding over convex surfaces, The contact procedure is modified
slightly, by keeping contact constraints on vertices that are no longer in contact but that have a
negative Lagrange multiplier. While this adds a small amount of contact stickiness (negligible
in the examples), it is found to be very beneficial for sliding on convex surfaces. After all
collision sites are updated, the linear system solve is performed and a new setC is generated
based on the new Lagrange multipliers.
3.5.4 Static Damping
In complex contact configurations, the solution can be prone to contact limit cycles, where
the shape does not converge, but jumps between two or more contact configurations, creating
unpleasant visual high-frequency vibrations on the mesh (Figure 3.10, left). This is addressed
by static damping. Initially, a line search was experimented to solve this issue; but it was
3.6 Solving the Linear System 36
Figure 3.10 Static damping alleviates contact limit cycles: Here, the shape of the hand of a
character was edited. Left: Six frames performed without static damping are superimposed with
transparency to show the high-frequency vibration on the thumb, highlighted in a red rectangle.
The thumb looks blurred because it is the superposition of several frames of a vibrating motion.
Right: Six frames, with static damping enabled (α = 0.5). Same (recorded) handle motion as in
the Left image. The thumb does not vibrate and the shape converges, therefore the superposition
is sharp.
found that the line search sometimes makes very slow progress and is not compatible with
interactivity. The idea of static damping is that once the target positionsp
t
is obtained while
under contact, the system does not attempt to update the mesh to that position, but instead
to (1−α)p
t
+αp
c
, where p
c
is the current vertex positions. The parameter 0≤α ≤ 1 is the
static damping parameter; forα = 0, it is the undamped solution and forα = 1 it completely
discards new shapes, stopping the shape update. The user can setα = 0.5 or other values after
the manipulation of the handles is over. it can also automatically start static damping each time
the user stops moving the user handle, and gradually increasesα from 0 to the desired value,
letting the contact “settle”. Figure 3.10 demonstrates how static damping can help resolving the
cycles.
3.6 Solving the Linear System
In Equation 3.1, the collision constraint (C,c) changes frequently during the user interaction.
If no collision exists, the system matrix in Equation 3.3 can be pre-factored. Unfortunately,
pre-factoring does not work well under contact. Each update of C requires re-factoring. One
3.6 Solving the Linear System 37
Figure 3.11 Comparison between using vs not using the additional thread to computeF
i
. Klara
cloth example.
potential way is to use an iterative solver such as the preconditioned conjugate gradient (PCG)
to avoid factoring. Although PCG only works on positive-definite systems, a filtered version
given by Tamstorf et al. [118] could work on positive-definite systems with linear constraints.
However, on the examples, it is found that iterative solving is much slower than using a direct
method. When lots of local bases are added inU to handle collisions,U becomes sparse and so
doesL. Therefore, the direct sparse solver in Eigen [45] is used in the end to solve Equation 3.1.
In order to maintain the speed of the direct method without the need to re-factor, The technique
of incremental solving is adopted. Suppose a symmetric matrix A is factored, then solving
⎡
⎢
⎢
⎢
⎢
⎣
A B
T
B
⎤
⎥
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎢
⎣
x
λ
⎤
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎣
b
d
⎤
⎥
⎥
⎥
⎥
⎦
(3.24)
can be performed as
⎡
⎢
⎢
⎢
⎢
⎣
x
λ
⎤
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎣
A
−1
−FDF
T
FD
DF
T
−D
⎤
⎥
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎢
⎣
b
d
⎤
⎥
⎥
⎥
⎥
⎦
, (3.25)
where D= S
−1
,S= BA
−1
B
T
,F = A
−1
B
T
. Such incremental solving is standard (see, e.g., Barbic
et al. [13] and Yeung et al. [140]). The dense solve on S is performed using LAPACK routines in
Intel MKL. Given n rows in A and k rows in B, the incremental solve requires k+1 pre-factored
n× n solves, and one k× k dense solve. For a large n, this is faster than to factor the sparse
(n+k)×(n+k) system at every iteration. However, the performance of the incremental solving
decreases as more rows are added to B.
The rest of the section describes the novel accelerations to speed up incremental solving. First
note that the fixed structure of A and B can be exploited. In this system, A=
˜
L,B=[
˜
S
T ˜
C
T
]
T
.
3.7 Results 38
Since
˜
L is a block-diagonal matrix, it performs
˜
L
−1 ˜
b by computing
˜
L
−1 ˜
b
[∶i]
,i= 0,1,2, and
stacking the results. The matrix
˜
S comes from user handle constraints and they are less likely
to change frequently. It stores
˜
F=
˜
L
−1 ˜
S
T
, and reuses it at subsequent iterations. Computing
˜
F can be accelerated by observing that
˜
S is also a block-diagonal matrix, consisting of SU.
It only needs to solve
˜
L
−1
(SU)
T
once, and stacks the results to build
˜
F. The matrix
˜
C equals
CU, and C has a structure indicated in Equations 3.21 and 3.23. For an external collision on
vertex i, it computes and storesF
i
=
˜
L
−1
(S
i
U)
T
. Then, solving for the constraint is simplified:
˜
L
−1 ˜
C
T
i
=[n
0
F
T
i
n
1
F
T
i
n
2
F
T
i
]
T
. Similarly, for a self-collision between vertex i and triangle
t, it computes and storesF
i
,F
t
0
,F
t
1
, andF
t
2
. Overall, givenU, it stores the result ofF
i
so that if
there is a contact occurring at vertex i in subsequent iterations, it can be reused. During contact,
one dedicated thread is used to compute F
i
, while the main thread performs the incremental
solve. Performance is analyzed in Figure 3.11. If the user slides the mesh along the contact
surface, newF
i
are added and the old entries can be discarded if memory is scarce.
3.7 Results
The subspace hierarchical method and the full method are tested on several examples. All
experiments were conducted on a workstation with Intel Xeon 2.3GHz 2×6 cores and 32
GB memory. Statistics for the examples are given in Table 3.1 and 3.2 and 3.3. For both
methods, the slowest frames occur during collisions. Because the application is interactive
shape modeling, the results are primarily focused on the maximum slowdown between the
full method and hierarchical method, as these slow-downs impact interactivity the most. The
average slowdowns are also reported. It can be seen that the hierarchical method is an order of
magnitude faster than the full method, both in largest and maximum slowdowns.
Example vtx tets avg cons min/max basis size
Armadillo 9,746 44,371 5.7 43 / 133
Klara (cloth) 11,718 40,399 34 100 / 885
Bunny 71,724 215,172 49 50 / 555
Chimpanzee 84,613 334,222 47 50 / 870
Huangshan pine 111,364 348,467 37 200 / 1049
Table 3.1 Example Complexity for #mesh vertices and #tetrahedra, #average collision con-
straints during contact, min/max basis size used during contact.
3.7 Results 39
Example at:pWv at:rhs at:col at:solve at:total at:sp mt:total mt:sp
Armadillo 13 % 8.4 % 58 % 6.1 % 3.5 ms 3.9× 30 ms 5.3×
Klara (cloth) 3.1 % 4.8 % 8.2 % 34 % 23 ms 5.2× 139 ms 27×
Bunny 12 % 2.6 % 37 % 20 % 19.9 ms 26× 160 ms 96×
Chimpanzee 4.2 % 16 % 29 % 15 % 88.3 ms 16× 246 ms 109×
Huangshan pine 9.2 % 17 % 22 % 10 % 100 ms 2.4× 437 ms 8.5×
Table 3.2 Runtime statistics for average computation time (at) broken down in terms of con-
structing positionsp= Wv (at:pWv), constructing right-hand sides (at:rhs), collision detection
(at:col), and solving the linear system (at:solve). The average total computation time per
iteration (at:total) and the average speedup compared to the full method (at:sp) are also given.
The table also reports the maximum computation time (mt) during the manipulation (mt:total)
and its speedup (mt:sp). The implementation includes multi-threading with Intel TBB; it par-
allelizes collision detection, rotation extraction, basis formation, and computing mesh vertex
positions from reduced coordinates. It also accelerates the full method with the same amount of
parallelism, for a fair comparison. No GPU computation is employed.
Example hier clust bases proj total
Armadillo 0.55 s 0.32 s 0.91 s 1.0 s 2.8 s
Klara (cloth) 0.83 s 2.5 s 0.87 s 2.4 s 6.6 s
Bunny 5.7 s 5.1 s 10 s 40 s 61 s
Chimpanzee 8.8 s — 24 s 55 s 88 s
Huangshan pine 21 s — 52 s 411 s 484 s
Table 3.3 Preprocessing Statistics for building the hierarchy (“hier”), clustering rotations
(“clust”), creating all basis vectors (“bases”), projectingL (“proj”). Column “total” gives the
total preprocessing time for each model in Table 3.2. For the Chimpanzee and Huangshan pine
examples, rotation clustering is not used because many clusters would be needed to accurately
model deformations on cylindrical-like shapes such as branches and limbs. Pre-processing such
clusters was relatively expensive and only gave marginal runtime speed improvements.
Figure 3.2 gives a comparison to the method of Harmon et al. [49]. Different from the
iterative system described in this chapter, Harmon’s method operates as a post-process. So
their method is compared by performing Harmon’s collision projection after a collision-free
solve in each iteration. The experiment (armadillo; Figure 3.2) shows that the hierarchical
method outperforms Harmon’s method in output quality. This is because the contact and elastic
energy are resolved together in one solve, producing smooth collision-free shapes. In contrast,
Harmon’s method ignores smoothness and only seeks the closest collision-free configuration.
In Figure 3.14, it is shown that the hierarchical method can work with other deformation
energies beyond ARAP. The one used in the figure is the volumetric strain energy (Section 5.1
of [20]) combined with Projective Dynamics [20]. To use this energy, it is required to compute
3.7 Results 40
Figure 3.12 Editing the bunny with self-contacts. Top row: the rest shape. Middle row:
edited shape with both ears in contact (left) and the same handle movement but without contact
handling (right). Bottom row: back view of the edited shape with (left) and without (right)
contact handling.
the hierarchical basis functions, and then project and solve the shape editing problem in this
basis, as usual.
In the bunny example (Figure 3.12), one ear of the bunny was pulled to collide against the back.
The bottom of the bunny is fixed via constraints. Then, the other ear is pulled to collide with
the first ear, demonstrating self-collisions. Another example involves edits on the armadillo
with self-collisions (Figure 3.13).
3.7 Results 41
Figure 3.13 Editing the armadillo with self-contacts. Left: the rest shape. Upper-right: the
edited shape with one hand in contact with the head. Lower-right: same handle movement but
without contact handling. Note that a finger penetrates into the head.
Figure 3.14 Comparison between the ARAP energy and Projective Dynamics. Although
the energies are different, the output visual differences are minor in this example.
In the Klara example, the user edited cloth in contact with an external object (Klara’s body;
Figure 3.15). At the top, the cloth is constrained to Klara’s chest. The cloth is composed of
11,718 vertices and 40,399 tetrahedra. Because the subspace is only defined on tetrahedral
meshes, the triangle mesh was extruded to create a closed, manifold mesh, and convert it into
a tetrahedral mesh using Tetgen [108]. The cloth model generates many collisions, which
are successfully resolved by the hierarchical method. The user adjusts the shape of the cloth,
similar, say, to silhouette sculpting in film post-production.
3.7 Results 42
Figure 3.15 Modeling cloth in contact with the body. The user adds two handles on the
cloth and adjusts the cloth shape while the subspace hierarchical algorithm properly respects
collisions with the body. Different levels of activated control points are visualized in green
(L0), yellow (L1) and brown (L2) in the right figure. The user handles are shown blue and the
collided manipulators are red. 11,518 cloth triangles, 32 contacts. Shape-editing runs at 43 FPS
on average.
The Huangshan pine triangle mesh is embedded into a tetrahedral mesh of 111,364 vertices
and 348,467 tetrahedra. In this challenging contact example, the user pulls several branches to
collide with other branches (Figure 3.3). The total number of distinct collision sites appearing
during the entire modeling session is ∼40. Static damping is added at the end of the session
to stabilize the results. The method is stable and makes it possible to perform the intended
modeling at a 8x speedup in maximum time, compared to full method (Table 3.2), despite
numerous contact sites and many activated basis functions.
Chapter 4
Adding Physics Using FEM Simulation
In this chapter, we present a complete animation system which adds contact handling and elastic
response using physically based Finite Element Method (FEM) simulation. The problem of
combining keyframe animation with physically based simulation is one of the most interesting
and important ones in computer animation, and has been examined in several publications in
recent years. This system adds simulated tetrahedral mesh motions onto existing tetrahedral
animation. It also enables control over the combination of edited animation from the input
and simulated result both in space and time. Given a triangle mesh animation from arbitrary
geometry input, which can be non-physical, crude or even incomplete, the system creates an
output triangle mesh animation that has been augmented with contact response and physical
secondary motions. The artist provides weights, specified using a minimal user interface, for
how much physically based simulation should be allowed to modify the animation in any region
of the model, and in time. The system then computes a physically based animation that is
constrained to the input animation to the amount prescribed by these weights. This permits
smoothly turning physics on and off over space and time, making it possible for the output to
strictly follow the input, to evolve purely based on physically based simulation, and anything
in between. The artist exerts this control by a semi-automatic scalar field painting on the
tetrahedral mesh. To achieve the tradeoff between physics and input, two blending methods
are presented: (1) using an implicit blending force during the simulation, and (2) blending the
shapes after the simulation. For both methods, it demonstrates how to incorporate collision
and self-collision response that simultaneously preserves the input motion. The end result is an
artistically directed simulation: it conforms to the input by a prescribed amount, but is enriched
4.1 Overview 44
with physics, secondary motion and collision response. The method operates in the full FEM
space without model reduction, simulating both global and local deformations. In addition to
providing good dynamics, it can also infer missing data in a triangle mesh, making it suitable
for completing deformed animation frames based on a rest configuration template.
Achieving such results requires a careful combination of several system components. These
components will be discussed and analyzed in this chapter, including converting triangle mesh
animations into animations of the simulation mesh, constrained simulation, proper blending of
physics and input, weights painting, and resolving collisions and self-collisions.
4.1 Overview
The input of the whole system is the triangle meshΓ and an animation of a submeshΓ
′
⊆Γ. The
output of this stage is an animation ofΓ that is identical to the edited animation from the first
stage in a user-selected subsetΓ
′′
⊆Γ
′
. Elsewhere onΓ
′
, it strikes a balance between physically-
balanced simulation and following the edited animation. This balance is user-controlled: it
can be spatially-varying and painted on the mesh by the artist. The animation in Γ∖Γ
′
is
reasonably and automatically extended fromΓ
′
using a physical process. If desired, the output
animation also obeys collisions. The algorithm is designed so that it can optionally preserve
the input detail that is beyond the resolution of the chosen tetrahedral mesh, everywhere onΓ
′
,
unlike typical physically based simulation that either requires a computationally slow detailed
tetrahedral mesh to produce rich detail, or else detail is destroyed using coarse meshes.
The pseudocode of the stage is given in Algorithm 1. First, a tetrahedral mesh is created to
represent the volume of the input shape (Section 4.2). Then, a fitting process is launched to
create a tetrahedral mesh animation that resembles the input polygon animation (Section 4.3).
After that, the tetrahedral submesh that follows the input motion is selected. It then runs
a physically based simulation with boundary constraints to produce physical motion on the
unfollowing parts of the mesh (Section 4.4), and optionally preserve the dynamic spatial high-
frequency detail (Section 4.5). The user can control how much physics is added to the animation
by painting weights on the tetrahedral mesh. The tradeoff between input animation and physics
is achieved either via implicitly integrated forces driving the simulation to the input animation,
or by kinematically blending the input with the simulation result (Section 4.6).
4.2 Tetrahedral Mesh Generation 45
ALGORITHM 1: Add Physics to Animation
input :input meshΓ, its submeshΓ
′
, sub mesh animation { ¯ q
i
}
output :input mesh animation {q
i
} with physics added
FunctionaddPhysics(Γ,Γ
′
,{ ¯ q
i
})
{Ω,C,C
′
}←createTetMesh(Γ,Γ
′
); ▷ C,C
′
are interpolation matrices toΓ,Γ
′
{ ¯ u
i
}←fitAnimation(Ω,Γ
′
,C
′
,{ ¯ q
i
})
assignΩ
′′
⊆Ω that follows input
if use implicit blending force then
initialize blending force f
w
according to weights W
else
f
w
← 0
{u
i
}←simulate(Ω,Ω
′′
,{ ¯ u
i
}, f
w
)
if use kinematic blending then
{u
i
}←blendMotions({u
i
},{ ¯ u
i
},W)
{q
i
}← C{u
i
}; ▷ interpolate u
i
to q
i
if preserve high-frequency detail then
{q
i
}←preserveDetail(Ω,C
′
, ¯ q
i
,{ ¯ u
i
},{u
i
})
return {q
i
}
end
4.2 Tetrahedral Mesh Generation
Given the input triangle meshΓ, the system creates a tetrahedral simulation meshΩ that encloses
the space occupied by the undeformed triangle mesh. Many existing tools such as Tetgen [108]
and NetGen [101] can generate good-shaped tetrahedral meshes given manifold surface triangle
meshes. In order to generate a quality tetrahedral mesh from an arbitrary, non-manifold input,
it first computes a signed distance field using the method described in [ 139]. Next, it uses an
iso-surface meshing algorithm [18] to generate a manifold triangular surface. Recent work by
Sacht et al. [98] can be used to coarsen the surface mesh with collision avoidance. Finally, a
tetrahedral mesh is built using constrained Delaunay tetrahedralization [106, 107].
By adjusting the iso-surface and tetrahedral mesher parameters, the user can obtain a tetrahedral
mesh at whatever resolution desired. The following stages in the system then adds edits and
physical effects that are resolved by the chosen tetrahedral mesh, but keeps the input high-
frequency spatial component of the animation unmodified. Note that this is different from
merely preserving the spatially high-frequency detail in the neutral pose: it does not preserve
4.2 Tetrahedral Mesh Generation 46
Figure 4.1 Preserving dynamic spatial high-frequency detail: Top: input animation of the
plant withering (aging), containing spatial high-frequency motion on the leaves. The model
and animation were downloaded from the Internet. Middle: physically based simulation loses
high-frequency detail. Bottom: the system preserves high-frequency detail. Tetrahedral mesh
can be coarse and is shown overlaid. The fitting process can be seen as extracting the spatial
low-frequency part of the motion and modifying it using physically based simulation. The
high-frequency motion can be added back depending on needs.
just the static geometric detail of the neutral pose (as has been done in many prior papers), but
also the input spatially high-frequency animation of the detailed geometry (Figure 4.1).
The resulting tetrahedral mesh encloses the space of the input mesh and has good quality.
However, extra effort is needed to accommodate a common practical situation: the undeformed
triangle mesh contains self-intersections (such as in the case of human lips in Figure 4.2 and
dragon horns in Figure 4.10), causing the mesh to “weld”.
4.3 Tetrahedral Mesh Fitting 47
Figure 4.2 Self-intersections in the input mesh (human face). Left: the lips collide in the
highlighted region. Right: an interior view of the self-intersection. The system accepts such
ill-formed inputs using the method described in Chapter 5.
Figure 4.3 Completing a tracked face using physics: Due to occlusions, vision tracking
systems cannot capture all parts of a human head (mouth cavity, inner lips, tongue, etc.). We
utilize physically based simulation to produce the deformations of the untracked parts. Left
(the input): tetrahedral mesh (18,391 tets), the tracked triangle mesh (49,924 triangles) and the
neutral complete triangle mesh of the head. Middle, right (the output): the deformed tracked
mesh and the complete mesh computed using our method. 960 frames, 1.35 sec per frame.
4.3 Tetrahedral Mesh Fitting
The input to fitting is the triangle mesh Γ and an animation of a submeshΓ
′
⊆Γ. The animation
ofΓ
′
consists of frames i= 0,...,T, each given by a displacement vector away from the neutral
pose ofΓ
′
. This input animation can be created using any means available: it can be keyframe-
animated, itself animated using some physical process, or it can be scanned from a real subject.
In several examples, it simply hasΓ
′
=Γ. TheΓ
′
⫋Γ case occurs, for example, with the human
face geometry, whereΓ refers to the entire human head, including the complete lips (both
inside and outside the mouth cavity) and the internal mouth structure (the “artist mesh”). The
computer vision tracking system, however, only manages to track a proper subsetΓ
′
that misses
4.3 Tetrahedral Mesh Fitting 48
the interior part of the lips, the internal mouth cavity and large peripheral regions of the face
(Figure 4.3).
For every input frame of the triangle meshΓ
′
, the system computes a tetrahedral mesh deforma-
tion that best conforms to it. This is done using an energy function that balances the match to
the input frame and the elastic strain energy of the tetrahedral mesh. Given the tracked triangle
mesh with t vertices and its deformation ¯ q
i
∈R
3t
at a frame i, the deformation of the tetrahedral
mesh ¯ u
i
∈R
3n
is obtained by minimizing
¯ u
i
= argmin
u
1
2
∥C
′
u− ¯ q
i
∥
2
Λ
+γΦ(u), (4.1)
where C
′
∈R
3t×3n
is the interpolation matrix that interpolates tetrahedral mesh deformation
toΓ
′
,Φ(u) is an internal elastic energy for deformation u to regularize the shape, andγ > 0
determines the tradeoff between the two terms.
Diagonal matrixΛ is used to weight the interpolation energy according to the surface area
of each vertex, ∥x∥
2
Λ
= x
T
Λx. The invertible St.Venant-Kirchhoff (StVK) elastic energy [55]
is chosen because it is robust, because it preserves volume under large stretches better than
the linear co-rotational material (Figure 4.4), and because its evaluation speed is similar to the
corotational linear FEM [113]. One limitation of such invertible models is that the energy is
not consistent in the extrapolated regions corresponding to extreme compression and inversion.
During experiments, there are no related stability problems encountered; but if needed, the
issue can be remedied using [115]. The interpolation matrix simply consists of the barycentric
coordinates of every triangle mesh vertex in its containing tetrahedron. One could instead adopt
other coordinates to perform the interpolation, such as [62] or [61]. The parameterγ is reported
in Table 4.1 and parameter tuning is discussed in Section 4.8.
For each frame i, given an initial guess, the tetrahedral mesh deformation is computed at
that frame using a process similar to numerical timestepping. The energy in Equation 4.1
is approximated using a second-order Taylor series in u. We then minimize this quadratic
function and iterate. The quadratic approximation involves the tangent stiffness matrix (second
derivative) of the StVK elastic energyΦ(u), which can be computed easily. The algorithm
computes the derivative of the fitting energy and set it to zero:
C
′T
Λ(C
′
u− ¯ q)+γF
int
(u)= 0, (4.2)
4.3 Tetrahedral Mesh Fitting 49
Figure 4.4 Comparison of invertible StVK to corotational linear FEM. The input is a
severely stretched cube triangle mesh. With the same amount of total fitting error, the shape
produced by the corotational linear FEM (right) has self-intersections along the cube edge,
whereas the invertible StVK (middle) has no such artifacts.
where F
int
is the internal elastic force. Given an initial guess u
0
for u, it approximates Equa-
tion 4.2 with:
C
′T
Λ(C
′
∆u+C
′
u
0
− ¯ q)+γ(K(u
0
)∆u+F
int
(u
0
))= 0. (4.3)
Here, K(u
0
) is the tangent stiffness matrix of the elastic energyΦ evaluated at u
0
. By solving
the above sparse system for∆u, the next iterative approximation for ¯ u
i
is obtained. The iteration
continues until the relative difference between the two consecutive iterations becomes small
enough (1% in the examples). For the first frame in the animation, it uses the rest shape as the
initial guess and apply a line search after each iteration to stabilize the optimization. For the
other frames, it uses the deformation of the previous frame as the initial guess.
Such Newton’s method is faster than a gradient-only optimization method such as the conjugate
gradient optimizer [95] applied directly to the energy in Equation 4.1, with no visible loss of
quality. There is an observed 30× speedup in the examples. Although the algorithm proceeds
frame by frame, it can still exploit parallelism since the bottleneck of the computation is the
evaluation of the second derivative ofΦ(u) and the linear system solve, which can be easily
parallelized. One can instead parallelize the solver across multiple frames. However, such
independent parallel solves present a difficulty for maintaining temporal coherence across the
fitted frames, whereas sequential implementation is less likely to create temporal artifacts. In
practice it gives good results without noticeable temporal artifacts.
4.4 Constrained Simulation 50
4.4 Constrained Simulation
The user first selects a subset tetrahedral mesh Ω
′′
ofΩ that will be constrained during the
simulation to follow the input. The animation ofΩ∖Ω
′′
is computed using a physically based
simulation where the vertices inΩ
′′
serve as boundary constraints. The equation of motion of a
FEM nonlinear deformable object is
M ¨ u+D(u)˙ u+ f
int
(u)= f
ext
(t), (4.4)
subject to u
c
= ¯ u
c
, (4.5)
where M is mass matrix, u is the vertex displacements ofΩ, f
int
(u) is the internal elastic force,
f
ext
(t) is the external force including gravity and contact forces, and D(u) is the damping
matrix. The vertex displacements ofΩ
′′
as input to this stage and during the simulation are
denoted by ¯ u
c
and u
c
, respectively. The same invertible StVK energy is used as in the fitting
process because of its robustness and fast evaluation. Rayleigh damping is adopted in all the
examples: D(u)=αM+βK(u), whereα andβ are Rayleigh damping parameters and K(u) is
the tangent stiffness matrix at u.
Here v= ˙ u, representing vertex velocities. Given the state (v
k
,u
k
) at frame k, The simulation
method performs one iteration of the implicit backward Euler integrator to get∆v so that the
velocity at next frame is v
k+1
= v
k
+∆v, and the displacement is u
k+1
= u
k
+∆tv
k+1
, where∆t is
the timestep.
In order to compute∆v, the linear system is partitioned into
⎡
⎢
⎢
⎢
⎢
⎣
A
f f
A
f c
A
c f
A
cc
⎤
⎥
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎢
⎣
∆v
f
∆v
c
⎤
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎣
b
f
b
c
⎤
⎥
⎥
⎥
⎥
⎦
, (4.6)
⎡
⎢
⎢
⎢
⎢
⎣
A
f f
A
f c
A
c f
A
cc
⎤
⎥
⎥
⎥
⎥
⎦
= A= M+∆tD(u
0
k
)+∆t
2
K(u
0
k
), (4.7)
⎡
⎢
⎢
⎢
⎢
⎣
b
f
b
c
⎤
⎥
⎥
⎥
⎥
⎦
= b=∆t(f
ext
− f
int
(u
0
k
)−(∆tK(u
0
k
)+D(u
0
k
))v
k
), (4.8)
where ∆v
c
and ∆v
f
correspond to the DOFs in Ω
′′
and Ω∖Ω
′′
, respectively, and u
0
k
is the
linearization state for internal elastic forces, discussed further below. The boundary constraint
∆v
c
=
¯
∆v
c
implies that the motion inΩ
′′
should follow the input motion. Therefore, we only
4.4 Constrained Simulation 51
need to solve
A
f f
∆v
f
= b
f
−A
f c
¯
∆v
c
. (4.9)
Intel PARDISO, a direct sparse solver is used for implementation. While a fully implicit Euler
involves multiple iterations of the Newton-Raphson method, in practice Equation 4.9 is solved
only once per timestep. In the experiments, for a fixed amount of computation cost, it was
better to simply subdivide the timestep, than to invest the same amount of computation time
into multiple Newton-Raphson iterations of a single timestep. A comparison is made where it
ran (A) one timestep of size∆t with five Newton iterations, versus (B) five timesteps of size
∆t/5 with one Newton iteration. It was found that the largest stable timestep under (B) is 5-10×
times larger than (A), plus (B) is higher quality (less damped) due to the smaller timestep.
Two choices for the evaluation state u
0
k
for internal forces were experimented with. The choice
can lead to a substantial visual influence on the output animation (Figure 4.5). One choice,
u
0
k
= u
k
+∆tv
k
, is used by Bergou [14], while the other is u
0
k
= u
k
. It is observed that the latter
produces stabler simulations, allowing very large timesteps with minimal tuning, at the cost of
more damped motion. The system gives a user control on the initial guess as u
0
k
= u
k
+δ∆tv
k
,
0≤δ ≤ 1, so that the artist can choose to prefer accuracy or stability given same amount of
computation cost.
To keep the integrator stable, the animation timestep is subdivided if necessary. One limitation of
this constrained simulation is that the boundary constraints onΩ
′′
only provide aC
0
deformation
continuity across the interface between Ω
′′
and Ω∖Ω
′′
. This is because linear tetrahedral
elements is used for simulation. Deformations are linear per element andC
0
across elements.
However, in order for non-C
1
continuity to be visible, shear is generally required, and solid
simulations tend to avoid shear. Continuity can also be improved with sufficient stiffness around
the interface, or with additional averaging or soft constraints around the interface. The user can
also use the blending methods described in Section 4.6 to blend the simulated motion around
the interface with the scripted motion to reduce possible discontinuity.
4.5 Preserving High-Frequency Detail 52
Figure 4.5 Comparison of different internal force evaluation strategies: A simple bar vi-
brates with one end fixed. The yellow shape was obtained using u
0
k
= u
k
+∆tv
k
(Bergou’s initial
guess). The gray shape was obtained using u
0
k
= u
k
(the other). The initial displacement and
velocity are the same in both methods. In the first row, ∆t= 0.003. Bergou’s initial guess exhibits
a less-damped motion. In the second row,∆t= 0.1. A position constraint is added to one vertex
(shown as a red dot). The other initial guess is stable under the large timestep, whereas the
Bergou’s is not.
4.5 Preserving High-Frequency Detail
After the simulation, the output is the deformation q
i
of Γ at frame i as q
i
= Cu
i
, where C
interpolates tetrahedral mesh deformation toΓ. If the input triangle mesh animation contains
spatial high-frequency motion, it will be lost during this interpolation step because the simulation
mesh is typically much coarser than the triangle mesh. Note that this only applies to dynamic
detail in the animation; any spatial detail present in the undeformed meshΓ is always preserved.
In this system, the dynamic detail (as opposed to only static) is also preserved as follows (see
also Figure 4.6). At each frame, for each tracked triangle mesh vertex, the difference between
the input vertex displacement and the displacement interpolated using the fitted tetrahedral
mesh deformations ¯ u
i
is stored. This difference is then rotated and added to the interpolated
displacement obtained using simulation mesh deformations u
i
, effectively reproducing the
4.6 Balancing Input Motion with Physics 53
high-frequency detail
q
′
i
= C
′
u
i
+R
i
(¯ q
i
−C
′
¯ u
i
), (4.10)
where q
′
i
∈R
3t
is the final displacement on the tracked mesh Γ
′
at frame i, R
i
∈R
3t×3t
is a
block diagonal matrix formed by 3×3 rotation matrices. Each rotation matrix is extracted by
polar decomposition from the relative deformation gradient between the fitted and simulated
tetrahedron containing each vertex.
This procedure has the property that the dynamic detail onΓ
′′
is always preserved. For any
triangle mesh vertex inΓ
′′
chosen to follow the input, its embedding tetrahedron should be
inΩ
′′
, which is constrained throughout the simulation. Consequently, the deformations of
this tetrahedron in the input shape and in the simulated shape are identical, at any frame. The
rotation matrix is thus always identity, causing q
′
i
= ¯ q
i
. After the transfer of the detail, the vertex
position is guaranteed to be identical to the input; i.e., the input is exactly preserved inΓ
′′
. The
result of this procedure is shown in Figure 4.1. An alternative method was considered before
where the tetrahedral deformation gradient is directly used to transform ¯ q
i
−C
′
¯ u
i
. While this
alternative method worked on most examples, it failed in cases where the input animation causes
large compressions in the fitted tetrahedral mesh shape (Figure 4.7). This produced nearly
singular deformation gradients on compressed tetrahedra, resulting in mesh “spikes” when
the fitting procedure causes the triangle mesh vertex to escape out of its original tetrahedron.
Altering barycentric coordinates for each individual frame does not work well either. It would
have difficulties if the tetrahedral mesh self-collides during the fitting, which leaves multiple
tetrahedra containing the same embedded vertex, or for vertices outside of the mesh with
poorly conditioned neighboring tetrahedra. The method does not suffer from these issues. One
limitation of using rotations is that they have slight discontinuities at tetrahedra boundaries. In
practice, the rotations of neighboring tetrahedra are similar, and no artifacts were observed.
4.6 Balancing Input Motion with Physics
In the previous section, the output animation comes either from physically based simulation,
or from the input directly. This section describes and compares two procedures that permit
the artist to control the balance between physics and input animation. Their advantages and
disadvantages are also discussed. The balance is prescribed by spatially-varying tetrahedral
4.6 Balancing Input Motion with Physics 54
Figure 4.6 Preserving dynamic high-frequency detail: The difference between the input
vertex displacement ¯ q
i
and the interpolated displacement ¯ q
∗
i
of the fitted triangle ¯ u
1
¯ u
2
¯ u
3
(left) is
rotated by the rotation component R of the deformation gradient between the two triangles (2D
example). It is then added to the interpolated displacement q
∗
i
of the simulated triangle u
1
u
2
u
3
(right) to give the final displacement q
i
.
Figure 4.7 High-frequency detail transfer failure produces spikes on the mesh (left) when
using tetrahedral deformation gradients to transfer detail from highly compressed tetrahedra, as
contrasted with the smooth shape (right) produced by transfer with the rotation component (used
in this thesis). An alternative method to preserve high-frequency detail is to alter barycentric
coordinates for the triangle mesh at each frame. This method also suffers from such a “spike”
failure.
mesh vertex weights w
i
≥ 0. How the artist can easily generate such weights is explained in
Section 4.6.3.
4.6.1 Implicit Blending Force
In the first method, an implicitly integrated external force is added to drive the mesh toward the
input, scaled by the weights
f
w
= W(¯ u−u), (4.11)
4.6 Balancing Input Motion with Physics 55
where W ∈R
3n×3n
is the diagonal matrix of weights w
i
for each vertex. Higher weights generate
motion closer to the input. In industry, such forces are often referred to as the kin(ematic)
springs. They can be integrated implicitly in Equation 4.9. Because Equation 4.11 is perfectly
linearized, the Jacobian of this force is exact. The implicit integration thus has an exact
prediction of the force. Besides, the presence of W increases only diagonal entries in the system
matrix, making it more well-conditioned. As a result, the integration of this force is highly
stable. Stable convergence to ¯ u has been experimentally observed as W grows to extremely
large values such as 10
30
, well beyond values that produce visually identical motion as ¯ u. It
is convenient for the artist to express weights on the interval [0,1]; the user interface re-maps
them to [0,∞] using the function f(x)=η(x
1/ζ
), where the artist can adjust parametersη > 0
and ζ > 1. Contact can be resolved during the simulation; no special handling is required.
For simplicity, the penalty-based method is used to handle contact, but any contact resolution
method could be used. The disadvantage of this method is that the relationship between the
weights and the input vs. physics tradeoff, while monotonic, is only indirect and unintuitive.
4.6.2 Kinematic Blending
The second approach is similar to as-rigid-as-possible interpolation [1, 54] that blends the
output physical shapes and the input shapes. It decomposes the deformation gradient for
each tetrahedron into a rotation and symmetric matrix, blend them separately according to the
weights, and solve a linear system to combine the deformation gradients of all tetrahedra as in
rotation-strain coordinate warping [110]. Weights used in this procedure are intuitive as they
can be assigned directly on [0,1]. The method permits easy post-simulation tuning of weights,
upon which one can resolve for the blended animation without re-simulating physics.
Such an approach, however, is not compatible with contact resolution. Therefore, the following
strategy is used to resolve contact. After kinematic blending is used after one simulation step to
produce a blended shape, if any collision is detected in the blended shape, penetration depths are
collected and penalty forces calculated with respect to the blended shape (Section 4.7). These
penalty forces cannot be applied to the simulation mesh directly, since the blended shape usually
has a different shape and orientation as the simulation mesh. Therefore, a per-tetrahedron
rotation is calculated, based on the deformation gradient of each tetrahedron of the blended
shape relative to the same tetrahedron in the simulation mesh. It then rotates the penalty forces
with these rotations to get the correct force directions before adding the contact forces to the
4.6 Balancing Input Motion with Physics 56
integrator. This strategy works reasonably well in the elephant example (see in Figure 4.8). A
more precise method could compute the gradient of the blending function (such as in [102]);
This is however regarded as too complicated and not pursued in designing the system. After
the two methods are implemented and compared (Figure 4.8), the implicit blending force is
preferred to the kinematic blending because it enables collision handling without any special
additional provisions, albeit at some loss of weight intuitiveness.
4.6.3 Generating Weights
This section explains how to easily determine suitable weights w
i
⩾ 0 for any vertex i of the
simulation meshΩ∖Ω
′′
, for blending purposes. Larger weights bring the mesh closer to the
input motion. An interface is provided to compute all the weights using a bi-Laplace solver.
The Laplacian matrix L is defined as
w
T
Lw= ∑
1≤i< j≤m
σ
i, j
(w
i
−w
j
)
2
, (4.12)
σ
i, j
=
⎧
⎪
⎪
⎨
⎪
⎪
⎩
1 if vertex i and j share an edge,
0 otherwise.
(4.13)
Here, w∈R
m×1
is the vector of the weights, and m is the number of vertices inΩ∖Ω
′′
. The
action of the tetrahedral mesh bi-Laplacian L
T
L is defined similarly as in [ 15, 141], except that
here Laplacian is defined on vertices instead of tetrahedra. The algorithm then minimizes
1
2
w
T
L
T
Lw, (4.14)
subject to Sw= ¯ w, (4.15)
where S is the selection matrix, and ¯ w are the user-defined weights on a few vertices. Because
the constraints are linear and the objective is quadratic, the minimization can be performed
using a single sparse linear system solve.
4.7 Collisions 57
4.7 Collisions
Collisions may occur and need to be handled in various parts of the system. Robust handling of
collisions is generally a difficult, but an essential task for any physically based simulator. This
section explains how to address collisions appearing in this stage.
4.7.1 Collision Mesh Creation
Since the system uses tetrahedral mesh to simulate the motion of the embedded triangle mesh,
it needs to detect collision on the embedded geometry. The input triangle mesh may contain ill-
formed triangles and non-manifold components not suitable for collision detection. Therefore,
the method of [139] that creates the iso-surface in Section 4.2 is re-used to create a manifold
collision surface mesh. The mesh resolution can be adjusted to trade accuracy for collision
speed.
4.7.2 Collision Detection and Response During Simulation
After a collision mesh is created, it is embedded into the simulation meshΩ using barycentric
coordinates. At runtime, collision detection is performed using a bounding volume hierar-
chy [43]. For each colliding collision mesh vertex c, its penalty force f
c
is calculated based on
the penetration depth. The penetration depth is computed in a manner similar to [23]. The force
f
c
is redistributed to the vertices of the tetrahedron that embeds c, weighted by its barycentric
coordinates [112],
f
i
= w
c
i
f
c
, (4.16)
where f
i
is the collision force on tetrahedral mesh vertex i, and w
c
i
is the barycentric weight of c
with respect to i.
4.8 Results 58
statistics elephant dragon bird head plant sumo
vtx 8,400 25,002 408 25,600 961 15,482
tri 16,636 50,000 7,878 49,924 1,600 30,642
tet-vtx 3,850 1,590 3,949 5,729 200 2,717
tet 14,942 5,387 11,846 18,391 413 7,767
free tet-vtx 1,621 1,401 2,414 – 199 1,642
free tet 5,354 4,801 6,915 – 410 4,599
frames 178 420 299 960 120 650
sim-frames 1,602 2,500 600 – 120 3,900
fitting 0.54s 0.11s 0.78s 1.35s 0.024s 0.24s
γ 10
−10
10
−9
10
−5
10
−8
10
−2
10
−9
sim 8.37s 0.68s 0.084s – 0.0058s 0.30s
Table 4.1 Simulation statistics. vtx, tri=#triangle mesh vertices and triangles; tet-vtx, tet=#tet
mesh vertices and tets; free tet-vtx, free tet=#tet mesh vertices and tets in the free regionΩ∖Ω
′′
;
frames=#graphical frames; sim-frames=#simulation frames, including intermediate frames that
were inserted for stable simulation; fitting, sim = time for fitting the tetrahedral deformations
and simulation for one graphical frame, respectively;γ: parameter used in fitting. Intel Xeon
2.3GHz 2×6 cores, 32 GB memory. In the elephant example, 95% of the simulation time is
spent processing collisions. The head example only uses frame fitting; hence, no simulation
time is reported.
4.8 Results
The system has been tested on several examples. Statistics of fitting and simulation are available
in Table 4.1.
Elephant: In the first example (Figure 4.8), it enriches a walking elephant motion with physics.
The motion was created by an artist in Maya, using rigging and keyframing. it adds secondary
motion to the ears, trunk, belly and tail, by putting them intoΓ∖Γ
′′
. In each of these regions,
it also adjusted a single Young’s modulus value to make the region more or less elastic. The
motion has 178 frames. It took three minutes to generate the tetrahedral mesh. Then, it iterated
seven times to fit the first several frames of the animation, searching for a good γ (Equation 4.1).
This process took eight minutes including user interaction and running the fitting. It took about
90 seconds to produce the entire fitted animation with the desired γ. Collision detection and
contact resolution were enabled in the simulation. Figure 4.8 shows the blended result, where it
blends between a physically based simulation and the input with different weights, using the
blending method of Equation 4.11. The result of using kinematic blending is also shown: it
produces a similar motion as implicit blending. Figure 4.14 demonstrates that it can also enrich
input artist animations with large physically based contact. Such contact can be useful when
4.8 Results 59
Figure 4.8 Adding physics to an elephant walk cycle. Left-most column: the fixed Ω
′′
and
freeΩ∖Ω
′′
region of the elephant tetrahedral mesh. Free regions include the trunk, ears, belly
and tail. Left four columns: Each column is a separate animation with a different blending
weight (Equation 4.11) between the artist input and physics. Implicit blending method was
used. From left to right, the weights are decreased from following the input 100%, to a 100%
physically enriched result. InΓ∖Γ
′′
, the method produces secondary motion according to the
weights, whereas inΓ
′′
, the output is identical to the input. Right-most column: kinematic
blending with 50% physics. It produces similar results as the implicit blending method, but
suffers from less intuitive tuning of collision handling, such as between the legs and belly.
the artist already designed the primary motion of a character, and then wants to add additional
(deformable or rigid) objects into the scene that should collide with the deformable character,
without substantially changing the original character motion.
Dragon: The second example (Figure 4.10) is an Asian dragon animation created by a surface
deformation method [114]. The dragon’s horns (nearly) self-collide with the head in the input
mesh, and standard methods produce welded tetrahedral meshes with incorrect topology. The
method described in Section 4.2 is used to create an undeformed tetrahedral mesh with correct
un-welded topology, despite self-collisions. It took three minutes to generate a (non-neutral)
self-intersection-free mesh in Maya. Then it took three minutes to create a tetrahedral mesh,
and one minute to produce the topologically correct neutral tetrahedral mesh. Next, it took six
iterations in four minutes of user time to fit the initial frames, and then 46 seconds to fit the
animation. The feet of the dragon are constrained to perform a rising action similar to the one
used in [19]. the entire dragon except the back feet is included intoΓ∖Γ
′′
. Implicit blending
forces are added on the head and front feet. It produces vivid motion on the horns, whiskers,
body and tail. By adjusting the weights temporally on the head and front feet, it makes them
4.8 Results 60
follow the input strictly during the initial part of the animation, and exhibit secondary motion
afterwards.
In Figure 4.11, a comparison to TRACKS [14] is made. The original TRACKS is designed
for cloth simulation. It is extended by building Petrov-Galerkin constraints between the input
triangle mesh animation and the embedded mesh animation driven by solid FEM simulation.
The triangle mesh displacements v in those constraintsG(v)= 0 are then substituted with Cu= v,
to create simulation constraints for the tetrahedral meshΩ. It is found that TRACKS constraints
suppress global, spatially low-frequency motion. The generated high-frequency motion is less
pronounced on solids than on cloth. In comparison, our method imposes no restrictions on
Γ∖Γ
′′
to allow secondary motion.
Bird: The third example adds physics to a keyframed bird animation (Figure 4.12). The model
and the animation were purchased online in a 3D model store. the legs were placed intoΓ
′′
, so
that the artist’s previous work on ground contact is preserved. Similarly, the ends of the hands
are placed intoΓ
′′
, so that the hand motion follows the input, which can be seen as a variant
of inverse kinematics. The rest of the hand including the arm feathers, however, is placed into
Γ∖Γ
′′
, which gives secondary motion. It adds physics to the head and eyes.
Sumo Wrestler: In the fourth example (Figure 4.9), the motion of a dancing sumo wrestler is
enriched with secondary motion on the belly, cheeks, derrière and thighs. the hands and legs
(except thighs) are placed intoΓ
′′
to follow the artist’s input. The neck, and parts of the back
and the hips are put intoΓ
′′
, confining the motion to reasonable poses. The lower half of the
head is put intoΓ∖Γ
′′
, to produce secondary motion on the cheeks. The rest of the head is in
Γ
′′
. It produces large deformations when the wrestler is dancing and jumping. The thighs and
cheeks exhibit jiggling motion as well.
Head: The fifth example completes the deformed shapes of a human head triangle mesh Γ,
based on the input motion of a submeshΓ
′
(Figure 4.3). The input motion was computed using
a proprietary optical flow-like tracking algorithm at a major computer animation company. Due
to occlusions, it is missing all the deformation detail in the inner lips and the mouth cavity. The
method is compared to a standard surface deformation method [114] (Figure 4.13). It is found
that the surface deformation method produces collisions, especially in the mouth cavity and
inside the eyes, as these regions contain extremely concave geometry. Surface deformation
methods also require a single connected surface mesh, whereas in the head model the tongue
and teeth are separate triangle shells; them are removed for this comparison. a tetrahedral mesh
4.8 Results 61
Figure 4.9 Adding physics to a keyframed sumo wrestler animation. Left-most column: the
fixed Ω
′′
and freeΩ∖Ω
′′
region of the sumo tetrahedral mesh. Free regions include the body,
thighs, cheeks and derrière. The top row on the right shows several frames of the input motion.
The bottom row on the right shows the result with secondary motion added.
Ω is created for the entire head using the procedure described in Section 4.7. Then tetrahedral
mesh deformations ¯ u
i
are fitted based on the input tracked triangle mesh motion. Since the
tracked data already contains facial dynamics, the simulation part of the pipeline is skipped, and
the complete meshΓ is interpolated using ¯ u
i
. The teeth are handled separately by extracting the
closest rigid transformation to the interpolated shapes because they are rigid. In this example,
No input high-frequency detail is added back as the input detail contains tracking noise, rather
than something that was intelligently designed by the artist (as in Figure 4.1). Actually, the
input data contains tracking imperfections, mostly along the perimeter of the lips, which our
method resolved. 12 facial animations were processed in the example, each containing about
80 frames. Fitting each of these animations took 1.8 minutes on average.
4.8 Results 62
Figure 4.10 Adding physics to a dragon animation. The fixed region consists of the two back
feet. The blending weights (third and fourth subfigures in the first row) change through time to
follow input at first, and add secondary motion afterwards. Note that the weights on the back
feet are ignored because they are inΓ
′′
. The second row shows the resulting animation where
the dragon rises up. For cinematic effect, a particle-effect fire animation is added. Secondary
motion is present in the horns, whiskers, body and tail.
Figure 4.11 The method described produces more global motion than TRACKS with 40
automatically-generated regions. Two comparisons are made at each frame: (1) the Euclidean
distance between the displacement u
k
∈R
3n
and the input frame ¯ u
k
∈R
3n
, and (2) the difference
in the x-axis position of a selected dragon vertex between the output and input.
Fitting Parameters: The meaning ofγ changes if one adjusts the material stiffness (Young’s
modulus). The default Young’s modulus for fitting is 10
7
N/m
2
. Some examples require het-
erogeneous materials, in which case the average Young’s modulus is kept in the same range.
Poisson’s ratio of 0.45 is used in all the examples.
4.8 Results 63
Figure 4.12 Adding physics to a keyframed bird animation. The top row shows the fixed
Ω
′′
and freeΩ∖Ω
′′
bird regions. Free regions are the head, belly and wings, except the wing
tips which are fixed to drive the wings, similar to inverse kinematics. The middle row shows
several frames of the input motion. The bottom row shows our result with subtle secondary
motion added. This example does not use blending.
A wide range of γ is used across our examples. A good γ depends not only on the size and
overall stiffness of the object, but also on the input animation quality. In the elephant example,
the input animation has quite large non-smooth deformations on the legs which generate huge
internal elastic potential. A very smallγ (10
−10
) has to be used to achieve a tight fit. In contrast,
the motion of the plant is very gentle and smooth, and even a largeγ can produce a good result.
Besides, a largeγ (10
−2
) helps filter out the local high-frequency dynamic details, which were
removed during fitting, and added back after the simulation.
4.8 Results 64
Figure 4.13 Comparison to as-rigid-as-possible energy: Left: the deformation produced by
our method. Right: the shape produced by as-rigid-as-possible energy [Sorkine and Alexa
2007]. In order to improve the as-rigid-as-possible result, the negative weights in the as-rigid-as-
possible energy are clamped; otherwise, the as-rigid-as-possible method causes the obtuse mesh
triangles to collapse. The comparison were performed by setting the positions of the tracked
meshΓ
′
as constraints forΓ and minimized the as-rigid-as-possible energy. Lacking interior
information, the as-rigid-as-possible energy resulted in a colliding shape at extremely concave
regions such as the mouth cavity and eyelids, whereas our method produces a good result.
Figure 4.14 Adding physically based contact to artist animations: The soft ball is launched
at the elephant. The ears, trunk, belly and the tail are enriched with physics just as in Figure 4.8.
The system performs collision detection between the elephant and the ball. Upon impact,
contact forces are computed and added to both objects, resulting in local contact deformations
on both objects, all the while the elephant still generally follows the original animation.
Chapter 5
Immersion of Self-Intersecting Solids and
Surfaces
In this chapter, we give an algorithm to analyze self-intersecting 3D polygon meshes and build
self-intersecting tetrahedral meshes to conform to the input topology. Self-intersecting, or
nearly self-intersecting, meshes are commonly found in 2D and 3D computer graphics practice.
Self-intersections occur, for example, in the process of artist manual work, as a by-product of
procedural methods for mesh generation, or due to modeling errors introduced by scanning
equipment. If the space bounded by such inputs is meshed naively, the resulting mesh joins
(“glues”) self-overlapping parts, precluding efficient further modeling and animation of the
underlying geometry, even if users have robust collision handling methods in modeling and
animation phrases of the workflow. Similarly, near self-intersections force the simulation
algorithm to employ an unnecessarily detailed mesh to separate the nearly self-intersecting
regions.
The method presented in this chapter addresses both of these challenges, by giving an algorithm
to generate an “un-glued” simulation mesh, of arbitrary user-chosen resolution, that properly
accounts for self-intersections and near self-intersections. In order to achieve this result, the
mathematical concept of immersion is studied, and a deterministic and constructive algorithm is
given to determine if the input self-intersecting triangle mesh is the boundary of an immersion.
For near self-intersections, a robust algorithm is given to properly duplicate mesh elements
and correctly embed the underlying geometry into the mesh element copies. Both the self-
intersections and near self-intersections are combined into one algorithm that permits successful
5.1 Topology of Self-Intersecting Meshes 66
meshing at arbitrary resolution. The method enables robust volumetric shape editing, physically
based simulation and animation, volumetric weight and geodesic distance computations on
self-intersecting inputs.
5.1 Topology of Self-Intersecting Meshes
The input to the method is a self-intersecting orientable triangle mesh M without boundary
(inR
3
), or a self-intersecting polygonal line M without boundary (i.e., closed) (inR
2
). By
a triangle mesh, it hereby means a collection of vertices with positions, and a list of integer
triples specifying the vertex indices for each triangle. While such a data-structure is common
in computer graphics and may even be taken for granted, it is important for the proofs in this
chapter to view M merely as a set-theoretic ordered pair (V,T) of vertices V and triangles T ; M
is not seen as a geometric object and should not be equated with the set of positions of all points
on M. The inputs M must have manifold connectivity in the usual sense, i.e., every edge must
be shared by exactly two triangles, and the triangles touching a vertex must form a continuous
fan. Informally, the method asks how to determine and volume-mesh the “interior volume”
bounded by M, in a way that respects self-intersections and does not “glue” self-intersecting
regions. In the chapter, this intuitive question is properly formulated mathematically using
topology. While it is common in computer graphics to use the informal term “manifold mesh”
to refer to any triangle mesh with manifold connectivity, self-intersecting meshes require
more precise mathematical treatment; otherwise, it is not possible to soundly formulate the
immersion algorithm or prove well-defined statements about it. In topology, a manifold without
boundary of dimension s is a topological space where every point has a neighborhood that
is homeomorphic toR
s
(see Appendix B for precise definitions). Self-intersecting meshes
with manifold connectivity (inputs to the immersion method) are not manifolds because the
intersection points have no such neighborhood.
Given a triangle mesh M with manifold connectivity, one can form a 2-manifold
ˆ
M which
intuitively corresponds to the concept of seeing M as a “topological mesh”, i.e., forming one
continuous surface whereby any self-intersections are properly ignored. Intuitively,
ˆ
M is formed
by “gluing” the triangles along the shared edges. Formally, one can form it by taking a disjoint
union (see Appendix B) of all triangles, and then applying the equivalence relation whereby
points of the same edge shared by two triangles are made equivalent. It is easy to verify that
ˆ
M is a manifold (proof is in Appendix B). This manifold is an “abstract” manifold, and, when
5.1 Topology of Self-Intersecting Meshes 67
M has self-intersections, cannot be directly equated with some manifold inR
d
. For any point
x on any triangle of M, denote byρ(x) the location of x inR
d
. It can be easily seen that the
mappingρ can be augmented to an immersion ˆ ρ of
ˆ
M ontoρ(M) (Appendix B). An immersion
is a continuous map between two topological spaces X and Y that is locally injective, i.e.,
for any x∈ X, there exists a neighborhood of x such that the restriction of the map onto this
neighborhood is injective (i.e., one-to-one).
Further assume that M intersects itself in a general way, i.e., the intersections are not degenerate.
Formally define non-degeneracy as follows. Define the collision set of x∈
ˆ
M as
γ(x)={y∈
ˆ
M; ˆ ρ(x)= ˆ ρ(y)}⊂
ˆ
M, (5.1)
i.e., points that share the same location inR
d
with x. The “self-intersecting set” of M assembles
all points x where there exists at least one more point with the same location,
Γ={x∈
ˆ
M;∣γ(x)∣≥ 2}⊂
ˆ
M. (5.2)
Also define the triple set
Γ
≥3
={x∈
ˆ
M;∣γ(x)∣≥ 3}⊂
ˆ
M. (5.3)
For d= 3, define M to be non-degenerate ifΓ is a union of zero or finite number of 1-dimensional
loops and Γ
≥3
consists of zero or more disjoint points. For d = 2, M is defined to be non-
degenerate ifΓ consists of zero or more disjoint points andΓ
≥3
is empty. For both d= 2,3,
degenerate inputs are prohibited where a surface touches itself at a single point, or along a
loop, but does not penetrate (exact definition is in Appendix B). For d= 2, an example of such
a degeneracy is a circle touching a square, and for d= 3, an example is a solid bowl placed
upside down on a solid box. Non-degenerate inputs correspond to the general way in which
surfaces intersect. The setΓ typically consists of pairs of loops on
ˆ
M. Note that any degenerate
input can be perturbed by an infinitesimal amount to remove the degeneracy. To summarize the
requirements on the input, valid input is defined to be a non-degenerate orientable triangle mesh
M (for d= 3; and a polygonal line for d= 2) without boundary and with manifold connectivity.
Given a valid input M, the algorithm determines if the surface immersion ˆ ρ can be extended
to a volume immersion for d= 3, and if the curve immersion ˆ ρ can be extended to a surface
immersion for d = 2. Formally, it answers the question of whether there exists a compact
5.1 Topology of Self-Intersecting Meshes 68
Figure 5.1 Basic steps of the immersion method are illustrated in this 2D example. Given the
input shape (in this 2D example, a closed polygonal line / curve inR
2
), the method extracts
the cells, patches and arcs (Section 5.2), then duplicates and connects cells to recover the
immersion of a two-dimensional disk such that the disk boundary immerses onto the input shape
(Section 5.3 and 5.4). Note that the largest cell self-touches, which is addressed in Section 5.2.2.
To produce the immersion of the disk, first triangulate the entire space covered by the input
shape, then create a submesh for each cell (Section 5.5). Some triangle vertices are duplicated
to avoid gluing parts of the cell together. For visualization purposes here, a cyan line indicates
that the two triangles joined by it are actually the same triangle. Finally, submeshes for the cells
are merged together according to the cell connectivity recovered by the algorithm (Section 5.5).
Figure 5.2 Examples of self-intersecting curves inR
2
and surfaces inR
3
which are not bound-
aries of immersions.
d-manifold
ˆ
S and an immersion ˆ σ from
ˆ
S into R
d
such that
ˆ
M is the boundary of
ˆ
S, and
the restriction of ˆ σ onto
ˆ
M is ˆ ρ. Compactness, intuitively, means that
ˆ
S contains all of its
limits (exact definition is in Appendix B), and is needed to avoid immersions from unbounded
manifolds, or manifolds
ˆ
S where one has artificially removed points from
ˆ
S; e.g., to rule out
5.2 Cell complex 69
Figure 5.3 Cells, patches and arcs in 2D. Left: input 2D curve. Right: 4 cells (open disks),
8 patches (open curves; Pi with different colors), four arcs (isolated points). Cells and their
B-patches are as follows: Cell 0: P0, P1, P2*, P3*; Cell 1: P2, P3, P4, P6; Cell 2: P4*, P5; Cell
3: P6*, P7. Here, * denotes that the B-patch orientation disagrees with the cell. Arcs and the
topological neighboring patches at those arcs are as follows: Arc 0: (P1,P2), (P1,P3); Arc 1:
(P0,P4), (P2,P5); Arc 2: (P0,P6), (P3,P7); Arc 3: (P4,P7), (P5,P6).
immersions from punctured disks or similar. Valid inputs M for which such a pair (
ˆ
S, ˆ σ) exists
are called volume-immersible. The algorithm also explicitly constructs the immersion if it
exists. Practically, in 3D, this means that one can construct a tetrahedral mesh that meshes
the space “occupied” by M and does not “glue” any self-intersections. A simple illustration
of the immersion method is shown in Figure 5.1. Note that not all closed, orientable and
non-degenerate triangle meshes are boundaries of valid immersions. For example, triangle
meshes shown in Figure 5.2 are not; intuitively, this is because these shapes require an inversion
of the space, which violates local injectivity. The immersion method is able to reject such
inputs.
5.2 Cell complex
A self-intersecting mesh cutsR
d
into multiple components, forming cells, patches and arcs.
These components are glued together “nicely” along their boundaries and form what is known
in topology as a cell complex [34]. These concepts are introduced below. They play a key role
in the immersion algorithm.
5.2 Cell complex 70
Figure 5.4 Cells, patches and arcs in 3D.
5.2.1 Cells, Patches and Arcs
Cells are the connected components ofR
d
∖ρ(M) inR
d
, except the infinite component and the
components with winding number 0. The latter are “interior voids” (cavities) inside objects,
and should not be meshed. Winding numbers will be discussed later in the chapter. Because
ρ(M) is a closed set, cells are open subsets ofR
d
. Topologically, cells are d-manifolds without
boundary (due to being open). For connected inputs M for d= 3, they are the interior of either
a topological solid 3-ball, or a topological solid torus with k≥ 1 handles. Note that the torus
case can occur even if the input mesh M is the boundary of a 3-ball (e.g., Cell 3 in Figure 5.4 is
a torus with k= 1 hole). For connected inputs M for d= 2, cells are 2-manifolds that are the
interior of a topological 2-disk. Justifications of these statements are in Appendix B. The closure
of a cell inR
d
may be “self-touching” (e.g., Cell 0 in Figure 5.3), which will be addressed
later. Patches are the connected components of ˆ ρ(
ˆ
M∖Γ) inR
d
. Topologically, patches are
2-manifolds without boundary (i.e., open disks with k≥ 0 holes) when d= 3, and 1-manifolds
without boundary (i.e., non-intersecting loops, or open non-intersecting curves) when d= 2.
Patches form the boundaries between cells. Arcs are the connected components of ˆ ρ(Γ∖Γ
≥3
).
Arcs are open subsets ofR
d
. They form the boundaries of the closures of the patches, and are
1-manifolds without boundary (i.e., open non-intersecting curves) for d= 3, and isolated points
for d= 2. The cell-patch-arc complex is illustrated in Figures 5.3 and 5.4.
Two patches P
1
and P
2
are called topological neighbors if they are adjacent as per the mesh
connectivity of M; formally, if ˆ ρ
−1
(P
1
)∩ ˆ ρ
−1
(P
2
) is a 1-manifold (d = 3), or a point (d = 2).
5.2 Cell complex 71
Figure 5.5 Finding patch geometric neighbors on self-touching cells. Cell 0 in Figure 5.3
self-touches at Arc 0. First identify self-touching arcs (there is only one in this example;
bottom-left), then find the two components that should be disconnected (bottom-right). Use
them to identify the geometric neighbor pairs at Arc 0 as (P2, P1) and (P1, P3).
Here, X denotes the topological closure operation (Appendix B). This means that topological
neighboring patches neighbor each other in
ˆ
M by a shared arc. For d= 3, this excludes cases
where patches neighbor only in a corner point in
ˆ
M. Figure 5.3 shows examples of topological
neighbors.
5.2.2 B-patches
Each cell is surrounded by one or more patches. During the discussions, it is convenient to
call the patches surrounding a cell the B-patches (boundary patches) of that cell. Two cells are
defined to be neighbors if they share a B-patch. Each B-patch has an orientation induced by the
orientation of M. For d= 3, the orientation of a B-patch agrees with its cell if its orientation is
outward from the cell, and disagrees if inward (Figure 5.3). The condition for d= 2 is that the
B-patch is oriented so that the cell interior is to the right.
For each cell, define two B-patches to be geometric neighbors with respect to that cell, as
follows. First remove any singularities of self-touching cells (such as Cell 0 in Figure 5.3 whose
5.2 Cell complex 72
B-patches P1, P2, P3 meet at Arc 0), as described in Section 5.2.2 and illustrated in Figure 5.5.
Then, two patches are geometric neighbors, by definition, if they are spatially adjacent in R
d
,
i.e., there exists an arc A such that P
1
∩P
2
= A. The purpose of this definition is that if, for a cell,
one stitches all pairs of geometrically neighboring B-patches along the common arc, one obtains
the complete and manifold boundary of the cell. Note that without self-touching singularity
removal, a non-manifold boundary would be produced. As seen in Figure 5.9, there can be cells
with only one B-patch. The immersion algorithm has to treat such a B-patch as a geometric
neighbor of itself (by definition).
Geometrically Neighboring B-Patches for Self-Touching Cells
Self-touching cells are detected as follows. For each arc, observe how the patches connect
to the arc, in some sufficiently small neighborhood of the arc. For valid inputs M, due to
non-degeneracy, there are locally four separate patch connections to the arc, defined as patch
junctions (see Figure 5.5, top-left). Formally, it is the connected components of the intersection
of a sufficiently small d-dimensional ball centered at any one arc point, with the union of
all patches. Note that patches and arcs are open sets by definition; they do not contain their
boundaries. Each of the four patch junctions is a subset of some patch. Some of the patch
junctions may be subsets of the same patch (e.g., two patch junctions are subsets of P1 in
Figure 5.5). A cell is self-touching at an arc if all four patch junctions are subsets of B-patches
of this cell. For convenience, the algorithm for handling self-touching cells is first described
for d= 2. To define geometric neighbors for a self-touching cell, observe that the cell around
the arc locally consists of two manifold regions (cyan and yellow in Figure 5.5, top-right). For
each of the four patch junctions, identify the last line segment on it, i.e., the one that touches
the arc. Denote these segments by e
i
, for i= 0,1,2,3 (see Figure 5.5, bottom-left). Let us orient
the line segments so that the interior of the cell is on the right. Pick e
0
to be one of the two line
segments whose direction points into the arc. Then find the line segment e
1
among the other
three line segments that has the smallest counterclockwise angle with e
0
. Then, e
0
and e
1
form
the boundary of one local manifold region of the cell, and the remaining two line segments
e
2
and e
3
form the boundary of the other manifold region (see Figure 5.5, bottom-right). The
B-patches containing e
0
and e
1
are thus assigned to be geometric neighbors at Arc 0; and same
for B-patches containing e
2
and e
3
. Note that the choice of e
0
is deterministic, but arbitrary; an
alternative algorithm would be to pick e
0
to point away from the arc and search clockwise. In
3D, an arc is an open curve, so one can arbitrarily pick one line segment on the arc, and find the
5.3 Immersion graph 73
triangle on each patch junction that shares the line segment. Then proceed in the same way as
with d= 2, except using dihedral angles between triangles as opposed to angles between line
segments.
5.2.3 Computing the Cell Complex
The immersion method computes the cells and patches robustly using exact arithmetic [142]
available inside libigl [59]. This procedure produces a “cut” version of M, where edges are
added to M along its self-intersections. It utilizes the exact arithmetic kernel of the CGAL
library [27]. The method also identifies all the arcs, the B-patches of all cells, the four patch
junctions at each arc, the topological and geometric neighbors of each patch, and whether their
orientation agrees with the cells. Patches are stored as triangle meshes, given by the indices of
the triangles of the “cut” version of M. Arcs are stored as edges of the “cut” mesh. The method
also computes the winding number [58] of each cell using libigl.
5.3 Immersion graph
Intuitively, immersions can be formed by “unwrapping” the volume “occupied” by M. During
this process, some cells of the cell complex (Section 5.2) will require multiple copies to form
the “unwrapped” d-dimensional manifold
ˆ
S. The core idea of the immersion algorithm is to
construct the manifold
ˆ
S by properly gluing together replicated copies of cells from the cell
complex. The information about which cell copy should be glued to which other cell copy and
where (across which patch) can be naturally encoded into a graph. The algorithm will therefore
operate on graphs whose nodes are duplicated cells of the cell complex of M, and edges denote
“gluing” of cell copies across shared B-patches (see Figure 5.1, top-right).
Note that the graph is in general a multigraph, i.e., two nodes may be joined by more than one
edge. This is because a cell may have more than one B-patch shared by the same neighboring
cell (e.g., Cells 0 and 1 share Patches 2 and 3 in Figure 5.3). Therefore, each graph edge needs
to record which B-patch it crosses. The requirements in this paragraph define the set of eligible
graphs. Additionally, due to the requirement that the input surfaceρ(M) must be the boundary
of the immersed volume ˆ σ(
ˆ
S) (i.e., the restriction of ˆ σ onto
ˆ
M is ˆ ρ), each patch must be on
the external geometric boundary of exactly one node. This node is said to “owns” this patch,
5.3 Immersion graph 74
Figure 5.6 2D examples of simple and non-simple immersions. Top: a disk-topology abstract
manifold
ˆ
S is immersed onto 2D. Its boundary
ˆ
M is immersed onto M. For cell C
i
with winding
number (WN) w, its pre-image, ˆ σ
−1
(C
i
) has w connected components (PI-CC). The restriction
of ˆ σ to each connected component in ˆ σ
−1
(C
i
) is a surjective immersion onto C
i
, and the number
of sheets for the surjective immersion is 1. Therefore, this is a simple immersion. Bottom:
a ring-topology abstract manifold
ˆ
S is immersed onto 2D. Its boundary
ˆ
M is immersed onto
a 2D multi-component curve, which creates three cells. There are two possible immersions
(each row, bottom-right). Both immersions wind the 2D ring twice and connect it back to itself.
Among the three cells, the interesting one is the ring-shaped cell C
1
. In both immersions, C
1
has a winding number of 2, but the number of connected components of its pre-image is only 1
(highlighted in red in the PI-CC column). The restriction of ˆ σ to ˆ σ
−1
(C
1
) (which has a sole
connected component) is a surjective immersion where the number of sheets is 2 (bottom-left).
Intuitively, this means that ˆ σ “wraps” ˆ σ
−1
(C
1
) onto C
1
twice. Therefore, ˆ σ is not a simple
immersion.
and all other nodes “declined” it. It is then the task of forming the graph and finding out which
node owns each patch.
5.3 Immersion graph 75
Figure 5.7 3D example of a non-simple immersion. Left: a torus mesh which winds around
the center axis twice, forming a 3D self-connecting cell. Intersecting mesh faces are colored
red. Right: the orange curve is the skeleton of the torus mesh.
5.3.1 Simple immersions
The next question to address is then how many copies of each cell are need to make in the
resulting immersion graph. This simple question turned out more difficult to answer than
expected. Mukherjee [84] hypothesized, without proof and in two dimensions only, that the
number of copies of a cell in
ˆ
S has to equal its winding number with respect to M. This is proven
to be true (Appendix B) for immersed disks (both for d= 2 and d= 3), but for general inputs
(both for d= 2 and d= 3), it is not true. The issue are self-connecting shapes (Figure 5.6,5.7)
that contain a cell whose two (or more) node copies connect to each other. In the presence of
such cells, the winding number hypothesis no longer holds.
Self-connecting cells are addressed as follows. First, observe that, if there exists a volume
immersion ˆ σ from some
ˆ
S ontoR
d
whose boundary is
ˆ
M, then the pre-image of each cell
ˆ σ
−1
(C)⊂
ˆ
S consists of a finite number of components (proof in Appendix B). The restriction
of ˆ σ to any componentS in
ˆ
S is a surjective immersion fromS ontoC, but is not necessarily
globally injective. It can be easily shown (Appendix B) that that this surjective immersion
is what is known in algebraic topology as a covering projection ofS ontoC. It follows that
the pre-image of each individual point x∈C has the same finite cardinality for all x∈C, called
the “number of sheets” (Appendix B). Visually, this means that the immersion wrapsS ontoC
number of sheet times, e.g., 2× in Figure 5.6. Guided by the application of generating unglued
tet meshes, the immersion ˆ σ is defined to be simple if the number of sheets for each cellC
is 1, i.e., all cell immersions are globally injective and therefore bijective. A valid input is
defined to be simply volume-immersible if it is volume-immersible with a simple immersion.
The geometric interpretation of simple immersions is that they disallow self-connecting cells, by
definition. Because modeling self-connecting cells is rarely the intent in computer graphics, it
5.3 Immersion graph 76
makes sense to develop the theory and algorithms for simple immersions. Below states the first
two theorems in this chapter, establishing that the winding number intuition holds for simple
immersions, and that immersions from disks are automatically simple, both for d= 2 and d= 3.
Theorem 1: If a valid input M is simply volume-immersible, then for each cellC of the cell
complex of M, the number of connected components of ˆ σ
−1
(C) equals the winding number
wind(x,M) of any point x∈C with respect to M. The winding number wind(x,M) is the same
for all points x∈C. This holds both for d= 2 and d= 3.
Theorem 2: If a valid input M is the boundary of an immersion from a disk inR
d
, then the
immersion is simple. This holds both for d= 2 and d= 3.
These theorems are proven using algebraic topology, in Appendix B. The following Corollary
proves Mukherjee’s intuition for immersions of disks, both for d= 2 and d= 3.
Corollary 3: For a valid input M that is the boundary of an immersion from a disk inR
d
, the
number of connected components of ˆ σ
−1
(C) equals the winding number ofC with respect to
M. This holds both for d= 2 and d= 3.
Proof: This follows directly from Theorems 1 and 2. ∎
5.3.2 Rules for Graph Edges
Now proceed to defining a set of rules that further restrict the edges of eligible graphs, and
the specific B-patch ownership assignment, such that the graphs and B-patch ownerships that
satisfy them correspond to volume immersions of M. The rules are also shown in Figure 5.8.
The section starts by the assumption that M is volume-immersible, and derives the rules as
a necessary consequence of this assumption. Later the converse is proven, namely that any
eligible graph that satisfies all the rules gives a simple volume-immersion.
Theorem 4: Suppose M is simply volume-immersible. Then, the following conditions must be
satisfied for the cell graph and B-patch ownership assignment:
1. A node has at most one edge across each B-patch.
2. A node cannot have an edge across a B-patch it owns, and must have an edge across one
it declines.
3. A node from cell C must decline a B-patch whose orientation does not agree with C.
5.3 Immersion graph 77
Figure 5.8 Immersion graph connectivity rules. Solid black curves denote “owned” patches,
dashed black curves denote “declined” patches, cyan lines are connections between nodes,
and blue arrows specify patch orientations. Checkmarks and red crosses denote that a rule is
followed or violated, respectively.
Figure 5.9 Examples of Rule 6: Top-left: one common case of Rule 6, where four patches
p,q,r,s and four nodes c,d,e, f from four cells join at the same arc. Bottom-left: a rare case of
Rule 6, where c and e are the same node and p and q are the same patch. Right: examples of
such situations. Red circles denote the four-patch cases. Red rectangles denote the three-patch
cases. Note that the rule does not apply to arcs c
0
and c
1
highlighted with blue circles in the
bottom-right figure, because they both have a surrounding empty cell that does not belong to
the cell complex of M.
4. If nodes c
1
,c
2
from cells C
1
,C
2
are connected across patch q, and p
i
is a B-patch of C
i
,
and geometric neighbor of q, for i= 1,2, and p
1
, p
2
are topological neighbors, then either
both c
i
own p
i
, or both decline p
i
. Conversely, if c
i
owns p
i
, for both i= 1,2, and patches
p
1
, p
2
are topological neighbors, then c
1
and c
2
must be connected across q.
5. If a node c from cell C owns a patch p
1
, and p
1
is a geometric neighbor of p
2
on C, then
c must decline p
2
.
5.3 Immersion graph 78
6. If four nodes c,d,e, f surround the same arc A, and x connects to y across a patch that
neighbors A, for (x,y) equaling (c,d) and (d,e) and (e, f), then c must connect to f
across a patch that neighbors A. This also applies to the case where c and e are the same
node, and/or d and f are the same node. More details are in Figure 5.9.
7. Each patch is owned by exactly one node.
Proof: If Rule 1 was violated, then a B-patch would be shared by at least three nodes. For any
point x on the shared B-patch, its neighborhood includes a part of interior of all three nodes,
which contradicts local injectivity of ˆ σ. Note that Rule 1 implies that a node from a cell that
has k B-patches has at most k edges. Rules 2, 3 and 4 follow from the requirement that
ˆ
M
is the boundary of
ˆ
S, and the restriction of ˆ σ onto
ˆ
M is ˆ ρ. Note that Rule 2 establishes an
equivalence between B-patch ownership and the absence of an edge connection. Rule 5 can
be proven similarly to Rule 1. Suppose p
1
and p
2
are both owned by c, then pick a point x
at the arc shared by p
1
and p
2
. Any neighborhood of x includes a part of interior of p
1
and
p
2
. However, as p
1
and p
2
are B-patches of the same cell and are geometric neighbors, they
cannot be topological neighbors. This neighborhood of x contradicts local injectivity of ˆ σ.
Rule 6 is proved as follows. In Figure 5.9 (top-left), node c connects to d across patch p, and d
connects to e across q, and e connects to f across r. If c does not connect to f across patch s,
then by Rule 2, c must be connected to another node f
2
at the same cell as f , and the same for
f. Then, since c has declined p, by Rule 4 f
2
must decline patch r and connect to another node
e
2
from the same cell as e. By Rule 4 again e
2
should connect to another node d
2
and d
2
to c
2
.
Repeating this process either requires an infinite number of nodes, which is a contradiction, or
requires c
k
to connect back to f. This last case implies that an immersion of a disk intoR
2
is
constructed such that the disk boundary winds around the immersed disk more than once. This
contradicts the total turning number Lemma given in Appendix B. Therefore, f must connect to
c across s. Similar reasoning can be applied to the other situations of Rule 6 (e.g., Figure 5.9,
bottom-left). ∎
5.3.3 Valid Graphs Give Immersions
As suggested in Figure 5.1, finding ˆ σ involves replicating cells and connecting them so that the
ˆ
M is the immersion boundary. Next a theorem is proved that makes it possible to construct an
immersion if the graph rules from Section 5.3.2 are satisfied.
5.4 Algorithm to discover immersions 79
Theorem 5: Let M be a valid input. Assume that there is a cell graph and a B-patch ownership
assignment that satisfy Rules 1-7 given in Section 5.3.2. Then, M is simply volume-immersible,
and the graph can be used to construct the immersion.
The proof of this theorem is very technical, with details in Appendix B. The key intuition is to
employ disjoint unions of the cells belonging to each graph node, and glue the cells according
to the edges of the immersion graph. One then uses Rules 1-7 of Section 5.3.2 to prove that
that there is a manifold neighborhood of every point x. In order to do so, one has to analyze
the different cases for the position of x (interior of cell, on an owned boundary, on a declined
boundary, etc.). Now the final immersion result is stated and proven:
Corollary 6: For a valid input M, M is simply volume-immersible if and only if there exists a
cell graph and a B-patch ownership assignment that satisfy Rules 1-7 of Section 5.3.2.
Proof: This follows directly from Theorems 4 and 5. ∎
5.4 Algorithm to discover immersions
A graph G and patch ownership information are defined to be compliant if they satisfy Rules 1-7.
The immersion algorithm (Algorithm 2) searches the space of eligible graphs and ownership
assignments to discover compliant ones. Figure 5.10 shows an example of algorithm execution.
Although in principle one could enumerate all eligible graphs and ownership assignments and
test each one, here an algorithm is given that explores the space of all eligible graphs much
faster. The algorithm builds the graph one edge at a time. At each step, it examines all cells
for possible connections that could be added, and adds it to a cell where the Rules 1-7 force
a single choice. If there is more than a single choice at each cell, it then makes an arbitrary
choice and pushes the other choices onto a stack for future traversal. In many cases in practice,
Rules 1-7 are powerful enough that there is always a forced choice to be made. However, in
more challenging cases such as those in Figure 2.1, there are two or more distinct immersions,
necessitating stack use to discover all immersions.
In addition to the “owned” and “declined” ownership information, a third state,“undecided”, is
used during the course of the algorithm. Initially, all relations are initialized to “undecided”; at
the end of the algorithm they must all be either “owned” or “declined”. In addition, a declined
B-patch p on node c is tagged as incomplete if there is no edge from c across p. Note that by
5.4 Algorithm to discover immersions 80
ALGORITHM 2: Test if M is simply volume-immersible; if yes, construct the immersion
graph and patch ownerships.
FunctionisVolumeImmersible(M)
generate cell complex for M ;
create nodes of graph G at each cell, #nodes = winding number of cell ;
edges(G)←∅ ; pick any node a, make it own any B-patch p;
S← geometric neighbors of p on a + available connections;
stack← empty ; push G, patch ownerships P, and S ontostack;
whilestack not empty do
pop G,P,S fromstack;
while exist incomplete B-patches in S do
for each incompl. B-patch at node n, cell C, B-patch p in S do
connectableNodes← 0 ; D← cell that neighbors C at p ;
if #isolated nodes of D in S≥ 1 then
connectableNodes← 1 // always connectable;
nodeToConnectTo← any isolated node of D in S ;
for each non-isolated node m of D in S do
status =connectNodes(n,m) ;
if status ==success then
reset the effect ofconnectNodes(n,m) ;
connectableNodes ++ ;nodeToConnectTo← m ;
end
ambiguous← (connectableNodes≥ 2 ) ;
ifambiguous then continue // next incomplete B-patch;
else break // connect non-ambiguously;
end
if ambiguous then push G,P, dropLastEntry(S) ontostack;
status =connectNodes(n,nodeToConnectTo) ;
S← current incomplete B-patches + available connections ;
if status ==failure then returnfailure;
end
end
return G and P // success ;
end
5.4 Algorithm to discover immersions 81
ALGORITHM 3: Additional functions for constructing the immersion graph.
FunctionconnectNodes(a, b, p) // connect nodes a and b across patch p
set a and b to decline p ;updateNode(A) ;updateNode(B) ;
for each non-connected node pair (c,d) in G do // apply Rule 4
if their owned patches q,r are topological neighbors then
find patch s between c and d that shares arc with q and r ;
connectNodes(c, d, s) ;
end
for each node set c,d,e, f around an arc where Rule 6 applies do
connect c and f by callingconnectNodes() ;
if any function called above failed then returnfailure;
returnsuccess;
end
FunctionupdateNode(a) // update patch ownership for a and its neighbors
for each node b connecting to a do
apply Rule 4 to update a’s patch ownership based on b ;
apply Rule 5 to decline patches of a ;
for each owned patch p of a do
for each node b having an undecided patch p do
set b to decline p ;updateNode(b) // apply Rule 7 ;
end
if a’s patch ownership has been updated in this function then
for each node b connecting to a doupdateNode(b) // apply Rule 4;
if Rule 3, 4, 5 broken at a, or a function above failed then returnfailure;
returnsuccess;
end
5.4 Algorithm to discover immersions 82
Figure 5.10 An example of the execution of the immersion algorithm to construct an
immersion. The input is the 2D curve shown in the bottom-right of Figure 5.9. Red, solid black
and dashed black curves denote “undecided”, “owned” and “declined” patches, respectively. At
each iteration, a blue line represents the first valid connection between two nodes, green lines
are the subsequent connections made by Rule 6 after the first connection, and cyan lines are
previous connections. This is a complete example showing all 11 iterations of the algorithm to
produce the output immersion graph (bottom-right).
Rule 2, there cannot be incomplete patches; hence these states are transient and must disappear
by the end of the algorithm.
To initialize, first pick any cell that has at least one B-patch whose orientation agrees with
the cell and let the node own that B-patch. At each iteration, first check whether the graph
is compliant by checking whether there are incomplete B-patches. If not compliant, try to
find a node to connect across an incomplete B-patch. For any candidate, use the function
connectNodes to check whether the connection satisfies Rules 1-7. This function first updates
the patch ownership of the two nodes, then recursively propagates the change to other nodes
as needed. Whenever any of the Rules becomes un-satisfiable during the process, the function
returns failure. If connectNodes returns success on connecting node a to b, and there is no
ambiguity where there is another node that can also be connected to a successfully to create a
different graph, then connect a to b. Continue building the graph until encountering an ambiguity
that cannot be avoided. Then, push the current graph and patch ownership on a stack, pick one
5.4 Algorithm to discover immersions 83
possible connection, and continue. When done with the algorithm, pop the graph and patch
ownership from stack, make another choice, and continue.
Theorem 7: For a valid input M, M is simply volume-immersible if and only if the immersion
algorithm terminates successfully, i.e., discovers a graph and patch ownership assignment that
satisfies all rules. If it terminates successfully, it also explicitly constructs the immersion.
Proof: Suppose the volume immersion exists. Then, by Corollary 6, there exists a cell graph
and patch ownership assignment with the stated properties. Because at each step, the algorithm
considers all possible steps that do not violate the rules, the algorithm tests all possible graphs.
Therefore, it will discover the compliant cell graph and patch ownership assignment, i.e.,
terminate successfully. By Theorem 5, the graph then gives the immersion. Suppose the
algorithm reports that there is no graph and ownership assignment satisfying the rules. Then, by
Corollary 6, M is not simply volume-immersible. ∎
5.4.1 Running time
For d= 2, it has been proven by Eppstein and Mumford [33] that the problem of determining
if the given planar curve is a boundary of an immersed surface, is NP-complete. The size of
the problem is measured as a number of intersections (arcs in 2D), or equivalently, patches or
cells. In Appendix B, it is proven that the problem is NP-complete also for d= 3, and therefore
the running time of the algorithm has to be theoretically non-polynomial (unless P= NP). In
practice, however, the running time is very small even for very complex examples with high
winding numbers and genus. It is always dominated by the meshing time (Table 5.2).
5.4.2 Extensions
The immersion algorithm can be extended to non-connected inputs M consisting of multiple
connected components. This can be achieved by running the algorithm as usual. When the
algorithm runs out of valid connections to make, but there are still undecided patches at nodes,
then pick one undecided patch and assign it to a node, and continue the algorithm.
The algorithm above can be extended to handle non-simple-immersions that have self-connecting
cells (Figure 5.6). Such immersions require connecting nodes of the graph to nodes that are
5.4 Algorithm to discover immersions 84
Figure 5.11 Resolving the double-loop torus. A cube is added to cut the self-connecting cell.
Left to right column-wise: the model (collided faces in red) and the tetrahedral mesh, the helper
cube pulled away to show the torus mesh, volumetric ARAP deformation exposing the topology,
the double-loop torus further deformed into a single-loop.
copies of the same cell. They have limited applicability in computer graphics practice, but can
introduce an arbitrary complex subgraphs on nodes from the same cell. Designing an algorithm
to connect those nodes is quite daunting. Instead, a simple workaround is adopted here: cut the
self-connecting cells. This is done by forming a box whereby one of the six faces intersects the
self-connecting cell. Modify the input to be M plus the box. Note that no CSG is performed
here; the input consists of two separate meshes, namely M and the box. Then run the immersion
algorithm. During the cell complex creation, the original self-connecting cell will disappear
because it was cut by a box face. The output tet mesh will have two components: one for M
which should be kept, and the other for the box mesh which to be discard (see Figure 5.11). If
there are still self-connecting cells remaining, add more boxes for more cuts. Always place
the box in such a way that it intersects the self-connecting cell along its longest axis. Because
the inputs are meshes with a finite number of line segments or triangles, It will, in a finite
number of cuts, reach a situation where the resulting cells are not self-connecting any more.
Note that by Corollary 3, inputs M that are topologically spheres are guaranteed not to produce
self-connecting cells and do not require this extension. Self-connecting cells have a very small
chance to exist. Even for inputs in the examples used in this chapter that are topologically tori
with k≥ 1 handles, there is no self-connecting cells encountered, unless it is created on purpose.
5.4 Algorithm to discover immersions 85
Figure 5.12 Building the region graph. Top: input geometry, the steps to build the region
graph, and the final region graph with the assigned +,− labels. Blue arrows are the triangle
normals. Black empty circles are region graph nodes. Bottom: nearly self-intersecting dragon.
Naive meshing glues the mouth and back and produces almost no motion in animation; whereas
the immersion method produces good dynamics with a coarse mesh.
Figure 5.13 Examples of geometry inside a tetrahedron. Top row: 3 pieces partition the tet
into 4 regions. The embedded triangle mesh is shown dashed. Blue arrows represent triangle
normals. Regions labeled+ are outside of the embedded triangle mesh. Regions labeled− are
inside of the triangle mesh, and are a part ofΩ. Note that the left and right image have identical
piece geometry, but opposite+,− labels, due to the opposite piece orientation. Bottom row:
corresponding region graphs. The nodes of the+ and− regions are denoted by empty blue and
solid yellow circles, respectively. Each orange arrow represents a directed edge corresponding
to a piece.
5.5 Cell Tetrahedralization 86
5.5 Cell Tetrahedralization
The boundary of the volume immersion matches the input triangle mesh. So far, no tet mesh
was needed in the previous sections; building the cell complex and the immersion algorithm
only require the input triangle mesh. This section addresses the situation where a given input
tet mesh embeds the input triangle mesh. Arbitrary tet meshes can be used, as long the volume
enclosed by the input triangle mesh is a subset of the volume covered by the tet mesh. This
makes the immersion method in this chapter versatile, as it permits to easily adjust the tet mesh
resolution. For example, such a tet mesh can be obtained by voxelizing, or by computing a
BCC lattice [81] of the input triangle mesh, and then flood-filling the interior with tets.
The difficulty here is how to “unglue” the given tet mesh to respect the immersion. For each
cell c, a submesh T
c
of T containing the tets of T that are intersecting c is generated. To avoid
gluing proximity features on the surface of c, a new method called the “virtual tets” algorithm
is used on T
c
(Section 5.5.1). Next, assign to each immersion graph node that is a copy of c a
unique copy of T
c
. Finally, connect the tet meshes based on the edges of the immersion graph,
as follows. In detail, suppose graph nodes a and b are connected by an edge across a B-patch p.
Then, fuse two identical tets from tet meshes of a and b if they both embed the same triangle of
p. Note that if the tet mesh covers a larger space than the input triangle mesh, the immersion
will not follow the boundary of the tet mesh.
5.5.1 Virtual tets
The input to the virtual tets algorithm is a closed manifold self-intersection-free triangle mesh,
and a tet mesh that covers the volumeΩ enclosed by the triangle mesh. This section gives a
novel fast algorithm to “virtualize” the tet mesh, i.e., duplicate tetrahedra that cover more than
one “local region” ofΩ. The method ensures each tet copy knows what “local region” ofΩ it is
embedding, then it connects the copied tetrahedra into one consistent global mesh according to
the connectivity ofΩ. The algorithm is similar in goal to [111], but is more than an order of
magnitude faster. Note that this idea of “duplicate and connect” appears twice in this chapter:
once in finding the volume immersion, and the second time here. Sifakis used CSG to cut the
triangle mesh with each tet to find out the local regions. However, the CSG operations are
slow, and are the bottleneck of virtualization for complex examples. The virtual tets method
uses the Sutherland-Hodgman tet-triangle clipping algorithm [116] in exact arithmetic, and
5.5 Cell Tetrahedralization 87
pseudo-normal tests to eliminate CSG, which boosts performance between 10− 30× in the
examples.
Duplicating Tetrahedra
The Sutherland-Hodgman algorithm [116] is applied to clip triangles against each tet. The
clipped triangles form connected components (called “pieces”) inside each tet. If there are no
triangles inside the tet, then this is an interior tet, and the entire tet volume is one local region in
Ω. Otherwise, the pieces partition the tet volume into disjoint regions. The regions are labeled
with a+ if they are outside ofΩ, and− if inside. The algorithm to label the regions, and identify
the pieces that form the boundary of each region, is given as follows.
Similarly to the immersion graph, a region graph is formed, where each region is a node and
each piece corresponds to an edge (Figure 5.13) joining two regions that have this piece as their
shared boundary. The region graph is connected and without cycles, i.e., a tree. The proof is
simple: if it incrementally adds pieces one by one into the tet, since each piece is a subset of a
manifold triangle mesh, it always partitions the tet volume into two volumes. Since the mesh is
self-intersection-free, a newly added piece will only subdivide one region. Therefore, if there
are k pieces in a tet, it will form k+ 1 regions. Since the triangles in each piece are oriented
consistently, one can define directions of the edges in the graph as directions from the region
inside the triangles to the outside.
The region graph is built incrementally as follows (see also Figure 5.12). First, initialize it by
placing only one piece into the tet, creating two region nodes and one piece edge. When adding
a new piece P, first identify which region will be subdivided by P. This is done by querying an
arbitrary vertex v on P against the boundary surface of every region, until discovering the one
that contains it. In Appendix C, it is proven that there is no need to explicitly form the region
boundary surface (which would require using CSG) to query this. It is sufficient to only perform
inside-outside queries on the boundary pieces of the region. It first uses the region graph to
identify the pieces that bound this region, then finds the closest site on each piece to v. Next,
perform the pseudo-normal test against the closest site [7] to identify if v is inside or outside
each piece. After finding the region ν that contains v (and therefore P), subdivideν intoν
+
and
ν
−
, split its node in the region graph into two nodes, and connect them to each other with the
piece edge of P. For a piece Q connectingν before P is added, move it to eitherν
+
orν
−
by
testing whether an arbitrary vertex w on Q is outside or inside P, using the same pseudo-normal
5.5 Cell Tetrahedralization 88
Figure 5.14 Duplication of tet vertices to embed a comb.
test as done on v. Since pieces do not intersect, a piece cannot be both inside and outside of P.
After all the pieces are added to the region graph, find all − regions by checking whether they
are inside their boundary pieces, i.e., they have only outward edges in the region graph. Then
duplicate this tet by the number of− regions and for each region, embed the triangles of its
boundary pieces into the tet for this region.
Connecting Duplicated Tetrahedra
The method connects a pair of neighboring tets if their embedded pieces share points on the
boundary. The boundaries of each piece are polygonal lines lying on tet faces. If the shared tet
face between the two tets does not have piece boundaries from either tet, and an interior point
on the face is considered inside for the pieces in both tets (again by pseudo-normal test), those
tets must also be connected. Points on faces of a tet with no pieces are considered inside in the
above test.
Once the method knows which tets should be connected, the assignment and duplication of
tet vertices is performed as explained in [121], Section 6. In short, each copied tet is initially
assigned a totally unique set of 4 vertices. If two tets are connected, merge the three pairs of
vertices on the shared tet face. The merging is done using a disjoint-set data structure. Such tet
vertex duplication produces correct topological connectivity of the tet mesh. This is illustrated
in Figure 5.14 whereby none of the teeth of the comb are welded to adjacent teeth.
5.6 Results 89
Table 5.1 Mesh complexity statistics. m.w.#: max winding number on the mesh.
Model vertices triangles genus cells patches arcs m.w.# #tets
head 3,020 6,036 0 3 3 2 2 566,515
helix 100,020 200,036 0 3 3 2 2 36,182
yeahright-on-ground 94,063 188,646 131 117 263 362 4 83,024
quintuple torus 15,211 30,438 5 314 886 1,700 7 84,905
tree 163,998 328,152 40 131 259 258 2 32,457
Sacht (small) 12,448 24,896 1 26 49 48 2 112,554
Sacht (large) 24,160 48,320 1 41 79 78 2 225,338
Table 5.2 Immersion algorithm timings on examples. t-cell, t-imm, t-mesh: times for gener-
ating the cell complex, immersion algorithm and meshing.
Model t-cell t-imm t-mesh t-total
head 0.34s 0.0090s 38.0s 38.4s
helix 56.4s 0.0003s 190s 246s
yeahright-on-ground 14.2s 0.12s 57.1s 71.4s
quintuple torus 11.6s 4.89s 35.9s 52.4s
tree 15.9s 0.19s 130s 146s
Sacht (small) 2.77s 0.0015s 19.7s 22.5s
Sacht (large) 6.01s 0.0045s 42.2s 48.2s
Details
For more detail and robustness considerations of the algorithm, please refer to Appendix C. The
tetrahedralization algorithm can be easily adapted to 2D to embed polygonal lines into triangles.
5.6 Results
The statistics of the immersion method is given in Table 5.1 and Table 5.2. (Intel Xeon 2.3GHz
2x6 cores, 32 GB memory). An illustration of the progress of the immersion algorithm is
shown in Figure 5.17. To help with implementation, the source code is released under a free
BSD license as a part of Vega FEM 4.0 [11]. The released code includes both self-intersecting
meshing (Sections 5.2, 5.3, 5.4) and nearly self-intersecting meshing (Section 5.5). The input tet
meshes in the examples are generated by red-green subdivision on tet grids [80], or voxelization
implemented in [11]. The immersion method is tested on multiple meshes. The upper and lower
lips in the head model in Figure 5.18 self-intersect. The immersion method was able to generate
an overlapping, but separated tetrahedral mesh despite the self-intersections, which makes it
5.6 Results 90
Figure 5.15 Yeahright-on-ground model. Top left: Yeahright dropped on the ground. Collided
faces are in red. Bottom left: generated tetrahedral mesh. Top right: the Yeahright-on-ground is
deformed using volumetric ARAP on the tetrahedral mesh. Bottom right: closer look on how
the immersion algorithm is able to untangle the intersecting features.
possible for FEM simulation to open the mouth. It can achieve this even with a very coarse
mesh that only has 7 tets across the length of the mouth. The helix model (Figure 5.19) consists
of 100 tightly self-overlapping loops. The method is able to mesh and animate it correctly.
The yeahright-on-ground model (Figure 5.15) is generated by modifying the “yeahright” model
from Keenan Crane’s 3D Model Repository [31]. Both the original model and the modified
model have genus 131. The modification is that the original model is dropped onto the ground
under gravity using a FEM simulation, and let it come to rest on a plane without any self-
collision handling. This resulted in a massively self-intersecting shape of genus 131. Here this
shape is called the “yeahright-on-ground”, and used as input in this example. The immersion
method was able to generate an unglued tetrahedral mesh for this severely self-intersecting
shape (Figure 5.15), which enabled us to correctly pull it apart using volumetric ARAP. The
high genus of the model causes zero-area faces in the forward flow, which causes Sacht’s
method to fail, whereas the immersion method succeeds.
5.6 Results 91
Figure 5.16 Quintuple Torus model. Top-left: side view of the Quintuple Torus. Collided
faces are in red. Bottom-left: top view with holes visible. Top-middle: Sacht’s implementation
became stuck in the reverse flow due to the complicated collisions, spending 8 hours without
progress until manually stopped. Bottom-middle: the immersion algorithm successfully gener-
ated a tetrahedral mesh to embed the shape. Top-right: simulating the output mesh using FEM.
Bottom-right: the deformed tetrahedral mesh.
Figure 5.17 Intermediate stages of building immersion graphs. The three rows correspond
to the examples of Sacht (small), yeahright-on-ground and quintuple torus. Leftmost column:
the first node added to the immersion graph. Rightmost column: the final node added.
In the quintuple torus example, 24 cylinders are attached onto a quintuple torus (genus 5) and
collided the cylinders (Figure 5.16), forming cells of a high winding number 7. It generated a
quality tetrahedral mesh that is capable of separating all of the collided cylinders.
5.6 Results 92
Figure 5.18 Unglued meshing of self-intersecting lips and mouth cavity: Column 1: side
view of the head and the self-intersecting mouth connected to an interior cavity. Such self-
intersecting modeling of the mouth cavity is very common in visual effects and games, to add a
plausible mouth appearance to characters. The self-intersecting mouth is a common technology
blocker for head simulation algorithms. Column 2: head triangle mesh (top) and the resulting
tetrahedral mesh (bottom). Column 3: zoom to the self-intersecting lips. Column 4: the
tetrahedral mesh “unglued” the input self-intersection, enabling the simulation to successfully
open the mouth.
In the tree example (Figure 5.20), the immersion method is applied to procedurally generated
geometry, and in the dragon example (Figure 5.12), the method processes nearly self-intersecting
geometry. In addition to direct modeling and animation, the immersion method also enables
topology-aware computation of volumetric weights such as bounded biharmonic weights [57]
(Figure 5.21), suitable for further use in shape deformation and simulation.
In Figure 5.22, The immersion method is compared to Sacht’s method. The immersion method
were successful on all of the non-inverted inputs of [97]. It can accommodate non-spherical
inputs (genus≥ 1), whereas Sacht’s method often fails on such inputs; for example, it fails to
produce any meaningful result on the yeahright-on-ground and quintuple torus models. The
torus example that Sacht’s work listed as problematic for their method (“Sacht large” example
in Table 5.1) is also conducted, and confirmed that Sacht’s method fails whereas the immersion
method produces a good result. Then a follow-up experiment decreased the difficulty of
5.6 Results 93
Figure 5.19 Self-intersection-aware tet meshing: The immersion method properly duplicated
and connected the tets to account for helix triangle mesh self-intersections. Left: input triangle
and tetrahedral mesh. Intersecting triangles are shown red. Right: the output tet mesh has been
successfully pulled apart using volumetric ARAP [114].
Figure 5.20 Tetrahedral meshing of self-intersecting procedurally generated geometry.
The triangle mesh for this tree (328,152 triangles; 2,634 branches) was generated using a
procedural modeling method without regard for branch intersections and self-intersections. The
immersion method successfully generated a simulation tet mesh that “unglues” all the branches.
Because our tets can be overlapping, our mesh can be relatively coarse (32,457) for a model
of this complexity, permitting faster FEM simulation. For comparison, naive meshing glues
branches and results in visibly suboptimal tree motion when simulated using FEM.
the example by removing some vertical columns (“Sacht small”), upon which Sacht’s work
succeeded, but, due to meshing in the unwrapped configuration, produces a very suboptimal tet
mesh compared to the world-space meshing (Figure 5.22, B). Furthermore, because Sacht’s
5.6 Results 94
Figure 5.21 The immersion method produces topology-aware weights. Bounded biharmonic
weights[57] are computed using the two handles “A” and “B” indicated; the figure (left) shows
the BBW weight function for handle “A”. BBWs computed on our volumetric mesh (left)
reflects the correct topology of the embedded surface (we re-use the Sacht’s test model [ 97]
from Figure 5.22). In contrast, naive voxelization (right) glues the cylinders together, resulting
in unsatisfactory weights.
Figure 5.22 Comparison to Sacht’s work: Part (A): Left to right: side and top view of the
model from [97] (intersections are in red), Sacht’s method result (failed in the middle of
the process), the tetrahedral mesh generated successfully by the immersion method, and the
tetrahedral mesh and its embedded mesh deformed by volumetric ARAP, demonstrating the
correct topology. Part (B): Left to right: side and top view of the “small” input model; Sacht’s
method produces a highly deformed shape that results in a low-quality tetrahedral mesh for the
original surface; the immersion method produces a good volumetric mesh.
5.6 Results 95
method requires repeated self-collision detection at each iteration, whereas the immersion
method only need it once in building the cell complex and once for tetrahedralization, The
immersion method is more than 100× faster than Sacht’s method.
The virtual tets method (Section 5.5) is compared with the CSG method of [111] on the
yeahright-on-ground model (Figure 5.15). Note that prior work [111] did not address self-
intersecting inputs (only nearly self-intersecting), but the prior work hereby is used to replace the
tetrahedralization stage. It turns out that the prior work [111] needed 2,545 seconds total, with
more then 98% of the time spent in running CSG operations to duplicate tetrahedra. The same
CSG routines ([142]) are used when implementing both the immersion method and the prior
work. The profiling shows that, in prior work, 97% of the total time is spent in exact-arithmetic
mesh booleans. The immersion method does not need mesh booleans for tetrahedralization,
and requires a lot fewer exact arithmetic operations (Table 5.3). Similar speedups are observed
on smaller examples.
Table 5.3 Timing comparison to Sifakis’s algorithm. The total time spent in four key CGAL
exact arithmetic routines when using CSG operations as in [111], versus the novel pseudo-
normal and exact arithmetic Sutherland-Hodgman algorithm. Yeah-right-on ground model.
“Meshing total” is the time to generate the tet mesh after the immersion algorithm has gener-
ated the graph, including CGAL and other operations specific to each method. The “grand
total” is the total time from loading mesh M to producing the output tet mesh; it equals the
immersion algorithm time (16 sec; same for both methods), plus the “meshing total” time. The
libigl::mesh_boolean times are for tetrahedralization in Sifakis’s method; the virtual tets method
does not need it for tetrahedralization. it is used in the immersion algorithm; where it occupies
12.5 sec of the running time.
Routine Sifakis 2007 [sec] Immersion [sec] Speedup
CGAL::Lazy_exact_Add 7.6 0.14 54×
CGAL::Lazy_exact_Sub 2.5 0.34 7.3×
CGAL::Lazy_exact_Mul 7.6 0.14 54×
CGAL::Lazy_rep_2 305 13.6 22.4×
libigl::mesh_boolean 2479 0 ∞
Meshing total 2529 55 50×
Grand total 2545 71 36×
Although the immersion theorems apply to shapes that are free of self-connecting shapes, it can
still accommodate such shapes (Figure 5.6) by cutting them with helper cubes (Figure 5.11).
This increases the number of cells of M from 5 only to 11, and converts the input into one without
self-connecting shapes. The method is then able to produce a valid separating tetrahedral mesh.
This tetrahedral mesh has two connected components: one for the cube (which is discarded),
5.6 Results 96
and one for the torus. In Figure 5.11, right, it is shown that it can then use the produced torus tet
mesh to deform the double-loop torus into a self-intersection-free configuration. The immersion
method is able to find all possible immersions (Figure 2.1). Unlike prior work [ 82], the input is
not restricted to a single connected component (see Figure 2.2), which means that the immersion
method accepts a broader input than the state of the art on 2D immersions.
Chapter 6
Conclusion
This thesis provided several methods to address collisions and large nonlinear elasticity in
computer modeling and animation pipelines. Given a polygon mesh, a tetrahedral mesh is
first created. Then in the modeling stage, our reduced hierarchical method provides interactive
shape editing along with local contact handling. Local bases from the hierarchy are added to
the system at runtime to resolve contact and users can control the tradeoff between accuracy
and speed. In the animation stage, we demonstrate how to fits a tetrahedral mesh animation
to the given polygon mesh animation, and then add collision handling and secondary motion
using physically based simulation with constraints. Two blending methods are also given to
adjust how much the new animation should follow the input.
Finally, a mathematical analysis of self-intersecting polygon meshes is presented to deal with
nearly self-intersecting, or intersecting, input meshes. By studying the topological concept
of immersions, we gave a deterministic and constructive algorithm to determine if the input
self-intersecting mesh is the boundary of an immersion. Our algorithm is able to construct
topologically-consistent self-intersecting tetrahedral meshes at arbitrary resolutions. For near
self-intersections, a robust method is introduced to properly duplicate mesh elements and
correctly embed the underlying geometry into the mesh element copies. With this algorithm,
even self-intersecting polygon meshes can be used as rest shapes for subsequent modeling and
animating stages with proper collision handling and physically based simulation.
6.1 Limitations and Future Work 98
6.1 Limitations and Future Work
The hierarchical method presented in the thesis is limited to vertex handles as user inputs.
Use of rigid regions, cages, skeletons and other control methods are left for future work. Our
energy operates on a three-dimensional solid tet mesh. It can be extended to shells and rods.
The method uses standard vertex weights for the ARAP energy. It would be possible to apply
different weights to model heterogeneous objects; the specifics are left for future work. Discrete
self collision detection is used on vertex-tet pairs, and therefore may miss collisions on thin
features. Collision detection on edge-edge pairs [119] would improve contact handling, at the
cost of additional computation. Continuous collision detection on triangles could be added at an
additional computational cost. GPU-based collision detection methods are also useful. While
the basis vectors can be dynamically added or removed at runtime, the method does not support
cutting or modifying the shape of the basis vectors at runtime. Changing the working basis at
runtime can create minor visual popping; this can be visually alleviated by temporally blending
the output shape. Much like all contact resolution schemes, the performance degrades when
large parts of the model undergo contact. In addition, collisions are generally less stable on
concave surfaces than on convex ones. Although the static damping improves stability, it adds
one extra dimension for parameter tuning and affects the shape globally. An interesting future
extension is local static damping for each unstable contact site.
The animation system described in this thesis has the limitation that the input animation must
be provided at every frame; completely missing frames would require solving space-time
optimization problems. The system also requires that at least a (small) part of the model has
to follow the input exactly. While the system is not tied to a particular interpolation method,
barycentric coordinates are only C
0
at the tetrahedral boundaries. Smoother interpolation
methods such as mean value coordinates [62] could be used instead. The system cannot
generate deformations smaller than the volumetric mesh resolution. Therefore, it can be hard to
create wrinkles and other small details. If the input triangle mesh animations are unreasonable
or chaotic, the fitting process may produce severely strained volumetric mesh animations, which
in turn will make the simulation unstable. Adding anatomical detail to our solid models would
improve animation realism.
Topologically-correct tetrahedralization is the most costly component in the immersion algo-
rithm, and it can be accelerated in the future with better engineering. The immersion method
is used only to build tetrahedral meshes, but the technique can be easily extended to other
6.1 Limitations and Future Work 99
elements such as hexahedral meshing. There are several remaining open challenges; namely,
how to detect and properly handle non-simple immersions such as multiple self-connecting tori,
immersions where the compactness assumption does not hold, and non-immersions such as
inversions. Challenges also include handling input meshes with boundaries, and constructing
immersions in higher dimensions. The surfaces of the resulting tetrahedral meshes do not match
the input triangle meshes. An interesting future work is to tetrahedralize a constrained boundary
that self-intersects.
Acknowledgements
This research was sponsored in part by NSF (CAREER-1055035, IIS-1422869), the Sloan
Foundation, the Okawa Foundation, USC Annenberg Fellowships to Yijing Li and Hongyi Xu,
Bosch Research and Adobe Research.
References
[1] Alexa, M., Cohen-Or, D., and Levin, D. (2000). As-rigid-as-possible shape interpolation.
In Proc. of ACM SIGGRAPH 2000, pages 157–164.
[2] Allard, J., Faure, F., Courtecuisse, H., Falipou, F., Duriez, C., and Kry, P. G. (2010). V olume
contact constraints at arbitrary resolution. In ACM SIGGRAPH 2010 papers, SIGGRAPH
’10, pages 82:1–82:10, New York, NY , USA. ACM.
[3] Alliez, P., Cohen-Steiner, D., Yvinec, M., and Desbrun, M. (2005). Variational tetrahedral
meshing. In ACM Transactions on Graphics (TOG), volume 24, pages 617–625. ACM.
[4] An, S. S., Kim, T., and James, D. L. (2008). Optimizing cubature for efficient integration
of subspace deformations. ACM Trans. on Graphics (SIGGRAPH Asia 2008), 27(5):165:1–
165:10.
[5] Attene, M. (2010). A lightweight approach to repairing digitized polygon meshes. The
Visual Computer, 26(11):1393–1406.
[6] Attene, M. (2014). Direct repair of self-intersecting meshes. Graphical Models, 76(6):658–
668.
[7] Bærentzen, J. and Aanæs, H. (2005). Signed distance computation using the angle weighted
pseudo-normal. IEEE Trans. on Visualization and Computer Graphics, 11(3):243–253.
[8] Baraff, D. and Witkin, A. P. (1998). Large Steps in Cloth Simulation. In Proc. of ACM
SIGGRAPH 98, pages 43–54.
[9] Barbiˇ c, J. and James, D. L. (2008). Six-dof haptic rendering of contact between geometri-
cally complex reduced deformable models. IEEE Transactions on Haptics, 1(1):39–52.
[10] Barbiˇ c, J., da Silva, M., and Popovi´ c, J. (2009). Deformable object animation using
reduced optimal control. ACM Trans. on Graphics, 28(3).
[11] Barbiˇ c, J., Li, Y ., Wang, B., and Zhao, D. (2018). Vega FEM Library 4.0.
http://www.jernejbarbic.com/vega.
[12] Barbiˇ c, J. and Popovi´ c, J. (2008). Real-time control of physically based simulations using
gentle forces. ACM Trans. on Graphics (SIGGRAPH Asia 2008), 27(5):163:1–163:10.
References 101
[13] Barbiˇ c, J., Sin, F., and Grinspun, E. (2012). Interactive editing of deformable simulations.
ACM Trans. on Graphics (SIGGRAPH 2012), 31(4):70:1–70:8.
[14] Bergou, M., Mathur, S., Wardetzky, M., and Grinspun, E. (2007). TRACKS: Toward
directable thin shells. ACM Trans. on Graphics (SIGGRAPH 2007), 26(3):50:1–50:10.
[15] Bickel, B., Baecher, M., Otaduy, M., Matusik, W., Pfister, H., and Gross, M. (2009).
Capture and modeling of non-linear heterogeneous soft tissue. ACM Trans. on Graphics
(SIGGRAPH 2009), 28(3):89:1–89:9.
[16] Boier-Martin, I., Zorin, D., and Bernardini, F. (2005). A survey of subdivision-based tools
for surface modeling. DIMACS Series in Discrete Math and Theoretical CS, 67:1.
[17] Boissonnat, J.-D. and Oudot, S. (2005a). Provably good sampling and meshing of surfaces.
Graphical Models, 67(5):405–451.
[18] Boissonnat, J.-D. and Oudot, S. (2005b). Provably good sampling and meshing of surfaces.
Graphical Models, 67(5):405–451.
[19] Botsch, M., Pauly, M., Gross, M., and Kobbelt, L. (2006). PriMo: Coupled Prisms for
Intuitive Surface Modeling. In Eurographics Symp. on Geometry Processing, pages 11–20.
[20] Bouaziz, S., Martin, S., Liu, T., Kavan, L., and Pauly, M. (2014). Projective dynamics:
fusing constraint projections for fast simulation. ACM Transactions on Graphics (TOG),
33(4):154.
[21] Brahana, H. R. (1921). Systems of circuits on two-dimensional manifolds. Annals of
Mathematics, pages 144–168.
[22] Brandt, C., Eisemann, E., and Hildebrandt, K. (2018). Hyper-reduced projective dynamics.
ACM Transactions on Graphics (TOG), 37(4):80.
[23] Bridson, R., Fedkiw, R., and Anderson, J. (2002). Robust Treatment of Collisions, Contact,
and Friction for Cloth Animation. ACM Trans. on Graphics, 21(3):594–603.
[24] Campen, M. and Kobbelt, L. (2010). Exact and robust (self-) intersections for polygonal
meshes. In Computer Graphics Forum, volume 29, pages 397–406.
[25] Capell, S., Burkhart, M., Curless, B., Duchamp, T., and Popovi´ c, Z. (2005). Physically
based rigging for deformable characters. In Symp. on Computer Animation (SCA), pages
301–310.
[26] Capell, S., Green, S., Curless, B., Duchamp, T., and Popovi´ c, Z. (2002). A Multiresolution
Framework for Dynamic Deformations. In Proc. of the Symp. on Comp. Animation 2002,
pages 41–48.
[27] CGAL (2018). Computational Geometry Algorithms Library.
[28] Chen, L., Huang, J., Sun, H., and Bao, H. (2010). Cage-based deformation transfer.
Computers & Graphics, 34(2):107–118.
References 102
[29] Choi, J. and Szymczak, A. (2009). Fitting solid meshes to animated surfaces using linear
elasticity. ACM Trans. on Graphics (TOG), 28(1):6:1–6:10.
[30] Coros, S., Martin, S., Thomaszewski, B., Schumacher, C., Sumner, R., and Gross, M.
(2012). Deformable objects alive! ACM Trans. on Graphics (SIGGRAPH 2012), 31(4):69.
[31] Crane, K. (2017). 3D Model Repository.
[32] Debunne, G., Desbrun, M., Cani, M.-P., and Barr, A. H. (2001). Dynamic Real-Time
Deformations Using Space & Time Adaptive Sampling. In Proc. of ACM SIGGRAPH 2001,
pages 31–36.
[33] Eppstein, D. and Mumford, E. (2009). Self-overlapping curves revisited. In Proc. of
ACM-SIAM Symp. on Discrete Algorithms, pages 160–169.
[34] Erickson, J. (2009). Computational topology; cell complexes. In Course Notes, University
of Illinois at Urbana-Champaign.
[35] Ericson, C. (2004). Real-time collision detection. CRC Press.
[36] Erleben, K. (2018). Methodology for assessing mesh-based contact point methods. ACM
Transactions on Graphics (TOG), 37(3):39.
[37] Faure, F., Gilles, B., Bousquet, G., and Pai, D. K. (2011). Sparse meshless models of
complex deformable solids. In ACM Trans. on Graphics (SIGGRAPH 2011), volume 30,
page 73.
[38] Frisch, D. (2010). Extending immersions into the sphere. arXiv preprint arXiv:1012.4923.
[39] Fröhlich, S. and Botsch, M. (2011). Example-driven deformations based on discrete shells.
In Computer Graphics Forum, volume 30, pages 2246–2257.
[40] Gain, J. E. and Dodgson, N. A. (2001). Preventing self-intersection under free-form
deformation. IEEE Transactions on Visualization and Computer Graphics, 7(4):289–298.
[41] Galoppo, N., Otaduy, M. A., Tekin, S., Gross, M., and Lin, M. C. (2007). Soft articulated
characters with fast contact handling. In Computer Graphics Forum, volume 26, pages
243–253. Wiley Online Library.
[42] Glickenstein, D. (2018). Covering space. Technical report, University of Arizona.
[43] Gottschalk, S., Lin, M. C., and Manocha, D. (1996). OBBTree: A Hierarchical Structure
for Rapid Interference Detection. In Proc. of ACM SIGGRAPH 96, pages 171–180.
[44] Grinspun, E., Krysl, P., and Schröder, P. (2002). CHARMS: A Simple Framework for
Adaptive Simulation. ACM Trans. on Graphics, 21(3):281–290.
[45] Guennebaud, G., Jacob, B., et al. (2010). Eigen v3. http://eigen.tuxfamily.org.
References 103
[46] Hahn, F., Martin, S., Thomaszewski, B., Sumner, R., Coros, S., and Gross, M. (2012).
Rig-space physics. ACM Transactions on Graphics (TOG), 31(4):72.
[47] Hahn, F., Thomaszewski, B., Coros, S., Sumner, R., and Gross, M. (2013). Efficient
simulation of secondary motion in rig-space. In Symp. on Computer Animation (SCA), pages
165–171.
[48] Hang Si (2011). TetGen: A Quality Tetrahedral Mesh Generator and a 3D Delaunay
Triangulator.
[49] Harmon, D., Panozzo, D., Sorkine, O., and Zorin, D. (2011). Interference aware geometric
modeling. ACM Trans. on Graphics (SIGGRAPH Asia 2011), 30(6):137:1–137:10.
[50] Harmon, D. and Zorin, D. (2013). Subspace integration with local deformations. ACM
Trans. on Graphics (TOG), 32(4):107.
[51] Hétroy, F., Rey, S., Andújar, C., Brunet, P., and Vinacua, A. (2011). Mesh repair with
user-friendly topology control. Computer-Aided Design, 43(1):101–113.
[52] Hildebrandt, K., Schulz, C., von Tycowicz, C., and Polthier, K. (2012). Interactive
spacetime control of deformable objects. ACM Trans. on Graphics (SIGGRAPH 2012),
31(4):71:1–71:8.
[53] Hu, Y ., Zhou, Q., Gao, X., Jacobson, A., Zorin, D., and Panozzo, D. (2018). Tetrahedral
meshing in the wild. ACM Transactions on Graphics (TOG), 37(4):60.
[54] Huang, J., Tong, Y ., Zhou, K., Bao, H., and Desbrun, M. (2011). Interactive shape
interpolation through controllable dynamic deformation. IEEE Trans. on Visualization and
Computer Graphics, 17(7):983–992.
[55] Irving, G., Teran, J., and Fedkiw, R. (2004). Invertible Finite Elements for Robust
Simulation of Large Deformation. In Symp. on Computer Animation (SCA), pages 131–140.
[56] Jacobson, A., Baran, I., Kavan, L., Popovi´ c, J., and Sorkine, O. (2012). Fast automatic
skinning transformations. ACM Trans. on Graphics (SIGGRAPH 2012), 31(4):77:1–77:10.
[57] Jacobson, A., Baran, I., Popovi´ c, J., and Sorkine, O. (2011). Bounded biharmonic weights
for real-time deformation. ACM Trans. Graph., 30(4):78:1–78:8.
[58] Jacobson, A., Kavan, L., and Sorkine-Hornung, O. (2013a). Robust inside-outside
segmentation using generalized winding numbers. ACM Transactions on Graphics (TOG),
32(4):33.
[59] Jacobson, A., Panozzo, D., Schüller, C., Diamanti, O., Zhou, Q., Pietroni, N., et al. (2013b).
libigl: A simple c++ geometry processing library.
[60] James, D. L. and Pai, D. K. (2003). Multiresolution Green’s Function Methods for Interac-
tive Simulation of Large-scale Elastostatic Objects. ACM Trans. on Graphics, 22(1):47–82.
References 104
[61] Joshi, P., Meyer, M., DeRose, T., Green, B., and Sanocki, T. (2007). Harmonic coordinates
for character articulation. In ACM Trans. on Graphics (SIGGRAPH 2007), volume 26, pages
71:1–71:9.
[62] Ju, T., Schaefer, S., and Warren, J. (2005). Mean value coordinates for closed triangular
meshes. In ACM Trans. on Graphics (SIGGRAPH 2005), volume 24, pages 561–566.
[63] Kass, M. and Anderson, J. (2008). Animating oscillatory motion with overlap: Wiggly
splines. ACM Trans. on Graphics (SIGGRAPH 2008), 27(3):28:1–28:8.
[64] Kaufman, D. M., Sueda, S., James, D. L., and Pai, D. K. (2008). Staggered Projections
for Frictional Contact in Multibody Systems. ACM Transactions on Graphics, 27(5):164:1–
164:11.
[65] Kim, J. and Pollard, N. S. (2011). Fast simulation of skeleton-driven deformable body
characters. ACM Trans. on Graphics (TOG), 30(5):121.
[66] Kim, T. and James, D. L. (2012). Physics-based character skinning using multidomain
subspace deformations. IEEE Trans. on Visualization and Computer Graphics, 18(8):1228–
1240.
[67] Kondo, R., Kanai, T., and ichi Anjyo, K. (2005). Directable animation of elastic objects.
In Symp. on Computer Animation (SCA), pages 127–134.
[68] Li, S., Huang, J., de Goes, F., Jin, X., Bao, H., and Desbrun, M. (2014). Space-time editing
of elastic motion through material optimization and reduction. ACM Trans. on Graphics
(SIGGRAPH 2014), 33(4):108:1–108:10.
[69] Li, S., Huang, J., Desbrun, M., and Jin, X. (2013). Interactive elastic motion editing
through space–time position constraints. Computer Animation and Virtual Worlds, 24(3-
4):409–417.
[70] Li, W. (2011). Detecting ambiguities in 3d polygons with self-intersecting projections. In
Computer-Aided Design and Computer Graphics CAD/Graphics, pages 11–16.
[71] Lin, M. C. and Otaduy, M. (2008). Haptic rendering: foundations, algorithms, and
applications. CRC Press.
[72] Liu, L., Yin, K., Wang, B., and Guo, B. (2013). Simulation and control of skeleton-driven
soft body characters. ACM Trans. on Graphics (SIGGRAPH Asia 2013), 32(6):215.
[73] Liu, T., Bouaziz, S., and Kavan, L. (2016). Towards real-time simulation of hyperelastic
materials. arXiv preprint arXiv:1604.07378.
[74] Lu, H., Xian, C., Li, G., and Zhang, Z. (2014). EC-CageR: error controllable cage reverse
for animated meshes. Computers & Graphics, 46:138—-148.
[75] Malgat, R., Gilles, B., Levin, D. I., Nesme, M., and Faure, F. (2015). Multifarious
hierarchies of mechanical models for artist assigned levels-of-detail. In Symp. on Computer
Animation (SCA).
References 105
[76] Manson, J. and Schaefer, S. (2011). Hierarchical deformation of locally rigid meshes. In
Computer Graphics Forum, volume 30, pages 2387–2396.
[77] McAdams, A., Zhu, Y ., Selle, A., Empey, M., Tamstorf, R., Teran, J., and Sifakis, E.
(2011). Efficient elasticity for character skinning with contact and collisions. ACM Trans.
on Graphics (SIGGRAPH 2011), 30(4):37:1–37:11.
[78] Mitchell, N., Aanjaneya, M., Setaluri, R., and Sifakis, E. (2015). Non-manifold level sets:
A multivalued implicit surface representation with applications to self-collision processing.
ACM Transactions on Graphics (TOG), 34(6):247.
[79] Molino, N., Bao, Z., and Fedkiw, R. (2004). A virtual node algorithm for changing mesh
topology during simulation. In Proc. of ACM SIGGRAPH 2004, pages 385–392.
[80] Molino, N., Bridson, R., and Fedkiw, R. (2003a). Tetrahedral mesh generation for
deformable bodies. In Proc. Symposium on Computer Animation.
[81] Molino, N., Bridson, R., Teran, J., and Fedkiw, R. (2003b). A crystalline, red green strategy
for meshing highly deformable objects with tetrahedra. In 12th Int. Meshing Roundtable,
pages 103–114.
[82] Mukherjee, U. (2014). Self-overlapping curves: Analysis and applications. Computer-
Aided Design, 46:227–232.
[83] Mukherjee, U. and Gopi, M. (2012). Tweening boundary curves of non-simple immersions
of a disk. In Proc. of Indian Conference on Computer Vision, Graphics and Image Processing,
page 36. ACM.
[84] Mukherjee, U., Gopi, M., and Rossignac, J. (2011). Immersion and embedding of self-
crossing loops. In Proc. of the Eurographics Symposium on Sketch-Based Interfaces and
Modeling, pages 31–38.
[85] Müller, M. and Chentanez, N. (2010). Wrinkle meshes. In Symp. on Computer Animation
(SCA), pages 85–92.
[86] Müller, M. and Gross, M. (2004). Interactive Virtual Materials. In Proc. of Graphics
Interface 2004, pages 239–246.
[87] Munkres, J. (2000). Topology. Featured Titles for Topology Series. Prentice Hall, Incor-
porated.
[88] Narain, R., Overby, M., and Brown, G. E. (2016). ADMM⊇ projective dynamics: fast
simulation of general constitutive models. In Symposium on Computer Animation, pages
21–28.
[89] Nesme, M., Kry, P. G., Jeˇ rábková, L., and Faure, F. (2009). Preserving topology and
elasticity for embedded deformable models. ACM Trans. on Graphics (SIGGRAPH 2009),
28(3):52:1–52:9.
References 106
[90] Nocedal, J. and Wright, S. (2006). Numerical optimization. Springer Science & Business
Media.
[91] Noe, K. Ø. and Sørensen, T. S. (2010). Solid mesh registration for radiotherapy treatment
planning. In Biomedical Simulation, volume 5958, pages 59–70. Springer.
[92] Otaduy, M. A., Germann, D., Redon, S., and Gross, M. (2007). Adaptive Deformations
with Fast Tight Bounds. In Symp. on Computer Animation (SCA), pages 181–190.
[93] Otaduy, M. A. and Lin, M. C. (2006). A modular haptic rendering algorithm for stable
and transparent 6-dof manipulation. IEEE Transactions on Robotics, 22(4):751–762.
[94] Paillé, G.-P., Poulin, P., and Lévy, B. (2013). Fitting polynomial volumes to surface meshes
with voronoï squared distance minimization. In Computer Graphics Forum, volume 32,
pages 103–112.
[95] Press, W., Teukolsky, S., Vetterling, W., and Flannery, B. (2007). Numerical recipes: The
art of scientific computing . Cambridge University Press, Cambridge, UK, third edition.
[96] Rohmer, D., Popa, T., Cani, M.-P., Hahmann, S., and Sheffer, A. (2010). Animation
wrinkling: augmenting coarse cloth simulations with realistic-looking wrinkles. In ACM
Trans. on Graphics (SIGGRAPH Asia 2010), volume 29, page 157.
[97] Sacht, L., Jacobson, A., Panozzo, D., Schüller, C., and Sorkine-Hornung, O. (2013).
Consistent volumetric discretizations inside self-intersecting surfaces. In Computer Graphics
Forum, volume 32, pages 147–156. Wiley Online Library.
[98] Sacht, L., V ouga, E., and Jacobson, A. (2015). Nested cages. ACM Trans. on Graphics
(SIGGRAPH Asia 2015), 34(6):170.
[99] Savoye, Y . and Franco, J.-S. (2010a). CageIK: dual-Laplacian cage-based inverse kine-
matics. In Articulated Motion and Deformable Objects, volume 6169, pages 280–289.
Springer.
[100] Savoye, Y . and Franco, J.-S. (2010b). Conversion of performance mesh animation into
cage-based animation. In ACM SIGGRAPH ASIA 2010 Posters, page 2. ACM.
[101] Schoberl, J. (1997). NETGEN - An advancing front 2D/3D-mesh generator based on
abstract rules. Comput.Visual.Sci, 1:41–52.
[102] Schulz, C., von Tycowicz, C., Seidel, H.-P., and Hildebrandt, K. (2014). Animating
deformable objects using sparse spacetime constraints. ACM Trans. on Graphics (SIGGRAPH
2014), 33(4):109:1–109:10.
[103] Seifert, H. (1931). Konstruktion drei dimensionaler geschlossener Räume. Hirzel.
[104] Shi, X., Zhou, K., Tong, Y ., Desbrun, M., Bao, H., and Guo, B. (2008). Example-based
dynamic skinning in real time. In ACM Transactions on Graphics (SIGGRAPH 2008),
volume 27, pages 29:1–29:8.
References 107
[105] Shor, P. W. and Van Wyk, C. J. (1989). Detecting and decomposing self-overlapping
curves. In Proc. of ACM Symp. on Computational geometry, pages 44–50.
[106] Si, H. (2008). Adaptive tetrahedral mesh generation by constrained Delaunay refinement.
International Journal for Numerical Methods in Engineering, 75:856–880.
[107] Si, H. and Gärtner, K. (2011). 3d boundary recovery by constrained delaunay tetrahedral-
ization. International Journal for Numerical Methods in Engineering, 85(11):1341–1364.
[108] Si, H. and TetGen, A. (2006). A quality tetrahedral mesh generator and three-dimensional
delaunay triangulator. Weierstrass Institute for Applied Analysis and Stochastic, Berlin,
Germany.
[109] Sifakis, E. (2007). Algorithmic Aspects of the Simulation and Control of Computer
Generated Human Anatomy Models. PhD thesis, Stanford University.
[110] Sifakis, E. and Barbic, J. (2012). FEM simulation of 3D deformable solids: a prac-
titioner’s guide to theory, discretization and model reduction. In ACM SIGGRAPH 2012
Courses, page 20.
[111] Sifakis, E., Der, K., and Fedkiw, R. (2007a). Arbitrary cutting of deformable tetrahedral-
ized objects. In Symp. on Computer Animation (SCA), pages 73–80.
[112] Sifakis, E., Shinar, T., Irving, G., and Fedkiw, R. (2007b). Hybrid simulation of de-
formable solids. In Symp. on Computer Animation (SCA), pages 81–90.
[113] Sin, F. S., Schroeder, D., and Barbiˇ c, J. (2013). Vega: Nonlinear fem deformable object
simulator. Computer Graphics Forum, 32(1):38–50.
[114] Sorkine, O. and Alexa, M. (2007). As-rigid-as-possible surface modeling. In Symp. on
Geometry processing, volume 4, pages 109–116.
[115] Stomakhin, A., Howes, R., Schroeder, C., and Teran, J. M. (2012). Energetically
consistent invertible elasticity. In Symp. on Computer Animation (SCA), pages 25–32.
[116] Sutherland, I. and Hodgman, G. (1974). Reentrant Polygon Clipping. Communications
of the ACM, 17:32–42.
[117] Talvas, A., Marchal, M., Duriez, C., and Otaduy, M. A. (2015). Aggregate constraints
for virtual manipulation with soft fingers. IEEE transactions on visualization and computer
graphics, 21(4):452–461.
[118] Tamstorf, R., Jones, T., and McCormick, S. F. (2015). Smoothed aggregation multigrid
for cloth simulation. ACM Transactions on Graphics (TOG), 34(6):245.
[119] Tang, M., Liu, Z., Tong, R., Manocha, D., et al. (2018). I-cloth: incremental collision
handling for gpu-based interactive cloth simulation. In SIGGRAPH Asia 2018 Technical
Papers, page 204. ACM.
References 108
[120] Teng, Y ., Meyer, M., DeRose, T., and Kim, T. (2015). Subspace condensation: full space
adaptivity for subspace deformations. ACM Transactions on Graphics (TOG), 34(4):76.
[121] Teran, J., Sifakis, E., Blemker, S., Ng-Thow-Hing, V ., Lau, C., and Fedkiw, R. (2005a).
Creating and simulating skeletal muscle from the visible human data set. IEEE Trans. on
Visualization and Computer Graphics, 11(3):317–328.
[122] Teran, J., Sifakis, E., Irving, G., and Fedkiw, R. (2005b). Robust Quasistatic Finite
Elements and Flesh Simulation. In Symp. on Computer Animation (SCA), pages 181–190.
[123] Terzopoulos, D., Platt, J., Barr, A., and Fleischer, K. (1987). Elastically Deformable
Models. Computer Graphics (Proc. of ACM SIGGRAPH 87), 21(4):205–214.
[124] Teschner, M., Heidelberger, B., Müller, M., Pomeranets, D., and Gross, M. (2003).
Optimized spatial hashing for collision detection of deformable objects. In Proc. Vision,
Modeling, and Visualization Conference, pages 47–54.
[125] Thiery, J.-M., Tierny, J., and Boubekeur, T. (2012). CageR: Cage-Based Reverse
Engineering of Animated 3D Shapes. In Computer Graphics Forum, volume 31, pages
2303–2316.
[126] Umetani, N., Kaufman, D. M., Igarashi, T., and Grinspun, E. (2011). Sensitive couture
for interactive garment modeling and editing. ACM Trans. on Graphics, 30(4):90.
[127] Vaillant, R., Barthe, L., Guennebaud, G., Cani, M.-P., Rohmer, D., Wyvill, B., Gourmel,
O., and Paulin, M. (2013). Implicit skinning: real-time skin deformation with contact
modeling. ACM Transactions on Graphics (TOG), 32(4):125.
[128] Van Kampen, E. R. (1933). On the connection between the fundamental groups of some
related spaces. American Journal of Mathematics, 55(1):261–267.
[129] Wang, H., Hecht, F., Ramamoorthi, R., and O’Brien, J. F. (2010). Example-based wrinkle
synthesis for clothing animation. In ACM Trans. on Graphics (TOG), volume 29, page 107.
ACM.
[130] Wang, Y ., Jacobson, A., Barbiˇ c, J., and Kavan, L. (2015). Linear subspace design for
real-time shape deformation. ACM Transactions on Graphics (TOG), 34(4):57.
[131] Wang, Y ., Jiang, C., Schroeder, C., and Teran, J. (2014). An adaptive virtual node
algorithm with robust mesh cutting. In Symposium on Computer Animation (SCA), pages
77–85.
[132] Weber, O. and Zorin, D. (2014). Locally injective parametrization with arbitrary fixed
boundaries. ACM Transactions on Graphics (TOG), 33(4):75.
[133] Whitney, H. (1937). On regular closed curves in the plane. Compositio Mathematica,
4:276–284.
References 109
[134] Whitney, H. (1944). The self-intersections of a smooth n-manifold in 2n-space. Ann. of
Math, 45(2):220–246.
[135] Winkler, T., Drieseberg, J., Alexa, M., and Hormann, K. (2010). Multi-scale geometry
interpolation. In Computer Graphics Forum, volume 29, pages 309–318.
[136] Wojtan, C., Thürey, N., Gross, M., and Turk, G. (2009). Deforming meshes that split and
merge. In ACM Transactions on Graphics (TOG), volume 28, page 76. ACM.
[137] Wojtan, C. and Turk, G. (2008). Fast viscoelastic behavior with thin features. ACM
Trans. on Graphics (TOG), 27(3):47.
[138] Wuhrer, S., Lang, J., Tekieh, M., and Shu, C. (2015). Finite element based tracking of
deforming surfaces. Graphical Models, 77:1–17.
[139] Xu, H. and Barbiˇ c, J. (2014). Signed distance fields for polygon soup meshes. In Proc.
of Graphics Interface, pages 35–41.
[140] Yeung, Y .-H., Crouch, J., and Pothen, A. (2016). Interactively cutting and constraining
vertices in meshes using augmented matrices. ACM Trans. on Graphics (TOG), 35(2):18.
[141] Zhang, H. (2004). Discrete combinatorial laplacian operators for digital geometry
processing. In Proceedings of SIAM Conference on Geometric Design and Computing, pages
575–592. Nashboro Press.
[142] Zhou, Q., Grinspun, E., Zorin, D., and Jacobson, A. (2016). Mesh arrangements for solid
geometry. ACM Transactions on Graphics (TOG), 35(4):39.
[143] Zhu, Y ., Sifakis, E., Teran, J., and Brandt, A. (2010). An efficient multigrid method for
the simulation of high-resolution elastic solids. ACM Trans. on Graphics (TOG), 29(2):16.
Appendix A
Fast Right-Hand Side Evaluation for
Reduced Shape Editing
Let N(i) be the 1-ring neighborhood of vertex i and R
i
the “best-fit” rotation on N(i). The
ARAP right-hand side can be expressed as
b=∑
i
B
i
R
T
i
, for B
i
= ∑
j∈N(i) w
i j
2
(S
i
−S
j
)
T
( ¯ p
i
− ¯ p
j
)
T
R
T
i
. (A.1)
Index i loops over all vertices, and S
i
∈R
1×n
selects the i-th vertex. For rotation clustering, we
partition the vertices using k-means. The right-hand side is then approximated by
b
c
=∑
k
B
c
k
R
T
k
=∑
k
(∑
i∈C
k
B
i
)R
T
k
. (A.2)
Index k loops over each cluster C
k
. Rotation R
k
is the rotation for cluster k. In the shape
subspace, we have
˜
b
c
=∑
k
˜
B
c
k
R
T
k
=∑
k
(U
T
B
c
k
)R
T
k
. (A.3)
We can pre-compute
˜
B
c
k
∈R
r×3
. At runtime, the complexity of evaluating
˜
b
s
is only O(rm
c
),
where m
c
is the number of clusters.
111
When using cubature, we first create some training shapes. Given a few vertices forming a
groupD of size m
d
, we wish to find
˜
B
d
k
so that the right-hand side computed via
˜
b
d
=∑
k∈D
˜
B
d
k
R
T
k
(A.4)
is on average the closest to the ground truth values
˜
b, evaluated at all the training shapes. Since
˜
b
d
is linear in each element in
˜
B
d
k
, a simple quadratic optimization can be used to minimize
the sum of the errors evaluated at all the training shapes. However, we should be aware that
this brute force optimization may not preserve linear precision, i.e., the method cannot recover
the rest shape ¯ p when all rotations are identity. The rotation clustering method preserves
this property naturally due to the way it is formulated. We resolve this problem by adding
a linear constraint to the cubature quadratic optimization,
˜
b(I
3
)=∑
k∈D
˜
B
d
k
I
3
, where I
3
is the
3× 3 identity matrix. Typically, m
c
and m
d
can be set to a size of r, reducing right-hand side
evaluation complexity to O(r
2
).
Appendix B
Immersion of Self-Intersecting Solids and
Surfaces: Supplementary Material 1:
Topology
B.1 Topological definitions
A mapping f ∶ X→ Y is injective (also called “one-to-one”) if for every x,y∈ X such that x≠ y,
we have f(x)≠ f(y). It is surjective (also called “onto”) if for each y∈ Y, there exists x∈ X
such that f(x)= y, i.e., f(X)= Y.
A topological space is a set X equipped with a topology. A topology is a collectionτ of subsets
of X satisfying the following properties: (1) the empty set and X belong toτ, (2) any union of
members ofτ still belongs toτ, and (3) the intersection of any finite number of members of τ
belongs toτ. A neighborhood of x∈ X is any set that contains an open set that contains x.
A homeomorphism between two topological spaces is a bijection that is continuous in both
directions. A mapping between two topological spaces is continuous if the pre-image of
each open set is continuous. Homeomorphic topological spaces are considered equivalent for
topological purposes.
B.1 Topological definitions 113
An immersion between two topological spaces X and Y is a continuous mapping f that is
locally injective, but not necessarily globally injective. I.e., for each x∈ X, there exists a
neighborhoodN of x such that the restriction of f ontoN is injective.
Topological closure of a subset A of a topological set X is the set A of elements of x∈ X with
the property that every neighborhood of x contains a point of A.
A topological space X is path-connected if, for each x,y∈ X, there is a continuous mapping
γ∶[0,1]→ X (a curve), such thatγ(0)= x andγ(1)= y.
A manifold of dimension d is a topological space with the property that every point x has a
neighborhood that is homeomorphic toR
d
(or, equivalently, to the open unit ball inR
d
), or a
halfspace ofR
d
(or, equivalently, to the open unit ball inR
d
intersected with the upper halfspace
ofR
d
). Intuitively, this means that the manifold is a space that is “locally Euclidean”. The
points that fall into the second case form the manifold boundary. In addition, one requires
that manifolds are Haussdorff spaces, i.e., distinct points have disjoint neighborhoods, and
second-countable, i.e., there exists a countable collectionU = {U
i
}
∞
i=1
of open sets such that
each open set is a union of a subfamily ofU. These additional conditions just ensure that
the manifold is a sufficiently “nice” space; all subsets of Euclidean spaces R
k
satisfy these
additional conditions. For d= 2, important examples of manifolds in this paper are the sphere
and torus with k≥ 1 handles, both of which have no boundary. For d= 3, in this paper, manifolds
have boundaries. Typical examples are solid balls, or solid tori with k≥ 1 handles.
A manifold M embedded into an Euclidean space is compact if, intuitively, there exists a
bounding box that completely contains it (i.e., the manifold is finite in size), and M is closed,
i.e., contains all of its limits. Formally, a topological space X is compact if for any infinite
family of open sets {U
α
}
α
, that covers X, there exists a finite subfamily that also covers X.
By the Whitney embedding theorem [134] and its extensions [87], every compact d−manifold
(even if non-smooth) can be embedded intoR
N
, for some properly high N, e.g., N= 2d+1 is
sufficient. Bijective continuous mapping from a compact topological space X onto a Haussdorff
topological space (such as anR
k
) are homeomorphisms.
B.1 Topological definitions 114
Disjoint union X∐Y of sets X and Y is a union of X and Y where we have ensured that the
elements of X∩Y are present twice. Formally, it can be formed as X∐Y =(X×{0})∪(Y×{1}),
where× denotes the Cartesian set product, and∪ is the usual set union.
Equivalence relation∼ on a set X is a relation with the properties that (1) x∼ x for each x∈ X,
(2) if x∼ y for some x,y∈ X, then also y∼ x, and (3) if x∼ y and y∼ z for some x,y,z∈ X, then
x∼ z. An equivalence class is a maximal set such that all elements are equivalent to one another.
A quotient space of a set X with respect to an equivalent relation∼ on X is the set of equivalence
classes of∼ on X.
A quotient topological space of a topological space X with respect to an equivalence relation
on X, is the quotient space of X with respect to∼, equipped with the quotient topology. The
quotient topology consists of all sets with an open pre-image under the projection map that
maps each element of X to its equivalence class. Intuitively speaking, a quotient topological
space is the result of identifying or "gluing together" certain points of a given topological space.
This is commonly done in order to construct new spaces from the existing ones.
Self-touching input triangle meshes are defined as follows. Let x∈
ˆ
M whereγ(x)≥ 2, i.e., x∈
Γ. Then, we call x to be a crossing intersection if there exists a (sufficiently small) neighborhood
N of ˆ ρ(x) such that ˆ ρ
−1
(N) consists ofγ(x) components in
ˆ
M, and the ˆ ρ-projections of any
two of these compoments cross each other inR
d
, as opposed to merely touching. Then, we
call the input mesh M self-touching if there exists at least one point x∈Γ that is not a crossing
intersection. The intent is to capture degenerate inputs where a surface touches itself at a single
point, or along a loop, but does not penetrate. For d= 2, an example of such a degeneracy is a
circle touching a square. For d= 3, an example is a solid bowl placed upside down on a solid
box. We note that self-touching input triangle meshes should not be confused with self-touching
cells. The former is a degenerate input that our algorithm does not permit, whereas the latter is
a common occurence in our algorithm, even with non-degenerate inputs.
B.1 Topological definitions 115
A covering projection between path-connected topological spaces X and Y is a continuous
surjective mapping p∶ X→ Y, such that for each y∈ Y, there exists a neighborhood U so that
p
−1
(U) is a union of disjoint open sets, each of which is mapped homeomorphically onto U.
The total turning number of a planar curve is the angle, divided by 2π, by which a passenger
riding down this curve will rotate when they make one loop around the curve. The total turning
number is a non-zero integer. It should not be confused with the winding number. The winding
number is a property of a point with respect to the curve. The total turning number is an inherent
property of the curve independent of any particular point.
A group is a set G equipped with a binary operation○, such that (1) x○y∈ G for each x,y∈ G,
and (2)○ is associative ((x○y)○z= x○(y○z) for all x,y,z∈ G), (3) x○e= e○x= x for a special
element e∈ G (the “unity”), and (4) for each x∈ G, there exists the inverse y so that x○y= y○x= e.
A homomorphism is a mapping between groups G and H that respects the group operation,
i.e., f(a○b)= f(a)○ f(b) for any a,b∈ G.
An isomorphism is a bijective homomorphism. Groups are isomorphic if there exists an iso-
morphism between them. Isomorphic groups are considered equivalent for algebraic purposes.
A subgroup H of a group G is a subset H of G which is closed under the group operation.
The free product G
1
∗G
2
∗...∗G
n
of groups G
1
,...,G
n
is the group of finite words (of arbitrary
length) a
j
1
i
1
a
j
2
i
2
...a
j
ℓ
i
ℓ
, where a
i
∈ G
i
,ℓ≥ 0, j
k
≥ 1, and i
k
are arbitrary integer indices between 1
and n.
The fundamental group π
1
(X) of a path-connected topological space X is the set of loops,
starting and ending at some basepoint b∈ X, under the equivalence relation that equates loops
that can be continuously morphed into one another. Group operation is path concatenation.
B.2 Topological proofs 116
B.2 Topological proofs
(Section 3)
ˆ
M is a manifold: If the equivalence class ˆ x∈
ˆ
M originates from x∈ M that is in the
interior of the triangle face (i.e., ˆ x={x}), then x trivially has a neighborhood homeomorphic
to the open ball inR
d−1
. If ˆ x originates from the interior of an edge between two triangles,
then, because M has manifold connectivity, there are exactly two triangle joining at that edge.
Therefore, it is possible to find a neighborhood homeomorphic to the open ball in R
d−1
. The
remaining case is that ˆ x originates from a vertex x of M. Then, take any triangle that has x as its
vertex. Because M has no boundary, there must be two adjacent triangles that share an edge
and vertex x with M. We can continue adding such adjacent triangles. Because the number of
triangles is finite, we must eventually form a complete 360 degree loop around x. Because, due
to the manifold connectivity assumption, the triangles touching x must form a continuous fan,
the patch of triangles looping around x is homeomorphic to a closed disk inR
d−1
, and therefore
there exists a neighborhood of ˆ x homeomorphic to the open ball inR
d−1
. The Haussdorff and
countable separability axioms are trivially satisfied because
ˆ
M is constructed from triangles in
the Euclidean space. ∎
(Section 3) Augmentation of ρ to an immersion ˆ ρ of
ˆ
M onto M∶ As stated in paper, we
assume that M is non-degenerate. In order for ρ to be augmentable, it needs to respect the
gluing equivalence relation of
ˆ
M, i.e., it maps all members of each equivalence class of
ˆ
M
into the same point inR
3
. This is true because all the “glued” equivalent points are co-located
inR
d
. We can therefore un-ambiguously define ˆ ρ(ˆ x)=ρ(x). Because ρ is continuous, ˆ ρ is
also continuous. It remains to be proven that ˆ ρ is locally injective. Given ˆ x∈
ˆ
M, we can form
the manifold neighborhood N of x in M, as in the proof of manifoldness of
ˆ
M. Because the
restriction of the mappingρ onto this neighborhood is injective, it follows that ˆ ρ is injective on
the augmentation
ˆ
N={[x];x∈ N} of N, where [x] is the equivalence class of x. ∎
(Section 4.1) Topology of cells: By construction, cells are open and connected. For connected
inputs, because M is not degenerate, the boundary B of the closure of each cell is a d− 1
manifold without boundary. B is also orientable, and a closed and bounded subset ofR
d
. For
d= 3, by the characterization theorem of 2-manifolds [21], B must topologically be either a
sphere, or a torus with k≥ 1 handles. For d= 2, it must be a topological circle. ∎
B.2 Topological proofs 117
(Section 5.1) Finiteness of cells: If there exists a volume immersion ˆ σ from a compact manifold
ˆ
S ontoR
d
whose boundary is
ˆ
M, then the pre-image of each cell ˆ σ
−1
(C)⊂
ˆ
S consists of a finite
number of components in
ˆ
S.
Proof: Suppose the pre-image of a cellC consists of an infinite number of components. Then,
construct a sequence y
i
∈
ˆ
S such that each y
i
is in a different connected component of ˆ σ
−1
(C).
Because
ˆ
S is compact, there exists a subsequence {x
i
}
i
of the sequence {y
i
}
i
such that x
i
has a limit x∈
ˆ
S. Suppose x is not a boundary point of the manifold
ˆ
S. Then, there exists a
neighborhoodN of x that is homeomorphic to the unit ball inR
d
, such that the restriction of ˆ σ
ontoN is a homeomorphism onto its image. This last property holds because ˆ σ is an immersion.
Because ˆ σ is continuous, ˆ σ(x) must be inC. Now, we need to consider several cases depending
on whether ˆ σ(x) lies in the interior ofC (i.e., insideC; note that by our definition, cells are
open sets), or on a patch, arc or corner point on the boundary ofC. Suppose ˆ σ(x)∈C. Then, it
is possible to select a neighborhoodN
′
⊂N of x such that ˆ σ(N
′
)⊂C, and the restriction of ˆ σ
ontoN
′
is a homeomorphism onto its image. However, this implies thatN
′
does not contain
any points that map toρ(M). Because x
i
→ x, the sequence will eventually enterN
′
, and so we
will find x
i
0
and x
i
0
+1
that are both inside ofN
′
. We can then join these two points with an arc
contained completely insideN
′
. None of the points on this arc map toρ(M), and therefore x
i
0
and x
i
0
+1
are in the same connected component of ˆ σ
−1
(C). This is a contradiction, proving the
original statement. The other cases where ˆ σ(x) lies on a patch, arc or corner point, or where x
is a boundary point of
ˆ
S, can be handled in the same way; we just have to consider “half-disk”,
“quarter-disk”, or “eight-disk” neighborhoodsN
′
of x. ∎
(Section 5.1) Restrictions of immersions to cells are covering projections. In addition, we
will simultaneously prove the following. Let p= ˆ σ
S
∶S→C be the restriction of the immersion
ˆ σ onto a pre-cellS ⊂
ˆ
S. Then, for any y∈C the cardinality of p
−1
(y) is finite and constant on C
(the “number of sheets”).
Proof: By the definition of S, p is surjective and continuous. Suppose the cardinality of p
−1
(y)
would be infinite for some y∈C. Then by the same compactness argument as in the proof of
the finiteness of cells (above), we can construct a sequence x
i
∈S such that p(x
i
)= y for all i,
and x
i
has a limit x∈S ⊂
ˆ
S. Because ˆ σ is an immersion, there exists a neighborhood of x such
that the restriction of ˆ σ is injective. However, because the sequence x
i
eventually has to enter
this neighborhood, injectivity is impossible, as p(x
i
)= y for all i. Therefore, the cardinality
of p
−1
(y) is finite for all y∈C. Let us now prove the existence of the neighborhood U from
B.2 Topological proofs 118
the definition of the covering space, and that the cardinality of p
−1
(y) is constant onC. Let
c(y)=∣p
−1
(y)∣ for y∈C. For each of the c(y) distinct elements z
i
∈S, such that p(z
i
)= y, there
exists a neighborhood Z
i
of z
i
that maps homeomorphically onto a neighborhood of y. We can
always shrink Z
i
so that they are pairwise disjoint. Because the number of z
i
is finite, we can set
U =⋂
i
p(Z
i
).
By construction, p
−1
(U) is a union of disjoint open sets, each of which is mapped homeomor-
phically onto U. It also follows that for any y
′
∈ U, c(y
′
)=∣p
−1
(y
′
)∣= c(y), i.e., c(y) is locally
constant on U. Any locally constant scalar function on a path-connected space is constant [87].
Because the cellC is path-connected, c(y) is constant on the entire cellC. ∎
(Section 5.1) Theorem 1: If a valid input M is simply volume-immersible, then for each cellC
of the cell complex of M, the number of connected components of ˆ σ
−1
(C) equals the winding
number wind(x,M) of any point x∈C with respect to M. The winding number wind(x,M) is
the same for all points x∈C.
Proof: Although we numbered this theorem as Theorem 1 for pedagogical reasons, it actually
depends on Theorems 4 and 5, whose proof must logically happen first. Note that Theorems 4
and 5, and Corollary 6 do not discuss or need the winding number, or this Theorem 1, so there
is no circular dependency. This is because Theorems 4 and 5 never specifically postulated any
number of copies of each cell in the immersion graph. So, we defer this proof until later on in
this document, when both Theorem 4 and 5 are proven.
(Section 5.1) Lemma (total turning number): Let i be an immersion of the unit disk inR
2
intoR
2
, such that its restriction on the unit circle is a surjective immersion onto a closed 2D
loopγ∶[0,1]→R
2
;γ(0)=γ(1). Then, the total turning number ofγ must be 1 or -1.
Proof: Because i is an immersion, there exists ε > 0 such that the the restriction of i on a
disk of radiusε centered at the origin is injective. Now, consider the following homotopy of
immersions:
H ∶ [0,1]×[0,1]→R
2
, (B.1)
H(s,t)=((1−s)+
sε
2
)e
2πt
, (B.2)
B.2 Topological proofs 119
Figure B.1 The generators ofπ
1
(C). In this case,C is a solid torus with 3 handles (k= 3), viewed
(on the left) under an orthographic projection from the top.
where we have identified R
2
with the complex plane. Because on theε-disk, i is a homeomorphic
onto its image, this is a homotopy of immersions between the immersion onto a closed 2D loop
γ (s= 0) and an immersion of circle onto S
1
(i.e., a topological circle; s= 1). By the Whitney-
Graustein theorem [133], two immersions are homotopic through a family of immersions if and
only if they have the same turning number. This means that the total turning number ofγ must
be 1 or -1. ∎
(Section 5.1) Uniqueness and existence of lifted paths: Let p∶ X→ Y be a covering space
projection. Letγ be a path in Y. Then, by picking an arbitrary starting point x∈ X, there exists a
unique path ˜ γ in X that starts at x and that “lifts”γ, i.e., p○ ˜ γ =γ. This is a standard result in
algebraic topology [42].
(Section 5.1) Theorem 2: If a valid input M is the boundary of an immersion from a disk in
R
d
, then the immersion is simple. This holds both for d= 2 and d= 3.
Proof: LetC be a cell of the cell complex of M, andS be its “pre-cell”, i.e., one of the connected
components of ˆ σ
−1
(S. Let p= ˆ σ
S
∶S→C be the restriction of the immersion ˆ σ ontoS ⊂
ˆ
S. By
the “Restrictions of immersions” lemma (above),S is a covering space overC. In order to prove
simple immersibility, we must prove that the number of sheets ofS overC is 1. Topologically,
cellC is the interior of a torus with k≥ 0 handles. This is because the boundary of the cell is
a connected orientable compact 2-manifold, and therefore, by the classification theorem of
2-manifolds [21], has to be homeomorphic to a torus with k≥ 0 handles (the case k= 0 denotes
the sphere). We can homeomorphically deform cellC into a canonical shape, whereby the walls
are planar at 90-degree angles (see Figure B.1), and the holes are perpendicular to the y-axis.
The spaceS can then be seen as “hovering” aboveC along the y-axis.
B.2 Topological proofs 120
For k= 0, the fundamental groupπ
1
(C) is {e}. For k> 0, the fundamental group ofC can be
computed using the Seifert-van Kampen theorem [103, 128] using the fundamental group of
the solid torus, which isZ. We therefore haveπ
1
(C)=Z∗...∗Z (free product; k times), i.e.,
π
1
(C) is the set of all words formed by k letters a
1
,...,a
k
, i.e., all finite words of the form
a
j
1
i
1
a
j
2
i
2
...a
j
ℓ
i
ℓ
. Each of the words corresponds to a loop which starts at some fixed base point b,
and loops around one of the holes exactly once; the finite words corresponds to all possible
loops inC (see Figure B.1). The projection p induces a homomorphism p
#
ofπ
1
(S) intoπ
1
(C),
whereby every loopγ inS with a basepoint in p
−1
(b) is mapped to p○γ. A well-known result
in algebraic topology is that the number of sheets equals the index of the subgroup p
#
(π
1
(S))
inπ
1
(C) [42].
We therefore need to prove that this index is 1. For k = 0, because π
1
(C) is trivial, so is
p
#
(π
1
(S)), and therefore, this index is 1. For k> 0, we use another standard result in algebraic
topology, which states that p
#
(π
1
(S)) is isomorphic to the subgroup of loops inπ
1
(C) which
have the property that they lift to loops inS [42]. We will now prove that every loop inπ
1
(C)
lifts to a loop inS. Then, it follows that p
#
(π
1
(S)) is isomorphic toπ
1
(C), and therefore, the
number of sheets equals 1. Because the generator loops a
1
,...,a
k
generateπ
1
(C), it is sufficient
to prove that they each lift to a loop inS. Suppose one of them does not (say, a
1
without loss
of generality). Let us place the coordinate system so that the origin is in the middle of the
handle of a
1
, and the y-axis is pointing through the handle. This means that the lifting of a
1
intoS is not a loop, but an arc joining two different points B
1
and B
2
from p
−1
(b)⊂S. We
can then perform the lifting again starting from B
2
, obtaining another arc, from B
2
to B
3
. This
process will eventually terminate with B
M
= B
1
, for M≥ 3, producing a loop inS that winds
itself around the y-axis once, and whose projection via immersion p intoC winds itself M−1
times around the y-axis. Because p is an immersion, this is not possible by the total turning
lemma (above). Hence, all the generator loops lift to a loop inS. ∎
(Section 5.2) Theorem 5: Let M be a valid input. Assume that there is a cell graph and
a B-patch ownership assignment that satisfy Rules 1-7 given in Section 5.1. Then, M is
volume-immersible, and the graph can be used to construct the immersion.
Proof: By gluing the cells according to the cell graph, we will form a d-manifold
ˆ
S, and an
immersion ˆ σ from
ˆ
S intoR
d
such that
ˆ
M is the boundary of
ˆ
S, and the restriction of ˆ σ onto
ˆ
M is ˆ ρ. We first define
ˆ
S, as follows. Note that
ˆ
S will be an “abstract” manifold, and if M has
self-intersections, not embedded intoR
d
. Number the cell graph nodes as 1,2,..., in arbitrary
B.2 Topological proofs 121
order. Then, define S as the disjoint union
S= ∐
i∈nodes(G) C(i)×{i}, (B.3)
where i runs over all the nodes of the cell graph, × is the set Cartesian product, C(i) is the
cell corresponding to node i, and is the topological closure operation. We then form the
set
ˆ
S as a quotient set of S with respect to a gluing relation whereby we identify points on
a closure of a patch shared by two connected nodes as equivalent. We define the mapping
σ∶ S→R
d
asσ(x,i)= x, where x∈C(i)∈R
d
. We then define our volume immersion ˆ σ∶
ˆ
S→R
d
,
by augmentingσ with respect to the gluing relation. Because each patch is owned by exactly
one node (Rule 7), we can embed
ˆ
M into
ˆ
S. The immersion ˆ σ trivially matches ˆ ρ on
ˆ
M.
We now prove that
ˆ
S is a d-manifold. Let x∈
ˆ
S. We consider the different cases. (1) If x is
in the interior of a node c, then since cells are d-manifolds, x trivially has a neighborhood
homeomorphic toR
d
. (2) If x is on the boundary of c, we have several subcases. (2.1) If x is on
a B-patch p that the node owns, we consider whether it is on the boundary of the patch or not.
(2.1.1) If it is in the interior of p, then x has a neighborhood homeomorphic to the halfspace of
R
d
due to manifoldness of cells. (2.1.2) Else, it is on the boundary of p, and therefore in the
closure of one or more arcs. Hence, x is shared with other topologically neighboring patches.
By Rule 4, the graph nodes owning the patches are connected. By Rule 5, no other B-patches
contribute to the neighborhood of x. So, there is a neighborhood of x homeomorphic to the
halfspace ofR
d
, bounded by the patches and emanating into the interior of the cells. (2.2) Else,
node c declined p. By Rule 2, there must be another node d connecting c across p. By Rule
1, there are no other nodes connected across p. (2.2.1) If x is in the interior of p, then there is
a neighborhood of x homeomorphic toR
d
straddling p and emanating into the interiors of c
and d. (2.2.2) Else, x is on the boundary of p. Since the boundary of a patch consists of arcs,
and each arc is shared by two pairs of topologically neighboring patches, there are four patch
connections surrounding x. If the four patch connections are from four different patches, we
are either in the situation in Figure 15 (in main paper), or (2.1.2) applies. By Rule 6, the nodes
surrounding the arcs whose closure contains x, are all connected. So, there is a neighborhood of
x homeomorphic toR
d
straddling patches and emanating into the interior of cells.
Finally, we prove that ˆ σ is an immersion. The mapping ˆ σ is continuous because the glued cells
are adjacent in space, and is trivially locally injective if x is in the interior of a node. If x is on
the boundary of a node, then local injectivity follows from the fact that the glued cells occupy
B.2 Topological proofs 122
distinct regions of space inR
d
. ∎
(Section 5.1) Theorem 1: If a valid input M is simply volume-immersible, then for each cellC
of the cell complex of M, the number of connected components of ˆ σ
−1
(C) equals the winding
number wind(x,M) of any point x∈C with respect to M. The winding number wind(x,M) is
the same for all points x∈C.
Proof: We now give our “deferred” proof of Theorem 1. By Corollary 6, there exists an
immersion graph G and a B-patch ownership assignment that satisfy Rules 1-7 given in Section
5.2. LetC
1
,...,C
n
be the cells of the cell complex of M. By the “Finiteness of Cells” lemma,
the number of connected components (the “pre-cells”) ofC
i
is finite; denote it by c
i
. Let x∈C
i
0
,
for some i
0
, 1≤ i
0
≤ n. Then, we have
wind(x,M)=
∫
y∈M
Ω(x,y)dS=
n
∑
i=1
c
i
∑
j=1
∫
y∈∂C
i
Ω(x,y)dS=
n
∑
i=1
c
i
∑
j=1
(if (x∈C
i
) then 1 else 0)=
(B.4)
=
n
∑
i=1
c
i
(if (x∈C
i
) then 1 else 0)= c
i
0
. (B.5)
where∂C
j
is the boundary of cellC
j
, and whereΩ(x,y) is the “winding number kernel”, namely
Ω(x,y)dS is the solid angle seen from x towards y, due to a small surface patch dS on M, as
shown in [58]. The second equality in (B.5) holds because of Rules 1, 2 and 7. Namely, every
boundary patch of a node is either shared by two nodes that are connected in G (in which case
the contributions by these two nodes on this patch cancel each other in the sum), or it is owned
by the node and therefore appear exactly once and forms a part of the boundary of M. ∎
(Section 6.1) Our 3D immersion problem is NP-complete: For d = 2, it has been proven
in [33] that the problem of determining if the given planar curve is a boundary of an immersed
surface, is NP-complete. The size of the problem is measured as a number of 2D intersections
(2D arcs), or equivalently, 2D patches or 2D cells. For d= 3, it is easiest to measure the problem
size as a number of patches. For d= 3, we can easily prove that the problem is also NP-complete.
Given a 2D problem, M
2D
, we can pick a 3D line that does not intersect the bounding box of
M
2D
and revolve M
2D
around it 360
○
, producing a 3D problem for a 3D input M. Due to the
geometry of our 3D surface of revolution, M
2D
is a boundary of an immersed surface for d= 2
if and only if M is a boundary of an immersed surface for d= 3. Therefore, the immersibility
B.2 Topological proofs 123
problem for d= 2 is polynomial-time reducible to the immersibility problem for d= 3. It now
follows from Eppstein’s result that the problem for d= 3 is NP-complete.
Appendix C
Immersion of Self-Intersecting Solids and
Surfaces: Supplementary Material 2:
Tetrahedralization
C.1 Inside-outside test without forming the region boundary
surface
In this section, we prove the statement from Section 7 of the main paper that we do not need
to form the region boundary surface (which would require using CSG) to perform inside-
outside region queries. The closest site theorem of Baerentzen [7] states that a point is located
outside of a closed manifold non-self-intersecting mesh if and only if the vector from the
closest site on the mesh to the query point has a positive dot product with the outward surface
pseudo-normal. Pseudo-normals are defined in [ 7]. We now generalize this test to non-closed
non-self-intersecting manifold meshes M (i.e., meshes with boundary). We say that a point p is
pseudo-outside of M, by definition, if the vector from the closest site on M to p has a positive
dot product with the outward pseudo-normal at the closest site.
Theorem: Let p be a point inside a tet, and let R be one of the regions (both+ and− regions
are ok) of this tet, as defined in Section 7 of the main paper. Then, p is inside the region if and
only if p is pseudo-inside the closest piece of R. This theorem also gives an algorithm to rapidly
C.2 Pseudo-normal on the piece boundary 125
perform the inside-outside test: by this theorem, it is sufficient to perform pseudo-tests against
all the pieces of a region. This means that the boundary mesh of R does not need to be explicitly
formed on the faces of the tet. We perform the point-vs-triangle distance calculations for the
pseudo-tests using exact arithmetic.
Proof: Suppose a query point p is inside a tetrahedron, but outside of a region. Then, we claim
that the closest site to p on the region has to be on a triangle on one of its pieces. If this was
not the case, the closest site is on the tetrahedron surface. Consider the line segment joining
this closest site and p. Because the tetrahedron is convex, the line segment is completely inside
the tetrahedron. If the segment does not intersect any pieces, this implies that the line segment
never leaves the region, which is a contradiction with p being outside of the region. Therefore,
it intersects one piece, but then this intersection is even closer to p than the original closest site,
which is a contradiction. Because the closest point is therefore on a boundary triangle, we can
now apply Baerentzen’s theorem, which proves that our algorithm will in this case correctly
identify that p is outside of the region.
Suppose the query point p is inside the region. Let P be the piece of this region which has the
closest point to p compared to the other pieces. Piece P cuts the tetrahedron into two parts,
e
0
,e
1
. With out loss of generality, we assume e
0
does not contain p, and therefore, e
1
contains
the region. Since p is outside e
0
, according to the proof in the last paragraph, the closest site on
e
0
to p must be on P and the pseudo-normal test on P to p should report outside for e
0
. Since
e
0
and the region share the piece P, P must have different orientations for e
0
and the region.
Therefore, the pseudo-normal test on P for the region should give the opposite result to the test
on P for e
0
, which proves that our algorithm correctly reports that p is inside the region. ∎
C.2 Pseudo-normal on the piece boundary
Although we do not need to form the complete region mesh, we still need the tet faces to
compute the pseudo-normal if the closest site to p on a piece is on the boundary of the piece,
i.e., it is on a tet face. The pseudo-normal that we compute must be identical to the one that
would be computed using a complete region mesh. We achieve this as follows. We have to
consider two separate cases.
C.2 Pseudo-normal on the piece boundary 126
Figure C.1 Handling closest sites at piece boundaries. Left: a piece inside a tetrahedron
T
0
T
1
T
2
T
3
, where T
3
is occluded by the piece. Vector N
1
is the normal of the tet face T
0
T
1
T
2
.
Vertices v
i
are on the boundary of a piece. Vector n is the normal of the triangle with the
boundary edge v
1
v
0
. Top-right: degeneracy cases for boundary edges. Bottom-right: tet face
angle cases for boundary vertices.
Boundary edges: We first consider the case where the pseudo-normal is located in the interior
of a boundary edge of a piece. Figure C.1 explains the procedure to find the correct pseudo-
normal. For the boundary edge v
1
v
0
in the figure (left), we first determine on which tet face this
edge lies. Then, we compute the pseudo-normal ˜ n as the average of the triangle normal n on the
piece neighboring v
1
v
0
, and tet face normal N∶ ˜ n= unit(n+N). As described by Baerentzen’s
work, the final inside-outside result is the sign of dot(p−c, ˜ n), where p is the query position
and c is the closest site on v
1
v
0
. Degeneracy may occur if n is close to−N, or p−c is almost
perpendicular to ˜ n. In either case, we reach the situation shown in Figure C.1, top-right. Looking
along the direction of the boundary edge v
0
v
1
, the tet face T
0
T
1
T
2
is seen as a line segment
T
0
T
1,2
, and the edge v
1
v
0
is seen as a point v
0,1
. When n is close to−N (left sub-figure), we
observe that the angle at v
0,1
of the triangle T
0
v
0,1
p is less or equal to the angle β, because
otherwise, the closest site will not be on v
1
v
0
. Sinceβ is entirely outside the region, we can
safely classify p as outside. When p−c is almost perpendicular to ˜ n, n must be close to N (right
sub-figure). In this case, the angle at v
0,1
of the triangle T
1,2
v
0,1
p is less or equal toβ, so we
can safely classify p as inside.
Boundary vertices: We now consider the case where we need the pseudo-normal on a
boundary vertex of a piece. For one boundary vertex v, we first find the location of this vertex
and its two neighboring boundary edges to find whether we are in the case of v
1
(interior of
C.3 Meshing self-touching cells 127
tet face) or v
4
(on tet edge), shown in Figure C.1, left. (1) In the case of v
1
, by Baerentzen’s
theorem, the pseudo-normal is ˜ n=∑
i
w
i
n
i
+ w
N
N, where w
i
is the weight for normal n
i
on
triangle i surrounding v
1
, and w
N
is the weight for the tet face normal N. Normal weights are
computed as the triangle angle at v
1
. To compute w
N
, we first compute the angle α between
edge v
1
v
0
and v
1
v
2
. Depending on the orientation of the piece, w
N
is eitherα or 2π−α. Our
intuition for choosing between these two options is that the faces contributing to the pseudo-
normal should have consistent orientations. Our specific procedure (illustrated in Figure C.1,
bottom-right) is as follows. Without loss of generality, assume the directions of edges v
2
v
1
and v
1
v
0
are consistent with the orientations of the piece triangles that they belong to. Then, if
cross(v
1
−v
2
,v
0
−v
1
) has an opposite sign to N, we deduce that w
N
must be 2π−α. Otherwise,
w
N
has to be α. (2) In the case of v
4
, we have ˜ n=∑
i
w
i
n
i
+ w
N
1
N
1
+ w
N
2
N
2
, where w
N
i
,N
i
are the weights, and tet face normals neighboring the tet edge containing v
4
. Let N
1
be the
tet face normal of T
0
T
1
T
2
and N
2
be the tet face normal of T
0
T
3
T
1
. Then, w
N
1
is the angle
between edge v
4
v
3
and edge v
4
T
i
, and w
N
2
is the angle between edge v
4
v
5
and edge v
4
T
i
,
where i is either 0 or 1. We again use the intuition that faces around v
4
should have consistent
orientations. Without loss of generality, assume that the directions of edges v
5
v
4
and v
4
v
3
are consistent with the orientation of the piece triangles that they belong to. We compute
dot(cross(v
4
− v
5
,v
3
− v
4
),T
0
− T
1
). If the resulting sign is positive, we select i= 1, otherwise
i= 0. We perform these computations in exact arithmetic.
C.3 Meshing self-touching cells
Special care must be taken when the tetrahedralization algorithm of Section 7 in the main paper
is used to mesh self-touching cells produced by our immersion algorithm. If a tetrahedron
includes two regions that touch each other at a self-touching edge of the cell, then during the
pseudo-normal test, we also check the distance between the closest site on the piece and the
query point. If the distance is zero under exact arithmetic, we return the “outside” result for the
pseudo-normal test.

Abstract (if available)

Abstract Although physically based deformations and contact are well-understood in science, their application to practical computer animation pipelines poses several challenges. In computer animation practice, characters, plants and inanimate objects are represented as triangle meshes, and these meshes easily and often inter-penetrate, both during modeling and animation. It is very difficult to design algorithms that resolve contact, yet simultaneously permit the artists to retain direction and control over the animations. ❧ This thesis introduces several methods to address intersections in various stages of computer animation pipelines and improve their physical realism, all the while retaining artistic control over the modeling and animation process. The first problem that we tackle is how to equip standard shape modeling algorithms with robust collision response handling. Shape deformation is a well-studied problem in computer graphics with numerous articles published on the topic over the years. However, to date, there has been very little work on adding contact handling to shape deformation modeling

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Plant substructuring and real-time simulation using model reduction

PDF

CUDA deformers for model reduction

PDF

Anatomically based human hand modeling and simulation

PDF

Closing the reality gap via simulation-based inference and control

PDF

Computational models for multidimensional annotations of affect

PDF

Feature-preserving simplification and sketch-based creation of 3D models

PDF

Data-driven 3D hair digitization

PDF

Multi-scale dynamic capture for high quality digital humans

PDF

Acquisition of human tissue elasticity properties using pressure sensors

PDF

3D urban modeling from city-scale aerial LiDAR data

PDF

Real-time simulation of hand anatomy using medical imaging

PDF

Correspondence-based analysis of 3D deformable shapes: methods and applications

PDF

3D deep learning for perception and modeling

PDF

Deep generative models for image translation

PDF

Deep representations for shapes, structures and motion

PDF

Effective data representations for deep human digitization

PDF

Particle simulations of colloid thruster beam

PDF

Recording, reconstructing, and relighting virtual humans

PDF

Scalable exact inference in probabilistic graphical models on multi-core platforms

PDF

Thermal modeling and control in mobile and server systems

Asset Metadata

Creator Li, Yijing (author)

Core Title Artistic control combined with contact and elasticity modeling in computer animation pipelines

School Andrew and Erna Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Computer Science

Publication Date 07/25/2019

Defense Date 06/07/2019

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag animation,collision,computer graphics,contact,directable simulation,FEM,hierarchical,Modeling,OAI-PMH Harvest,physically based modeling

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Barbic, Jernej (committee chair), Kuo, C.-C. Jay (committee member), Nakano, Aiichiro (committee member)

Creator Email strivinggene@gmail.com,yijingl@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-193994

Unique identifier UC11663590

Identifier etd-LiYijing-7632.pdf (filename),usctheses-c89-193994 (legacy record id)

Legacy Identifier etd-LiYijing-7632.pdf

Dmrecord 193994

Document Type Dissertation

Format application/pdf (imt)

Rights Li, Yijing

Type texts

Source 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 a...

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