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.

Screen Shot 2020-03-26 at 7.17.49 PM

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.

    vector<double> create_scalar_filter(int radius, double sigma) {
        vector<double> filter;
        for (int x = -radius; x <= radius; x++) {
            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;
        for (int x = -radius; x <= radius; x++) {
            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;
    }

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.