Spherical Harmonics Lighting的代码实现(基于OpenGL)

1、二维空间的勒让德多项式

勒让德多项式定义在[-1,1]范围内，其递归式是





下面这个函数的参数是给定的x，给定的l和m，其中l必须是正整数，而且m在[-l,l]范围内。

//勒让德多项式计算方法
double ALPStd(float x,int l,int m)
{
if(l==m)
//doubleFactorial(x)计算x!!
return (pow(-1,m)*doubleFactorial(2*m-1)*pow(sqrt(l-x*x),m));
if(l==m+1)
return (x*(2*m+1)*ALPStd(x,m,m));
return ((x*(2*l-1)*ALPStd(x,l-1,m)-(l+m-1)*ALPStd(x,l-2,m))/(l-m));

}

int doubleFactorial(int x)
{
int result;
if(x==0 || x==-1)
return 1;
result=x;
while((x-=2)>0)
result*=x;
return result;
}

下面是利用勒让德多项式进行绘制

void display()
{
glColor3f(r,v,b);
glBegin(GL_LINE_STRIP);
for(i=0;i<=samples,++i)
{
glVertex2f(x,ALPStd(x,l,m));
x+=step;
}
glEnd();
}

2、三维空间的Spherical Harmonics

在三维空间下绘制SH基函数，使用的是球坐标。









下面将会引进几个函数

//计算整数的阶乘，即x!=x*(x-1)*(x-2)*...*2*1
int factorial(int);
//给定带宽l和m，计算阶乘因子K
double evaluateK(int,int);
//给定带宽l和m，以及球坐标的两个角度，计算SH基函数的值
double evaluateSH(int,int,double,double);

int factorial(int x)
{
int result;
if(x==0 || x==-1)
return 1;
result=x;
while((x-=1)>0)
result*=x;

return result;
}

double evaluateK(int l,int m)
{
double result;
result=(2.0*l+1.1)*factorial(l-m)/(4*PI*factorial(l+m));
result=sqrt(result);
return result;
}

double evaluateSH(int l,int m,double theta,double phi)
{
double SH=0.0;
if(m==0)
SH=evaluateK(l,0)*ALPStd(cos(theta),l,0);
else if(m>0)
SH=sqrt(2)*evaluateK(l,m)*cos(m*phi)*ALPStd(cos(theta),l,m);
else
SH=sqrt(2)*evaluateK(l,-m)*sin(-m*phi)*ALPStd(cos(theta),l,-m);

return SH;
}

绘制SH的函数传递theta和phi值循环调用evaluateSH()，计算返回值，之后使用四边形绘制网格。

3、使用SH基函数重构复合函数

我们已经创建了SH的基函数，现在要使用它们来重构一个给定的复合函数(也就是求出因子c)，这个复合函数将会被用来描述光线在半球上的分布。





















下面将会引进两个函数。

//计算spherical stratified sampling
sphericalStratifiedSampling();
//将输入的函数投影到SH基上，即计算因子
SHProjectionSphericalFunction();

接下来还要定义一些数据结构。首先定义一个结构SHVector3d，容纳向量或顶点的笛卡尔坐标；同样定义一个结构，容纳分层取样(stratified sampling)的信息,即单位球坐标,样品的三维坐标，和SH因子。

typedef struct
{
double x,y,z;
}SHVector3d;

typedef struct
{
SHVector3d sph;
SHVector3d vec;
double *coeff;
}SHSample;

void sphericalStratifiedSampling(SHSample* samples,int sqrtNumSamples,
int nBands)
{
//Indexes
int a,b,index,l,m;

//Loop Index
int i=0;

//Inverse of the number of samples
double invNumSamples=1.0/sqrtNumSamples;

//Cartesian and spherical coordinates
double x,y,theta,phi;

//Loop on width of grid
for(a=0;a<sqrtNumSamples,++a)
{
//Loop on height of grid
for(b=0;b<sqrtNumSamples;++b)
{
//jitter center of current cell
x=(a+rnd())*invNumSamples;
y=(b+rnd())*invNumSamples;

//Calculate corresponding spherical angles
theta=2.0*acos(sqrt(1.0-x));
phi=2.0*PI*y;

//Sample spherical coordinates
samples[i].sph.x=theta;
samples[i].shp.y=phi;
samples[i].shp.z=1.0;

//Sample normal
samples[i].vec.x=sin(theta)*cos(phi);
samples[i].vec.y=sin(theta)*sin(phi);
samples[i].vec.z=cos(theta);

//Calculate SH coefficients of current sample
for(l=0;l<nBands;++l)
{
for(m=-l;m<=l;++m)
{
index=l*(l+1)+m;
samples[i].coeff[index]=evaluateSH(l,m,theta,phi);
}
}
++i;
}
}
}

//复合函数的定义
typedef double (*SphericalFunction)(double theta,double phi);

void SHProjectSphericalFunction(SphericalFunction mySPFunc,
SHSample* samples,double *result,int numSamples,int numCoeffs)
{
//Weighting factor
double dWeight=4.0*PI;
double factor=dWeight/numSamples;

//Spherical angles
double theta,phi;
//Loop indexes
int i,n;

//Loop through all samples
for(i=0;i<numSamples;++i)
{
//Get current sample spherical coordinates
theta=samples[i].sph.x;
phi=samples[i].sph.y;

//Update SH coefficients
for(n=0;n<numCoeffs;++n)
{
result
+=mySPFunc(theta,phi)*samples[i].coeff
;
}
}
for(n=0;n<numCoeffs;++n)
result
*=factor;
}