Commit f1c3d4ddbe054fcc883fcbcc680dae0d406f3453
1 parent
5795a81c
added spherical harmonics evaluation header file
Showing
1 changed file
with
72 additions
and
0 deletions
Show diff stats
| 1 | +/* Solving the associated legendre polynomials iteratively. Most of this code was taken from | ||
| 2 | + * the work described by Robin Green from Sony Computer Entertainment America: | ||
| 3 | + * | ||
| 4 | + * R. Green, Spherical Harmonic Lighting: The Gritty Details, Sony Computer Entertainment America, 2003 | ||
| 5 | + * http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf | ||
| 6 | + */ | ||
| 7 | + | ||
| 8 | + | ||
| 9 | +/// This function computes the factorial n! | ||
| 10 | +/// where n! = n * (n-1) * (n-2) * (n-3) * ...... * 2 * 1 | ||
| 11 | +unsigned long factorial(unsigned int n) | ||
| 12 | +{ | ||
| 13 | + if (n == 0) //the standard definition of a factorial specifies that 0! = 1 | ||
| 14 | + return 1; | ||
| 15 | + return n * factorial(n - 1); //recursively call this function | ||
| 16 | +} | ||
| 17 | + | ||
| 18 | +/// Evaluate an Associated Legendre Polynomial P(l,m,x) at x | ||
| 19 | +/// This implements evaluation of the associated Legendre polynomials using | ||
| 20 | +/// a recursive method with the following rules: | ||
| 21 | +/// | ||
| 22 | +/// 1) (l - m)P(l, m) = x(2l-1)P(l-1, m) - (l + m - 1)P(l-2, m) | ||
| 23 | +/// 2) P(m, m) = (-1)^m (2m - 1)!! (1 - x^2)^{m/2} | ||
| 24 | +/// 3) P(m+1, m) = x (2m + 1)P(m, m) | ||
| 25 | +double P(int l,int m,double x) | ||
| 26 | +{ | ||
| 27 | + | ||
| 28 | + double pmm = 1.0; //first iteration, P(0, 0) = 1 | ||
| 29 | + if(m>0) { //solve for positive values of m | ||
| 30 | + double somx2 = sqrt((1.0-x)*(1.0+x)); | ||
| 31 | + double fact = 1.0; | ||
| 32 | + for(int i=1; i<=m; i++) { | ||
| 33 | + pmm *= (-fact) * somx2; | ||
| 34 | + fact += 2.0; | ||
| 35 | + } | ||
| 36 | + } | ||
| 37 | + if(l==m) return pmm; //if the degree equals the order, the solution is found above | ||
| 38 | + double pmmp1 = x * (2.0*m+1.0) * pmm; | ||
| 39 | + if(l==m+1) return pmmp1; | ||
| 40 | + double pll = 0.0; | ||
| 41 | + for(int ll=m+2; ll<=l; ++ll) { | ||
| 42 | + pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m); | ||
| 43 | + pmm = pmmp1; | ||
| 44 | + pmmp1 = pll; | ||
| 45 | + } | ||
| 46 | + return pll; | ||
| 47 | +} | ||
| 48 | + | ||
| 49 | +/// Calculate the scaling factor for a normalized Legendre polynomial | ||
| 50 | +/// l is the degree, ranging from [0..N] | ||
| 51 | +/// m is the order in the range [-l..l] | ||
| 52 | +double K(int l, int m) | ||
| 53 | +{ | ||
| 54 | + static const double K_PI = 3.14159; | ||
| 55 | + // renormalisation constant for SH function | ||
| 56 | + double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*K_PI*factorial(l+m)); | ||
| 57 | + return sqrt(temp); | ||
| 58 | +} | ||
| 59 | + | ||
| 60 | +/// return a point sample of a Spherical Harmonic basis function | ||
| 61 | +/// l is the degree, range [0..N] | ||
| 62 | +/// m is the order in the range [-l..l] | ||
| 63 | +/// theta in the range [0..Pi] | ||
| 64 | +/// phi in the range [0..2*Pi] | ||
| 65 | +double SH(int l, int m, double theta, double phi){ | ||
| 66 | + | ||
| 67 | + const double sqrt2 = sqrt(2.0); | ||
| 68 | + if(m==0) return K(l,0)*P(l,m,cos(theta)); | ||
| 69 | + else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta)); | ||
| 70 | + else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta)); | ||
| 71 | +} | ||
| 72 | + |