// This code is in the public domain -- castanyo@yahoo.es #ifndef NV_MATH_SPHERICALHARMONIC_H #define NV_MATH_SPHERICALHARMONIC_H #include "nvmath.h" #include // memcpy namespace nv { class Vector3; class Matrix; NVMATH_API float legendrePolynomial( int l, int m, float x ) NV_CONST; NVMATH_API float shBasis( int l, int m, float theta, float phi ) NV_CONST; NVMATH_API float shBasis( int l, int m, const Vector3 & v ) NV_CONST; NVMATH_API float hshBasis( int l, int m, float theta, float phi ) NV_CONST; NVMATH_API float hshBasis( int l, int m, const Vector3 & v ) NV_CONST; class Sh; float dot(const Sh & a, const Sh & b) NV_CONST; /// Spherical harmonic class. class Sh { friend class Sh2; friend class ShMatrix; public: /// Construct a spherical harmonic of the given order. Sh(int o) : m_order(o) { m_elemArray = new float[basisNum()]; } /// Copy constructor. Sh(const Sh & sh) : m_order(sh.order()) { m_elemArray = new float[basisNum()]; memcpy(m_elemArray, sh.m_elemArray, sizeof(float) * basisNum()); } /// Destructor. ~Sh() { delete [] m_elemArray; m_elemArray = NULL; } /// Get number of bands. static int bandNum(int m_order) { return m_order + 1; } /// Get number of sh basis. static int basisNum(int m_order) { return (m_order + 1) * (m_order + 1); } /// Get the index for the given coefficients. static int index( int l, int m ) { return l * l + l + m; } /// Get sh order. int order() const { return m_order; } /// Get sh order. int bandNum() const { return bandNum(m_order); } /// Get sh order. int basisNum() const { return basisNum(m_order); } /// Get sh coefficient indexed by l,m. float elem( int l, int m ) const { return m_elemArray[index(l, m)]; } /// Get sh coefficient indexed by l,m. float & elem( int l, int m ) { return m_elemArray[index(l, m)]; } /// Get sh coefficient indexed by i. float elemAt( int i ) const { return m_elemArray[i]; } /// Get sh coefficient indexed by i. float & elemAt( int i ) { return m_elemArray[i]; } /// Reset the sh coefficients. void reset() { for( int i = 0; i < basisNum(); i++ ) { m_elemArray[i] = 0.0f; } } /// Copy spherical harmonic. void operator= ( const Sh & sh ) { nvDebugCheck(order() <= sh.order()); for(int i = 0; i < basisNum(); i++) { m_elemArray[i] = sh.m_elemArray[i]; } } /// Add spherical harmonics. void operator+= ( const Sh & sh ) { nvDebugCheck(order() == sh.order()); for(int i = 0; i < basisNum(); i++) { m_elemArray[i] += sh.m_elemArray[i]; } } /// Substract spherical harmonics. void operator-= ( const Sh & sh ) { nvDebugCheck(order() == sh.order()); for(int i = 0; i < basisNum(); i++) { m_elemArray[i] -= sh.m_elemArray[i]; } } // Not exactly convolution, nor product. void operator*= ( const Sh & sh ) { nvDebugCheck(order() == sh.order()); for(int i = 0; i < basisNum(); i++) { m_elemArray[i] *= sh.m_elemArray[i]; } } /// Scale spherical harmonics. void operator*= ( float f ) { for(int i = 0; i < basisNum(); i++) { m_elemArray[i] *= f; } } /// Add scaled spherical harmonics. void addScaled( const Sh & sh, float f ) { nvDebugCheck(order() == sh.order()); for(int i = 0; i < basisNum(); i++) { m_elemArray[i] += sh.m_elemArray[i] * f; } } /*/// Add a weighted sample to the sh coefficients. void AddSample( const Vec3 & dir, const Color3f & color, float w=1.0f ) { for(int l = 0; l <= order; l++) { for(int m = -l; m <= l; m++) { Color3f & elem = GetElem(l, m); elem.Mad( elem, color, w * shBasis(l, m, dir) ); } } }*/ /// Evaluate void eval(const Vector3 & dir) { for(int l = 0; l <= m_order; l++) { for(int m = -l; m <= l; m++) { elem(l, m) = shBasis(l, m, dir); } } } /// Evaluate the spherical harmonic function. float sample(const Vector3 & dir) const { Sh sh(order()); sh.eval(dir); return dot(sh, *this); } protected: const int m_order; float * m_elemArray; }; /// Compute dot product of the spherical harmonics. inline float dot(const Sh & a, const Sh & b) { nvDebugCheck(a.order() == b.order()); float sum = 0; for( int i = 0; i < Sh::basisNum(a.order()); i++ ) { sum += a.elemAt(i) * b.elemAt(i); } return sum; } /// Second order spherical harmonic. class Sh2 : public Sh { public: /// Constructor. Sh2() : Sh(2) {} /// Copy constructor. Sh2(const Sh2 & sh) : Sh(sh) {} /// Spherical harmonic resulting from projecting the clamped cosine transfer function to the SH basis. void cosineTransfer() { const float c1 = 0.282095f; // K(0, 0) const float c2 = 0.488603f; // K(1, 0) const float c3 = 1.092548f; // sqrt(15.0f / PI) / 2.0f = K(2, -2) const float c4 = 0.315392f; // sqrt(5.0f / PI) / 4.0f) = K(2, 0) const float c5 = 0.546274f; // sqrt(15.0f / PI) / 4.0f) = K(2, 2) const float normalization = PI * 16.0f / 17.0f; const float const1 = c1 * normalization * 1.0f; const float const2 = c2 * normalization * (2.0f / 3.0f); const float const3 = c3 * normalization * (1.0f / 4.0f); const float const4 = c4 * normalization * (1.0f / 4.0f); const float const5 = c5 * normalization * (1.0f / 4.0f); m_elemArray[0] = const1; m_elemArray[1] = -const2; m_elemArray[2] = const2; m_elemArray[3] = -const2; m_elemArray[4] = const3; m_elemArray[5] = -const3; m_elemArray[6] = const4; m_elemArray[7] = -const3; m_elemArray[8] = const5; } }; /// Spherical harmonic matrix. class ShMatrix { public: /// Create an identity matrix of the given order. ShMatrix(int o = 2) : m_order(o), m_identity(true) { nvCheck(m_order > 0); m_e = new float[size()]; m_band = new float *[bandNum()]; setupBands(); } /// Destroy and free matrix elements. ~ShMatrix() { delete m_e; delete m_band; } /// Set identity matrix. void setIdentity() { m_identity = true; } /// Return true if this is an identity matrix, false in other case. bool isIdentity() const { return m_identity; } /// Get number of bands of this matrix. int bandNum() const { return m_order+1; } /// Get total number of elements in the matrix. int size() const { int size = 0; for (int i = 0; i < bandNum(); i++) { size += square(i * 2 + 1); } return size; } /// Get element at the given raw index. float element(int idx) const { return m_e[idx]; } /// Get element at the given with the given indices. float & element(int b, int x, int y) { nvDebugCheck(b >= 0); nvDebugCheck(b < bandNum()); return m_band[b][(b + y) * (b * 2 + 1) + (b + x)]; } /// Get element at the given with the given indices. float element(int b, int x, int y) const { nvDebugCheck(b >= 0); nvDebugCheck(b < bandNum()); return m_band[b][(b + y) * (b * 2 + 1) + (b + x)]; } /// Copy matrix. void copy(const ShMatrix & m) { nvDebugCheck(m_order == m.m_order); memcpy(m_e, m.m_e, size() * sizeof(float)); } /// Rotate the given coefficients. /*void transform( const Sh & restrict source, Sh * restrict dest ) const { nvCheck( &source != dest ); // Make sure there's no aliasing. nvCheck( dest->m_order <= m_order ); nvCheck( m_order <= source.m_order ); if (m_identity) { *dest = source; return; } // Loop through each band. for (int l = 0; l <= dest->m_order; l++) { for (int mo = -l; mo <= l; mo++) { Color3f rgb = Color3f::Black; for( int mi = -l; mi <= l; mi++ ) { rgb.Mad( rgb, source.elem(l, mi), elem(l, mo, mi) ); } dest->elem(l, mo) = rgb; } } }*/ NVMATH_API void multiply( const ShMatrix &A, const ShMatrix &B ); NVMATH_API void rotation( const Matrix & m ); NVMATH_API void rotation( int axis, float angles ); NVMATH_API void print(); private: // @@ These could be static indices precomputed only once. /// Setup the band pointers. void setupBands() { int size = 0; for( int i = 0; i < bandNum(); i++ ) { m_band[i] = &m_e[size]; size += square(i * 2 + 1); } } private: // Matrix order. const int m_order; // Identity flag for quick transform. bool m_identity; // Array of elements. float * m_e; // Band pointers. float ** m_band; }; } // nv namespace #endif // NV_MATH_SPHERICALHARMONIC_H