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)
- Geometric
- Kinematic
- Physical
- Behavioral
- 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
- 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
- First estimate of slope: $\displaystyle k_1 = f(x_0,y_0)$
- Prediction of tangent at midpoint: $\displaystyle k_2 = f(x_0 + \frac{h}{2}, y_0 +\frac{h}{2} * k_1)$
- Correction of prediction: $\displaystyle k_3 = f(x_0 + \frac{h}{2}, y_0 + \frac{h}{2} * k_2)$
- Prediction of slope at full step $\displaystyle k_4 = f(x_0 + h, y_0 + h*k_3)$
- 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
- External forces
- Euler Step
- Position
- Velocity
- Rotation
- Angular velocity
- 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
- External forces
- Euler Step
- Position
- Velocity
- Rotation
- Angular velocity
- Inverse inertia tensor
- Angular momentum
- 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
- Computation of grid / kd-tree
- Density & color-field
- Forces
- Pressure
- Viscosity
- Surface tension
- 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
- 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
- Vertex hashing: {cell -> list of vertices in cell}
- Tetraheda hashing
- tetraheda hashed to all cells touching bounding box
- 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
- 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
- Integration
- Collision detection
- Obtaining list of colliding object pairs
- Usage of accelerated data structure
- Minimal face-corner collisions
- List building over time
- Timestep
- Gather single deepest collections
- Timestep
- Collision verification
- Check of $v_{\text{rel}}$
- 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
- 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