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[n]+=mySPFunc(theta,phi)*samples[i].coeff[n];
          }        
      }
      for(n=0;n<numCoeffs;++n)
         result[n]*=factor;
}