Created
May 7, 2019 09:48
Revisions
-
codewings created this gist
May 7, 2019 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,390 @@ # 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() """