You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
246 lines
6.3 KiB
C++
246 lines
6.3 KiB
C++
// This code is in the public domain -- castanyo@yahoo.es
|
|
|
|
#include "SphericalHarmonic.h"
|
|
|
|
#include "Vector.h"
|
|
|
|
using namespace nv;
|
|
|
|
|
|
namespace
|
|
{
|
|
|
|
// Basic integer factorial.
|
|
inline static int factorial( int v )
|
|
{
|
|
const static int fac_table[] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800 };
|
|
|
|
if(v <= 11){
|
|
return fac_table[v];
|
|
}
|
|
|
|
int result = v;
|
|
while (--v > 0) {
|
|
result *= v;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
|
|
// Double factorial.
|
|
// Defined as: n!! = n*(n - 2)*(n - 4)..., n!!(0,-1) = 1.
|
|
inline static int doubleFactorial( int x )
|
|
{
|
|
if (x == 0 || x == -1) {
|
|
return 1;
|
|
}
|
|
|
|
int result = x;
|
|
while ((x -= 2) > 0) {
|
|
result *= x;
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
/// Normalization constant for spherical harmonic.
|
|
/// @param l is the band.
|
|
/// @param m is the argument, in the range [0, m]
|
|
inline static float K( int l, int m )
|
|
{
|
|
nvDebugCheck( m >= 0 );
|
|
return sqrtf(((2 * l + 1) * factorial(l - m)) / (4 * PI * factorial(l + m)));
|
|
}
|
|
|
|
/// Normalization constant for hemispherical harmonic.
|
|
inline static float HK( int l, int m )
|
|
{
|
|
nvDebugCheck( m >= 0 );
|
|
return sqrtf(((2 * l + 1) * factorial(l - m)) / (2 * PI * factorial(l + m)));
|
|
}
|
|
|
|
/// Evaluate Legendre polynomial. */
|
|
static float legendre( int l, int m, float x )
|
|
{
|
|
// piDebugCheck( m >= 0 );
|
|
// piDebugCheck( m <= l );
|
|
// piDebugCheck( fabs(x) <= 1 );
|
|
|
|
// Rule 2 needs no previous results
|
|
if (l == m) {
|
|
return powf(-1.0f, float(m)) * doubleFactorial(2 * m - 1) * powf(1 - x*x, 0.5f * m);
|
|
}
|
|
|
|
// Rule 3 requires the result for the same argument of the previous band
|
|
if (l == m + 1) {
|
|
return x * (2 * m + 1) * legendrePolynomial(m, m, x);
|
|
}
|
|
|
|
// Main reccurence used by rule 1 that uses result of the same argument from
|
|
// the previous two bands
|
|
return (x * (2 * l - 1) * legendrePolynomial(l - 1, m, x) - (l + m - 1) * legendrePolynomial(l - 2, m, x)) / (l - m);
|
|
}
|
|
|
|
|
|
template <int l, int m> float legendre(float x);
|
|
|
|
template <> float legendre<0, 0>(float ) {
|
|
return 1;
|
|
}
|
|
|
|
template <> float legendre<1, 0>(float x) {
|
|
return x;
|
|
}
|
|
template <> float legendre<1, 1>(float x) {
|
|
return -sqrtf(1 - x * x);
|
|
}
|
|
|
|
template <> float legendre<2, 0>(float x) {
|
|
return -0.5f + (3 * x * x) / 2;
|
|
}
|
|
template <> float legendre<2, 1>(float x) {
|
|
return -3 * x * sqrtf(1 - x * x);
|
|
}
|
|
template <> float legendre<2, 2>(float x) {
|
|
return -3 * (-1 + x * x);
|
|
}
|
|
|
|
template <> float legendre<3, 0>(float x) {
|
|
return -(3 * x) / 2 + (5 * x * x * x) / 2;
|
|
}
|
|
template <> float legendre<3, 1>(float x) {
|
|
return -3 * sqrtf(1 - x * x) / 2 * (-1 + 5 * x * x);
|
|
}
|
|
template <> float legendre<3, 2>(float x) {
|
|
return -15 * (-x + x * x * x);
|
|
}
|
|
template <> float legendre<3, 3>(float x) {
|
|
return -15 * powf(1 - x * x, 1.5f);
|
|
}
|
|
|
|
template <> float legendre<4, 0>(float x) {
|
|
return 0.125f * (3.0f - 30.0f * x * x + 35.0f * x * x * x * x);
|
|
}
|
|
template <> float legendre<4, 1>(float x) {
|
|
return -2.5f * x * sqrtf(1.0f - x * x) * (7.0f * x * x - 3.0f);
|
|
}
|
|
template <> float legendre<4, 2>(float x) {
|
|
return -7.5f * (1.0f - 8.0f * x * x + 7.0f * x * x * x * x);
|
|
}
|
|
template <> float legendre<4, 3>(float x) {
|
|
return -105.0f * x * powf(1 - x * x, 1.5f);
|
|
}
|
|
template <> float legendre<4, 4>(float x) {
|
|
return 105.0f * (x * x - 1.0f) * (x * x - 1.0f);
|
|
}
|
|
|
|
} // namespace
|
|
|
|
|
|
float nv::legendrePolynomial(int l, int m, float x)
|
|
{
|
|
switch(l)
|
|
{
|
|
case 0:
|
|
return legendre<0, 0>(x);
|
|
case 1:
|
|
if(m == 0) return legendre<1, 0>(x);
|
|
return legendre<1, 1>(x);
|
|
case 2:
|
|
if(m == 0) return legendre<2, 0>(x);
|
|
else if(m == 1) return legendre<2, 1>(x);
|
|
return legendre<2, 2>(x);
|
|
case 3:
|
|
if(m == 0) return legendre<3, 0>(x);
|
|
else if(m == 1) return legendre<3, 1>(x);
|
|
else if(m == 2) return legendre<3, 2>(x);
|
|
return legendre<3, 3>(x);
|
|
case 4:
|
|
if(m == 0) return legendre<4, 0>(x);
|
|
else if(m == 1) return legendre<4, 1>(x);
|
|
else if(m == 2) return legendre<4, 2>(x);
|
|
else if(m == 3) return legendre<4, 3>(x);
|
|
else return legendre<4, 4>(x);
|
|
}
|
|
|
|
// Fallback to the expensive version.
|
|
return legendre(l, m, x);
|
|
}
|
|
|
|
|
|
/**
|
|
* Evaluate the spherical harmonic function for the given angles.
|
|
* @param l is the band.
|
|
* @param m is the argument, in the range [-l,l]
|
|
* @param theta is the altitude, in the range [0, PI]
|
|
* @param phi is the azimuth, in the range [0, 2*PI]
|
|
*/
|
|
float nv::shBasis( int l, int m, float theta, float phi )
|
|
{
|
|
if( m == 0 ) {
|
|
// K(l, 0) = sqrt((2*l+1)/(4*PI))
|
|
return sqrtf((2 * l + 1) / (4 * PI)) * legendrePolynomial(l, 0, cosf(theta));
|
|
}
|
|
else if( m > 0 ) {
|
|
return sqrtf(2.0f) * K(l, m) * cosf(m * phi) * legendrePolynomial(l, m, cosf(theta));
|
|
}
|
|
else {
|
|
return sqrtf(2.0f) * K(l, -m) * sinf(-m * phi) * legendrePolynomial(l, -m, cosf(theta));
|
|
}
|
|
}
|
|
|
|
|
|
/**
|
|
* Real spherical harmonic function of an unit vector. Uses the following
|
|
* equalities to call the angular function:
|
|
* x = sin(theta)*cos(phi)
|
|
* y = sin(theta)*sin(phi)
|
|
* z = cos(theta)
|
|
*/
|
|
float nv::shBasis( int l, int m, Vector3::Arg v )
|
|
{
|
|
float theta = acosf(v.z);
|
|
float phi = atan2f(v.y, v.x);
|
|
return shBasis( l, m, theta, phi );
|
|
}
|
|
|
|
|
|
/**
|
|
* Evaluate the hemispherical harmonic function for the given angles.
|
|
* @param l is the band.
|
|
* @param m is the argument, in the range [-l,l]
|
|
* @param theta is the altitude, in the range [0, PI/2]
|
|
* @param phi is the azimuth, in the range [0, 2*PI]
|
|
*/
|
|
float nv::hshBasis( int l, int m, float theta, float phi )
|
|
{
|
|
if( m == 0 ) {
|
|
// HK(l, 0) = sqrt((2*l+1)/(2*PI))
|
|
return sqrtf((2 * l + 1) / (2 * PI)) * legendrePolynomial(l, 0, 2*cosf(theta)-1);
|
|
}
|
|
else if( m > 0 ) {
|
|
return sqrtf(2.0f) * HK(l, m) * cosf(m * phi) * legendrePolynomial(l, m, 2*cosf(theta)-1);
|
|
}
|
|
else {
|
|
return sqrtf(2.0f) * HK(l, -m) * sinf(-m * phi) * legendrePolynomial(l, -m, 2*cosf(theta)-1);
|
|
}
|
|
}
|
|
|
|
|
|
/**
|
|
* Real hemispherical harmonic function of an unit vector. Uses the following
|
|
* equalities to call the angular function:
|
|
* x = sin(theta)*cos(phi)
|
|
* y = sin(theta)*sin(phi)
|
|
* z = cos(theta)
|
|
*/
|
|
float nv::hshBasis( int l, int m, Vector3::Arg v )
|
|
{
|
|
float theta = acosf(v.z);
|
|
float phi = atan2f(v.y, v.x);
|
|
return hshBasis( l, m, theta, phi );
|
|
}
|
|
|
|
|
|
|