Python Implementation of Spherical Harmonics Stratified Sampling

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:

#!/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