Port from google/spherical-harmonics [https://github.com/google/spherical-harmonics]

PythonSH.py

# porting from https:#github.com/google/spherical-harmonics for blender [WIP]
import math, mathutils
from random import random, uniform


PI = math.pi
TWO_PI = 2.0 * PI
FACTORIAL_CACHE = [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000]
DOUBLE_FACTORIAL_CACHE = [1, 1, 2, 3, 8, 15, 48, 105, 384, 945, 3840, 10395, 46080, 135135, 645120, 2027025]
FACTORIAL_CACHESIZE = 16


# Hardcoded spherical harmonic functions for low orders (l is first number
# and m is second number (sign encoded as preceeding 'p' or 'n')).
#
# As polynomials they are evaluated more efficiently in cartesian coordinates,
# assuming that @d is unit. This is not verified for efficiency.
def HardcodedSH00(d):
    # 0.5 * sqrt(1/pi)
    return 0.282095

def HardcodedSH1n1(d):
    # -sqrt(3/(4pi)) * y
    return -0.488603 * d.y

def HardcodedSH10(d):
    # sqrt(3/(4pi)) * z
    return 0.488603 * d.z

def HardcodedSH1p1(d):
    # -sqrt(3/(4pi)) * x
    return -0.488603 * d.x

def HardcodedSH2n2(d):
    # 0.5 * sqrt(15/pi) * x * y
    return 1.092548 * d.x * d.y

def HardcodedSH2n1(d):
    # -0.5 * sqrt(15/pi) * y * z
    return -1.092548 * d.y * d.z

def HardcodedSH20(d):
    # 0.25 * sqrt(5/pi) * (-x^2-y^2+2z^2)
    return 0.315392 * (-d.x * d.x - d.y * d.y + 2.0 * d.z * d.z)

def HardcodedSH2p1(d):
    # -0.5 * sqrt(15/pi) * x * z
    return -1.092548 * d.x * d.z

def HardcodedSH2p2(d):
    # 0.25 * sqrt(15/pi) * (x^2 - y^2)
    return 0.546274 * (d.x * d.x - d.y * d.y)

def HardcodedSH3n3(d):
    # -0.25 * sqrt(35/(2pi)) * y * (3x^2 - y^2)
    return -0.590044 * d.y * (3.0 * d.x * d.x - d.y * d.y)

def HardcodedSH3n2(d):
    # 0.5 * sqrt(105/pi) * x * y * z
    return 2.890611 * d.x * d.y * d.z

def HardcodedSH3n1(d):
    # -0.25 * sqrt(21/(2pi)) * y * (4z^2-x^2-y^2)
    return -0.457046 * d.y * (4.0 * d.z * d.z - d.x * d.x - d.y * d.y)

def HardcodedSH30(d):
    # 0.25 * sqrt(7/pi) * z * (2z^2 - 3x^2 - 3y^2)
    return 0.373176 * d.z * (2.0 * d.z * d.z - 3.0 * d.x * d.x - 3.0 * d.y * d.y)

def HardcodedSH3p1(d):
    # -0.25 * sqrt(21/(2pi)) * x * (4z^2-x^2-y^2)
    return -0.457046 * d.x * (4.0 * d.z * d.z - d.x * d.x - d.y * d.y)

def HardcodedSH3p2(d):
    # 0.25 * sqrt(105/pi) * z * (x^2 - y^2)
    return 1.445306 * d.z * (d.x * d.x - d.y * d.y)

def HardcodedSH3p3(d):
    # -0.25 * sqrt(35/(2pi)) * x * (x^2-3y^2)
    return -0.590044 * d.x * (d.x * d.x - 3.0 * d.y * d.y)

def HardcodedSH4n4(d):
    # 0.75 * sqrt(35/pi) * x * y * (x^2-y^2)
    return 2.503343 * d.x * d.y * (d.x * d.x - d.y * d.y)

def HardcodedSH4n3(d):
    # -0.75 * sqrt(35/(2pi)) * y * z * (3x^2-y^2)
    return -1.770131 * d.y * d.z * (3.0 * d.x * d.x - d.y * d.y)

def HardcodedSH4n2(d):
    # 0.75 * sqrt(5/pi) * x * y * (7z^2-1)
    return 0.946175 * d.x * d.y * (7.0 * d.z * d.z - 1.0)

def HardcodedSH4n1(d):
    # -0.75 * sqrt(5/(2pi)) * y * z * (7z^2-3)
    return -0.669047 * d.y * d.z * (7.0 * d.z * d.z - 3.0)

def HardcodedSH40(d):
    # 3/16 * sqrt(1/pi) * (35z^4-30z^2+3)
    z2 = d.z * d.z
    return 0.105786 * (35.0 * z2 * z2 - 30.0 * z2 + 3.0)

def HardcodedSH4p1(d):
    # -0.75 * sqrt(5/(2pi)) * x * z * (7z^2-3)
    return -0.669047 * d.x * d.z * (7.0 * d.z * d.z - 3.0)

def HardcodedSH4p2(d):
    # 3/8 * sqrt(5/pi) * (x^2 - y^2) * (7z^2 - 1)
    return 0.473087 * (d.x * d.x - d.y * d.y) * (7.0 * d.z * d.z - 1.0)

def HardcodedSH4p3(d):
    # -0.75 * sqrt(35/(2pi)) * x * z * (x^2 - 3y^2)
    return -1.770131 * d.x * d.z * (d.x * d.x - 3.0 * d.y * d.y)

def HardcodedSH4p4(d):
    # 3/16*sqrt(35/pi) * (x^2 * (x^2 - 3y^2) - y^2 * (3x^2 - y^2))
    x2 = d.x * d.x
    y2 = d.y * d.y
    return 0.625836 * (x2 * (x2 - 3.0 * y2) - y2 * (3.0 * x2 - y2))


# Return floating mod x % m.
def FastFMod(x, m):
    return x - (m * math.floor(x / m))

#
def Abs(x):
    return int(math.fabs(x))

# Get the total number of coefficients for a function represented by
# all spherical harmonic basis of degree <= @order (it is a point of
# confusion that the order of an SH refers to its degree and not the order).
def SHCoeffCount(order):
    return (order + 1) * (order + 1)

# Get the one dimensional index associated with a particular degree @l
# and order @m. This is the index that can be used to access the Coeffs
# returned by SHSolver.
def SHIndex(l, m):
    return l * (l + 1) + m

# Compute the factorial for an integer @x. It is assumed x is at least 0.
# This implementation precomputes the results for low values of x, in which
# case this is a constant time lookup.
#
# The vast majority of SH evaluations will hit these precomputed values.
def Factorial(x):
    if x < FACTORIAL_CACHESIZE:
        return float(FACTORIAL_CACHE[x])
    else:
        s = 1.0
        for n in range(2, x+1):
            s *= n

        return s

# Compute the double factorial for an integer @x. This assumes x is at least
# 0.  This implementation precomputes the results for low values of x, in
# which case this is a constant time lookup.
#
# The vast majority of SH evaluations will hit these precomputed values.
# See http:#mathworld.wolfram.com/DoubleFactorial.html
def DoubleFactorial(x):
    if x < FACTORIAL_CACHESIZE:
        return float(DOUBLE_FACTORIAL_CACHE[x])
    else:
        s = 1.0
        n = x
        while n > 1.0:
            s *= n
            n -= 2.0

        return s

# Evaluate the associated Legendre polynomial of degree @l and order @m at
# coordinate @x. The inputs must satisfy:
# 1. l >= 0
# 2. 0 <= m <= l
# 3. -1 <= x <= 1
# See http:#en.wikipedia.org/wiki/Associated_Legendre_polynomials
#
# This implementation is based off the approach described in [1],
# instead of computing Pml(x) directly, Pmm(x) is computed. Pmm can be
# lifted to Pmm+1 recursively until Pml is found
def EvalLegendrePolynomial(l, m, x):
    # Compute Pmm(x) = (-1)^m(2m - 1)!!(1 - x^2)^(m/2), where !! is the double factorial.
    pmm = 1.0
    # P00 is defined as 1.0, do don't evaluate Pmm unless we know m > 0
    if m > 0:
        sign = 1 if m % 2 == 0 else -1
        pmm = sign * DoubleFactorial(2 * m - 1) * math.pow(1 - x * x, m / 2.0)

    if l == m:
        # Pml is the same as Pmm so there's no lifting to higher bands needed
        return pmm

    # Compute Pmm+1(x) = x(2m + 1)Pmm(x)
    pmm1 = x * (2 * m + 1) * pmm
    if l == m + 1:
        # Pml is the same as Pmm+1 so we are done as well
        return pmm1

    # Use the last two computed bands to lift up to the next band until l is
    # reached, using the recurrence relationship:
    # Pml(x) = (x(2l - 1)Pml-1 - (l + m - 1)Pml-2) / (l - m)
    for n in range(m+2, l+1):
        pmn = (x * (2 * n - 1) * pmm1 - (n + m - 1) * pmm) / (n - m)
        pmm = pmm1
        pmm1 = pmn

    # Pmm1 at the end of the above loop is equal to Pml
    return pmm1

#
def Spherical2Cartesian(phi, theta):
    r = math.sin(theta)
    return [r * math.cos(phi), r * math.sin(phi), math.cos(theta)]

#
def Cartesian2Spherical(dir):
    assert math.fabs(dir.length - 1.0) < 1e-5
    # Explicitly clamp the z coordinate so that numeric errors don't cause it
    # to fall just outside of acos' domain.
    theta = math.acos(max(min(dir.z, 1.0), -1.0))
    # We don't need to divide dir.y or dir.x by sin(theta) since they are
    # both scaled by it and atan2 will handle it appropriately.
    phi = math.atan2(dir.y, dir.x)

    return theta, phi

#
def ImageXtoPhi(x, width):
    # The directions are measured from the center of the pixel, so add 0.5
    # to convert from integer pixel indices to float pixel coordinates.
    return TWO_PI * (x + 0.5) / width

#
def ImageYtoTheta(y, height):
    return PI * (y + 0.5) / height

#
def SphericaltoImageXY(phi, theta, width, height):
    # Allow theta to repeat and map to 0 to pi. However, to account for cases
    # where y goes beyond the normal 0 to pi range, phi may need to be adjusted.
    theta = max(min(FastFMod(theta, TWO_PI), TWO_PI), 0.0)
    if theta > PI:
        # theta is out of bounds. Effectively, theta has rotated past the pole
        # so after adjusting theta to be in range, rotating phi by pi forms an
        # equivalent direction.
        theta = TWO_PI - theta  # now theta is between 0 and pi
        phi += PI

    # Allow phi to repeat and map to the normal 0 to 2pi range.
    # Clamp and map after adjusting theta in case theta was forced to update phi.
    phi = max(min(FastFMod(phi, TWO_PI), TWO_PI), 0.0)

    # Now phi is in [0, 2pi] and theta is in [0, pi] so it's simple to inverse
    # the linear equations in ImageCoordsToSphericalCoords, although there's no
    # -0.5 because we're returning floating point coordinates and so don't need
    # to center the pixel.
    return (width * phi / TWO_PI), (height * theta / PI)

# According to Archimedes' Theorem, when we project the surface of a cylinder 
# to a sphere, the area preserved. So, we can distribute the points on the cylinder, 
# then map those points onto a sphere. This is much easier approach that handling 
# the solid angle (ϕ) dependency.
def UniformSphereDistributing(alpha, beta):
    theta = 2.0 * PI * alpha
    z = beta
    x = math.sqrt(1.0-z * z) * math.cos(theta)
    y = math.sqrt(1.0-z * z) * math.sin(theta)
    return [x,y,z]

#
def EvalSHValue(l, m, phi, theta):
    assert l >= 0
    assert -l <= m and m <= l

    kml = math.sqrt((2.0 * l + 1) * Factorial(l - Abs(m)) / (4.0 * PI * Factorial(l + Abs(m))))
    
    if m > 0:
        return math.sqrt(2.0) * kml * math.cos(m * phi) * EvalLegendrePolynomial(l, m, math.cos(theta))
    else:
        if m < 0:
            return math.sqrt(2.0) * kml * math.sin(-m * phi) * EvalLegendrePolynomial(l, -m, math.cos(theta))
        else:
            return kml * EvalLegendrePolynomial(l, 0, math.cos(theta))

# 
def EvalSH(l, m, phi, theta):
    return EvalSHValue(l, m, phi, theta), SHIndex(l, m)

#
def ProjectFunctionToSH(order, fn, samples):
    assert order >= 0
    assert samples > 0

    # This is the approach demonstrated in [1] and is useful for arbitrary
    # functions on the sphere that are represented analytically.
    side = int(math.floor(math.sqrt(samples)))
    coeff = SHCoeffCount(order) * [0.0]

    for t in range(side):
        for p in range(side):
            alpha = (t + uniform(0.0, 1.0)) / side
            beta = (p + uniform(0.0, 1.0)) / side

            # See http:#www.bogotobogo.com/Algorithms/uniform_distribution_sphere.php
            phi = TWO_PI * beta
            theta = math.acos(2.0 * alpha - 1.0)

            # evaluate the analytic function for the current spherical coords
            f = fn(phi, theta)

            # evaluate the SH basis functions up to band O, scale them by the
            # function's value and accumulate them over all generated samples
            for l in range(order+1):
                for m in range(-l, l+1):
                    sh, index = EvalSH(l, m, phi, theta)
                    coeff[index] += f * sh

    # scale by the probability of a particular sample, which is
    # 4pi/sample_side^2. 4pi for the surface area of a unit sphere, and
    # 1/sample_side^2 for the number of samples drawn uniformly.
    weight = 4.0 * PI / (side * side)
    for i in range(len(coeff)):
        coeff[i] *= weight

    return coeff

#
def ProjectEnvImageToSH(order, fn, width, height):
    assert order >= 0

    # An environment map projection is three different spherical functions, one
    # for each color channel. The projection integrals are estimated by
    # iterating over every pixel within the image.
    pixel_area = (TWO_PI / width) * (PI / height)
    coeffs = SHCoeffCount(order) *[[0,0,0]]

    for t in range(height):
        theta = ImageYtoTheta(t, height)
        # The differential area of each pixel in the map is constant across a
        # row. Must scale the pixel_area by sin(theta) to account for the
        # stretching that occurs at the poles with this parameterization.
        weight = pixel_area * math.sin(theta)
        for p in range(width):
            phi = ImageXtoPhi(p, width)
            color = fn(p, t)

            for l in range(order+1):
                for m in range(-l, l+1):
                    sh, index = EvalSH(l, m, phi, theta)
                    coeffs[index][0] += sh * weight * color[0]
                    coeffs[index][1] += sh * weight * color[1]
                    coeffs[index][2] += sh * weight * color[2]
    return coeffs

#
def EvalSHSum(order, coeffs, phi, theta):
    assert SHCoeffCount(order) == len(coeffs)

    sum = 0.0
    for l in range(order + 1):
        for m in range(-l, l+1):
            sh, index = EvalSH(l, m, phi, theta)
            sum += sh * coeffs[index]

    return sum

# simple test
"""
#
def TEST_ProjectFunction_Fn(phi, theta):
    coeffs = [-1.028, 0.779, -0.275, 0.601, -0.256, 1.891, -1.658, -0.370, -0.772]
    return EvalSHSum(2, coeffs, phi, theta)

#
def TEST_ProjectFunction():
    # The expected coefficients used to define the analytic spherical function
    expected_coeffs = [-1.028, 0.779, -0.275, 0.601, -0.256, 1.891, -1.658, -0.370, -0.772]

    coeffs = ProjectFunctionToSH(2, TEST_ProjectFunction_Fn, 5000)

    for i in range(9):
        if math.fabs(coeffs[i] - expected_coeffs[i]) > 1e-5:
            print('Failed ' + str(coeffs[i]) + ' ' + str(expected_coeffs[i]))

TEST_ProjectFunction()
"""