sh.h
2.38 KB
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
/* Solving the associated Legendre polynomials iteratively. Most of this code was taken from
* the work described by Robin Green from Sony Computer Entertainment America:
*
* R. Green, Spherical Harmonic Lighting: The Gritty Details, Sony Computer Entertainment America, 2003
* http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf
*/
/// This function computes the factorial n!
/// where n! = n * (n-1) * (n-2) * (n-3) * ...... * 2 * 1
unsigned long factorial(unsigned int n)
{
if (n == 0) //the standard definition of a factorial specifies that 0! = 1
return 1;
return n * factorial(n - 1); //recursively call this function
}
/// Evaluate an Associated Legendre Polynomial P(l,m,x) at x
/// This implements evaluation of the associated Legendre polynomials using
/// a recursive method with the following rules:
///
/// 1) (l - m)P(l, m) = x(2l-1)P(l-1, m) - (l + m - 1)P(l-2, m)
/// 2) P(m, m) = (-1)^m (2m - 1)!! (1 - x^2)^{m/2}
/// 3) P(m+1, m) = x (2m + 1)P(m, m)
double P(int l,int m,double x)
{
double pmm = 1.0; //first iteration, P(0, 0) = 1
if(m>0) { //solve for positive values of m
double somx2 = sqrt((1.0-x)*(1.0+x));
double fact = 1.0;
for(int i=1; i<=m; i++) {
pmm *= (-fact) * somx2;
fact += 2.0;
}
}
if(l==m) return pmm; //if the degree equals the order, the solution is found above
double pmmp1 = x * (2.0*m+1.0) * pmm;
if(l==m+1) return pmmp1;
double pll = 0.0;
for(int ll=m+2; ll<=l; ++ll) {
pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
pmm = pmmp1;
pmmp1 = pll;
}
return pll;
}
/// Calculate the scaling factor for a normalized Legendre polynomial
/// l is the degree, ranging from [0..N]
/// m is the order in the range [-l..l]
double K(int l, int m)
{
static const double K_PI = 3.14159;
// renormalisation constant for SH function
double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*K_PI*factorial(l+m));
return sqrt(temp);
}
/// return a point sample of a Spherical Harmonic basis function
/// l is the degree, range [0..N]
/// m is the order in the range [-l..l]
/// theta in the range [0..Pi]
/// phi in the range [0..2*Pi]
double SH(int l, int m, double theta, double phi){
const double sqrt2 = sqrt(2.0);
if(m==0) return K(l,0)*P(l,m,cos(theta));
else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
}