My idea is to have few concentric spherical layers of equidistant points like this:
where each layer has its own transformation matrix. At start all the layers have the same Transformation matrix and after each frame (or simulation timer tick) the matrices are passed (like cyclic ring buffer) so the lowest radius layer has actual matrix, the next one has matrix from previous frame and so on. Its basically geometric version of motion blur...
However my attempts to ray tracing this in fragment shader (as merged particles) based on these:
Hit a wall with accuracy and other math edge cases and rounding problems leading to ugly artifacts which to debug (if even possible) will take forever and porting to 64 bit will solve just a part of the problem...
So I decided to try to create geometry data for this on per frame basis which might be doable in geometry shader (later on). So first I tried to do this on CPU side (C++ and old GL api for now just for simplicity and as a proof of concept)
So I got list of equidistant points on unit sphere surface (static points) which I use to render the layers. So each point is transformed into line strip where the point is just scaled by layer radius and transformed by the matrix of layer it belong to. Here C++ code for this:
//---------------------------------------------------------------------------
// just some vector and matrix math needed for this
void vectorf_mul(float *c,float *a,float b) { for(int i=0;i<3;i++) c[i]=a[i]*b; } // c[3] = a[3]*b
void matrixf_mul_vector(float *c,float *a,float *b) // c[3] = a[16]*b[3]
{
float q[3];
q[0]=(a[ 0]*b[0])+(a[ 4]*b[1])+(a[ 8]*b[2])+(a[12]);
q[1]=(a[ 1]*b[0])+(a[ 5]*b[1])+(a[ 9]*b[2])+(a[13]);
q[2]=(a[ 2]*b[0])+(a[ 6]*b[1])+(a[10]*b[2])+(a[14]);
for(int i=0;i<3;i++) c[i]=q[i];
}
//---------------------------------------------------------------------------
// equidistant sphere points
const int sphere_na=20; // number of slices (latitude)
float sphere_pnt[sphere_na*sphere_na*6]; // equidistant sphere points
int sphere_n=0; // 3 * number of points
// create list of equidistant sphere points r=1 center=(0,0,0)
void sphere_init()
{
int ia,ib,na=sphere_na,nb;
float x,y,z,r;
float a,b,da,db;
da=M_PI/float(na-1); // latitude angle step
sphere_n=0;
for (a=-0.5*M_PI,ia=0;ia<na;ia++,a+=da) // slice sphere to circles in xy planes
{
r=cos(a); // radius of actual circle in xy plane
z=sin(a); // height of actual circle in xy plane
nb=ceil(2.0*M_PI*r/da);
if ((ia==0)||(ia==na-1)) { nb=1; db=0.0; } // handle edge cases
db=2.0*M_PI/float(nb); // longitude angle step
for (b=0.0,ib=0;ib<nb;ib++,b+=db) // cut circle to vertexes
{
x=r*cos(b); // compute x,y of vertex
y=r*sin(b);
sphere_pnt[sphere_n]=x; sphere_n++;
sphere_pnt[sphere_n]=y; sphere_n++;
sphere_pnt[sphere_n]=z; sphere_n++;
}
}
}
// render sphere as lines from center to surface
void sphere_draw()
{
int i;
glBegin(GL_LINES);
for (i=0;i<sphere_n;i+=3)
{
glColor3f(0.0,0.0,0.0); glVertex3f(0.0,0.0,0.0);
glColor3f(0.5,0.6,0.7); glVertex3fv(sphere_pnt+i);
}
glEnd();
}
//---------------------------------------------------------------------------
// puff ball
const int puff_n=8; // number of layers
float *puff_m[puff_n]={NULL}; // transform matrix for each layer
float puff_matrices[puff_n*16]; // transform matrix for each layer
float puff_col[puff_n][3]; // color for each layer
// render sphere as spicules dived to layers
void puff_draw(float r0,float r1)
{
int i,j;
float p[3],*p0,r,dr=(r1-r0)/float(puff_n);
glColor3f(0.5,0.6,0.7);
for (i=0;i<sphere_n;i+=3)
{
p0=sphere_pnt+i;
glBegin(GL_LINE_STRIP);
for (r=r0,j=0;j<puff_n;j++,r+=dr)
{
vectorf_mul(p,p0,r);
matrixf_mul_vector(p,puff_m[j],p);
glColor3fv(puff_col[j]);
glVertex3fv(p);
}
glEnd();
}
}
// update matrices
void puff_update()
{
int i;
float *p,t;
if (puff_m[0]==NULL) // in first pass all the matrices are actual
{
for (i=0;i<puff_n;i++)
{
puff_m[i]=puff_matrices+(i<<4);
glGetFloatv(GL_MODELVIEW_MATRIX,puff_m[i]);
t=1.0-(float(i)/float(puff_n-1));
puff_col[i][0]=0.1+0.3*t;
puff_col[i][1]=0.2+0.5*t;
puff_col[i][2]=0.2+0.6*t;
}
return;
}
p=puff_m[puff_n-1]; // insert actual matrix to ring FIFO buffer
for (i=puff_n-1;i>0;i--) puff_m[i]=puff_m[i-1];
puff_m[0]=p;
glGetFloatv(GL_MODELVIEW_MATRIX,puff_m[0]);
}
//---------------------------------------------------------------------------
And usage:
// init
sphere_init();
// render
glEnable(GL_DEPTH_TEST);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(0.0,0.0,-5.0); // this is just putting ball before my perspective camera
glRotatef(animt,0.0,0.3,0.5); // animt is just changing angle updated in timer
puff_update(); // this will update the matrices in layers (must be called before first draw)
glMatrixMode(GL_MODELVIEW); // however render itself is using unit matrix !!!
glLoadIdentity();
puff_draw(0.5,1.0); // render the hairs with radiuses from 0.5 up to 1.0
And here preview (for low enough hair count to actually see the geometry):
As you can see the hairs bend in correct direction (they are affected also by speed and this should be affected also by translation not just rotation as a bonus) so as a proof of concept this approach works. If you want them to bend more just lower the difference between r0,r1 radiuses.
As you can see this should be possible to port to GLSL geometry shader easily. The CPU just pass the unit sphere points once and geometry will emit a line strip for each point in the same way. You just need to have the matrices as uniforms.
Here another preview with changing rotation and lower r0,r1 difference:
Aside of tweaking the number of points sphere_na and layers puff_n this can be added to improve things:
- distribute the layers between
r0,r1 non linearly to achieve different look/behavior
- add sphere surface (
r0 so you do not see the hairs behind it)
- textures
- motion blur
- variable hair width (so at base its more wide)
- light model (maybe just for the surface and not for the hairs itself)
I think this might be used for improvement (curvature and variable width):
[Edit1] some tweaks and geometry shader
I succesfully ported this to shaders and added some stuff like colors, line width and inner sphere surface. Here the full VCL/C++/OpenGL/GLSL code:
VCL app window (just ignore the VCL stuff)
//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
float animx=0.0;
float animy=0.0;
bool _glsl=true;
//---------------------------------------------------------------------------
// just some vector and matrix math needed for this
void vectorf_mul(float *c,float *a,float b) { for(int i=0;i<3;i++) c[i]=a[i]*b; } // c[3] = a[3]*b
void matrixf_mul_vector(float *c,float *a,float *b) // c[3] = a[16]*b[3]
{
float q[3];
q[0]=(a[ 0]*b[0])+(a[ 4]*b[1])+(a[ 8]*b[2])+(a[12]);
q[1]=(a[ 1]*b[0])+(a[ 5]*b[1])+(a[ 9]*b[2])+(a[13]);
q[2]=(a[ 2]*b[0])+(a[ 6]*b[1])+(a[10]*b[2])+(a[14]);
for(int i=0;i<3;i++) c[i]=q[i];
}
//---------------------------------------------------------------------------
// equidistant sphere points
const int sphere_na=30; // number of slices (latitude)
float sphere_pnt[sphere_na*sphere_na*4]; // equidistant sphere points (approx size <= regular grid size / 1.5)
int sphere_n=0; // 3 * number of points
int sphere_slc[sphere_na+1]; // start index of each slice
// create list of equidistant sphere points r=1 center=(0,0,0)
void sphere_init()
{
int ia,ib,na=sphere_na,nb;
float x,y,z,r;
float a,b,da,db;
da=M_PI/float(na-1); // latitude angle step
sphere_n=0;
for (a=-0.5*M_PI,ia=0;ia<na;ia++,a+=da) // slice sphere to circles in xy planes
{
sphere_slc[ia]=sphere_n;
r=cos(a); // radius of actual circle in xy plane
z=sin(a); // height of actual circle in xy plane
nb=ceil(2.0*M_PI*r/da);
if ((ia==0)||(ia==na-1)) { nb=1; db=0.0; } // handle edge cases
db=2.0*M_PI/float(nb); // longitude angle step
for (b=0.0,ib=0;ib<nb;ib++,b+=db) // cut circle to vertexes
{
x=r*cos(b); // compute x,y of vertex
y=r*sin(b);
sphere_pnt[sphere_n]=x; sphere_n++;
sphere_pnt[sphere_n]=y; sphere_n++;
sphere_pnt[sphere_n]=z; sphere_n++;
}
}
sphere_slc[na]=sphere_n;
}
void sphere_gl_draw(float r)
{
int i,i0,i1,i2,n,n0,n1,j,k;
glMatrixMode(GL_MODELVIEW);
glPushMatrix();
glScalef(r,r,r);
glBegin(GL_TRIANGLE_STRIP);
for (k=1;k<sphere_na;k++)
{
i0=sphere_slc[k-1];
i1=sphere_slc[k+0];
i2=sphere_slc[k+1];
n0=(i1-i0)/3;
n1=(i2-i1)/3;
n=n0; if (n<n1) n=n1;
for (i=0;i<n;i++)
{
j=i0+(3*((i*n0)/n));
glNormal3fv(sphere_pnt+j);
glVertex3fv(sphere_pnt+j);
j=i1+(3*((i*n1)/n));
glNormal3fv(sphere_pnt+j);
glVertex3fv(sphere_pnt+j);
}
}
glEnd();
glPopMatrix();
}
//---------------------------------------------------------------------------
// puff ball
float r0=0.80,r1=1.0; // hair radiuses
const int puff_n=16; // number of layers (the more the more its bendy like)
float *puff_m[puff_n]={NULL}; // transform matrix for each layer
float puff_matrices[puff_n*16]; // transform matrix for each layer
float puff_col[puff_n*3]; // color for each layer
// render sphere as spicules dived to layers
void puff_gl_draw(float r0,float r1)
{
int i,j;
float p[3],*p0,r,dr=r1-r0;
float t,dt=1.0/float(puff_n);
glColor3f(0.5,0.6,0.7);
for (i=0;i<sphere_n;i+=3)
{
p0=sphere_pnt+i;
glBegin(GL_LINE_STRIP);
for (t=0.0,j=0;j<puff_n;j++,t+=dt)
{
r=r0+t*dr;
vectorf_mul(p,p0,r);
matrixf_mul_vector(p,puff_m[j],p);
glColor3fv(puff_col+j+j+j);
glVertex3fv(p);
}
glEnd();
}
}
void puff_glsl_draw(float r0,float r1)
{
int i;
glBegin(GL_POINTS);
for (i=0;i<sphere_n;i+=3)
{
glVertex3fv(sphere_pnt+i);
}
glEnd();
}
// update matrices
void puff_update()
{
int i;
float *p,t;
if (puff_m[0]==NULL) // in first pass all the matrices are actual
{
for (i=0;i<puff_n;i++)
{
puff_m[i]=puff_matrices+(i<<4);
glGetFloatv(GL_MODELVIEW_MATRIX,puff_m[i]);
t=float(i)/float(puff_n-1);
puff_col[i+i+i+0]=0.1+0.1*t;
puff_col[i+i+i+1]=0.1+0.2*t;
puff_col[i+i+i+2]=0.1+0.4*t;
}
return;
}
p=puff_m[puff_n-1]; // insert actual matrix to ring FIFO buffer
for (i=puff_n-1;i>0;i--) puff_m[i]=puff_m[i-1];
puff_m[0]=p;
glGetFloatv(GL_MODELVIEW_MATRIX,puff_m[0]);
}
//---------------------------------------------------------------------------
void gl_draw()
{
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glEnable(GL_DEPTH_TEST);
// animate matrix
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(0.0,0.0,-10.0*r1);
glRotatef(animx,1.0,0.0,0.0);
glRotatef(animy,0.0,1.0,0.0);
puff_update();
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
glEnable(GL_COLOR_MATERIAL);
glColor3f(0.4,0.2,0.1);
sphere_gl_draw(r0);
glDisable(GL_LIGHTING);
glDisable(GL_COLOR_MATERIAL);
glLineWidth(3.0);
if (_glsl)
{
GLint id;
glUseProgram(prog_id);
id=glGetUniformLocation(prog_id,"r0"); glUniform1f(id,r0);
id=glGetUniformLocation(prog_id,"r1"); glUniform1f(id,r1);
id=glGetUniformLocation(prog_id,"n" ); glUniform1i(id,puff_n);
// pass ring FIFO of matrices in their order
// can be improved by passing start index instead of reordering
int i,j,k;
float m[16*puff_n];
for (k=0,i=0;i<puff_n;i++)
for (j=0;j<16;j++,k++)
m[k]=puff_m[i][j];
id=glGetUniformLocation(prog_id,"mv"); glUniformMatrix4fv(id,puff_n,false,m);
glGetFloatv(GL_PROJECTION_MATRIX,m);
id=glGetUniformLocation(prog_id,"mp"); glUniformMatrix4fv(id,1,false,m);
id=glGetUniformLocation(prog_id,"c" ); glUniform3fv(id,puff_n,puff_col);
puff_glsl_draw(r0,r1);
glUseProgram(0);
}
else{
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
puff_gl_draw(r0,r1);
}
glLineWidth(1.0);
glFlush();
SwapBuffers(hdc);
}
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{
gl_init(Handle);
int hnd,siz; char vertex[4096],geom[4096],fragment[4096];
hnd=FileOpen("puff.glsl_vert",fmOpenRead); siz=FileSeek(hnd,0,2); FileSeek(hnd,0,0); FileRead(hnd,vertex ,siz); vertex [siz]=0; FileClose(hnd);
hnd=FileOpen("puff.glsl_geom",fmOpenRead); siz=FileSeek(hnd,0,2); FileSeek(hnd,0,0); FileRead(hnd,geom ,siz); geom [siz]=0; FileClose(hnd);
hnd=FileOpen("puff.glsl_frag",fmOpenRead); siz=FileSeek(hnd,0,2); FileSeek(hnd,0,0); FileRead(hnd,fragment,siz); fragment[siz]=0; FileClose(hnd);
glsl_init(vertex,geom,fragment);
// hnd=FileCreate("GLSL.txt"); FileWrite(hnd,glsl_log,glsl_logs); FileClose(hnd);
int i0,i;
mm_log->Lines->Clear();
for (i=i0=0;i<glsl_logs;i++)
if ((glsl_log[i]==13)||(glsl_log[i]==10))
{
glsl_log[i]=0;
mm_log->Lines->Add(glsl_log+i0);
glsl_log[i]=13;
for (;((glsl_log[i]==13)||(glsl_log[i]==10))&&(i<glsl_logs);i++);
i0=i;
}
if (i0<glsl_logs) mm_log->Lines->Add(glsl_log+i0);
sphere_init();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
{
gl_exit();
glsl_exit();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
{
gl_resize(ClientWidth,ClientHeight-mm_log->Height);
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
{
gl_draw();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Timer1Timer(TObject *Sender)
{
gl_draw();
animx+=1.5; if (animx>=360.0) animx=-360.0;
animy+=2.5; if (animy>=360.0) animy=-360.0;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormMouseDown(TObject *Sender, TMouseButton Button, TShiftState Shift, int X, int Y)
{
_glsl=!_glsl;
}
//---------------------------------------------------------------------------
The gl_simple.h can be found in here:
However the version posted there is old (the one I used has added geometry shader support...) Here the shaders:
// Vertex
#version 400 core
in vec3 pos;
void main()
{
gl_Position=vec4(pos,1.0);
}
// Geometry
#version 400 core
layout(points) in;
layout(line_strip, max_vertices = 32) out;
uniform int n; // hair layers
uniform float r0; // min radius for hairs
uniform float r1; // max radius for hairs
uniform mat4 mp; // global projection
uniform mat4 mv[32]; // hair modelview per layer
uniform vec3 c[32]; // color per layer
out vec3 col;
void main()
{
int i;
vec3 p;
float r,dr=(r1-r0)/float(n-1);
// input just single point
p=gl_in[0].gl_Position.xyz;
// emit line strip
for (r=r0,i=0;i<n;i++,r+=dr)
{
col=c[i];
gl_Position=mp*mv[i]*vec4(p.xyz*r,1.0);
EmitVertex();
}
EndPrimitive();
}
// Fragment
#version 400 core
in vec3 col;
out vec4 fcol;
void main()
{
fcol=vec4(col,1.0);
}
The VCL app is just simple window with one 40ms update timer and memo for GLSL logs so ignore tha and or mimic the VCL events in your environment. The _glsl just determines if shaders are used or not. This is far from optimized there is a lot of stuff to improve like using VAO instead glBegin/glEnd, recoding the matrix ring buffer so it does not need to reorder, etc ...
And finally new preview:
