Animations with Particles in Computer Graphics

Lets say we are asked to make an animation of throwing a basketball. How can we do?

• Option 1: Manual Animation
  • Manually move the basketball according to a certain curve
    • Doesn’t look realistic
• Option 2: Physically Based Animation
  • Numerically solve the governing equations to move the basketball
    • New animations automatically look realistic

Particles

• Moving mass points in 2D or 3D space

Particle Systems

• Simulation of many particles at once
• Particles can be:
  • Independent
  • Dependent (Interact with each other)

Particles

The Particle class is about each particle,

and the ParticleSystem is about a collection of particles.

Each particle is described by:

• Position $\vec{x}(t)$
• Velocity $\vec{v}(t)$
• Accumulated force $\vec{F}(t)$
• Mass $m$

where $t$ is the time.

Given the Accumulated force and the mass, we want to find the position.

Forces

• Gravitational $\vec{F}(t)=m \vec{g} \rightarrow \frac{\mathrm{d}^2 \vec{x}(t)}{\mathrm{d} t^2}=\frac{m \vec{g}}{m}=\vec{g}$
• Viscous drag $\vec{F}(t)=-C \vec{v}(t)$
• Spring force $\vec{F}(t)=C\left(\vec{x}_o-\vec{x}(t)\right)$

In general, we cannot mathematically solve the equations of motion

• We numerically estimate them using Time Discretization.

• We cannot use continuous time on a computer

• We consider discrete time instead

  • $t \approx t_n = n\Delta t$
  • where $\Delta t$ is called “time step”

Using Verlet Integration:

• Has less numerical issues comparing to Euler’s method
• Based on the central difference approximation

$$\frac{\vec{x}_{n+1}-2 \vec{x}_n+\vec{x}_{n-1}}{\Delta t^2} \approx \frac{\vec{F}_n}{m}$$

$$\vec{x}_{n+1} \approx \vec{x}_n+\left(\vec{x}_n-\vec{x}_{n-1}\right)+\Delta t^2 \frac{\vec{F}_n}{m}$$

newPosition = position + (position - prevPosition) + deltaT*deltaT*globalGravity;

Collisions and constraints

Constraints:

• The things to restrict the motion of each particles
  • e.g. Walls, obstacles, other particles etc.
    • Collision against a wall
    • Bead on a wire

We can use Force-based or Position-based to solve the problem.

• Force-based constraints
  • Exact solutions
  • Sometimes difficult to apply
  • In Implementation:
    • Check if there is any collision
    • then solve for constraint forces
• Position-based constraints
  • Easy and general
  • May not correspond to the real physics
  • In Implementation:
    • Check if there is any collision
    • then remove collisions by moving points

Collision against a Wall

Force-based

General recipe:

• Reflect the velocity (i.e., adjust the force)
• Push the position to resolve the collision
• Reflect the velocity by the collision boundary (e.g., the wall surface)

Position-based

General recipe:

• Check if there is any collision
• Reflect the previous position
• Push both the current and previous positions to resolve the collision

Example: a wall at $x = K$ where $K$ is a constant

We dont want x pass through K that $x < K$ .

Visualization from.

$K = 0$ in this case.

Check for collision:

condition: old position didnt pass through the wall and the new position pass through the wall

$$\vec{x}_n . x>K \text { and } \vec{x}_{n+1} . x<K$$

is the condition “old position didnt pass through the wall” necessary?

Collision response:

Update the old position with (new position - old position) and update new position with the value at wall

$$\begin{aligned} \vec{x}_n . x & \leftarrow C_R(-\vec{x}_n . x+\vec{x}_{n+1} . x) \\ \vec{x}_{n+1} . x & \leftarrow K \end{aligned}$$

where . means the attribute of the vector in code here

$C_R$ can be incorporated by scaling the distance

Bead on a wire

Force-based

Solve for constraint force

• Constraint on position
• Constraint on kinetic energy

Position-based

• Project the position back to the circle wire

$$\vec{x}_{n+1}=\vec{c}+r \frac{\vec{x}_{n+1}-\vec{c}}{\left\|\vec{x}_{n+1}-\vec{c}\right\|}$$

where $\vec{c}$ is the position of center, $r$ is the radius.

It is basically denominator is basically normalizing the thing.

Multi-body Dynamics

Interacting particles for more advanced effects

Collision Processing

Key aspect of any dynamics simulation

• Collision detection
• Collision response

Collision Detection

Detect if there is any interpenetration occurred

Lets consider two circle with center $\vec{x_i}$ and $\vec{x_j}$

Collision can be detected if $\| \vec{x_i} - \vec{x_j} \| \leq r_i + r_j$

static bool collisionDetection(const Particle part_i, const Particle part_j, const float sphereRadius_i, const float sphereRadius_j){
    float centerDistance = length(part_i.position - part_j.position);
    float radiusSum = sphereRadius_i + sphereRadius_j;
    if(centerDistance <= radiusSum){return true;}
    return false;
}

Simple Pseudocode for simultaneous collisions:

for some number of times {
  for all the pairs of objects (i, j) {
  	if (collision(i, j)) { 
    resolveCollision(i, j)
  	} 
  }
}
# Do this in each time step

Note that checking all the combinations detections is costly.

There are some methods to make the detection efficient.

• Basically using acceleration data structure

Collision Response: Force-based

• New velocities based on laws of conservation

Energy Conservation does by definition

Collision Response: Position-based

• Just push circles away from each other lmao

Note that the overlap distance = $(r_i + r_j) - \| \vec{x_i} - \vec{x_j} \|$

Then push each other

Energy Conservation in Position-based Response

• Modify the previous locations according to the conservation of energy

Gravitational fields

Newton’s law of gravity

$$\vec{F}_{i \rightarrow j} = Gm_im_j\frac{\vec{x}_j - \vec{x}_i}{\| \vec{x}_j - \vec{x}_i \|^3}$$

where $G$ is the gravitational constant (Not $g$)

Note that $\| \vec{x}_j - \vec{x}_i \|^3$ may become a very small value. Always add a small value to the distance.

static float3 gravityField(const Particle part_i, const Particle part_j) {
    float G = 9.8f;
    float mass_i = 0.1f, mass_j = 0.1f;
    float eps = 0.001f;
    float3 distance = normalize(part_j.position - part_i.position + eps);
    float3 F_ij = (G*mass_i*mass_j) * distance/ (length(distance) * length(distance) * length(distance));
    return F_ij;
}

^Note that it really depends on how you gonna tweak the parameters.

Pseudocode of Basic implementation (Do this in each time step before time integration):

• Double loop over all the particles
  • For each particle, add the gravitational force by all the others
• ^Repeat for all particles
• O(N^2) for N particles

for all the particles i { 
  	accumulatedForce[i] = 0
  
    for all the other particles j { 
      F_ij = gravityField(i, j) 
      accumulatedForce[i] += F_ij
    } 
}

My implementation:

void step() {
  float3 accumulatedForce[globalNumParticles];
  for (int i = 0; i < globalNumParticles; i++) {
          accumulatedForce[i] = float3(0.0f);
          for (int j = 0; j < globalNumParticles; j++) {
              if(j == i) continue;
              float3 F_ij = gravityField(particles[i], particles[j]);
              accumulatedForce[i] += F_ij;
          }
  }
  for (int i = 0; i < globalNumParticles; i++) {
      particles[i].step(accumulatedForce[i]);
  }
  updateMesh();
}