Spring-Mass Cloth simulation

Overview

A piece of cloth can be considered as a grid of particles and several springs that connect different particles. Springs act as a constraint on distance: applying a force to push away objects that are too close, or applying tension between objects that are too far. For this project, I created cloth that consists of mass nodes connected by springs and dampers and explore the effect of different stiffness and damping settings on simulation speed, stability and cloth appearance.

Duration: 3 weeks

Role: Programmer

Team: Solo

Platform: OpenGL

Techniques

Particle

class Particle{
public:
// ...

Vector3f m_Position;
Vector3f m_Velocity;
Vector3f m_ForceAccumulator;
float m_Mass;
};

Spring Force

class SpringForce : public Force {
// ...
public:
virtual void apply_force();

private:
Particle *const m_p1; // particle1
Particle *const m_p2; // particle2
double const m_dist;  // rest length
double const m_ks, m_kd; // spring constants
};

To apply force to the two particles, we can use the equation $f_1=-[k_s(||p_1-p_2||-l_0)+k_d\cdot\frac{(v_1-v_2)\cdot(p_1-p_2)}{||p_1-p_2||}]\cdot\frac{p_1-p_2}{||p_1-p_2||}\\ f_2 = -f_1$ which is implemented by following code:

// Compute and apply forcce to the two particles
void SpringForce::apply_force() {
float norm = L2Norm(m_p1->m_Position - m_p2->m_Position);
float product = DotProduct(
m_p1->m_Velocity - m_p2->m_Velocity,
m_p1->m_Position - m_p2->m_Position
);

Vector3f force = -(m_ks * (norm - m_dist) + m_kd * (product / norm))
* (m_p1->m_Position - m_p2->m_Position) / norm;

m_p1->m_ForceAccumulator += force;
m_p2->m_ForceAccumulator += -force;
}

Cloth

First, a piece of cloth should be formed by creating a grid of particles.

Cloth::Cloth() {
// ...
// Create a grid of particles to form a piece of cloth
for (int i = 0; i < width; i++) {
for (int j = 0; j < width; j++) {
Vector3f pos = center;
pos.x += i * offset;
pos.z += j * offset;
pVector.push_back(new Particle(pos));
}
}
}

Then, particles should be connected to each other by spring-dampers. There are three types of spring-dampers that are shown below:

• stretch springs (also known as structural springs)
• shear springs
• bend springs

All these types can be easily implemented by code:

// Stretch springs
for (int i = 0; i < width; i++) {
for  (int j = 0; j < width; j++) {
Particle *current = pVector[i * width + j];
if (i < width - 1) {
Particle *bottom = pVector[(i + 1) * width + j];
fVector.push_back(new SpringForce(current, bottom, dist, ks, kd));
}

if (j < width - 1) {
Particle *right = pVector[i * width + j + 1];
fVector.push_back(new SpringForce(current, right, dist, ks, kd));
}
}
}

// Shear springs
for (int i = 0; i < width; i++) {
for (int j = 0; j < width; j++) {
Particle *current = pVector[i * width + j];
if (i < width - 1 && j < width - 1) {
Particle *right_bottom = pVector[(i + 1) * width + j + 1];
fVector.push_back(new SpringForce(current, right_bottom, dist * sqrt(2), ks, kd));
}

if (i < width - 1 && j > 0) {
Particle *left_bottom = pVector[(i + 1) * width + j - 1];
fVector.push_back(new SpringForce(current, left_bottom, dist * sqrt(2), ks, kd));
}
}
}

// Bend springs
for (int i = 0; i < width; i++) {
for (int j = 0; j < width; j++) {
Particle *current = pVector[i * width + j];
if (i < width - 2) {
Particle *bottom = pVector[(i + 2) * width + j];
fVector.push_back(new SpringForce(current, bottom, dist * 2, ks, kd));
}

if (j < width - 2) {
Particle *right = pVector[i * width + j + 2];
fVector.push_back(new SpringForce(current, right, dist * 2, ks, kd));
}
}
}   

The actual spring structures in run time are shown as follows:

Finally, the cloth is simulated and updated frame by frame using explicit euler method.

void Cloth::simulation_step() {
// Clear force accumulators for all the particles
// And add gravity to all of them
for (int i = 0; i < pVector.size(); i++) {
pVector[i]->clear_force();
pVector[i]->apply_force();
}

// Apply sprinf forces to all the particles
for (int i = 0; i < fVector.size(); i++) {
fVector[i]->apply_force();
}

// Implement hard constraints forces
unsigned int fixedParticleIndex[] = {
0,
width - 1,
(width - 1) * width,
width * width - 1,
};

for (int i = 0; i < sizeof(fixedParticleIndex) / 4; i++) {
pVector[fixedParticleIndex[i]]->clear_force();
}

// Handle integration
euler_step();
}

• Handle contact with the environment and self-collisions
• Explore controlled cloth

Integrators

Overview

Intergration is one key step of simulation. There are different types of integration techniques that have different effect on simulation speed and stability. For this project, I applied different integrators on the spring-mass cloth simulation implemented before and compared the difference.

Duration: 3 weeks

Role: Programmer

Team: Solo

Platform: OpenGL

Techniques

Explicit Euler

void Cloth::euler_step() {
for(int i = 0; i < pVector.size(); i++) {
pVector[i]->m_Position += dt * pVector[i]->m_Velocity;
pVector[i]->m_Velocity = DAMP * pVector[i]->m_Velocity
+ dt * pVector[i]->m_ForceAccumulator / pVector[i]->m_Mass;
}
}

Semi-implicit Euler / Symplectic Euler

void Cloth::symplectic_step() {
for(int i = 0; i < pVector.size(); i++) {
pVector[i]->m_Velocity = DAMP * pVector[i]->m_Velocity + dt * pVector[i]->m_ForceAccumulator / pVector[i]->m_Mass;
pVector[i]->m_Position += dt * pVector[i]->m_Velocity;
}
}

Verlet

void Cloth::verlet_step() {
for(int i = 0; i < pVector.size(); i++) {
Vector3f temp = make_vector<float>(pVector[i]->m_Position);
pVector[i]->m_Position = pVector[i]->m_Position + DAMP * (pVector[i]->m_Position - pVector[i]->m_PreviousPos) + (pVector[i]->m_ForceAccumulator / pVector[i]->m_Mass) * dt * dt;
pVector[i]->m_PreviousPos  = temp;
}
}

RK4

https://gafferongames.com/post/integration_basics/

• Implement implicit Euler method

Motion Capture Editing

Overview

Very detailed motion data can be obtained from motion capture, after which motion edting and filtering is common and important. In this project, I created three kinds of edits of there walk cycles in ASF/AMC format, making exaggerting, damping and cartoon effects.

Duration: 3 weeks

Role: Programmer

Team: Solo

Platform: OpenGL

Techniques

File parsing

The AMCViewer already provides all the functions that parse data from the ASF/AMC files. ASF file contains basic information of the skeleton, including mass, unit length, skeleton hierarchy, etc.

// mass, unit length, unit of angle
:units
mass 1.0
length 1.0
angle deg
// info of root and the first element of the bone list
// dof, position, orientation, etc
:root
order TX TY TZ RX RY RZ
axis XYZ
position 0 0 0
orientation 0 0 0
:bonedata
begin
id 1
name lhipjoint
direction 0.678412 -0.646833 0.348374
length 2.5146
axis 0 0 0  XYZ
end
...
// hierarchy
// the first one is the parent joint, the others are children joints that are attached to it
:hierarchy
begin
root lhipjoint rhipjoint lowerback
lhipjoint lfemur
lfemur ltibia
...
end

Data retreiving

The Library::init() function scans the given directory for AMC/ASF files and loads them into data, which is wrapped into classes like Library::Motion, Character::Pose, Character::Angles. The basic data retreiving process is shown as the following.

for (int index = 0; index < Library::motion_count(); index++) {
// Retreive all the motion from the file
Motion & motion = Library::motion_nonconst(index);
for (int j = 0; j < motion.frames(); j++) {
Character::Pose pose;
Character::Angles angles;

// Retreive pose and angle data
motion.get_pose(j, pose);
pose.to_angles(angles);

// Calculate new angles by applying some motion filters
Character::Angles new_angles = Character::Angles(angles);
// ...

// Set new values to the angles
motion.set_angles(j, new_angles);
}
}

Motion Filter

I implemented three kinds of motion filters - scalar filter, Gaussian filter and cartoon filter (which seems to have some problem and doesn’t work).

• Scalar

Scalar exaggeration can be considered as a filter where only the middle number is non-zero. As shown in the video, scalar filter usually exaggerate the motion when the amplitude is greater than 1.

vector<double> create_scalar_filter(int radius, double sigma) {
vector<double> filter;
if (x == 0) filter.push_back(1.25);
else filter.push_back(0);
}
return filter;
}
• Gaussian filter

Just like Gaussian filter in the image processing area, Gaussian motion filter provides a damping effect on the motion and sigma controls the degree of the effect. For example, when sigma is 5, visible damping effect can be seen; when it is 50, the whole model almost turns to the original pose.

• Cartoon filter

I haven’t figured out yet why this filter works very much like the Gaussian filter. I guess I misunderstand the use of LOG filter. The current filter is simply the second derivative of the Gaussian filter.

vector<double> create_cartoon_filter(int radius, double sigma) {
vector<double> filter;
double g = 1 / (M_PI * pow(sigma, 4)) * (1 - x * x / (2 * sigma * sigma)) * pow(M_E, -x * x / (2 * sigma * sigma));
filter.push_back(g);
}

return filter;
}

CCD IK and FABRIK

Overview

CCD IK and FABRIK are two approaches to implement IK. For this mini project, I created a Unity project where CCD and FABRIK are applied to a very simple skeleton. The main purpose is to go through the basic process of implementing IK.

Duration: 3 weeks

Role: Programmer

Team: Solo

Platform: Unity

Techniques

Unity

I considered Unity as the first choice because it provides very detailed and established interfaces for manipulating position and rotation of objects and implementing interactive interface.

For example, by directly controlling transform.position field of an gameObject, I can modify positions of joints and by mutiplying on transform.rotation , which is a Quanternion, rotations can be applied to bones.

public void SetPositionAndRotation (Vector3 tipOrigin, Vector3 tipTarget, Vector3 rootTarget) {
bone.transform.position = (tipTarget + rootTarget) / 2;

Quaternion rotation = Quaternion.FromToRotation (
tipOrigin - rootTarget,
tipTarget - rootTarget
);
bone.transform.rotation = bone.transform.rotation * rotation;
}

As for user input, although Unity has a class Input handling common inputs like keyboard and mouse, I simply utilized Unity Scene View.

Also, Unity provides a hierarchy of objects, so that I don’t need to calculate position for every part of a bone. Instead, I can simply calculate the position and rotation of the parent object and detail parts will change automatically.