How to make a PIXAR movie: Deformable Bodies
Daniel Lagler∗
Technische Universität München
Abbildung 1: Deformierbare Körper in Computeranimationsfilmen (Toy Story 3, Ice Age 3, Oben)
Zusammenfassung
Computeranimierte Filme gibt es jetzt schon seit einigen Jahren.
Die ersten Filme waren noch sehr einfach gestaltet, doch mit der
raschen Entwicklung der Prozessoren ließen sich immer realistischere Szenen berechnen. Aber mit dem steigenden Realismus erhöhte sich auch der Modellierungsaufwand, weshalb es für einige Anwendungen sinnvoll geworden ist, von manueller Steuerung
zu automatischer Simulation überzugehen. Um deformierbare Körper darzustellen, werden aus diesem Grund meistens physikalische
Simulationen eingesetzt. Sie werden oft mithilfe von Feder-Masse
Systemen oder der Finite Elemente Methode berechnet. Diese beiden Techniken werden näher betrachtet.
Stichwörter: Deformierbare Körper, Feder-Masse Systeme, Finite
Elemente Methode
1
Einleitung
Die Idee, Körper, die ihre Form oder ihr Volumen ändern, in
Computer-Animationsfilmen einzusetzen, existiert schon seit geraumer Zeit [Lasseter 1987]. Allerdings bedurfte es mehr, als nur
dem Wissen alleine, diese realistisch umzusetzen. Früher wurden
eher kinematische Aspekte von Bewegungen betrachtet, heutzutage
verlegt sich der Schwerpunkt immer mehr auf physikalische Simulationen. Deformierbare Körper stellen im Bereich der Computergrafik physikalische Körper dar, die sich unter Krafteinwirkung verändern. In Computeranimationsfilmen werden sie mittlerweile für
Kleidung ([Goldenthal et al. 2007], [Waggoner and Baraff 2007]),
Haare ([Petrovic et al. 2008], [Selle et al. 2008]), Gesicht, Seile und
andere elastische Körper verwendet. Allerdings erfordert die Simulation von solchen Körpern umfangreiche Kenntnisse im Gebiet der
Dynamik, einem Teilbereich der technischen Mechanik. Sie stellt
etliche physikalische Zusammenhänge bereit, welche als Grundlagen für heute verwendete Simulationsverfahren dienen. Zwei dieser
Techniken, Feder-Masse Systeme und deformierbare Körper unter
Einsatz der Finite Elemente Methode, werden in den nachfolgenden
Kapiteln vorgestellt.
∗ e-mail:
[email protected]
2
Feder-Masse Systeme
Feder-Masse Systeme bestehen aus einer beliebigen Anzahl an
Punktmassen, die über elastische Federn miteinander verbunden
sind [Möllenhoff 2010].
ks
xj
xi
L
mi
mj
Abbildung 2: Feder-Masse System
In Abbildung 2 sei ein einfaches Feder-Masse System dargestellt.
Das System besteht nur aus zwei Punktmassen (Positionen xi , xj ∈
R3 , Massen mi , mj ∈ R+ ) und einer Feder (Ruhelänge L ∈ R+ ,
Federkonstante ks ∈ R+ ), welche die beiden verbindet. Der Verbindungsvektor der Punkte ergibt sich aus dij = xj − xi . Damit
resultiert aus dem Hooke’schen Gesetz eine Berechnungsvorschrift
für die Kraft fij , die Punktmasse j auf i auswirkt [Müller 2008].
fij =
dij
· ks · (|dij | − Lij )
|dij |
(1)
Daraus folgt, dass die Kräfte fij und fji betragsmäßig gleich sein
müssen und die Richtung entgegengesetzt ist fij = −fji .
vi
xj
xi
mi
vj
ks
kd
L
mj
Abbildung 3: Feder-Masse-Dämpfer System
Allerdings sind solche Federn für Anwendungen eher unpraktikabel, da sie, wenn sie einmal zum Schwingen angeregt wurden,
unendlich lang weiterschwingen. Deswegen werden in der Praxis
häufig gedämpfte Federn eingesetzt. Der Grad der Dämpfung wird
über die Dämpfungskonstante kd ∈ R+ geregelt. Mithilfe der Geschwindigkeiten vi , vj ∈ R3 der Punktmassen lässt sich die Dämpfungskraft errechnen. Diese kann nun zu der Federkraft einfach aufaddiert werden, um die Gesamtkraft zu bekommen.
dij
dij
· ks · (|dij | − Lij ) + kd · (vj − vi ) ·
(2)
fij =
|dij |
|dij |
Die Erweiterung auf ein komplexeres System mit beliebig vielen
Federn und N ∈ N Punktmassen ist einfach. Für einen Massenpunkt i müssen nur alle Kräfte resultierend aus den angrenzenden
Federn aufsummiert werden. Zusätzlich können hier auch noch externe Kräfte Fiext berücksichtigt werden.
Fi =
N
X
fij + Fiext
Um schließlich von der Kraft Fi zu der Beschleunigung ai der
Punktmasse zu kommen, wird das zweite Newton’sche Gesetz herangezogen F = m·a, wodurch die Beschleunigung für Punktmasse
i folgendermaßen ermittelt werden kann.
Fi
mi
(4)
Durch Umformung der gewöhnlichen Differentialgleichung ai =
v˙i = x¨i lässt sich die Geschwindigkeit für einen Zeitpunkt t + ∆t
dann folgendermaßen (5) bestimmen. Für die Position kann eine
ähnliche Vorschrift (6) aufgestellt werden.
Z t+∆t
vi (t + ∆t) = vi (t) +
ai (τ ) dτ
(5)
t
Z
t+∆t
xi (t + ∆t) = xi (t) +
vi (τ ) dτ
(6)
t
Da diese Integrale allerdings oft nur sehr schwer bis unmöglich
analytisch auszuwerten sind, wird die Lösung mittels Methoden
der Numerik approximiert. Ein einfaches numerisches Integrationsverfahren stellt der explizite Euler dar. Dabei wird die Zeit ∆t in
Schritte mit der Länge ∆τ zerlegt, über diese approximativ eine
konstante Beschleunigung wirkt, was zu dem anschließend gezeigten Iterationsverfahren (7) führt. Der Zusammenhang zwischen Position und Geschwindigkeit lässt sich wieder analog (8) ermitteln.
vit+1 = vit + ∆τ · ati
(7)
xt+1
= xti + ∆τ · vit
i
(8)
Die Numerik bietet noch jede Menge andere Integrationsverfahren,
die alle ihre Vor- und Nachteile bezüglich Rechenaufwand, Genauigkeit, Komplexität und Stabilität haben. Beispiele dafür sind impliziter Euler, Runge-Kutta und Verlet-Integration.
2.2
3D
Abbildung 4: Gitter-Typen
3
Zeitintegration
ai =
2D
(3)
j=1
2.1
1D
Unterschiedliche Arten von Gittern
Abbildung 4 zeigt verschiedene Möglichkeiten für die Erstellung
von Gittern. Bei eindimensionalen Strukturen, einer Kette aus
Punktmassen und Federn, können die Federn praktisch frei um die
Massenpunkte rotieren und so zusammenklappen oder sich beliebig
lang um ihre Achse drehen. Dieses Verhalten ist meist unerwünscht.
Deswegen werden zusätzliche Federn zur Versteifung benötigt, um
eine Biege- und Eindrehresistenz zu gewährleisten. Bei zweidimensionalen Strukturen, einem Netz aus Punktmassen, verhält es sich
ähnlich. Zusätzliche Versteifungen werden benötigt, um Überfaltungen zu vermeiden. Es gibt jedoch etliche verschiedene Ansätze,
diese Federn zu wählen. Im Allgemeinen wird dies anwendungsspezifisch entschieden.
Elastizitätstheorie
Für die Simulation von deformierbaren Objekten mithilfe der Finite Elemente Methode werden Zusammenhänge aus der Elastizitätstheorie verwendet. Bei gegebener Ausgangsstellung Ω ⊂ R3
kann das deformierte Modell durch die Verschiebungsfunktion
u(x) : Ω → R3 beschrieben werden. Die Verschiebungsfunktion ordnet jedem Punkt x ∈ Ω einen Vektor zu, der die Differenz
von Ausgangs- und deformierter Stellung angibt. Die deformierten
Punkte ergeben sich dann durch {x + u(x)|x ∈ Ω}. Die potentielle Energie eines Systems im statischen Gleichgewicht kann durch
folgende Gleichung beschrieben werden [Hauth 2004]:
1
Π=
2
Z
Z
T
ε σ dx −
Ω
Z
T
f T u ds = min
g u dx −
Ω
(9)
∂Ω
Π gibt dabei die potentielle Energie als Funktional der Verschiebung u an. Der erste Term repräsentiert die elastische Energie, die
der Körper speichert. Die beiden anderen Terme beschreiben die
Arbeit, die durch Wirkung von Volumenkräften g und Oberflächenkräften f aufgebracht wird. Das System befindet sich im Gleichgewicht, wenn die Verschiebung u bei gegebenen externen Kräften g
und f die potentielle Energie Π minimiert. So ergibt sich die tatsächliche Lösung durch Minimierung von Gleichung 9. und σ
repräsentieren dabei 6-dimensionale Vektoren, welche die Komponenten des symmetrischen Dehnungstensors ε und symmetrischen
Spannungstensors Σ enthalten.
= (1 , . . . , 6 )T = (ε11 , ε22 , ε33 , ε12 , ε13 , ε23 )T
σ = (σ1 , . . . , σ6 )T = (Σ11 , Σ22 , Σ33 , 2Σ12 , 2Σ13 , 2Σ23 )T
Der Green-Lagrange-Dehnungstensor beschreibt die Beziehung
zwischen Deformation und Verschiebung. Diese Beziehung ist im
Allgemeinen nicht-linear. Anschaulich gibt der Dehnungstensor
die relativen Verschiebungen der Flächen eines infinitesimal kleinen achsenorientierten Würfels in allen Raumrichtungen an [Lagler
2009].
ij =
1
2
∂ui
∂uj
+
∂xj
∂xi
+
3
1 X ∂uk ∂uk
2
∂xi ∂xj
(10)
k=1
Der Spannungstensor gibt analog die Kräfte auf die Würfelflächen
für die einzelnen Raumrichtungen an (Normalkräfte, Scherkräfte).
Für isotrope, elastische Körper gilt für die Beziehung zwischen
Dehnungs- und Spannungstensor das Hooke’sche Gesetz:
Σ=λ
3
X
i=1
!
ii
· I3,3 + 2µ
(11)
Für die Vektoren und σ kann das Hooke’sche Gesetz auf eine
Multiplikation mit der Matrix D überführt werden σ = D · ε.

λ + 2µ
 λ

 λ
D=


λ
λ + 2µ
λ

λ
λ
λ + 2µ






4µ
4µ
(12)
λ und µ werden Lamé-Koeffizienten genannt. Diese beiden können
mithilfe des Elastizitätsmoduls E und der Querkontraktion (Poissonzahl) ν bestimmt werden.
Eν
(1 + ν) (1 − 2ν)
E
µ=
2 (1 + ν)
(13)
d
∆l
∆d
Kraft
Abbildung 5: Querkontraktion (Poissonzahl)
4
Nd−1 (x)
(14)
l
∆d
d
∆l
l
Mithilfe der sogenannten shape-Matrix kann diese Interpolation
e
als Matrix-Multiplikation
geschrieben werden f (x) = Φ(x)f .
e
T
T
f = f0 , · · · , fd−1 gibt die Linearisierung der Vektoren fi eines Finiten Elements an.
T

N0 (x)
N0 (x)



N0 (x) 




..
..
..

(16)
Φ(x) = 


.
.
.


Nd−1 (x)



N
(x)
d−1
Das Elastizitätsmodul ist ein Maß für die Steifigkeit eines Materials. Die Querkontraktion gibt das Verhältnis der relativen Dickenänderung ∆d
zur relativen Längenänderung ∆l
an.
d
l
ν=−
(15)
i
4µ
λ=
innerhalb eines Finiten Elements ausgewertet werden.
X
u (x) =
Ni (x) · ui
Tetraeder werden häufig für Finite-Elemente Simulationen verwendet. Um die Koeffizienten cij , 0 ≤ i, j ≤ d der shape-Funktion zu
bestimmen, wird aus den Interpolationsbedingungen
1
i=j
Ni (vj ) =
(17)
0
i 6= j
ein lineares Gleichungssystem aufgestellt. In den Koeffizienten cij
sind selbst die shape-Funktionen für die komplexeren Finiten Elemente linear. Deshalb können diese Bedingungen als Gleichungssystem mit S ∈ Rd×d formuliert werden. ei ist hier der i-te Einheitsvektor. Die Koeffizienten cij können damit effizient über die
Inverse von S bestimmt werden.
Sci = ei 0 ≤ i ≤ d
cij = S −1 ji
Finite Elemente Methode
Finite Elemente sind diskrete Elemente mit einer festen Interpolationsfunktion als Basiseigenschaft. Im Gegensatz zu Finite Differenzen Methoden, die Werte nur an diskreten Abtastpunkten betrachten, berücksichtigen Finite Elemente Methoden das Kontinuum. Deswegen erreichen Finite Elemente Methoden im Allgemeinen eine bessere Genauigkeit. Heutzutage werden Finite Elemente Methoden in den unterschiedlichsten Gebieten verwendet, beispielsweise in der Strukturanalyse, für akustische Berechnungen, in
der Wärmeleitung, zur Flüssigkeitssimulation und zur Berechnung
von elektrischen und magnetischen Feldern [Georgii 2007].
4.1
(18)
(19)
Lineare Approximation der Dehnung
Im Folgenden wird gezeigt, dass durch lineare Approximation des
Dehnungstensors dieser aus einer einfachen Matrixmultiplikation
mit den Vertex-Verschiebungen hervorgeht. Gemäß 20 können die
Verschiebungsvektoren innerhalb des Finiten Elements aus den Verschiebungsvektoren an den Eckpunkten Ui ∈ R3 berechnet werden
[Georgii 2007].
u(x) =
d−1
X
Ni (x)Ui
(20)
i=0
Tetraeder
Serendipity Tetraeder
Hexaeder
u(x) wird nun in den linearisierten Dehnungstensor eingesetzt.
1 ∂ui
∂uj
ij =
+
(21)
2 ∂xj
∂xi
Die einzelnen Komponenten des Dehnungstensors sind somit lineare Terme in Ui . Dadurch lässt sich als Matrixmultiplikation
ε = B · ue
Abbildung 6: Finite Elemente
Finite Elemente können dabei verschiedene Formen (Tetraeder, Serendipity-Tetraeder, Hexaeder, ...) annehmen. Ein Finites Element
mit d Freiheitsgraden, wird durch seine Eckpunkte vj , 0 ≤ j < d
beschrieben. Basierend auf dieser Form können dann sogenannte
shape-Funktionen Ni : R3 → R definiert werden, die eine Interpolation innerhalb des Finiten Elements beschreiben. Angenommen
es seien Vektoren fi ∈ R3 für jeden Eckpunkt (Vertex) des Finiten Elements gegeben, dann kann f ∈ R3 für jeden Punkt x ∈ R3
e
(22)
T
(U0T , . . . , Ud−1
)T
schreiben, wobei u =
die Linearisierung des
Vektors Ui eines Finiten Elements darstellt. Die Matrix B ∈ R6×3d
enthält die partiellen Ableitungen der shape-Funktion Ni (x). Lösungen u(x) müssen die potentielle Energie minimieren. Mit dem
für Finite Elemente konstruierten Ansatz folgt die Bestimmungs∂Π
gleichung aus ∂u
e = 0.
Der elastische Anteil von 9 berechnet sich wie folgt:
Z
Z
∂ 1
T
T
ε Dεdx =
B DBdx ue = K e · ue
∂ue 2 Ω
Ω
(23)
Die Steifigkeitsmatrizen für die einzelnen Finiten Elemente werden
anschließend unter Berücksichtigung der Nachbarschaftsbeziehungen zu einer globalen Steifigkeitsmatrix zusammengesetzt, was zu
folgendem Gleichungssystem führt:
F =K ·x
4.2
(24)
Externe Kräfte
Dynamik
Bisher wurde das Elastizitätsproblem als statisch angenommen. Bei
dynamischen Vorgängen spielen aber Dämpfung und Massenträgheit eine nicht vernachlässigbare Rolle. Zeitliche Ableitungen müssen also berücksichtigt werden. Dieser Zusammenhang wird durch
die Lagrange’sche Bewegungsgleichung beschrieben.
F = K · u + C · u˙ + M · u
¨
(26)
M wird Massenmatrix und C Dämpfungsmatrix genannt. Die Massenmatrix für ein Finites Element M e lässt sich folgendermaßen
bestimmen:
Z
Me =
ρΦT (x) Φ (x) dx
(27)
Ω
Φ stellt dabei die shape-Matrix und ρ die Dichte des Finiten Elements dar, die innerhalb des Finiten Elements variieren kann. Für
die Dämpfung gibt es unterschiedliche Ansätze. Die Dämpfung
kann proportional zur Massenmatrix C e = αM e gewählt werden. α gibt die Dämpfungskonstante an. Eine allgemeinere Variante der Dämpfung (Rayleigh Dämpfung) kann durch einen zusätzlichen Summanden, der proportional (Konstante β) zur Steifigkeitsmatrix K e ist, C e = αM e + βK e erreicht werden. Die
globale Dämpfungs- und Massenmatrix können analog zur Steifigkeitsmatrix aus den einzelnen Matrizen für die Finiten Elemente
zusammengesetzt werden. Die in der Lagrange’schen Bewegungsgleichung enthaltenen Ableitungen nach der Zeit werden oft mithilfe von Finite-Differenzen-Verfahren approximiert, wodurch eine
einfachere Lösung der Bewegungsgleichung ermöglicht wird.
5
BARAFF , D., W ITKIN , A., UND K ASS , M. 2003. Untangling
Cloth. In Proceedings of SIGGRAPH 2003, Annual Conference
Proceedings, 862–870.
G EORGII , J. 2007. Real-Time Simulation and Visualization of
Deformable Objects. Dissertation, Technische Universität München.
Drei verschiedene Arten von externen Kräften können berücksichtigt werden: Volumenkräfte, Oberflächenkräfte und Punktkräfte.
Gravitation, als Beispiel für eine Volumenkraft, kann durch Gradi∂
entenbildung ∂u
e aus dem zweiten Term des Potentials berechnet
werden:
Z
Z
∂
T
e
g
Φ
(x)
u
dx
=
g T Φ (x) dx
(25)
∂ue Ω
Ω
4.3
Literatur
Ausblick
Die Möglichkeiten für den Einsatz von deformierbaren Körpern in
Computeranimationsfilmen sind praktisch unbegrenzt. Immer realistischere Szenen erfordern es, physikalische Simulationsmethoden stetig weiterzuentwickeln. Stabilität spielt ebenfalls eine entscheidende Rolle, ob eine Technik Verwendung findet oder nicht.
Um den Realismus und die Effizienz zu erhöhen, beschäftigen sich
zahlreiche Forschungsgruppen mit Volumenerhaltung [Irving et al.
2008], Kollisionserkennung und -reaktion ([Baraff et al. 2003],
[Govindaraju et al. 2006]), Parallelisierung und speziellen Implementierungen, die einen Teil der Berechnungen auf die GPU auslagern. Einen weiteren interessanten Schwerpunkt stellt die Simulation von Zustandsübergängen (fest, flüssig, gasförmig) [Müller et al.
2004] und die Zerstörung von Körpern dar. Bis physikalische Simulationen allerdings in allen denkbaren Anwendungen einen Grad
der Realitätsnähe erreichen, sodass das menschliche Auge nicht
mehr Realität von Simulation unterscheiden kann, wird es jedoch
noch einige Zeit dauern.
G EORGII , J., 2008. Simulation und Animation.
G IBSON , S., UND M IRTICH , B., 1997. A Survey of Deformable
Modeling in Computer Graphics.
G OLDENTHAL , R., H ARMON , D., FATTAL , R., B ERCOVIER , M.,
UND G RINSPUN , E., 2007. Efficient Simulation of Inextensible
Cloth.
G OVINDARAJU , N., K NOTT, D., JAIN , N., K ABUL , I., TAM STORF, R., G AYLE , R., L IN , M., UND M ANOCHA , D., 2006.
Interactive Collision Detection between Deformable Models
using Chromatic Decomposition.
H AUTH , M. 2004. Visual Simulation of Deformable Models. Dissertation, Eberhard-Karls-Universität Tübingen.
I RVING , G., S CHROEDER , C., UND F EDKIW, R., 2008. Volume
Conserving Finite Element Simulations of Deformable Models.
K AUFMANN , P., M ARTIN , S., B OTSCH , M., UND G ROSS , M.
2008. Flexible Simulation of Deformable Models Using Discontinuous Galerkin FEMs. In Proceedings of Eurographics 2008,
Symposium on Computer Animation.
L AGLER , D. 2009. Skelett-basierte Deformationen. Bachelorarbeit, Technische Universität München.
L ASSETER , J. 1987. Principles of Traditional Animation applied to
3D Computer Animation. In Proceedings of SIGGRAPH 1987,
Computer Graphics 21(4), 35–44.
M EZGER , J. 2008. Simulation and Animation of Deformable Bodies. Dissertation, Eberhard-Karls-Universität Tübingen.
M ÖLLENHOFF , T. 2010. Deformable Objects. Seminararbeit,
Technische Universität München.
M ÜLLER , M., K EISER , R., N EALEN , A., PAULY, M., G ROSS ,
M., UND A LEXA , M. 2004. Point Based Animation of Elastic,
Plastic and Melting Objects. In Proceedings of Eurographics
2004, Symposium on Computer Animation, 141–151.
M ÜLLER , M., 2008. Real Time Physics.
N EALEN , A., M ÜLLER , M., K EISER , R., B OXERMAN , E., UND
C ARLSON , M. 2005. Physically Based Deformable Models in
Computer Graphics. In Proceedings of Eurographics 2005, State
of The Art Report.
P ETROVIC , L., H ENNE , M., UND A NDERSON , J. 2008. Volumetric Methods for Simulation and Rendering of Hair. Pixar
Technical Memo.
RYU , D. 2009. 500 Million and Counting: Hair Rendering on
Ratatouille. Pixar Technical Memo.
S ELLE , A., L ENTINE , M., UND F EDKIW, R. 2008. A Mass Spring
Model for Hair Simulation. In Proceedings of SIGGRAPH 2008,
Annual Conference Proceedings.
WAGGONER , C., UND BARAFF , D. 2007. Virtual Tailoring for Ratatouille: Clothing the Fattest Man in the World. Pixar Technical
Memo.