www.pudn.com > ZFXMath-latest.zip > SphericalHarmonics.h
/// \file
///
/// \if DE
/// @brief Basis-SH Klassen
/// \else
/// @brief Basic SH Classes
/// \endif
#ifndef _ZFXMATH_INCLUDE_SH_H_
#define _ZFXMATH_INCLUDE_SH_H_
namespace ZFXMath
{
/// \if DE
/// @brief enthält Basis SH Klassen
/// \else
/// @brief contains Basic SH Classes
///
/// SphericalHarmonics::TSamples
/// SphericalHarmonics::TSphericalFunction
/// SphericalHarmonics::TCoeffs
/// SphericalHarmonics::Project
///
/// Based on Spherical Harmonic Rotation
/// Don Williamson, July 2004
/// http://www.donw.co.uk
/// \endif
namespace SphericalHarmonics
{
/// \if DE
/// @brief Punkt auf eienr Einheitssphäre
/// \else
/// @brief Point on a unit Sphere
/// \endif
template
struct TSample
{
// Default constructor
TSample() { };
// Construct from spherical co-ordinates
TSample( const PrecisionType u, const PrecisionType v ) :
theta(u),
phi(v)
{
x = PrecisionType( sin(theta) * cos(phi) );
y = PrecisionType( sin(theta) * sin(phi) );
z = PrecisionType( cos(theta) );
}
/// Spherical co-ordinates
PrecisionType theta, phi;
/// Cartesian co-ordinates
PrecisionType x, y, z;
};
/// \if DE
/// @brief Spherical Function Interface
///
/// Dessen Implementation ist der Weg eine beliebige Beleuchtungsumgebung in SH-Koeffizienten zu komprimieren.
/// \else
/// @brief Spherical Function Interface
///
/// Its implementation is the way to compress an arbitrary lighting situation into SH Coefficients.
/// \endif
template
struct TSphericalFunction
{
/// \if DE
/// Ermittelt den Funktionswert am gegebenen Punkt
/// \else
/// Evaluates the function at the given sample point
/// \endif
virtual FuncValueType EvalFunction(const TSample& s) const = 0;
};
/// \if DE
/// @brief beinhaltet eine SH-kompimierte sphärische Funktion
///
/// T = PrecisionType/FuncValueType
/// \else
/// @brief Contains a SH comressed spherical function
///
/// T = PrecisionType/FuncValueType
/// \endif
template
class TCoeffs
{
public:
/// \if DE
/// @brief Speicherzuteiler für Koeffizienten
/// \else
/// @brief Allocator for coefficients
/// \endif
template
class TAllocator
{
public:
TAllocator( int size )
{
m_Size = size;
m_pMemory = new TAlloc[size];
m_pNextMem = m_pMemory;
}
~TAllocator()
{
delete[] m_pMemory;
}
TAlloc* GetMem( int amount )
{
m_Size -= amount;
assert( m_Size >= 0 );
TAlloc* pMem = m_pNextMem;
m_pNextMem += amount;
return pMem;
}
private:
TAlloc* m_pMemory;
TAlloc* m_pNextMem;
int m_Size;
};
typedef TAllocator Allocator;
TCoeffs() : m_Coeffs(0), m_NbBands(0), m_NbCoeffs(0) { }
TCoeffs(const int nb_bands) :
m_Coeffs(0),
m_NbBands(nb_bands),
m_NbCoeffs(m_NbBands * m_NbBands),
m_FromAllocator(false)
{
// Allocate the coeff list
m_Coeffs = new T[m_NbCoeffs];
}
TCoeffs(const int nb_bands, Allocator* pAlloc) :
m_Coeffs(0),
m_NbCoeffs(nb_bands * nb_bands),
m_FromAllocator(true)
{
// Allocate the coeff list
m_Coeffs = pAlloc->GetMem( m_NbCoeffs );
}
~TCoeffs()
{
// Stack allocator cleans itself up
if (!m_FromAllocator)
delete [] m_Coeffs;
}
TCoeffs(const TCoeffs& other) :
m_Coeffs(0),
m_NbBands(other.m_NbBands),
m_NbCoeffs(other.m_NbCoeffs),
m_FromAllocator(other.m_FromAllocator)
{
// Just copy the pointer in this case
if (m_FromAllocator)
m_Coeffs = other.m_Coeffs;
else
{
// Make a copy when not using the stack allocator
m_Coeffs = new T[m_NbCoeffs];
for (int i = 0; i < m_NbCoeffs; i++)
m_Coeffs[i] = other.m_Coeffs[i];
}
}
// Const/non-const accessors for band/arg
T operator () (const int l, const int m) const
{
return (m_Coeffs[Check(l * (l + 1) + m)]);
}
T& operator () (const int l, const int m)
{
return (m_Coeffs[Check(l * (l + 1) + m)]);
}
// Const/non-const accessors by index
T operator () (const int i) const
{
return (m_Coeffs[Check(i)]);
}
T& operator () (const int i)
{
return (m_Coeffs[Check(i)]);
}
int GetNbBands(void) const
{
return (m_NbBands);
}
int GetSize(void) const
{
return (m_NbCoeffs);
}
private:
inline int Check(const int index) const
{
// Check bounds in debug build
assert(index >= 0 && index < m_NbCoeffs);
return (index);
}
// List of SH coefficients
T* m_Coeffs;
// Number of bands used
int m_NbBands;
// Number of coefficients (=bands^2)
int m_NbCoeffs;
// Allocated from a special-purpose allocator?
bool m_FromAllocator;
};
/// \if DE
///
/// \else
// http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.html
// Project a spherical function onto the spherical harmonic bases
//
// a. Spherical Function to project
// b. List of sample locations on the sphere
// c. List of y(l,m,theta,phi) for each sample
// d. Number of samples on the sphere
// e. Destination projection
/// \endif
template
void Project(const TSphericalFunction& func,
const TSample* samples,
const TCoeffs* coeffs,
const int nb_samples,
TCoeffs& dest )
{
int nb_coeffs = dest.GetSize();
// Check input
assert(samples && coeffs);
assert(nb_samples >= 0);
assert(coeffs[0].GetSize() == nb_coeffs);
// Clear out sums
for (int i = 0; i < nb_coeffs; i++)
dest(i) = 0;
for (int i = 0; i < nb_samples; i++)
{
// Take the sample at this point on the sphere
const TSample& s = samples[i];
const TCoeffs& c = coeffs[i];
FuncValueType sample = func.EvalFunction(s);
// Sum the projection of this sample onto each SH basis
for (int j = 0; j < nb_coeffs; j++)
dest(j) += sample * FuncValueType( c(j) );
}
// Divide each coefficient by the number of samples and multiply by weights
for (int i = 0; i < nb_coeffs; i++)
dest(i) = dest(i) * FuncValueType( PrecisionType( (4 * PI) / nb_samples ) );
}
//---------------------------------------------------------------------------------------
// Evaluate Associated Legendre Polynomial
template
PrecisionType P(const int l, const int m, const PrecisionType x)
{
// Start with P(0,0) at 1
PrecisionType pmm = 1;
// First calculate P(m,m) since that is the only rule that requires results
// from previous bands
// No need to check for m>0 since SH function always gives positive m
// Precalculate (1 - x^2)^0.5
PrecisionType somx2 = (PrecisionType)sqrt(1 - x * x);
// This calculates P(m,m). There are three terms in rule 2 that are being iteratively multiplied:
//
// 0: -1^m
// 1: (2m-1)!!
// 2: (1-x^2)^(m/2)
//
// Term 2 has been partly precalculated and the iterative multiplication by itself m times
// completes the term.
// The result of 2m-1 is always odd so the double factorial calculation multiplies every odd
// number below 2m-1 together. So, term 3 is calculated using the 'fact' variable.
PrecisionType fact = 1;
for (int i = 1; i <= m; i++)
{
pmm *= -1 * fact * somx2;
fact += 2;
}
// No need to go any further, rule 2 is satisfied
if (l == m)
return (pmm);
// Since m
PrecisionType y(const int l, const int m, const PrecisionType theta, const PrecisionType phi)
{
if (0 == m)
return (PrecisionType(K(l, 0) * P(l, 0, (PrecisionType)cos(theta))));
if (m > 0)
return PrecisionType(sqrt(2.0) * K(l, m) * cos(m * phi) * P(l, m, (PrecisionType)cos(theta)));
// m < 0, m is negated in call to K
return PrecisionType(sqrt(2.0) * K(l, -m) * sin(-m * phi) * P(l, -m, (PrecisionType)cos(theta)));
}
// Re-normalisation constant for SH function
template
PrecisionType K(const int l, const int m)
{
struct FactorialGen
{
FactorialGen(const int count) :
factorials(0),
nb_factorials(count)
{
factorials = new int[nb_factorials];
// Build a factorial LUT
factorials[0] = 1;
for (int i = 1; i < nb_factorials; i++)
factorials[i] = i * factorials[i - 1];
}
~FactorialGen(void)
{
delete [] factorials;
}
int operator () (const int i) const
{
// Check input for debug builds
assert(i >= 0 && i < nb_factorials);
return (factorials[i]);
}
int* factorials;
int nb_factorials;
};
// Generate lut the first time this function is called
// Note that this method doesn't interact well with explicit init/shutdown memory managers
// Additionally, 32-bits is not enough to store any higher bands (7 max) !!!
static FactorialGen factorial(14);
// Note that |m| is not used here as the SH function always passes positive m
return PrecisionType(sqrt(((2 * l + 1) * factorial(l - m)) / (4 * PI * factorial(l + m))));
}
}; // namespace SphericalHarmonics
}; // namespace ZFXMath
#endif //_ZFXMATH_INCLUDE_SH_H_