From f1c3d4ddbe054fcc883fcbcc680dae0d406f3453 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Tue, 11 Aug 2015 09:49:32 -0500 Subject: [PATCH] added spherical harmonics evaluation header file --- programming/random/cpp/sh.h | 72 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+), 0 deletions(-) create mode 100644 programming/random/cpp/sh.h diff --git a/programming/random/cpp/sh.h b/programming/random/cpp/sh.h new file mode 100644 index 0000000..9160697 --- /dev/null +++ b/programming/random/cpp/sh.h @@ -0,0 +1,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)); +} + -- libgit2 0.21.4