## 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();
}
```

### Follow up

- 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/

### Follow up

- 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.

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.

### Follow up

- Bad cartoon filter effects

## 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.

#### CCD IK

#### FABRIK

### Follow up

Jacobian IK

CCD IK and FABRIK are in some way similar. The main idea of both two approaches is to interatively modify the joints until reaching the error threshold. On the other hand, Jacobian methods are relatively based on math and calculation. With the help of an established math library, Jacobian IK is also not difficult to implement.

Maya

Implementing IK in Maya by making some tools might be very meaningful.