CS 184 Project 4 Writeup

Tianchen Liu 3034681144, Riley Peterlinz 3036754368

Link:
https://cal-cs184-student.github.io/project-webpages-sp23-CardiacMangoes/proj4/index.html

Overview

In Part 1 we create a row-major order mesh using point masses and springs (three different types of springs). In Part 2 we first compute total force acting on each point mass from springs and other global factors before using Verlet integration to compute new point mass positions after each timestep. In Part 3 and 4 we handle collision of the point masses between other objects (spheres, planes) as well as other point masses (self-collision).In Part 5 we use GLSL to produce different kinds of lighting shaders.

Part 1

Implementation

Here, we have the basic structure of the points and springs being initialized in row-major order. Structural and Shearing constraints are along the axes with Shearing constraints being 2x as long as strucutural. Shearing constraints are diagonal.

The initial masses depends on the orientation of the plane and checks if a mass is pinned or not.

double x,y,z;

float dw = width / (num_width_points - 1);
float dh = height / (num_height_points - 1);

// masses
for (int j = 0; j < num_height_points; j++) {
    for (int i = 0; i < num_width_points; i++) {

       x = dw * i;

       // xz plane
       if (orientation == HORIZONTAL) {
           z = dh * j;
           y = 1;
       }

       // xy plane
       if (orientation == VERTICAL) {
           y = dh * j;
           z = rand() - (0.5 * RAND_MAX);
           z /= RAND_MAX;
           z /= 1000;
       }

       bool pinned_bool = false;
       for (auto v : pinned) {
         if (i == v[0] && j == v[1]) {
           pinned_bool = true;
         }
       }
       point_masses.emplace_back(Vector3D(x, y, z), pinned_bool);
    }
}

Initializing the springs has conditions for whether you are on the edge of the plane. For example, if you are at the top edge, you cannot place a STRUCTURAL spring on the mass above it as it does not exist!

for (int j = 0; j < num_height_points; j++) {
        for (int i = 0; i < num_width_points; i++) {

             target = &point_masses[num_width_points * i + j];

             // STRUCTURAL
             if (j != 0) {
                 left = &point_masses[num_width_points * i + (j - 1)];
                 springs.emplace_back(target, left, STRUCTURAL);
             }
             if (i != 0) {
                 above = &point_masses[num_width_points * (i - 1) + j];
                 springs.emplace_back(target, above, STRUCTURAL);
             }

             // SHEARING
             if (i != 0 && j != 0) {
                 above_left = &point_masses[num_width_points * (i - 1) + (j - 1)];
                 springs.emplace_back(target, above_left, SHEARING);
             }
             if (i < (num_width_points - 1) && j != 0) {
                 above_right = &point_masses[num_width_points * (i + 1) + (j - 1)];
                 springs.emplace_back(target, above_right, SHEARING);
             }

             // BENDING
             if (j > 1) {
                 left_two = &point_masses[num_width_points * i + (j - 2)];
                 springs.emplace_back(target, left_two, BENDING);
             }
             if (i > 1) {
                 above_two = &point_masses[num_width_points * (i - 2) + j];
                 springs.emplace_back(target, above_two, BENDING);
             }

         }
     }

Results

No Shearing Constraints

Just Shearing Contsraints

All Constraints

Part 2

Implemenation

We first iterate through the point masses to find the force due to external accelerations

for (int i = 0; i < point_masses.size(); i++){
        PointMass *pm = &point_masses[i];
        for (int j = 0; j < external_accelerations.size(); j++) {
            pm->forces += mass * external_accelerations[j];
        }
    }

Then, for each spring we calculate the net force on each of its point masses by considering what type of constraints are on the point mass. These are internal forces.

for (int i = 0; i < springs.size(); i++) {
    Spring* spring = &springs[i];

    float k_s = 0;
    double d = (spring->pm_a->position - spring->pm_b->position).norm() - spring->rest_length;

    if (cp->enable_structural_constraints) k_s += cp->ks;
    if (cp->enable_shearing_constraints) k_s += cp->ks;
    if (cp->enable_bending_constraints) k_s += cp->ks * 0.2;

    spring->pm_a->forces += k_s * d * (spring->pm_b->position - spring->pm_a->position).unit();
    spring->pm_b->forces += k_s * d * (spring->pm_a->position - spring->pm_b->position).unit();
}

Next, we use Verlet Integration to calculate the new position given the forces on the point mass

for (int i = 0; i < point_masses.size(); i++){
      PointMass *pm = &point_masses[i];
      if (pm->pinned == false){
          Vector3D prev = pm->last_position;
          Vector3D curr = pm->position;
          pm->position = curr + ((1-(cp->damping/100.0)) * (curr - prev)) + ((pm->forces / mass) * pow(delta_t, 2));
          pm->last_position = curr;
      }
      pm->forces = Vector3D(0);
  }

Finally, we constrain the changes to be such that the spring does not change in length more than 10% per timestep [Provot 1995]

for (int i = 0; i < springs.size(); i++) {
  Spring *spring = &springs[i];
  PointMass *pm_a = spring->pm_a;
  PointMass *pm_b = spring->pm_b;

  double spring_length = (pm_a->position - pm_b->position).norm();
  double extra_length = spring_length - (spring->rest_length * 1.1);

  if (extra_length > 0) {
      if (pm_a->pinned) {
          pm_b->position += extra_length * (pm_a->position - pm_b->position).unit();
      } else if (pm_b->pinned) {
          pm_a->position += extra_length * (pm_b->position - pm_a->position).unit();
      } else if (!pm_a->pinned && !pm_b->pinned) {
          pm_b->position += extra_length / 2 * (pm_a->position - pm_b->position).unit();
          pm_a->position += extra_length / 2 * (pm_b->position - pm_a->position).unit();
      }
  }
}

Results

After simulation, the state with default parameters

ks = 5000 N/m
density = 15 g/cm^2
dampening = 0.2%

ks

ks is the measure of how tight or loose a spring is. Small ks makes a spring loose while large ks makes the spring tight.

ks = 10 N/m
density = 15 g/cm^2
dampening = 0.2%

ks = 50000 N/m
density = 15 g/cm^2
dampening = 0.2%

density

density affects the weight of our masses, so it will make the cloth tighter when density is low and will cause more sagging when density is high.

ks = 5000 N/m
density = 1 g/cm^2
dampening = 0.2%

ks = 5000 N/m
density = 50 g/cm^2
dampening = 0.2%

dampening

Unlike ks and density, dampening affects only the motion of the cloth, not its final resting state. With low dampening, more energy is conserved and thus the cloth will seem “bouncier”, like a rubber sheet in the air. With high dampening the sheet of cloth falls slower and has less “fidgety” behavior.

ks = 5000 N/m
density = 50 g/cm^2
dampening = 0%

ks = 5000 N/m
density = 50 g/cm^2
dampening = 1%

Part 3

Implementation

To handle collisions with object privitives, we simply loop over each point mass and each primitive to check their collisions

for (auto &obj : *collision_objects) {
    for (auto &pm : point_masses) {
      obj->collide(pm);
    }
  }

Sphere

Implementation

Collisions are radially outward thus the tangent vector is in the direction away from the origin at all times

float distance = (pm.position - origin).norm();
  if (distance <= radius) {
    auto dir = (pm.position - origin).unit();
    auto tangent_point = origin + dir * radius;
    auto correction_vector = (tangent_point - pm.last_position);
    pm.position = pm.last_position + correction_vector * (1 - friction);
  }

Result

We vary ks to notice as ks grows, the cloth falls on the sphere more rigidly

ks = 500 N/m
density = 15 g/cm^2
dampening = 0.2%

ks = 5000 N/m
density = 15 g/cm^2
dampening = 0.2%

ks = 50000 N/m
density = 15 g/cm^2
dampening = 0.2%

Plane

Implementation

For collision, we simply look at the normal vs the direction of our point mass compared with the plane and make sure to correct the point to above the plane slightly on the side it is entering from.

Vector3D dir = pm.position - point;

float costheta = dot(normal, dir);

if (costheta < 0) {
  float offset = SURFACE_OFFSET - costheta;
  Vector3D correction = (pm.position + normal * offset) - pm.last_position;
  pm.position = pm.last_position + (1-friction) * correction;
}

Result

A plane of cloth lying peacefully at rest on a plane

Part 4

Implementation

First, we have to build a spatial map for efficient positional encoding for self-collisions. We first map a small box of the cloth to a hash like so:

float Cloth::hash_position(Vector3D pos) {
  float w = 3 * width / num_width_points;
  float h = 3 * height / num_height_points;
  float t = max(w,h);

  float x = floor(pos.x / w);
  float y = floor(pos.y / h);
  float z = floor(pos.z / t);

  // raise values to powers in order to keep unique
  return x + pow(y, 2) + pow(z, 3);
}

We take the floor to get discrete values, then power and sum them to spreead the hash out so we are less likely to get multiple hits in a single box.

Nowe, we actually build the spatial map using our hashing. This basically hashes the position which is a key in a hashmap to find close points.

void Cloth::build_spatial_map() {
  for (const auto &entry : map) {
    delete(entry.second);
  }
  map.clear();
    
  for (int i = 0; i < point_masses.size(); i++){
      PointMass* pm = &point_masses[i];
      // Key is not present
      float hash = hash_position(pm->position);
      if (map.find(hash) == map.end())
          map[hash] = new vector<PointMass*>();
      map[hash]->push_back(pm);
  }

}

The actual self-collision will take close point masses (that are not itself) and then if the points are less than 2 * thickness away from each other, we can then seperate them so they are at least 2 * thickness apart. We average over all the close points corrections, dividing by the simulation steps to keep the simulation temporaly consistent.

void Cloth::self_collide(PointMass &pm, double simulation_steps) {
  // TODO (Part 4): Handle self-collision for a given point mass.
  float hash = hash_position(pm.position);

  // correction initialization
  Vector3D correction(0);
  int correction_count = 0; // for averaging

  auto close_points = *(map[hash]);

  // separation between pm and close_pm
  for (PointMass* close_pm : close_points){
      if (&*close_pm == &pm) {
      } else{
          Vector3D separation = (close_pm->position - pm.position);
          if (separation.norm() < 2 * thickness) {
              correction += separation - 2 * thickness * separation.unit();
              correction_count += 1;
          }
      }

  }

  if (correction_count > 0){
      correction /= correction_count;
      correction /= simulation_steps;
      pm.position += correction;
  }

}

In simulate we just call self_collide on all the point masses

build_spatial_map();
for (int i = 0; i < point_masses.size() ; i++) {
    PointMass* pm = &point_masses[i];
    self_collide(*pm, simulation_steps);
}

Results

Falling Sequence

density vs ks

Lower density causes the cloth to be flatter when it falls since it is light and can spread away from itself given self-collision

Higher density causes it to clump onto itself since it is heavier now and weighs itself down so self-collision doesn’t have as large of an effect

The effect of ks is inverse to density

Lower ks causes the cloth to clump onto itself more since the structure becomes loose

Higher ks causes the cloth to flatten more since it has more rigid structure

Part 5

Explain in your own words what is a shader program and how vertex and fragment shaders work together to create lighting and material effects.

Shaders are isolated programs that run in parallel on GPU, executing sections of the graphics pipeline, which makes it faster than the CPU. Good for processing parallelizable stuff like information about each vertex and each triangle formed by three vertices.

Vertex shaders change positions, tangents/normals of each vertex, and where the vertex should be in the screen space
Fragment shaders takes in info from the vertex shaders and output the color for each area segment formed by the vertices. Together using the vertex shaders in the .vert files and the fragment shaders in the .frag files we get lighting and shading effects.

Explain the Blinn-Phong shading model in your own words.

Blinn-Phong shading allows us to create highlights and shading with ambient light and diffuse reflections and specular highlights. It achieves this in three parts in the equation.

First term is ambient lighting, second is diffuse reflection, and third is spectural highlights.

Show a screenshot of your Blinn-Phong shader outputting only the ambient component, a screen shot only outputting the diffuse component, a screen shot only outputting the specular component, and one using the entire Blinn-Phong model.

Only ambient component, ka = 0.5, Ia = [1, 1, 1]

Only diffuse component, kd = 0.5

Only specular component, ks = 0.9, p = 10

Entire Blinn-Phong model, parameters listed above:

Show a screenshot of your texture mapping shader using your own custom texture by modifying the textures in /textures/.

This is Thom Yorke

Show a screenshot of bump mapping on the cloth and on the sphere.

The cloth:

The sphere:

Show a screenshot of displacement mapping on the sphere.

Use the same texture for both renders. You can either provide your own texture or use one of the ones in the textures directory, BUT choose one that’s not the default texture_2.png. Compare the two approaches and resulting renders in your own words.

Bump mapping changes the “look” of the render to make it “look” like the texture but it’s actually the same shape (the vertices doesn’t change places). Whereas in displacement mapping, the vertex position is actually displaced to a different position. As you can see, in the rendering of displacement mapping, the sphere is not a strictly perfect sphere anymore.

Compare how your the two shaders react to the sphere by changing the sphere mesh’s coarseness by using -o 16 -a 16 and then -o 128 -a 128.

The less the coarse the mapping is the less detailed the render is. This is especially true in displacement shading more than bump shading. In the following four pictures we can see that there’s really not that much of a difference in bump shading with the two different parameters. But in displacment shading the with param=16 it doesn’t look circular at all, and when param=128 we can see the ups and downs of the bricks.

param = 16:

param = 128:

Show a screenshot of your mirror shader on the cloth and on the sphere.