## October 31, 2010

### Python Implementation of Spherical Harmonics Stratified Sampling

Filed under: Uncategorized — admin @ 7:32 pm

Its Sunday afternoon and its time to write some code, this is pretty much lifted verbatim from Robin Green’s 2003 Paper: Spherical Harmonic Lighting the Gritty Details

So for the next trick is to put this as a large table into a shader so I can bake out a set of 9 cooefficients (4 bands) for Shadowed Diffuse Transfer, these I will store per sample in a point cloud to be looked up during a shading stage.

Then after that I can implement the Image Based Lights, then I can do really quick Image Based lighting using Spherical Harmonics

Lots of fun!

I made a big mess of this but luckily Markus Kransler was able to fix it up:

Here is the amended code:

```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 #!/usr/bin/env python   class SHSample(): sph=(0.0,0.0) vec=(0.0,0.0,0.0) coeff={} pass   def P(l,m,x): import math #Associated Legendre Polynomial P(l,m,x) at x pmm = 1.0 if m > 0: somx2=math.sqrt(1.0-(x*x)) fact = 1.0 for i in xrange(1,m+1,1): pmm *= (-fact)*somx2 fact += 2.0 if l == m: return pmm   pmmp1 = x * ((2.0*m)+1.0)*pmm   if l == m+1: return pmmp1   plm = 0.0   for ll in xrange(m+2,l+1,1): plm = ((2.0*ll - 1.0) * x * pmmp1 - (ll + m - 1.0) * pmm) / (ll - m); pmm = pmmp1 pmmp1 = plm   return plm   def K(l,m): import math temp = float((((2.0*l)+1.0)*math.factorial(l-m))/(4.0*math.pi*math.factorial(l+m))) return math.sqrt(temp)   def SH(l,m,theta,phi): import math sqrt2 = math.sqrt(2.0) if m==0: return K(l,0)*P(l,0,math.cos(theta)) elif m > 0: return sqrt2*K(l,m)*math.cos(m*phi)*P(l,m,math.cos(theta)) else: return sqrt2*K(l,-m)*math.sin(-m*phi)*P(l,-m,math.cos(theta))     def setupSamples(sqrtNumSamples=64,numBands=4): import random,math counter = 0 oneOverN = 1.0/float(sqrtNumSamples) samples = [SHSample() for i in range(sqrtNumSamples*sqrtNumSamples)]   for i in range(sqrtNumSamples): for j in range(sqrtNumSamples): x = (i+ random.random())*oneOverN theta = 2.0*math.acos(math.sqrt(1-x))   y = (j+ random.random())*oneOverN phi = 2.0*math.pi*y   samples[counter].sph=(theta,phi)   vec = (math.sin(theta)*math.cos(phi),\ math.sin(theta)*math.sin(phi),\ math.cos(theta)) samples[counter].vec = vec   tmpDict = {} for l in range(numBands): for m in xrange(-l,l+1,1): index = l*(l+1)+m sh= SH(l,m,theta,phi) tmpDict[index]=sh samples[counter].coeff=tmpDict counter +=1 return samples     for i in setupSamples(): print i.coeff```