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



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

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

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

1 2 3 4 5 6 7 8 9 10 11

//勒让德多项式计算方法 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));         }

1 2 3 4 5 6 7 8 9 10

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

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

1 2 3 4 5 6 7 8 9 10 11

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基函数，使用的是球坐标。

下面将会引进几个函数

1 2 3 4 5 6

//计算整数的阶乘，即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 );

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

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)，这个复合函数将会被用来描述光线在半球上的分布。

下面将会引进两个函数。

1 2 3 4

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

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

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; }