Archived

This topic is now archived and is closed to further replies.

mechamecha

help using jakobsen method to simulate cloth

Recommended Posts

mechamecha    122
hello ive read through the following paper at http://www.ioi.dk/Homepages/thomasj/publications/gdc2001.htm and im still a bit confused with jakobsen''s approach to simulate cloth. anyone have or know of any source code demonstrating jakobsen''s method? from my understanding, cloth is simulated by creating a mesh of triangles with each edge represented by a "spring". As the springs contract and expand, the overall effect lets the cloth bend and wrap around other objects. but doesnt the jakobsen method relax all the edges to its restlength every frame? so as a result, your cloth of springs would be rendered as a "stiff" like board....or am i missing something? i appolgize for the crude manner in which this is written...too tried to think straight right now.

Share this post


Link to post
Share on other sites
oliii    2196
yes, it gives a pretty stiff mesh. however, you can relax it with various means. by reducing the numbers of iterations, and softening the constraint (by no forcing the particles to be at restlength apart, but rather, push the particles slowly towards the desired distance).


void SatisfyConstraint()
{
if (!m_pxParticle[0] || !m_pxParticle[1])
return;
//-----------------------------------------------------------

// cache stuff

//-----------------------------------------------------------

float restlength = m_fRestLength;
float restlength2 = restlength*restlength;
float invmass1 = m_pxParticle[0]->GetInvMass();
float invmass2 = m_pxParticle[1]->GetInvMass();
float invmass = invmass1+invmass2;

if (invmass < 0.00001f)
return;

//-----------------------------------------------------------

// relative position of particles

//-----------------------------------------------------------

Vector& x1 = m_pxParticle[0]->GetTargetPos();
Vector& x2 = m_pxParticle[1]->GetTargetPos();
Vector delta = x2 - x1;

float diff;
float delta2 = delta*delta;

//-----------------------------------------------------------

// Using Square root approximation

// calcualte the difference from the current distance to the ideal distance

//-----------------------------------------------------------

diff = restlength2/(delta2 + restlength2)-0.5f;
diff *= -2.0f;

//-----------------------------------------------------------

// move particles so their distance = ideal distance

//-----------------------------------------------------------

float sloppyness = 0.5f;
delta *= diff / invmass;
x1 += delta*(invmass1) * sloppyness;
x2 -= delta*(invmass2) * sloppyness;
}


but nothing prevents you from using a spring equation instead of a stiff constraint equation if you want.

Share this post


Link to post
Share on other sites
anti_alia    127
Also I find that if you make a mistake in the constraint section (ie iterate over the same constraint n times before moving to the next!! - yeah exactly like oliii suggested but with waisting heaps of cpu!) you can get some pretty nice results. I found that using nice heavy particles helped too..

Edit: Oh! and also found that putting cross bars made it move nicer too.. i had just the triangles at first.. kept squeing to one side.. also if you get rope working (jsut a string of particles and constraints) then its very easy to go from there to cloth. (just add a dimention)


A.

[edited by - anti_alia on March 30, 2004 5:59:10 PM]

Share this post


Link to post
Share on other sites
oliii    2196
yeah, it''s simple stable, nice looking, relatively efficient, you can bind the cloth to walls very easily, break links, tear the cloth... the problem, like in most of these situations, is to combine the rendering to the phsyics, which isn''t hard, but the vertex normal calculations can be interesting, and need careful thoughts about making it as fast as possible, because it''s gonna cost!

Share this post


Link to post
Share on other sites
mechamecha    122
thanks for the response guys...
oliii:
my code looks a lot like your code, except you have a "slopyness" value..is this what you mean by softening the restraint?

anti_ala:
could you further explain what you mean by putting "cross bars"?

thanks

Share this post


Link to post
Share on other sites
anti_alia    127
sorry my bad... it depends on how you generate your mesh.. but basically if you have 4 particles in a square - then the ''cross bars'' would be from p1->p4 and p2->p3 .. the other mistake i made first time was not giving them (p1->p4 and p2->p3) the proper length ie the long side of a right angled triangle.

(prehaps you dont even need these?)

silly mistakes i know.. oh and I also used a slightly springy constraint. (so that the cloth could strech etc)

I havent had much time to muck around with the cloth cause i got distracted by other things.. (damn work!).. and I dont actually need it for what im working on.. its kinda perdy tho.

A.

Share this post


Link to post
Share on other sites
oliii    2196
''sloppiness''. yeah, that''s one wayt of softening the constraint. I used it for soft bodies. although you have to be careful, things can blow up more easily.

Share this post


Link to post
Share on other sites
mechamecha    122
thanks for the input guys.
so ive been messing around a bit w/ vertlet integration scheme and constraint sys to see how it simulates cloths.

one simulation im working on is rendering a scene w/ a sphere and a cloth composed of a mesh of triangles....each vertex is a particle and the edges connecting them are constrained to a length, just like jokobsen suggests in his paper.
im trying to get the cloth to fall on the sphere and slightly conform to its shape just like a real cloth would...but im not getting the results i want. basicaly, as soon as the cloth begins to wrap itself around the sphere, it reflects off and collapses/folds itself togethor and bounces all over my scene....

im using the source code oliii has provided to satisfy my constraints. when an intersection btw the particle on the cloth and the surface of the sphere occur, i put the particle on the surface of the sphere. i would post my exact code, but i dont know how to put it in its own scroll box.....

im wondering if theres something wrong in my implementation or understanding of vertlets? Or perhaps using vertlets doesnt work well in this kind of simulation?

im think that its the latter..that this approach doesnt work b/c the constraints imposed on the mesh are not allowing the mesh to conform to the shape of the sphere. maybe i need to allow edges that intersect the sphere to relax to a different length?

any help would be appreciated..
thanks.


Share this post


Link to post
Share on other sites
oliii    2196
it should be fine, really.

here''s some code I have.
It''s very, frictionless



#include <errno.h>
#include <iostream.h>
#include <math.h>
#include <stdio.h>
#include <GL/glut.h>


static int window_w = 600, window_h = 400;
static int mouse_x=0, mouse_y=0;
static int mouse_dx=0, mouse_dy=0;



inline void ValidateFloat(float fNum)
{
UINT32 Num = *(UINT32 *)&fNum;

Num = (Num >> 23) & 255;

if (Num == 255)
cerr <<"Invalid floating point number"<<endl;
}

inline float FABS(float x) { return (float) fabs(x); }

inline float SquareRoot(float x) { return (float) sqrt(x); }

inline float Random(float x) { return (rand() / (float) RAND_MAX) * (x); }

inline float Random(float min, float max)
{
return Random(max - min) + min;
}

struct Vector
{
float x,y,z;

inline Vector(void){}; //no initialisation

inline Vector(float Ix,float Iy,float Iz):x(Ix),y(Iy),z(Iz) {} //instantiate members


inline float operator * (const Vector &pOtherVec)const
{
return ( (x*pOtherVec.x) + (y*pOtherVec.y) + (z*pOtherVec.z) );
}

inline float GetLengthSquared(void)const { return (*this) * (*this);}
inline float GetLength(void)const
{
float Ret = SquareRoot( GetLengthSquared() );
ValidateFloat(Ret);
return Ret;
}

float Normalise(void); //normalise and return length


inline Vector &operator *=(const float Scalar)
{
x *= Scalar; y *= Scalar; z *= Scalar;

return *this;
}

inline Vector &operator /=(const float Scalar)
{
*this *= (1.0f / Scalar);

return *this;
}

inline float Dot(const Vector& V) const
{
return (*this) * V;
}

inline Vector Cross(const Vector& V) const
{
Vector T;
return T.Cross(*this, V);
}

inline Vector &Cross(const Vector &In1,const Vector &In2)
{
float Tempx,Tempy;

Tempx = (In1.y * In2.z) - (In1.z * In2.y);
Tempy = (In1.z * In2.x) - (In1.x * In2.z);
z = (In1.x * In2.y) - (In1.y * In2.x);

x = Tempx;
y = Tempy;

return *this;
}

inline Vector operator & (const Vector &Vn2)
{
Vector Temp;

return Temp.Cross(*this, Vn2);
}

inline void SetMinimum(const Vector &Comparison)
{
if (x > Comparison.x) x = Comparison.x;
if (y > Comparison.y) y = Comparison.y;
if (z > Comparison.z) z = Comparison.z;
}

inline void SetMaximum(const Vector &Comparison)
{
if (x < Comparison.x) x = Comparison.x;
if (y < Comparison.y) y = Comparison.y;
if (z < Comparison.z) z = Comparison.z;
}

inline Vector &operator +=(const Vector &Other)
{
x += Other.x; y += Other.y; z += Other.z;

return *this;
}

inline Vector &operator -=(const Vector &Other)
{
x -= Other.x; y -= Other.y; z -= Other.z;

return *this;
}

inline Vector operator - (const Vector &Other) const
{
Vector Temp = *this;

return Temp -= Other;
}

inline Vector operator + (const Vector &Other) const
{
Vector Temp = *this;

return Temp += Other;
}

inline Vector& Clamp(const Vector& Min, const Vector& Max)
{
if (x < Min.x) x = Min.x; else if (x > Max.x) x = Max.x;
if (y < Min.y) y = Min.y; else if (y > Max.y) y = Max.y;
if (z < Min.z) z = Min.z; else if (z > Max.z) z = Max.z;

return *this;
}

inline Vector operator * (float k) const
{
Vector Temp = *this;

return Temp *= k;
}

inline Vector operator ^ (const Vector& V) const
{
return Cross(V);
}

inline Vector operator / (float k) const
{
Vector Temp(*this);
return Temp /= k;
}

inline Vector& operator ^=(const Vector& V)
{
*this = this->Cross(V);

return *this;
}

inline Vector &AddScaledVector(float Scalar,const Vector &Vec)
{
x += Scalar * Vec.x;
y += Scalar * Vec.y;
z += Scalar * Vec.z;

return *this;
}

inline Vector &SetToScaledVector(float Scalar,const Vector &Vec)
{
x = Scalar * Vec.x;
y = Scalar * Vec.y;
z = Scalar * Vec.z;

return *this;
}

inline BOOL operator != (const Vector &Other)
{
const float Epsilon = 0.0001f;

if ((FABS(Other.x - x) > Epsilon) ||
(FABS(Other.y - y) > Epsilon) ||
(FABS(Other.z - z) > Epsilon))
{
return TRUE;
}

return FALSE;
}

inline BOOL operator == (const Vector &Other) const
{
const float Epsilon = 0.0001f;

if ((FABS(Other.x - x) > Epsilon) ||
(FABS(Other.y - y) > Epsilon) ||
(FABS(Other.z - z) > Epsilon))
{
return FALSE;
}

return TRUE;
}

static Vector Random(const Vector& Min, const Vector& Max)
{
return Vector(::Random(Min.x, Max.x), ::Random(Min.y, Max.y), ::Random(Min.z, Max.z));
}

static float DistanceSquared(const Vector& A, const Vector& B) // returns the distance between the two given vectors, (where we assume the vetors are being used as positions)

{
Vector Delta = A;
Delta -= B;
return Delta.GetLengthSquared();
}

inline void Clear(void) {x = y = z = 0.0f;}

#ifdef _DEBUG
void Print(void) const;

inline void Validate(void) const {ValidateFloat(x); ValidateFloat(y); ValidateFloat(z);}
#endif _DEBUG
};



template <const int width, const int height>
class CSheet
{
public:
void AddConstrainVolume(const Vector& Cen, float Rad)
{
m_Cen[m_Num] = Cen;
m_Rad[m_Num] = Rad;
m_Num++;
}

CSheet(float XSize=400.0f, float YSize=400.0f)
{
if ((width & (width-1)) || (height & (height-1)))
{
printf("Cloth size is not power of 2\n.");
}

float lx = XSize / (float) width;
float ly = XSize / (float) height;
m_lx2 = lx*lx;
m_ly2 = ly*ly;

for(int i = 0; i < width; i++)
{
for(int j = 0; j < height; j++)
{
m_Pos [i][j] = Vector(i * lx, 0, j * ly);
m_Pos [i][j]+= Vector(0, YSize, 0);
m_Accel [i][j] = Vector(0, 0, 0);
m_OldPos[i][j] = m_Pos[i][j];
}
}
m_Num=0;
}

void Update(void)
{
Verlet();
Constrain();
ResetForces();
}

void Render(void)
{
//-------------------------------------------------

// render points

//-------------------------------------------------

glPointSize (3.0f);
glColor3f (1.0f, 0.3f, 0.3f);
glBegin (GL_POINTS);

Vector* px = &m_Pos[0][0];

for(int i = width*height; i > 0; i --)
{
glVertex3fv((GLfloat*)px);
px++;
}
glEnd();

int w = (width-1);
int h = w * height;

//-------------------------------------------------

// render lines

//-------------------------------------------------

glLineWidth (1.0f);
glColor3f (1.0f, 1.0f, 1.0f);
glBegin (GL_LINES);

for (i = 0; i < width*height; i ++)
{
if (i & w)
{
glVertex3fv((GLfloat*) px);
glVertex3fv((GLfloat*) (px-1));
}

if (i >= width)
{
glVertex3fv((GLfloat*) px);
glVertex3fv((GLfloat*) (px-width));
}

px++;
}

glEnd();

//-------------------------------------------------

// render volumes

//-------------------------------------------------

glColor4f(0.3f, 1.0f, 0.3f, 0.3f);

for(i = 0; i < m_Num; i ++)
{
glPushMatrix();
glTranslatef(m_Cen[i].x, m_Cen[i].y, m_Cen[i].z);
glutSolidSphere(m_Rad[i], 16, 16);
glPopMatrix();
}

}

private:
Vector m_Pos [width][height];
Vector m_OldPos[width][height];
Vector m_Accel [width][height];
float m_lx2;
float m_ly2;

Vector m_Cen[32];
float m_Rad[32];
int m_Num;


void Verlet(void)
{
Vector* px = &m_Pos [0][0];
Vector* pox = &m_OldPos[0][0];
Vector* pa = &m_Accel [0][0];

const float dt = 1.0f / 30.0f;
const float dt2 = dt*dt;

for(int i = width*height; i > 0; i --)
{
Vector Temp = *px;

px->x += (px->x - pox->x) + pa->x * dt2;
px->y += (px->y - pox->y) + pa->y * dt2;
px->z += (px->z - pox->z) + pa->z * dt2;

*pox = Temp;
px++;
pox++;
pa++;
}
}

void ResetForces(void)
{
Vector* pa = &m_Accel [0][0];

for(int i = width*height; i > 0; i --)
{
pa->x = 0.0f;
pa->y = -10.0f;
pa->z = 0.0f;

pa++;
}
}

void Constrain(void)
{
WorldConstrain();

for(int i = 0; i < 10; i ++)
{
StiffConstrain();
}
}

void WorldConstrain(void)
{
Vector* px = &m_Pos[0][0];
for(int i = width*height; i > 0; i --)
{
VolumeConstrain(px);
px++;
}
}

void VolumeConstrain(Vector* px)
{
Vector *pCen = m_Cen;
float *pRad = m_Rad;

for(int i = m_Num; i > 0; i --)
{
VolumeConstrain(*px, *pCen, *pRad);
pCen++;
pRad++;
}
}

void VolumeConstrain(Vector& x, const Vector& Cen, float Rad)
{
Vector D = (x - Cen);
float d2 = D * D;
float r2 = Rad * Rad;
if (d2 > r2)
return;

D /= (float) sqrt(d2);
x = Cen + D * Rad;
}


void StiffConstrain(void)
{
Vector *px = &m_Pos[0][0];

int w = (width-1);

for (int i = 0; i < width*height; i ++)
{
if (i & w)
Constrain(px, px-1, m_lx2);

if (i >= width)
Constrain(px, px-width, m_ly2);

px++;
}
}

void Constrain(Vector* px1, Vector* px2, float restlength2)
{
Vector delta(px2->x - px1->x, px2->y - px1->y, px2->z - px1->z);

float delta2 = delta*delta;

float diff = restlength2 / (delta2 + restlength2) - 0.5f;
delta *= diff;

*px1 -= delta * 0.5f;
*px2 += delta * 0.5f;
}
};

typedef CSheet<16, 16> CCloth;

CCloth Cloth;

void ShowStroke(char *message)
{
char *string = message;
glScaled(0.00833333, 0.00833333, 1.0); //// scaling to 1 / 120.0 (the size of the font).

glPushMatrix();
for(; *string; string++)
{
switch(*string)
{
case 10:
glPopMatrix();
glTranslated(0.0, -140.0, 0.0);
glPushMatrix();
break;
case 9:
glTranslated(420.0, 0.0, 0.0);
break;
case 0 :
return;
break;
default:
glutStrokeCharacter(GLUT_STROKE_ROMAN, *string);
break;
}
}
glPopMatrix();
}

static void init (void);
static void reshape (int w, int h);
static void Motion (int x, int y);
static void PassiveMotion (int x, int y);
static void Motion (int x, int y);
static void render (void);
static void timer (int t);
static void keyboard (unsigned char key, int x, int y);

void main(int argc, char** argv)
{
glutInit(&argc, argv);

init();

glutMainLoop();
}

void init(void)
{
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
glutInitWindowSize(window_w, window_h);
glutInitWindowPosition(100, 100);
glutCreateWindow("OpenGLViewer");
glutReshapeFunc (reshape);
glutDisplayFunc (render);
glutTimerFunc(30, timer, 0);
glutPassiveMotionFunc (PassiveMotion);
glutMotionFunc (Motion);
glutKeyboardFunc (keyboard);
glClearColor (0.3f, 0.3f, 0.3f, 1.0f);

glPointSize(3.0f);

glShadeModel(GL_SMOOTH);
glFrontFace(GL_CCW);
glPolygonMode(GL_FRONT, GL_FILL);
glPolygonMode(GL_BACK , GL_FILL);
glDisable(GL_CULL_FACE);

glEnable (GL_BLEND);
glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

}

void timer(int t)
{
render();

glutTimerFunc(10, timer, 0);
}

void keyboard(unsigned char key, int x, int y)
{
Cloth = CCloth(window_w, window_h);

Cloth.AddConstrainVolume(Vector(window_w/2, -window_h/2, window_w/2), window_w / 2);
Cloth.AddConstrainVolume(Vector(window_w/2, 0, window_w/2), window_w / 4);
Vector Min(-0, -window_h, 0);
Vector Max(window_w, 0, window_w);
float rad(window_w / 2);
for(int i = 0; i < 4; i ++)
{
Cloth.AddConstrainVolume(Vector::Random(Min, Max), Random(rad));
}
}

void reshape(int w, int h)
{
window_w = w;
window_h = h;

Cloth = CCloth(window_w, window_h);

Cloth.AddConstrainVolume(Vector(window_w/2, -window_h/2, window_w/2), window_w / 2);
Cloth.AddConstrainVolume(Vector(window_w/2, 0, window_w/2), window_w / 4);
Vector Min(-0, -window_h, 0);
Vector Max(window_w, 0, window_w);
float rad(window_w / 2);
for(int i = 0; i < 4; i ++)
{
Cloth.AddConstrainVolume(Vector::Random(Min, Max), Random(rad));
}
}

void glutBitmapString(void* font, char* str)
{
char* chr = str;

while (*(chr))
glutBitmapCharacter(GLUT_BITMAP_9_BY_15, *(chr++));
}


void render(void)
{
glClear(GL_COLOR_BUFFER_BIT);

glViewport(0, 0, (GLsizei) window_w, (GLsizei) window_h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();

gluPerspective(90.0f, window_w / (float) window_h, 0.0f, 1000.0f);

static float a=0;
static float d =1.0f;

a += mouse_dx / (float) (window_w) * 4.0f;
d += mouse_dy / (float) (window_h) * 0.5f;

float xpos = cos(a) * (window_w * d) + window_w/2;
float zpos = sin(a) * (window_w * d) + window_w/2;


glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(xpos , window_h/2 , zpos,
window_w/2, window_h/2, window_w/2,
0, 1, 0);
Cloth.Update();
Cloth.Render();

glutSwapBuffers();
}

void Motion(int x, int y)
{
mouse_dx=x - mouse_x;
mouse_dy=y - mouse_y;

mouse_x = x;
mouse_y = y;
}


void PassiveMotion(int x, int y)
{
mouse_x = x;
mouse_y = y;
mouse_dx=0;
mouse_dy=0;
}

Share this post


Link to post
Share on other sites
mechamecha    122
oliii,
thanks so much man!
after comparing your code to mine, i figured out what i was doing wrong. apparently, my collision detection between the particle of cloth and sphere was waaay off. It seems your code constrains the motion of each particle to the surface of the sphere upon collision, where i was trying to calculate the point of collision and constrain the particle to that point. as result, after the next vertelet calcuation, that particle would reflect in the opposite direction, causing the whole mesh to bounce off the sphere.

i think i am begining to better understand this technique and how constraints should be defined.

i appreciate all your help!
thanks again

Share this post


Link to post
Share on other sites