Hair simulation

A while ago I implemented the hair simulation algorithm which is described in the paper Fast Simulation of Inextensible Hair and Fur by M.Müller, T.Y.Kim and N.Chentanez, using position based simulation. The simulation runs in real time can and handle quite some hairs (on CPU).

FTL.cpp

#include <glr/GL.h>
#include "FTL.h"
 
using namespace gl;
 
namespace ftl {
 
  Particle::Particle(Vec3 p, float m) 
    :position(p)
    ,tmp_position(p)
    ,enabled(true)
  {
 
    if(m < 0.001) {
      m = 0.001;
    }
 
    mass = m;
    inv_mass = 1.0 / mass;
  }
 
  // ------------------------
 
  FTL::FTL()
    :len(10)
  {
  }
 
  void FTL::setup(int num, float d) {
    float dim = 50;
    len = d;
    Vec3 pos(300, 300, 0);
    float mass = rx_random(1.0f, 10.0f);
    for(int i = 0; i < num; ++i) {
      Particle* p = new Particle(pos, mass); //rx_random(1.0f, 10.0f));
      particles.push_back(p);
      pos.y += d;
    }
 
    particles[0]->enabled = false;
  }
 
  void FTL::addForce(Vec3 f) {
    for(std::vector<Particle*>::iterator it = particles.begin(); it != particles.end(); ++it) {
      Particle* p = *it;
      if(p->enabled) {
        p->forces += f;
      }
    }
  }
 
  void FTL::update() {
    float dt = 1.0f/20.0f;
 
    // update velocities
    for(std::vector<Particle*>::iterator it = particles.begin(); it != particles.end(); ++it) {
      Particle* p = *it;
      if(!p->enabled) {
        p->tmp_position = p->position;
        continue;
      }
      p->velocity = p->velocity + dt * (p->forces * p->inv_mass);
      p->tmp_position += (p->velocity * dt);
      p->forces = 0;
      p->velocity *= 0.99;
    }    
 
    // solve constraints
    Vec3 dir;
    Vec3 curr_pos;
    for(size_t i = 1; i < particles.size(); ++i) {
      Particle* pa = particles[i - 1];
      Particle* pb = particles[i];
      curr_pos = pb->tmp_position;
      dir = pb->tmp_position - pa->tmp_position;
      dir.normalize();
      pb->tmp_position = pa->tmp_position + dir * len;
      pb->d = curr_pos - pb->tmp_position; //  - curr_pos;
    }    
 
    for(size_t i = 1; i < particles.size(); ++i) {
      Particle* pa = particles[i-1];
      Particle* pb = particles[i];
      if(!pa->enabled) {
        continue;
      }
      pa->velocity = ((pa->tmp_position - pa->position) / dt) + 0.9 *  (pb->d / dt);
      pa->position = pa->tmp_position;
    }
 
    Particle* last = particles.back();
    last->position = last->tmp_position;
 
  }
 
  void FTL::draw() {
 
    glLineWidth(0.1f);
    glr_begin(GL_LINE_STRIP);
    for(std::vector<Particle*>::iterator it = particles.begin(); it != particles.end(); ++it) {
      Particle* p = *it;
 
      if(!p->enabled) {
        glr_color(217, 41, 41);
      }
      else {
        glr_color(251, 251, 251);
      }
 
      glr_vertex(p->position);
    }
    glr_end();
  }
 
}

FTL.h

#ifndef ROXLU_FTL_H
#define ROXLU_FTL_H
 
#include <roxlu/Roxlu.h>
 
namespace ftl {
 
  struct Particle {
    Particle(Vec3 position, float m);
    Vec3 position;
    Vec3 tmp_position;
    Vec3 forces;
    Vec3 velocity;
    Vec3 d;
    float mass;
    float inv_mass;
    bool enabled;
 
  };
 
  class FTL {
  public:
    FTL();
    void setup(int num, float d);
    void addForce(Vec3 f);
    void update();
    void draw();
  public:
    float len;
    std::vector<Particle*> particles;
    Vec3 color;
  };
 
}
#endif