Create a gist now

Instantly share code, notes, and snippets.

@dwilliamson /SH.h
Last active Sep 8, 2016

What would you like to do?
Shared SH Library for C++/HLSL/OpenCL/CUDA using pycgen to generate/document optimal transforms
//
// References
//
// [1] Sloan, P.-P., Luna, B., and Snyder, J. 2005.
// Local, deformable precomputedradiance transfer.
// ACM Trans. Graph. 24, 3, 1216–1223
//
// [2] Stupid Spherical Harmonics (SH) Tricks
// Peter-Pike Sloan, Microsoft Corporation
//
#ifndef INCLUDED_SPHERICAL_HARMONICS_H
#define INCLUDED_SPHERICAL_HARMONICS_H
#include <cbpp/Math.h>
#if defined(__cplusplus)
namespace math
{
#endif
/*$pycgen
from math import pi
from math import sqrt
k_2sqrtpi = 2 * sqrt(pi)
k_4sqrtpi = 4 * sqrt(pi)
# Constants from [2], Appendix A2 Polynomial Forms of SH Basis
y = [
1.0 / k_2sqrtpi,
-sqrt(3) / k_2sqrtpi,
sqrt(3) / k_2sqrtpi,
-sqrt(3) / k_2sqrtpi,
sqrt(15) / k_2sqrtpi,
-sqrt(15) / k_2sqrtpi,
sqrt(5) / k_4sqrtpi,
-sqrt(15) / k_2sqrtpi,
sqrt(15) / k_4sqrtpi
]
*/
// Return single DC-term for order 0 SH
cmp_device_fn float SH_Y0()
{
/*$pycgen
code = """
return {y[0]}f;
"""
EmitFmt(code)
*/
//$pycgen-begin
return 0.282094791774f;
//$pycgen-end
}
// Calculate order 1 SH co-efficients along the given direction
cmp_device_fn void SH_Y1(float3 n, cmp_out float sh[4])
{
/*$pycgen
code = """
sh[0] = {y[0]}f;
sh[1] = {y[1]}f * n.y;
sh[2] = {y[2]}f * n.z;
sh[3] = {y[3]}f * n.x;
"""
EmitFmt(code)
*/
//$pycgen-begin
sh[0] = 0.282094791774f;
sh[1] = -0.488602511903f * n.y;
sh[2] = 0.488602511903f * n.z;
sh[3] = -0.488602511903f * n.x;
//$pycgen-end
}
// Calculate order 2 SH co-efficients along the given direction
cmp_device_fn void SH_Y2(float3 n, cmp_out float sh[9])
{
/*$pycgen
code += """
sh[4] = {y[4]}f * n.x * n.y;
sh[5] = {y[5]}f * n.y * n.z;
sh[6] = {y[6]}f * (3 * n.z * n.z - 1);
sh[7] = {y[7]}f * n.x * n.z;
sh[8] = {y[8]}f * (n.x * n.x - n.y * n.y);
"""
EmitFmt(code)
*/
//$pycgen-begin
sh[0] = 0.282094791774f;
sh[1] = -0.488602511903f * n.y;
sh[2] = 0.488602511903f * n.z;
sh[3] = -0.488602511903f * n.x;
sh[4] = 1.09254843059f * n.x * n.y;
sh[5] = -1.09254843059f * n.y * n.z;
sh[6] = 0.315391565253f * (3 * n.z * n.z - 1);
sh[7] = -1.09254843059f * n.x * n.z;
sh[8] = 0.546274215296f * (n.x * n.x - n.y * n.y);
//$pycgen-end
}
// Returns SH co-efficients at lm(0,0), lm(1,0) and lm(2, 0) for cosine lobe along positive z-axis.
// As the cosine lobe is circularly symmetric about the z-axis, all other SH co-efficients are zero.
// This allows the cosine lobe to be specified as ZH co-efficients.
cmp_device_fn float3 ZH_CosineLobe2()
{
// Projection of cosine lobe onto SH upper hemisphere
// Mathematica source from Sébastien Lagarde
// http://www.gamedev.net/topic/631913-cosine-lobe-in-spherical-harmonics/#entry4985124
/*
RealSH[L_, M_, \[Theta]_, \[Phi]_] := \[Sqrt]((2 L + 1) / (4 \[Pi]) (L - Abs[M])! / (L + Abs[M])!) *
LegendreP[L, Abs[M], Cos[\[Theta]]] *
If[M < 0, \[Sqrt]2 Cos[M.\[Phi]],
If[M > 0, \[Sqrt]2 Sin[-M.\[Phi]], 1]]
DiffuseTransfer[\[Theta]_, \[Phi]_] := Max[0, Cos[\[Theta]]]
ProjectOntoSH[F_, L_, M_] := \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(\[Pi]/2\)]\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(2 \[Pi]\)]\((RealSH[L,
M, \[Theta], \[Phi]] F[\[Theta], \[Phi]] Sin[\[Theta]])\) \
\[DifferentialD]\[Phi] \[DifferentialD]\[Theta]\)\)
Table[FullSimplify[ProjectOntoSH[DiffuseTransfer, L, 0], Exponent > 0], {L, 0, 3}]
*/
/*$pycgen
zh_cos = [
sqrt(pi) / 2,
sqrt(pi / 3),
sqrt(5 * pi) / 8
]
code = """
return float3_make({zh_cos[0]}f, {zh_cos[1]}f, {zh_cos[2]}f);
"""
EmitFmt(code)
*/
//$pycgen-begin
return float3_make(0.886226925453f, 1.02332670795f, 0.495415912201f);
//$pycgen-end
}
// Returns SH co-efficients for spherical light source that subtends an angle 'a',
// along the positive z-axis. Stored as zonal harmonics.
cmp_device_fn float3 ZH_SphericalCap2(float sin_a, float cos_a)
{
const float s = sin_a;
const float c = cos_a;
// Derived from [2], Appendix A3 ZH Coefficients for Spherical Light Source
// Can be derived in Mathematica by swapping DiffuseTransfer above for:
//
// SphericalCap[\[Theta]_, \[Phi]_, a_] := \[Theta] < a
float3 zh;
/*$pycgen
zh_sphc_k = [
-sqrt(pi),
0.5 * sqrt(3) * sqrt(pi),
-0.5 * sqrt(5) * sqrt(pi)
]
zh_sphc_t = [
"(-1 + c)",
"(s * s)",
"c * (-1 + c) * (c + 1)"
]
code = """
zh.x = {zh_sphc_k[0]}f * {zh_sphc_t[0]};
zh.y = {zh_sphc_k[1]}f * {zh_sphc_t[1]};
zh.z = {zh_sphc_k[2]}f * {zh_sphc_t[2]};
"""
EmitFmt(code)
*/
//$pycgen-begin
zh.x = -1.77245385091f * (-1 + c);
zh.y = 1.53499006192f * (s * s);
zh.z = -1.9816636488f * c * (-1 + c) * (c + 1);
//$pycgen-end
return zh;
}
// Given the SH function f and ZH function g, convolution of the two is defined as:
//
// f.g(s) = sum[lm] y_lm(s).k.g_l.f_lm (3)
//
// where k is the ZH convolution co-efficient.
//
// This can also be used to convert a ZH function to SH along an arbitrary axis, s:
//
// g*_lm = ylm(s).k.g_l (6)
//
// These equations are referenced from [1]
cmp_device_fn float ZH_ConvolutionCoefficient(float l)
{
/*$pycgen
k_4pi = 4 * pi
ZH_Convolve = lambda l: sqrt(k_4pi / (2 * l + 1))
# Complete equation (6) as much as possible before the direction gets applied
def ZH_ToSHScale(zh):
sh = [0] * 9
sh[0] = y[0] * zh[0] * ZH_Convolve(0)
sh[1] = y[1] * zh[1] * ZH_Convolve(1)
sh[2] = y[2] * zh[1] * ZH_Convolve(1)
sh[3] = y[3] * zh[1] * ZH_Convolve(1)
sh[4] = y[4] * zh[2] * ZH_Convolve(2)
sh[5] = y[5] * zh[2] * ZH_Convolve(2)
sh[6] = y[6] * zh[2] * ZH_Convolve(2)
sh[7] = y[7] * zh[2] * ZH_Convolve(2)
sh[8] = y[8] * zh[2] * ZH_Convolve(2)
return sh
EmitFmt("return sqrtf({k_4pi}f / (2l + 1));")
*/
//$pycgen-begin
return sqrtf(12.5663706144f / (2l + 1));
//$pycgen-end
}
cmp_device_fn float SH_CosineLobe0()
{
/*$pycgen
sh_cos = ZH_ToSHScale(zh_cos)
EmitFmt("return {sh_cos[0]}f;")
*/
//$pycgen-begin
return 0.886226925453f;
//$pycgen-end
}
cmp_device_fn void SH_CosineLobe1(float3 n, cmp_out float sh[4])
{
/*$pycgen
code = """
sh[0] = {sh_cos[0]}f;
sh[1] = {sh_cos[1]}f * n.y;
sh[2] = {sh_cos[2]}f * n.z;
sh[3] = {sh_cos[3]}f * n.x;
"""
EmitFmt(code)
*/
//$pycgen-begin
sh[0] = 0.886226925453f;
sh[1] = -1.02332670795f * n.y;
sh[2] = 1.02332670795f * n.z;
sh[3] = -1.02332670795f * n.x;
//$pycgen-end
}
cmp_device_fn void SH_CosineLobe2(float3 n, cmp_out float sh[9])
{
/*$pycgen
code += """
sh[4] = {sh_cos[4]}f * n.x * n.y;
sh[5] = {sh_cos[5]}f * n.y * n.z;
sh[6] = {sh_cos[6]}f * (3 * n.z * n.z - 1);
sh[7] = {sh_cos[7]}f * n.x * n.z;
sh[8] = {sh_cos[8]}f * (n.x * n.x - n.y * n.y);
"""
EmitFmt(code)
*/
//$pycgen-begin
sh[0] = 0.886226925453f;
sh[1] = -1.02332670795f * n.y;
sh[2] = 1.02332670795f * n.z;
sh[3] = -1.02332670795f * n.x;
sh[4] = 0.85808553081f * n.x * n.y;
sh[5] = -0.85808553081f * n.y * n.z;
sh[6] = 0.2477079561f * (3 * n.z * n.z - 1);
sh[7] = -0.85808553081f * n.x * n.z;
sh[8] = 0.429042765405f * (n.x * n.x - n.y * n.y);
//$pycgen-end
}
cmp_device_fn float SH_SphericalCap0()
{
/*$pycgen
sh_sphc = ZH_ToSHScale(zh_sphc_k)
EmitFmt("return {sh_sphc[0]}f;")
*/
//$pycgen-begin
return -1.77245385091f;
//$pycgen-end
}
cmp_device_fn void SH_SphericalCap1(float3 n, float sin_a, float cos_a, cmp_out float sh[4])
{
const float s = sin_a;
const float c = cos_a;
/*$pycgen
code = """
float zh_sphc_t0 = {zh_sphc_t[0]};
float zh_sphc_t1 = {zh_sphc_t[1]};
sh[0] = {sh_sphc[0]}f * zh_sphc_t0;
sh[1] = {sh_sphc[1]}f * n.y * zh_sphc_t1;
sh[2] = {sh_sphc[2]}f * n.z * zh_sphc_t1;
sh[3] = {sh_sphc[3]}f * n.x * zh_sphc_t1;
"""
EmitFmt(code)
*/
//$pycgen-begin
float zh_sphc_t0 = (-1 + c);
float zh_sphc_t1 = (s * s);
sh[0] = -1.77245385091f * zh_sphc_t0;
sh[1] = -1.53499006192f * n.y * zh_sphc_t1;
sh[2] = 1.53499006192f * n.z * zh_sphc_t1;
sh[3] = -1.53499006192f * n.x * zh_sphc_t1;
//$pycgen-end
}
cmp_device_fn void SH_SphericalCap2(float3 n, float sin_a, float cos_a, cmp_out float sh[9])
{
const float s = sin_a;
const float c = cos_a;
/*$pycgen
code += """
float zh_sphc_t2 = {zh_sphc_t[2]};
sh[4] = {sh_sphc[4]}f * n.x * n.y * zh_sphc_t2;
sh[5] = {sh_sphc[5]}f * n.y * n.z * zh_sphc_t2;
sh[6] = {sh_sphc[6]}f * (3 * n.z * n.z - 1) * zh_sphc_t2;
sh[7] = {sh_sphc[7]}f * n.x * n.z * zh_sphc_t2;
sh[8] = {sh_sphc[8]}f * (n.x * n.x - n.y * n.y) * zh_sphc_t2;
"""
EmitFmt(code)
*/
//$pycgen-begin
float zh_sphc_t0 = (-1 + c);
float zh_sphc_t1 = (s * s);
sh[0] = -1.77245385091f * zh_sphc_t0;
sh[1] = -1.53499006192f * n.y * zh_sphc_t1;
sh[2] = 1.53499006192f * n.z * zh_sphc_t1;
sh[3] = -1.53499006192f * n.x * zh_sphc_t1;
float zh_sphc_t2 = c * (-1 + c) * (c + 1);
sh[4] = -3.43234212324f * n.x * n.y * zh_sphc_t2;
sh[5] = 3.43234212324f * n.y * n.z * zh_sphc_t2;
sh[6] = -0.990831824402f * (3 * n.z * n.z - 1) * zh_sphc_t2;
sh[7] = 3.43234212324f * n.x * n.z * zh_sphc_t2;
sh[8] = -1.71617106162f * (n.x * n.x - n.y * n.y) * zh_sphc_t2;
//$pycgen-end
}
cmp_device_fn void SH_Scale2(cmp_out float sh_result[9], float scale)
{
sh_result[0] *= scale;
sh_result[1] *= scale;
sh_result[2] *= scale;
sh_result[3] *= scale;
sh_result[4] *= scale;
sh_result[5] *= scale;
sh_result[6] *= scale;
sh_result[7] *= scale;
sh_result[8] *= scale;
}
cmp_device_fn void SH_AddScaled2(cmp_out float sh_result[9], float sh[9], float scale)
{
sh_result[0] += sh[0] * scale;
sh_result[1] += sh[1] * scale;
sh_result[2] += sh[2] * scale;
sh_result[3] += sh[3] * scale;
sh_result[4] += sh[4] * scale;
sh_result[5] += sh[5] * scale;
sh_result[6] += sh[6] * scale;
sh_result[7] += sh[7] * scale;
sh_result[8] += sh[8] * scale;
}
cmp_device_fn float SH_Convolve1(float a[4], float b[4])
{
return
a[0] * b[0] +
a[1] * b[1] +
a[2] * b[2] +
a[3] * b[3];
}
cmp_device_fn float SH_Convolve2(float a[9], float b[9])
{
return
a[0] * b[0] +
a[1] * b[1] +
a[2] * b[2] +
a[3] * b[3] +
a[4] * b[4] +
a[5] * b[5] +
a[6] * b[6] +
a[7] * b[7] +
a[8] * b[8];
}
#if defined(__cplusplus)
}
#endif
#endif // INCLUDED_SPHERICAL_HARMONICS_H
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment