國
立
交
通
大
學
多媒體工程研究所
碩
士
論
文
即 時 頭 髮 模 擬 使 用 之 有 效 率 單 線 動 力 學
Efficient Strand Dynamics for Real-Time Hair Simulation
研 究 生：王上文
指導教授：林文杰 教授
Efficient Strand Dynamics for Real-Time Hair
Simulation
Student: Shang-Wen Wang
Advisor: Dr. Wen-Chieh Lin
Institute of Multimedia Engineering
College of Computer Science
National Chiao Tung University
ABSTRACT
Simulating hair motion is difficult due to their high-stiffness property, which usu-ally causes numerical stability problem and needs the numerical simulation be performed at a small time step. We observed that a force applied to one part of a strand will transport to other parts of the strand where the force can contribute the deformation of the strand. We model this force transport process and simplify the special Cosserat Theory to simulate inextensible hair strands. Our experi-ments show that the proposed approach can simulate hair dynamics robustly and efficiently such that hair simulation with user interaction can be executed in real time.
Acknowledgments
I would like to thank my advisor, Dr. Wen-Chieh Lin for his guidance. Thanks to my colleagues in the CG Group.: Fang-Shuo Li, Kuei-Li Fang, Chin-Hsian Chang, Hsin-Hsiao Lin, Kuang-Wei Fu, Ying-Tsung Li, Yueh-Tse Chen, Zhi-Han Yan, Ming-Shang Hsu, Jia-Ru Lin, Ying-Shou Lan, Yi-Jheng Huang, Wei-Chien Cheng,and Hsi-Chou Huang for their assistance and discussion.
Contents
1 Introduction 1 2 Related Work 3
2.1 Spring-Damper Model . . . 3
2.2 Rigid-body Chain Model . . . 4
2.3 Cosserat Model . . . 4
3 Strand Simulation 5 3.1 The Cosserat Model and Special Cosserat Theory . . . 5
3.2 Representing the Rod State in Reduced Coordinate . . . 7
3.3 Force Transport . . . 8
3.3.1 Tangential Force . . . 9
3.3.2 Normal Force . . . 9
3.4 Discretization of Strand Dynamics Equation of the strand . . . . 12
3.5 Dynamics Integration . . . 14
4 Hair Simulation 15 4.1 Hair Clustering . . . 15
4.2 Guide Hair Interpolation with GPU . . . 16
4.3 Collision Detection and Response . . . 16
5 Experiment and Result 18 5.1 Single Strand Test . . . 18
5.1.1 Statics Test . . . 18
5.1.2 Dynamics Test . . . 19 iii
5.1.3 Curliness and Twist Test . . . 20 5.1.4 Performance . . . 21 5.2 Hair Simulation . . . 22 5.2.1 Result . . . 22 5.2.2 Performance . . . 27 6 Conclusion 30 7 Appendix 32 iv
List of Figures
3.1 Rod Representation . . . 6 3.2 Tangential Force and Normal Force: The component of Fin
paral-lel to d3 is the tangential force, and the other is the normal force. . 10
3.3 Transport of Tangential Force, the component parallel to d3 of
force in segment i will be transported to proximal segment i − 1. . 10 3.4 Normal Force and Torque: d1 and d2 are the axes of the local
coordinate on centerline. FN is the normal force, and tf is the
local force transformed from torque. FN is a normalized vector of
FN. The neutralized part is max(kFNk, −tf ·FkF N
Nk)FN . . . 11
5.1 Statics Test: Young’s modulus is set as 10MPa, and shear modulus is 5MPa. The error of our force transport may increase as the segment numbers increases. . . 19 5.4 The Curliness Phenomena by twist: After we applied a twist torque
to the whole strand and self-collision, we have a phenomena that the strand entwined. . . 20 5.5 The strands entwined by rotating constrained point. . . 21 5.6 Under Wind: hair deformed by a constant wind force. . . 23 5.7 Simulated Hair Motion with Head Shaking: We simulate the head
as a rotating rigid body. Then we get a dynamics of head shaking. 24 5.8 Hair Moving: (a)(b)(c) the result that head moves vertically. (d)(e)(f)the
result that head moves horizontally. . . 25 5.9 The difference between different stiffness in moving: (a) stiff (b)
soft . . . 26
5.10 The difference of hair after head shaking between different inter-action parameter (a)with D = 2 × 10−5 (b) with D = 2 × 10−6. Other interaction parameters are set to K = 2 × 10−6, Pinner =
3mm. The result with high viscous force is more cluttered. . . 26 5.2 Dynamics Result: Our method is blue. The CORDE is green. We
can observe our result of dynamics is more stiff. This is because we assume the rest segments will move together in the simplifica-tion of Special Cosserat Theory. There is a big difference, but it is still acceptable in the hair simulation. . . 28 5.3 Curly Strand Dynamics . . . 29
List of Tables
C H A P T E R
1
Introduction
Realistic hair has been studied for years, including modeling, simulating, and rendering. Simulating hair strands is hard because of their high stiffness property. Simulating such hair usually surfers from numerical stability problem. It needs a very small time step to simulate high stiffness of hair.
The focus of this thesis is the transport of force in a hair strand with one fixed end and one free end. Finding where the force contributes is helpful to the stabil-ity of strand simulation, especially in high stiffness model, because force will not take so many time steps to transport. In experiment, our method indeed gives fast stability and good performance in hair simulation and make it be real-time.
For high-stiffness strands, we observed that forces applied to one segment will transport to root segment by segment to find where to contribute deformation. In a zero stiffness model, force is transported along the tangent of the strand, which is called the axial force. In a high stiffness model, force is not only transported by the axial force but also the stiffness torque which generate by material properties. The speed of transport through energy method is bounded by simulation time step and spatial step,though the speed is fast in high stiffness material. Too slow force transport also reduces stability of simulation. This is the bottleneck of traditional energy methods, like the spring-damping model. Finding the Jacobian matrix to map a local force to global may help, but update the matrix is a huge overhead.
2
We observed that hair strands is a constrained object rather than a free-deform object.The difference of hair shape at each time step is small. That means the motion of hair strand is more similar to a rigid body than a free rigid body chain. Therefore we assume that there will be no large local change in hair shape between two time step, and hair strands tend to have a global motion. Shape and stiffness torques will be similar between two time step. This make it able to compute the distribution of forces from its current shape and stiffness torques.
We then use the the Special Cosserat Theory[1] to model the curliness of strands; However solving the Special Cosserat Theory is hard and computation-ally expensive[2], which is not suitable for real time simulation. We did some simplification to make it suitable for our force distribution.
C H A P T E R
2
Related Work
There have been much work on hair simulation, including hairstyle modeling, hair dynamics, and hair rendering. It is beyond the scope of this thesis to survey all literature here. As our goal is to simulate hair motion, we will only review the related work on the physics-based model of single hair strand motion in this chapter. For a thorough survey on hair simulation, please see [3].
2.1
Spring-Damper Model
The spring-damper model is a popular model for simulating elastic objects for its simplicity. In this model, a hair strand is often modeled as a chain of connected particles[4][5]. This method is fast and easy to implement, but is not suitable for simulating high-stiffness objects, like hair. Since the stiff dynamics usually re-quire using small time step, hair motion simulation is usually a bottleneck of the whole system. Additionally,simulating torsion or curliness of hair is also a hard problem. Although Selle et al.[6] simulates torsion and curliness successfully, their solution for spring stretch still makes it lack of physical accuracy.
2.2 Rigid-body Chain Model 4
2.2
Rigid-body Chain Model
In order to simulate high stiffness of hair, the Rigid-body Chain model is more appreciated. The Rigid-body Chain removes stretching DOFs of strand, makes each strand inextensible[7][8]. The Rigid-body Chain model needs more compu-tation than the Spring-Damper Model, but still an O(n) algorithm[9]. Although Rigid-body Chain solve the high stiffness problem of the Spring-Damper model, it still has problem in simulating curliness.
2.3
Cosserat Model
Cosserat Rod is first introduced to the computer graphics field by Pai [10] in 2002. The kinematic state of a Cosserat rod is generated by integrating the spatial speed and spatial angular speed over the rod, which can preserve local coordinate at every point of the whole strand. Therefore, Cosserat Rod model can model tor-sion and tortor-sion stiffness. Pai used the Cosserat model to get only the statics of a whole strand, not the dynamics.
Bertails et al. (Super Helices)[11] simulated curly hair dynamics using the Cosserat Rod. They assumed a curly hair strand can be discretized into several segments. This method preserves torsion of curly hair and improves computa-tional efficiency with reduced DOFs. When the hair length is long, external force is applied or hair strands interact with other objects, curly hair needs to be simu-lated using more DOFs. The computational cost grows rapidly as O(n2) where n
is the DOFs of the hair strand.
Spillmann et al. also proposed a 1-D elastic model [12] based on Cosserat Rod. They suggested the strand is composed of many finite elements. The dynamics of the model can be solved in O(n). Unlike the Super-Helices, the method proposed by Spillmann et al. simulated each control point on the strand. A drawback of their method is that the stability of the method will be affected by increasing num-ber of segments when we model high stiffness of hair.
C H A P T E R
3
Strand Simulation
In this chapter, we describe our strand simulation method. First, we briefly review the Cosserat Model and Special Cosserat Theory. Second, we describe our reduced coordinate that makes the conversion between torque and force pos-sible and reduces the computational time. Third, we describe our forces transport method. Then we rewrite the Special Cosserat Theory to obtain the dynamics equation. Finally we solve this equation numerically to simulate hair motion.
3.1
The Cosserat Model and Special Cosserat
The-ory
Cosserat Rod was first introduced to the computer graphics field by Pai [10] in 2002. The dynamics of the rod is governed by an integration of positional variation and orientational variation1over the rod:
d30 = u × di i = 1 to 3, (3.1)
where di is the local axes on the rod as shown in Fig 3.1, and u is the Darboux
vector. In order to model the unstretchedness, the tangent along the center line of 1Notation: We use f0to denote spatial derivative along the center line of a strand, f0 = ∂f /∂s,
and ˙f = ∂f /∂t to denote time derivative.
3.1 The Cosserat Model and Special Cosserat Theory 6
Figure 3.1: Rod Representation
the rod is:
r0 = d3 (3.2)
where r is the center line position of the strand. If we integrate r0 along the center line, we can get the position of the Cosserat Rod (Strand):
r(l) = Z l
0
d3(s)ds + r(0) (3.3)
where l is the distance from the root along the center line and r(0) is the position of the root of the strand.
Because the Cosserat Strand can model torsion and curliness of strands, we use this model as our hair strand model.
The Special Cosserat Theory uses two equations to describe the dynamics of the strand. One is for position:
n0+ df = ρA¨r (3.4) where n0 is the constraint force generated by the material law, df is an external force applied to an infinitesimal segment, ρ is the linear density, and A is the area
3.2 Representing the Rod State in Reduced Coordinate 7 of the cross section. This is basically Newton’s Second Law. Because n is af-fected by the shape of strand and torque, n usually does not have a simple formula to compute[2]. The other equation of the Special Cosserat Theory describing the orientation of the rod is:
m0+ r0× n + dτ = ρI ˙ω (3.5) where m0 is the torques of stiffness generated by the material law, dτ is an ex-ternal torque applied on infinitesimal segment, I is the moment of inertia of an infinitesimal segment,and ω is the angular velocity.
In the Special Cosserat Theory, there are two unknowns, n and ˙ω)in the equa-tion. The resulting acceleration is affected by the force and torque simultaneously. From Equation 3.4, we know that a force that do not contribute to the linear acceleration of the strand will become an internal force. Furthermore, we know that the component of a force that does not induce any effective torque at a point of the rod will also become an internal force. Based on these two findings, we compute the effects of any external force applied to the rod by propagating the force segment by segment. We consider this problem as a force transport problem by finding where the force contributes to effective torque. However, different from the Special Cosserat Theory, this force will be the summation of many segments’s force. We need to consider the affected mass point when we compute dynamics. After computing transport forces, we use an approximated equation of Equation 3.5 to compute strand dynamics.
3.2
Representing the Rod State in Reduced
Coordi-nate
In order to represent the 3D position and 4D quaternion orientation of an in-finitesimal segment on a rod and their derivatives, a 14-dimension vector is needed to represent the state of the strand. However, based on the Cosserat Model’s def-inition, we can compute the position and velocity of each point on the centerline using Equation 3.3.
Also, from the rigid body dynamics we know: ˙
3.3 Force Transport 8 where ω is the angular velocity.
Therefore the speed on center line, which is the temporal differentiation of equation 3.3, is ˙r(l) = Z l 0 ˙ d3ds = Z l 0 (ω × d3)ds (3.6)
The temporal derivative of ˙d3 is
¨
d3 = ( ˙ω × d3+ ω × ˙d3) = ( ˙ω × d3+ ω × (ω × d3))
so the linear acceleration on the center line is ¨ r(l) = Z l 0 ¨ d3ds = Z l 0 ( ˙ω × d3+ ω × (ω × d3))ds (3.7)
which means the position of a segment on the rod can be obtained from the ori-entation of the segment if we add the centrifugal force to the external force. The centrifugal force will be described in the appendix. Therefore, we only need an 8-dimensional vector to represent the state of the rod.
3.3
Force Transport
From the classical mechanics, we know that a force f applies to a point with distance lsto the rotation center can be converted into a torque τ :
lsd3× F = τ = I ˙ω (3.8)
This means a non-twisting force can be transformed into a bending torque. There-fore, we can find where the force results a torque by checking shapes and torque on one strand. This equation makes our force transport possible.
Beginning from the tip, we find where force becomes a torque segment by seg-ment using Equation 3.8. In a strand segseg-ment, force can be decomposed into two components,tangential force and normal force, as shown in Fig 3.2. We observed that there are two situations that force will be transported to proximal segment:
1. Tangential Force
A force component, which is a tangential force, results no effective torque. It will be transported to next segment.
3.3 Force Transport 9 2. Normal Force which results no ineffective Force
A force component, which neutralized by a applied torque, will result no effective torque. It will be transported to next segment.
We will describe these two cases below.
3.3.1
Tangential Force
It is obvious that a tangential force does not induce any torque on the segment, and is just transported to other segments as shown in Fig 3.3. If a tangential force is transported to a segment that can induce a torque to the segment, we remove it from the transported force and set is as a torque to bend the segment.
The net force Fin at position r(l) is the sum of the external force applied to
the position r(l) and the tangential force transported from the position r(l + ds). Fin(l) = (fex(l) + FT(l + ds)) (3.9)
where fex is the external force, and FT(l) = (Fin(l) · d3(l))d3(l) is the tangential
force at r(l). Fin(l) is the summation of external force and transported axial force.
Thus, we can obtain the normal force FN(l) at r(l).
FN(l) = Fin(l) − FT(l) (3.10)
The normal force will result a torque at r(l).
3.3.2
Normal Force
Recall our assumption in Section 3.1 that the component of force which can not result a effective torque will become a internal force. When torques involves, we need to compute how much the effective torque resulted by FN. We call the
component of force which results a torque as fef f ective,
fef f ective(l) = FN − Fneutralized(l)
Fneutralized is the component of force neutralized by stiffness torque, which can
be written as:
3.3 Force Transport 10
Figure 3.2: Tangential Force and Normal Force: The component of Finparallel to
d3 is the tangential force, and the other is the normal force.
Figure 3.3: Transport of Tangential Force, the component parallel to d3 of force
3.3 Force Transport 11
Figure 3.4: Normal Force and Torque: d1 and d2 are the axes of the local
co-ordinate on centerline. FN is the normal force, and tf is the local force
trans-formed from torque. FN is a normalized vector of FN. The neutralized part is
3.4 Discretization of Strand Dynamics Equation of the strand 12 where c is a constant and ν is a normalized vector of FN(l).
In order to compute c, we transform the total torque τ applied here to a force to find FN’s contribution.
tf = T × d3/ls (3.12)
where lsis the length of segment. Then the length of FN neutralized by tf is:
c = 0, if (projLength(FN(l), tf (l)) > 0 min(−projLength(FN(l), tf (l)), kFN(l)k) otherwise (3.13) projLength is: projLength(FN(l), tf (l)) = tf (l) · FN(l) kFN(l)k (3.14)
The force neutralized by torque is also shown in Fig 3.4. We need to rewrite Equation 3.9 for Fneutralized:
Fin(l) = f (l) + FT(l + ds) + Fneutralized(l + ds) (3.15)
After integrating from tip to root, we get a series of FN, which is the force
make the strand deform.
3.4
Discretization of Strand Dynamics Equation of
the strand
We describe how we discretize the Special Cosserat Theory in this section. As we only concern the orientation of the rod, we do not need to discretize Equation 3.4. Recall that the orientation equation of the Special Cosserat Theory is:
m0+ r0× n + dτ = ρI ˙ω If we discrete it, we get
3.4 Discretization of Strand Dynamics Equation of the strand 13
∆mi+ (d3,i× ni) · ∆s + τi = ρI ˙ωi· ∆s (3.16)
Equation 3.16 is only for single segment case. In practice, the torque applied to a segment will result accelerations of other segments and decrease the accelera-tion of itself. The resulting acceleraaccelera-tion is d3× n. Since n does not have a simple
formula , we approximate n by considering segments that are accelerated by τi.
First we need to know how many segments will be accelerated by the torque. We observe that in a high-stiffness strand, the local deformation caused by other segments’ acceleration will be canceled by stiffness very fast. This is similar to a rotating rigid body.
So we consider the number of segments being accelerated by the force. We assume that the accelerated segments are located at the distal side of the segment. The total mass will be x·mass which x is the number of the accelerated segments. If we consider the force contributed to acceleration of x · mass as (x · mass)¨r. Then the dynamics equation will become:
(d3· ls) × (F − (x · mass)¨r) = ρI ˙ω
where F = τ × d3 is the bending force transformed from torque, and ls is the
length of segment. Then we write ¨r as angular acceleration: d3× F = ρI ˙ω + (x · mass)(d3× ω × d3)(ls)2
If we consider it in bending direction (torque on d1or d2):
d3× F = (ρI + (x · mass) · (ls)2) ˙ω
Therefore the bending moment of inertia will be I +(x·mass)·(ls)2. This equation
is the same with Equation 3.16 where n is approximated by −(x · mass)¨r · ls.
3.5 Dynamics Integration 14
3.5
Dynamics Integration
From above section, we know the approximation of Equation 3.5 will have this form:
∂ ˙Iωi
∂t = τtotal,i
where τtotalis
τtotal,i= FT,i × d3+ τistif f ness+ τ external
i + τ
damping
i (3.17)
We use Spillmann’s method to compute the stiffness torqueτstif f ness. The detail
is describe in the appendix.
τdamping is
τdamping = −Dq· ω (3.18)
where Dqis the damping constant. Then we can solve ˙ωi:
˙ ωi = τtotal,i− ωi× Iωi (3.19) Then integrate ω: ωt+hi = I−1(τtotal,i−ωt i×Iωi)h + ω t i (3.20)
Here we use a method similar as Spillmann’s semi-implicit method[12] but remove dimensions of position.
First, transform ω to quaternion velocity ˙q: ˙ qi = 1 2Qi 0 w (3.21) where Q is the quaternion matrix,and q is the orientation of the segment i. Then we can integrate it on quaternion:
ˆ qit+h= 1 2Q t i 0 wt+hi h + qit (3.22) qit+h= qˆ t+h i kˆqit+hk (3.23) After this integration, we have the simulation result of hair strands.
C H A P T E R
4
Hair Simulation
A whole hairstyle usually has tens of thousands of hair strands. We describe here how we apply strand dynamics model to hair simulation, and implement interaction between hair to hair and hair to other objects.
4.1
Hair Clustering
In order to simulate the motion of styled hair which have been designed, we cluster modeled hair strands to find relations between guide hairs and hair strand geometry.
We choose some guided strands from hair model first. Then we compute av-erage distances E between a guide hair g and a strand geometry h as an error function: E(g, h) = 1 N N X i=1 minDis(g, vi) (4.1)
where minDis is the minimal distance between g and vi. vi is the i-th vertex of
h, and n is the number of vertices of h. Then we use the K-means clustering algo-rithm in textbook [13] to find the cluster of hair model. There are two guide hair for each strand.
4.2 Guide Hair Interpolation with GPU 16
4.2
Guide Hair Interpolation with GPU
We use the vertex shader to accelerate our hair interpolation. From the Section 4.1 , we have got the relation of guide hair and hair geometry. In order to interpo-late vertex between two guide hairs, we need the weights between guide hair for each vertex of geometry. We compute the weight as:
w(g, h) = 1 E(g, h)
f oreachg
X 1
E(g, h) (4.2) where g is a guide hair, and h is a strand geometry.
After we compute weight of each vertex of hair geometry for corresponding guide hair using Equation 4.2, we put them into texture coordinate of hair ge-ometry. We also put the corresponding guide hair index into texture coordinate. The above tasks are completed in off-line. At each render step in runtime, we put the vertex data of guide hair which we simulated in strand simulation part into texture’s rgb channel. Then we interpolate vertex between multiple guide hair in vertex shader, using the method in the GPU Gems for skinning.
4.3
Collision Detection and Response
We use the interaction method proposed by Choe [14] to handle collision de-tection and collision response, which handles collision dede-tection and response with two cylinders. Their method uses the outer cylinder to handle damping ef-fects between hair cluster and use the inner cylinder to avoid deep intersection between hair cluster. Different from their method, we use the 1D Swept Sphere Volume(SSV) as our bounding volume instead of cylinder for the concern of per-formance. We use a user-define inner SSV radius and a outer SSV radius from hair cluster. The collision detection algorithm for SSV is described in the text-book by Ericson[15].
The damping force and resist force are describes in Equation 4.3 and 4.4.
4.3 Collision Detection and Response 17 where Fs is the resist force to push two clusters away, b is the bounding volume,
K is the spring parameter, the Pinner is the penetration depth of inner SSV, which
is generated from collision detection, and the rest is the force direction. r is the position of SSV.
Fviscous(bi, bj) = −D ∗ Pouter∗ ( ˙r(bi) − ˙r(bj)) (4.4)
where Fviscousis the viscous force, D is the damping parameter of hair interaction,
the Pouteris the penetration depth of outer SSV, and the rest is the relative velocity
C H A P T E R
5
Experiment and Result
We compute the stiffness torque from Young’s Modulus by using the method proposed by Spillmann[12]. The detail of Young’s Modulus and stiffness torque is described in the appendix.
5.1
Single Strand Test
5.1.1
Statics Test
In order to check the accuracy of statics, we applied gravity to strands and record its final shape. The static result is compared with CORDE, which is shown as Fig 5.1. We set parameters of strand property as:
Radius 5 mm
Density 1300 kg/m3 ρ
The result shows that the statics is almost the same with the CORDE. The force transport we use results a little force lost. The deform becomes smaller with more segments. But from the result, we know the lost of force is very small.
5.1 Single Strand Test 19
(a) CORDE (b) Our method with 10 segs
(c) Our method with 50 segs (d) Our method with 100 segs
Figure 5.1: Statics Test: Young’s modulus is set as 10MPa, and shear modulus is 5MPa. The error of our force transport may increase as the segment numbers increases.
5.1.2
Dynamics Test
Here we verify our dynamics model by releasing strand from a force, which is compared with CORDE. The result is shown in Fig 5.2. The CORDE result is col-ored by green, and our method is blue. There is a big difference between CORDE and our method. Remember our assumption in Section 3.4.1 that the deform of shape without a force will be absorbed fast in high stiffness object, which is one of properties of high stiffness object. As our assumption, we have the result of dynamics which looks like a high stiffness object. The curly strand dynamics is also shown in Fig 5.3
5.1 Single Strand Test 20
5.1.3
Curliness and Twist Test
In order to make sure our model has the advantage of the Cosserat Model that can deal with curliness. We use 50 segments to simulate a 1m strand. The Young’s and shear modulus are 200KPa. The result is shown as Fig 5.4. The average FPS of simulation is 47.6.
Also we do a twist experiment with many strands. We control 16 strands with constrained tip to make them entwined. The average FPS is 20. 5.1.3
(a)
Figure 5.4: The Curliness Phenomena by twist: After we applied a twist torque to the whole strand and self-collision, we have a phenomena that the strand entwined.
5.1 Single Strand Test 21
(a) (b)
(c) (d)
Figure 5.5: The strands entwined by rotating constrained point.
5.1.4
Performance
Our test environment is Intel Core Duo 2 6700 + GeForce 8800
The computing time of simulation and self-collision, and the stability test com-pared to CORDE are shown below. The stretching Young’s Modulus in CORDE is set to 1MPa. We test the stability by applying gravity and release.
The stability of our method is better than CORDE. This is because we remove the stretching DOF and approximate of the Special Cosserat Theory. Removing the stretching DOF make it more stable when segment number increases and our assumption removes a local deformation by force which may make it unstable.
5.2 Hair Simulation 22 Computational time in one step
Segment # Simulation Time(ms) CORDE Simulation Time Self-Collision Time(ms) 50 0.051 0.058 0.235
100 0.087 0.108 0.270 200 0.170 0.22 0.472 500 0.443 0.548 2.049
The stability of single strand with length 5m
Segment # Young’s Modulus (MPa) Max Time Step(ms) CORDE Max Time Step 50 1 7.0 3.1 50 10 2.4 1.6 50 100 0.7 0.5 100 100 0.3 0.2 200 100 0.15 0.1
5.2
Hair Simulation
5.2.1
Result
We use our strand model to simulate some animation of hair. Since there are not many self contact in hair simulation, we turn off the self-collision to get better performance. The Young’s Modulus of the material parameters are set to 1MPa. The hair under wind is as Fig 5.6. The wind is modeled as a constant wind multi-plied by a sin function. The wind force is set to 10−4kg.
5.2 Hair Simulation 23
(a) (b)
(c)
Figure 5.6: Under Wind: hair deformed by a constant wind force.
The dynamics results from head shaking and moving are as Fig 5.7 and Fig 5.8. The head dynamics is modeled by rigid body dynamics.
5.2 Hair Simulation 24
(a) (b)
(c) (d)
Figure 5.7: Simulated Hair Motion with Head Shaking: We simulate the head as a rotating rigid body. Then we get a dynamics of head shaking.
5.2 Hair Simulation 25
(a) (b)
(c) (d)
(e) (f)
Figure 5.8: Hair Moving: (a)(b)(c) the result that head moves vertically. (d)(e)(f)the result that head moves horizontally.
5.2 Hair Simulation 26 The difference motion of different Young’s Modulus is shown as Fig 5.9. We can see the stiff hair style has less deformation when head moves.
(a) (b)
Figure 5.9: The difference between different stiffness in moving: (a) stiff (b) soft The difference resulted by D in hair interaction is presented in Fig 5.10.
(a) (b)
Figure 5.10: The difference of hair after head shaking between different interac-tion parameter (a)with D = 2 × 10−5 (b) with D = 2 × 10−6. Other interaction parameters are set to K = 2 × 10−6, Pinner = 3mm. The result with high viscous
5.2 Hair Simulation 27
5.2.2
Performance
We does above demo with simulating 100 guide hair of about 1000 segments. The time step of simulation and collision detection is 3ms. The average simula-tion time is 0.64ms, the average rendering time is 1.4ms, and the average collision detection and response time is 2.3ms. The average FPS is 23.2
5.2 Hair Simulation 28
(a) (b) (c) (d)
(e) (f) (g) (h)
(i) (j) (k) (l)
(m) (n) (o) (p)
Figure 5.2: Dynamics Result: Our method is blue. The CORDE is green. We can observe our result of dynamics is more stiff. This is because we assume the rest segments will move together in the simplification of Special Cosserat Theory. There is a big difference, but it is still acceptable in the hair simulation.
5.2 Hair Simulation 29
(a) (b)
(c) (d)
C H A P T E R
6
Conclusion
We have proposed a dynamic model for simulating high stiffness model more efficiently. Our method is faster, more stable than traditional method, and also can model the curliness phenomena of strand. With our proposed method, we make the simulation of hair with full interaction able to be performed in real time. But there are still some problems needed to be solved.
Limitation
Our method is based on a shape analysis. It means the shape of simulated strand is needed to be continuous in time and space. It needs enough stiffness to keep strand continuous. In the case which we consider about, the shape is con-strained strongly. But in low stiffness case, the shape may be very discontinuous in time. In this case, it will result some obvious artifacts in dynamics.
Future Work
We can observe from the statics experiment that the force transport method has still some force lost in each integration. As the segment number increases, the error increases. Besides, from the performance result, we can observe that the
31 computational time for handling interaction is still the bottleneck in hair simula-tion. We would like to solve the hair interaction problem in the future.
C H A P T E R
7
Appendix
Centrifugal Force
Before searching axial force, we consider the internal force generate by angu-lar speed from Equation 3.7. We need to add it into the external force to compute
˙
ω. If we multiply a mass on Equation 3.7’s both side, we get: F = mass · ¨rl = mass ·
Z l
0
( ˙ω × d3+ ω × (ω × d3))ds
where ω × (ω × d3) is the centripetal force. What we want to find is the relation
between force and angular acceleration, so we move it to the left side: F − mass · Z l 0 ω × (ω × d3)ds = mass · Z l 0 ( ˙ω × d3)ds (7.1)
The left side is total force.
Stiffness
The potential energy of quaternion proposed by the CORDE: Vk = li 2Kkk(Bk(qj + qj+1) · 1 lj (qj+1− qj) − ˆuk)2 (7.2) 32
33 where K11= K22= E πr2 4 , K33= G πr2 2
where E is the Young’s Modulus for bending stiffness, and G is the Shear Modulus for torsional stiffness.
We can do symbolic differentiation to get the force τq on the quaternion:
τq = 3 X k=1 ∂ ∂qj Vk where ∂ ∂qj Vk = li 2Kkk ∂ ∂qj Bk(qj+ qj+1) · 1 lj (qj+1− qj) − ˆuk 2 = liKkk· Bk lj ∂ ∂qj (qj + qj+1) (qj+1− qj) · Bk(qj + qj+1) · 1 lj (qj+1− qj) − ˆuk = liKkk· Bk lj ((qj+1− qj) − (qj + qj+1)) · Bk(qj+ qj+1) · 1 lj (qj+1− qj) − ˆuk = liKkk Bk lj (−2qj) Bk(qj + qj+1) · 1 lj (qj+1− qj) − ˆuk = −2liKkk Bk lj qj· Bk(qj + qj+1) · 1 lj (qj+1− qj) − ˆuk
But the mass to quaternion is a dense matrix,which is hard to inverse. We use [16] method to transform a quaternion force to a torque.
0 τ = 1 2Q Tτ q (7.3)
Bibliography
[1] Stuart S. Antman. Nonlinear Problems of Elasticity. Springer New York, 1995.
[2] Lay T. Weber A. Sobottka, G. Stable integration of the dynamic cosserat equations with application to hair modeling. In International Conferences in Central Europe on Computer Graphics, Visualization and Computer Vision, 2008.
[3] Kelly Ward, Florence Bertails, Tae-Yong Kim, Stephen R. Marschner, Marie-Paule Cani, and Ming C. Lin. A survey on hair modeling: Styling, simulation, and rendering. IEEE Transactions on Visualization and Com-puter Graphics, 13(2):213–234, 2007.
[4] E. Plante, M.P. Cani, and P. Poulin. A Layered Wisp Model for Simulating Interactions inside Long Hair. In Computer Animation and Simulation 2001: Proceedings of the Eurographics Workshop in Manchester, UK, September 2-3, 2001. Springer, 2001.
[5] F. Bertails, TY Kim, MP Cani, and U. Neumann. Adaptive Wisp Tree: a multiresolution control structure for simulating dynamic clustering in hair motion. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics sym-posium on Computer animation, pages 207–213. Eurographics Association Aire-la-Ville, Switzerland, Switzerland, 2003.
[6] Lentine M. G. Selle, A. and R. Fedkiw. A mass spring model for hair simu-lation. In SIGGRAPH ’08: ACM SIGGRAPH 2008 Papers, 2008.
Bibliography 35 [7] S. Hadap and N. Magnenat-Thalmann. Modeling Dynamic Hair as a Contin-uum. In Computer Graphics Forum, volume 20, pages 329–338. Blackwell Synergy, 2001.
[8] J.T. Chang, J. Jin, and Y. Yu. A practical model for hair mutual interactions. In Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 73–80. ACM New York, NY, USA, 2002. [9] R. Featherstone. Robot Dynamics Algorithms. Kluwer Academic Publishers,
1987.
[10] Pai D.K. Strands: Interactive simulation of thin solids using cosserat models. Computer Graphics Forum, 21:347–352(6), September 2002.
[11] Florence Bertails, Basile Audoly, Marie-Paule Cani, Bernard Querleux, Fr´ed´eric Leroy, and Jean-Luc L´evˆeque. Super-helices for predicting the dy-namics of natural hair. In SIGGRAPH ’06: ACM SIGGRAPH 2006 Papers, pages 1180–1187, New York, NY, USA, 2006. ACM.
[12] J. Spillmann and M. Teschner. Corde: Cosserat rod elements for the dy-namic simulation of one-dimensional elastic objects. In SCA ’07: Proceed-ings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 63–72, Aire-la-Ville, Switzerland, Switzerland, 2007. Eu-rographics Association.
[13] S. Theodoridis and K. Koutroumbas. Pattern Recognition. Academic Press, 2006.
[14] B. Choe, M.G. Choi, and H.S. Ko. Simulating complex hair with robust col-lision handling. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 153–160. ACM New York, NY, USA, 2005.
[15] C. Ericson. Real-Time Collision Detection. Morgan Kaufmann, 2005. [16] AL Schwab and JP Meijaard. HOW TO DRAW EULER ANGLES AND
UTILIZE EULER PARAMETERS. In Proc. IDETC/CIE 2006, ASME 2006 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, 2006.
