Physikalische Grundlagen für Computerspiele

Maximilian Stark (mail@dakror.de), WS2018

Stand: 17.02.2019

Context

Real-time physics

• Interactive -> Simulation
• Accuracy vs. Speed
• Stability
• Applications
  • Training
  • Education
  • Entertainment

Modeling pyramid (bottom up)

1. Geometric
2. Kinematic
3. Physical
4. Behavioral
5. Cognitive

Animation techniques

• Conventional
  • keyframe interpolation
  • Full control
  • Tedious
  • No interaction
• Motion capture
  • recording of markers on actor
  • Real performance
  • Loads of data
  • Bound to realism
  • No real interaction
• Physics-based
  • fully interactive
  • Computationally expensive

Goals

• Immersion
• Gameplay
• Exploration
• Asset creation reduction

Mass spring systems

• Simulation of mesh structure by mass points with springs
• Hooke's law
  • Spring deformation force proportional to deformation distance
  • $\displaystyle F = -k(l-L)$, $\displaystyle k$ Stiffness, $\displaystyle l$ current length, $\displaystyle L$, original length
• Equations of motion (Newton 2nd law): $\displaystyle m_i*a_i = F_i^{\text{int}}(x_i(t)) + F_i^{\text{ext}}(t) - \underbrace{\gamma v_i(t)}_{\text{Damping}}$
• Implementation
  • Euler step
  • Generic force generators
    • linked to single mass point
    • evaluation & accumulation of force
    • + Clean interface
    • - Two spring instances per spring needed
• Simulation
  • Analytic solution
    • ODE: $\displaystyle x''=-\frac{k}{m}x$
    • Harmonic oscillation with damping: $\displaystyle x(t) = e^{-\gamma t}d\cos(\sqrt{\frac{k}{m}}t)$
    • - Only simplified setup
  • Numerical integration
    • From Taylor Series
    • - Explosions
  • Leapfrog
    • 1 force evaluation, 2nd order accuracy
    • Interleaved update of velocity and position
    • $\displaystyle v(t+\frac{h}{2}) = v(t-\frac{h}{2}) + ha(t)$
      $\displaystyle x(t+h) = x(t) + hv(t+\frac{h}{2})$
    • Initialisation of velocity at $\displaystyle \frac{h}{2}$

Numerical integration

• Reason
  • $a(t) = v'(t) = x''(t)$
• Methods
  • Explicit euler
    • $\displaystyle y_{n+1} = y_n + hf(x_n,y_n)$
    • Step in direction of current slope
  • Heun
    • $\displaystyle \widetilde{y} = y_0 + hf(x_0,y_0)$
      $y_1 = y_0 + \frac{h}{2}(f(x_0,y_0) + f(x_1,\widetilde{y}))$
    • Step with avg. tangent from current and from euler step → 2nd order
    • "Predictor-corrector" method
  • Midpoint
    • $\displaystyle \widetilde{y} = y_0 + \frac{h}{2}f(x_0,y_0)$
      $y_1 = y_0 + hf(x_0 + h/2,\widetilde{y})$
    • derivative at half step
    • Runge Kutta 2nd order
  • Runge Kutta
    1. First estimate of slope: $\displaystyle k_1 = f(x_0,y_0)$
    2. Prediction of tangent at midpoint: $\displaystyle k_2 = f(x_0 + \frac{h}{2}, y_0 +\frac{h}{2} * k_1)$
    3. Correction of prediction: $\displaystyle k_3 = f(x_0 + \frac{h}{2}, y_0 + \frac{h}{2} * k_2)$
    4. Prediction of slope at full step $\displaystyle k_4 = f(x_0 + h, y_0 + h*k_3)$
    5. Weighted step taken: $\displaystyle y_1 = y_0 + \frac{h}{6} (k_1 + 2k_2+2k_3+k_4)$
  • Velocity verlet
    • Higher order interpolation instead of multiple evaultions of derivative
    • $\displaystyle y_{n+1} = y_n + hy_n' + \frac{h^2}{2}y_n''+\mathcal{O}(h^3)$
    • $\displaystyle y_{n+1} = y_n' + h \frac{y_n'' + y_{n+1}''}{2} + \mathcal{O}(h^3)$
    • data vs. calculation tradeoff ($y_n''$ stored)
  • Implicit integration
    • Approximation of derivative at next time step instead of current
    • $\displaystyle y_{n+1} = y_n + hf(x_{n+1}, y_{n+1})$
    • Solution of (linear) $\displaystyle Ay_{n+1}=b$
    • + Unconditionally stable
    • - High complexity
    • - Uncontrollability when non-linear
    • - Boring aesthetics
• Accuracy
  • n step 1st-order much worse than 1 step n-order
  • Speed vs accuracy → Visual validity
  • Error bounded by n-th derivative in interval: $\displaystyle 0 \leq e \leq \frac{h^n}{n}y^n$

Rigid Bodies

• Center of mass $\displaystyle x_{cm} = \frac{\sum_i m_i x_i}{\sum_i m_i}$
• Angular velocity $\displaystyle w$
• Angular momentum $\displaystyle L = \sum_i (x_i \times m_i v_i)$
• Torque $\displaystyle q = \frac{\mathrm{d}L}{\mathrm{d}t} = \sum_i (x_i \times f_i)$: Rotational force
• Inertia tensor $\displaystyle i$: Resistance to rotation

2D

• Inertia tensor scalar: $\displaystyle i = \sum_n (m_n x_n \cdot x_n)$
• Change of angular velocity: $\displaystyle w(t+h) = w(t) + h\frac{q}{i}$
• Change of rotation: $\displaystyle r(t+h) = r(t) + hw(t)$
• Simulation
  1. External forces
  2. Euler Step
     1. Position
     2. Velocity
     3. Rotation
     4. Angular velocity
  3. World position
• Collision detection
  • Relative velocity: $\displaystyle v_{rel} = n \cdot (v_A - v_B)$
    • $\displaystyle v_{rel} < 0$: Collision
    • $\displaystyle v_{rel} > 0$: Separation
    • $\displaystyle v_{rel} = 0$: Sliding
• Collision response
  • using forces
    • unpredictable
    • full rigidity → infinite force
  • Impulses
    • Velocity change without time
    • $\displaystyle J = m\Delta v$
    • linear motion: $\displaystyle J_{lin} = \frac{-(1+c)v_{rel}\cdot n}{\frac{1}{M_a}+\frac{1}{M_b}}$
    • full motion: $\displaystyle J = \frac{-(1+c)v_{rel}\cdot n}{\frac{1}{M_a}+\frac{1}{M_b} + \frac{(x_a \times n)^2}{I_a} + \frac{(x_b \times n)^2}{I_b}}$
    • $\displaystyle v_{rel} < 0$ → $\displaystyle J >= 0$ scalar
    • $\displaystyle v_a' = v_a + \frac{Jn}{M_a}$
    • $\displaystyle v_b' = v_b - \frac{Jn}{M_b}$
    • $\displaystyle w_a' = w_a + \frac{(x_a \times Jn)}{I_a}$
    • $\displaystyle w_b' = w_b + \frac{(x_b \times Jn)}{I_b}$
  • Bounciness: coefficient of restitution
    • $\displaystyle c=1$: fully elastic
    • $\displaystyle c=0$: fully plastic
    • Loss of kinetic energy for $\displaystyle c<1$ → deformation of body

      Rotations

• Fixed Angles
  • Rotation about fixed axes
  • Reverse order: $\displaystyle R_{x,y,z} = R_z R_y R_x$
• Euler Angles
  • Axes fixed to object
  • Yaw (z), Pitch (x), Roll (y)
• Quaternions
  • $\displaystyle q = (\cos(\frac{\phi}{2}), n\sin(\frac{\phi}{2}))$, angle $\displaystyle \phi$, axis $\displaystyle n$
  • Rotation of point $\displaystyle p$: $\displaystyle q(0,p)\widetilde{q}=p'$ with $\displaystyle \widetilde{q} = (s, -v)$
  • $\displaystyle q_1 \cdot q_2 = (s_1, v_1) \cdot (s_2, v_2) = (s_1 s_2 - v_1 \cdot v_2, ~ s_1 v_2 + s_2 * v_1 + v_1 \times v_2)$

3D

• Inertia Tensor -$3 \times 3$ matrix, Computation using mass-weighted co-variance matrix of body coordinate positions
  $I = \mathrm{Id} ~ \mathrm{trace}(C)-C; ~ C = \sum_n m_n x_n x_n^\intercal; ~ C_{j,k} = \sum_n m_n, x_{n,j} x_{n,k}$
  • Dependant on current orientation
  • Generally working with inverse
  • Computation from initial values: $\displaystyle I_\text{current} = \mathrm{Rot}_r ~I_0~ \mathrm{Rot}_r^{-1} = \mathrm{Rot}_r ~I_0~ \mathrm{Rot}_r^\intercal$
• Angular motion
  • velocity not constant, momentum constant
  • velocity dependant on rotational axis
  • $\displaystyle L=Iw$
  • $\displaystyle w(t+h) = I^{-1} ~ L(t+h)$
• Orientation
  • Derivative of quaternion: $\displaystyle \frac{\mathrm{d}r}{\mathrm{d}t} = \frac{1}{2}\binom{0}{w}$
  • Integration using $\displaystyle r' = r + \frac{h}{2} \binom{0}{w} r$
  • Rotation matrix integration using shear matrix
• Simulation
  1. External forces
  2. Euler Step
     1. Position
     2. Velocity
     3. Rotation
     4. Angular velocity
     5. Inverse inertia tensor
     6. Angular momentum
  3. World position
• Inapplicability of leapfrog: Forces dependant on position
• Impulse
  • $\displaystyle J = \frac{-(1+c)v_{rel}\cdot n}{\frac{1}{M_a}+\frac{1}{M_b} + \left[ (I_a^{-1}(x_a \times n)) \times x_a + (I_b^{-1}(x_b \times n)) \times x_b \right] \cdot n}$
  • Updates
    • $\displaystyle v_a' = v_a + \frac{Jn}{M_a}$
    • $\displaystyle v_b' = v_b - \frac{Jn}{M_b}$
    • $\displaystyle L_a' = L_a + (x_a \times Jn)$
    • $\displaystyle L_b' = L_b - (x_b \times Jn)$

Fluid Simulation

Navier-Stokes equations, based on Euler

• Velocity $\displaystyle u$
• Pressure $\displaystyle p$
• Density $\displaystyle \rho$
• Conservation of momentum: $\displaystyle \text{advection} = \text{pressure} + \text{viscosity} + g$
• Conservation of mass
  • Divergence: Existance of sink / source

Representations

• Eulerian: grid based property field
  • + Precision
  • + Ease of calculation
  • - Problems with continuous motion
  • Conservation of mass: mass per cell, sum constant
• Lagrangian: particle like

Smoothed Particle Hydrodynamics (SPH)

• Solution of NS equations using particles
• Euler integration of position & velocity
• Goal: Calculation of forces for NS-equation approximation
• Smooth particles
  • Mass as continuous field
  • Smoothing kernel $\displaystyle W$ for interpolation
  • Attribute evaluation: $\displaystyle a_W(x) = \sum_i a_i \frac{m}{\rho_i}W(x-x_i)$
  • Gradient evaluation: $\displaystyle \nabla a_W(x) = \sum_i a_i \frac{m}{\rho_i} \nabla W(x-x_i)$
• Forces
  • Pressure
    • Conversion of density using equation of state
      $p_i = k((\frac{\rho_i}{\rho_0})^7 - 1); ~ k$ Stiffness
    • $\displaystyle f_\text{pressure} (x_j) = - \sum_i \frac{p_i + p_j}{2} \frac{m}{\rho_i} \nabla W$
  • Viscosity
    • $v$ viscosity, internal friction
    • $\displaystyle f_\text{visc}(x_j) = - \sum_i (v_i - v_j) \frac{m}{\rho_i} \nabla^2 W$
  • Color field
    • $\displaystyle c(x_j) = - \sum_i \frac{m}{\rho_i} W$
  • Surface tension
    • $\sigma$
    • $c$
    • $\displaystyle f_\text{st} = -\sigma \nabla^2 c \frac{n}{|n|}$
• Simulation
  1. Computation of grid / kd-tree
  2. Density & color-field
  3. Forces
     1. Pressure
     2. Viscosity
     3. Surface tension
  4. Integration of position & velocity (Euler)
• Issues
  • Collision with obstacles: Clumping
  • No flow orthogonal to walls
  • High pressure forces from clusters
  • Negative pressure: single pressure smaller than of rest
  • Explicit integration: stability problems
  • (Too) many interacting parameters & systems
• Improvement
  • Position based method
    • Velocity based on position update
    • Predictions
    • No iterative force computation

Smoothing kernels

• radius $\displaystyle d$
• Functions
  • Cubic splines, distance to center $\displaystyle q$
  • Polynomial functions
• Key properties
  • Fast evaluation
  • Computation of derivatives
• Improvements
  • Different kernels for different forces
  • Kernel with non-negative 2nd derivative
  • Spiky kernel: high values near origin

Collision detection & response

Broad phase

Culling

• Sweep and Prune
  • SAT with AABBs
  • $\mathcal{O}(n\log n + k)$: sorting ($\mathcal{O}(n\log n)$) + reporting intersection: $\mathcal{O}(k)$
  • Sorting with temporal coherence (insertion sort): $\mathcal{O}(n)$
• Cullide
  • Visibility based
  • No collision if fully in front of them (z-buffer pixel count)
  • $\mathcal{O}(n)$
  • Profits from hardware acceleration

Bounding Volume Hierarchies

• Bounding Volumes (asc quality, desc complexity)
  • Spheres
  • AABBs
  • OBBs
    • axes: Eigenvectors of covariance matrix $C$
    • $\displaystyle \mu = \frac{1}{n}\sum_i x_i; ~ C_{j,k} = \frac{1}{n} \sum_i (x_{i,j}-\mu_j)(x_{i,k} - \mu_k)$
    • SAT using 15 axes: 3x face normals of A, 3x face normals of B, 3*3 edge pairs
  • k-discrete oriented polytopes (k-DOP)
    • k fixed directions
    • k max values
    • SAT on k axes
  • Convex hulls
• Hierarchies
  • Subdivision of volumes to generate hierarchy
  • Bottom-up: from triangles to hierarchy by grouping
  • Top-down: recursive splitting
  • Insertion
    • Iterative addition of primitives
    • Graph building by proximity grouping
    • + Dynamic updates
    • - Quality dependant on insertion order
• Collision tests
  • SAT
  • Voronoi marching
  • GJK algorithm
• Deformable objects
  • Sphere tree
  • Update of leaves only for deformed primitives
  • Information propagation to root
  • + Efficient
  • - Suboptimal BV

Spatial partitioning

• Data structures
  • Uniform grid
  • Quadtree
  • kd-tree
  • BSP-tree
    • Binary space partitioning tree
    • Generalization of kd-tree
    • Creation
      • Selection of plane with arbitrary normal
      • Sorting of geometry into the two spaces
      • Splitting of polygons intersection selected plane
      • Recursive division
      • - time consuming
    • Useful for efficient rendering
• Spatial hashing for self-collisions
  • Heavily twisted deformable objects
  • Volumetric meshes
    • tetrahedal mesh
  • Hash function $H(i,j,k) = (ip_1 \oplus jp_2 \oplus kp_3) ~ \mathrm{mod} ~ n$
    • $i,j,k$: cell coordinates
    • $p_1,p_2,p_3$: large primes
    • $n$: hash table size
  • Cell based partition & hashing
    • Implicit uniform grid
      1. Vertex hashing: {cell -> list of vertices in cell}
      2. Tetraheda hashing
         • tetraheda hashed to all cells touching bounding box
      3. Vertices and tetraheda in same cells tested for intersection
• Vertex-Tetrahedon collision test
  • Barycentric coordinate based
    • All points $P$ in tetrahedon: $P = (b_1, b_2, b_3); ~ b_1+b_2+b_3 = 1; ~ b_i = a_i / (a_1+a_2+a_3)$
    • Intersection if: $\mathrm{coord}(A) = \mathrm{vol}(PBCD) / \mathrm{vol}(ABCD) \in (0,1) ~ \forall$ vertices
  • Oriented faces
    • $n$ face normal
    • $n \cdot (P-A) > 0 ~ \forall$ faces
• Signed-distance functions
  • Sampling of distance to object with dense grid
  • Indication of inside (neg.) / outside (pos.) by sign of value
  • part of level set framework (Niveaumengen)
    • Embedding of surface into higher dimensional representation
      • flat surface -> volumetric function
    • Surface given by level set $\Gamma = \{ x | \phi (x)=0\}; ~ \phi: \mathbb{R}^3\mapsto \mathbb{R}$
    • Implicit modification of surface
  • Calculatation of Closest Surface Point
    • Gradient descent instead of triangle intersection tests
    • Gradient always pointing away from object (Normalization)
  • Data Acquisition
    • Precomputation
    • Triangle rasterization
    • Ray intersections
    • Data from previous stages
  • Pros/Cons
    • + Fast checking
    • + Single precomputation
    • - Lots of memory
    • Improvements using adaptive grid structure

Narrow Phase

Contact & Friction

• Stationary objects
  • Balance of forces
  • Gravity canceled by reaction force from contact
  • → Resting contacts
• Multiple impulses
  • Change of total velocity by every impulse
  • Final relative velocity must be separating:$\Delta v_i + \sum_k J_k \geq 0$
  • Solution in single system of equations: No unique solution
  • → Constraint: Impulses always positive and minimal
  • Quadratic programming problem
  • Linear system matrix as representation of collision influence graph
• Micro Collisions
  • Resting contacts as small collisions
  • Removal of gravity acceleration from last time step
  • Detection of tiny velocities, decrease of restitution coefficient
  • → Cancels out bounce back
• Impulse without gravity
• Friction
  • static: material parameter, solid ground
  • dynamic: two moving, touching objects
  • Impulse with friction: apply impulse & friction (along tangent)

Position correction

• Approaches
  • Teleport to surface of partner
  • Backtrace velocity & orientation
• Account for mass
• Minimal movement
• Based on penetration depth instead of $v_{\text{rel}}$
• Orientation change: possibility of new collisions
  • minimize change
  • e.g proportional to rel. collision point on object
• Address friction

Full RB simulation

1. Integration
2. Collision detection
   • Obtaining list of colliding object pairs
   • Usage of accelerated data structure
   • Minimal face-corner collisions
   • List building over time
     1. Timestep
        • Gather single deepest collections
     2. Timestep
        • Collision verification
        • Check of $v_{\text{rel}}$
3. Penetration resolution
   • Sorted list, descending
   • Impulse based PosRot change
   • In theory: re-check for collisions
   • Usually finite time / attempts for resolution
     • Iterative Convergence
     • Priority based
4. Velocity resolution
   • Sorted list, ascending $v_{\text{rel}}$
   • Recomputation of velocity for all objects colliding with either partner
   • Collision tolerance
     • Ignoring of small collisions
     • Minimal penetration depth & velocity
     • Only visual validity

Sleeping

• Performance boost
• Sleep state basis
  • kinetic energy
  • velocity
  • recent history (EMA: exponential moving average)
• Still participation in CD with 0 velocity
• Waking up
  • Collision with active object
  • Hysteresis (threshold range)

Complex cloth simulation

• Mass spring system with springs for stretching & bending
• Thickness $c$
• Difficulty of tangling prevention
• Collision detection
  • AABB Hierarchy
    • Bounding box for movement
    • Including old & new position
    • Thickness
    • Safety margin
  • Visual validity
  • Fully plastic collisions
  • Position correction: relaxation (factor 0.25)
  • Penetration depth = thickness → Collision ignored, friction-esque
• Heun integration
• Smoothed rendering
  • Visual collision checking based on vertex displacement "velocity"
  • Temporal incoherence