From c591c5f8b4c40d0a9ee37b95bb740a4bb483a8e4 Mon Sep 17 00:00:00 2001 From: Ignacio Date: Thu, 31 Jan 2019 18:06:02 -0800 Subject: [PATCH] Compute spherical harmonics from cube maps. Work in progress. --- src/nvmath/SphericalHarmonic.h | 132 +++++++++++++++++++-------------- src/nvtt/CubeSurface.cpp | 70 +++++++++++++++-- src/nvtt/CubeSurface.h | 3 + src/nvtt/nvtt.h | 4 + 4 files changed, 144 insertions(+), 65 deletions(-) diff --git a/src/nvmath/SphericalHarmonic.h b/src/nvmath/SphericalHarmonic.h index b7f5774..f2eb741 100644 --- a/src/nvmath/SphericalHarmonic.h +++ b/src/nvmath/SphericalHarmonic.h @@ -4,6 +4,7 @@ #define NV_MATH_SPHERICALHARMONIC_H #include "nvmath.h" +#include "Vector.h" #include // memcpy @@ -31,33 +32,33 @@ namespace nv public: /// Construct a spherical harmonic of the given order. - Sh(int o) : m_order(o) + Sh(int o) : order(o) { - m_elemArray = new float[basisNum()]; + coef = new float[basisNum()]; } /// Copy constructor. - Sh(const Sh & sh) : m_order(sh.order()) + Sh(const Sh & sh) : order(sh.order) { - m_elemArray = new float[basisNum()]; - memcpy(m_elemArray, sh.m_elemArray, sizeof(float) * basisNum()); + coef = new float[basisNum()]; + memcpy(coef, sh.coef, sizeof(float) * basisNum()); } /// Destructor. ~Sh() { - delete [] m_elemArray; - m_elemArray = NULL; + delete [] coef; + coef = NULL; } /// Get number of bands. - static int bandNum(int m_order) { - return m_order + 1; + static int bandNum(int order) { + return order + 1; } /// Get number of sh basis. - static int basisNum(int m_order) { - return (m_order + 1) * (m_order + 1); + static int basisNum(int order) { + return (order + 1) * (order + 1); } /// Get the index for the given coefficients. @@ -65,46 +66,40 @@ namespace nv return l * l + l + m; } - /// Get sh order. - int order() const - { - return m_order; - } - /// Get sh order. int bandNum() const { - return bandNum(m_order); + return bandNum(order); } /// Get sh order. int basisNum() const { - return basisNum(m_order); + return basisNum(order); } /// Get sh coefficient indexed by l,m. float elem( int l, int m ) const { - return m_elemArray[index(l, m)]; + return coef[index(l, m)]; } /// Get sh coefficient indexed by l,m. float & elem( int l, int m ) { - return m_elemArray[index(l, m)]; + return coef[index(l, m)]; } /// Get sh coefficient indexed by i. float elemAt( int i ) const { - return m_elemArray[i]; + return coef[i]; } /// Get sh coefficient indexed by i. float & elemAt( int i ) { - return m_elemArray[i]; + return coef[i]; } @@ -112,47 +107,47 @@ namespace nv void reset() { for( int i = 0; i < basisNum(); i++ ) { - m_elemArray[i] = 0.0f; + coef[i] = 0.0f; } } /// Copy spherical harmonic. void operator= ( const Sh & sh ) { - nvDebugCheck(order() <= sh.order()); + nvDebugCheck(order <= sh.order); for(int i = 0; i < basisNum(); i++) { - m_elemArray[i] = sh.m_elemArray[i]; + coef[i] = sh.coef[i]; } } /// Add spherical harmonics. void operator+= ( const Sh & sh ) { - nvDebugCheck(order() == sh.order()); + nvDebugCheck(order == sh.order); for(int i = 0; i < basisNum(); i++) { - m_elemArray[i] += sh.m_elemArray[i]; + coef[i] += sh.coef[i]; } } /// Substract spherical harmonics. void operator-= ( const Sh & sh ) { - nvDebugCheck(order() == sh.order()); + nvDebugCheck(order == sh.order); for(int i = 0; i < basisNum(); i++) { - m_elemArray[i] -= sh.m_elemArray[i]; + coef[i] -= sh.coef[i]; } } // Not exactly convolution, nor product. void operator*= ( const Sh & sh ) { - nvDebugCheck(order() == sh.order()); + nvDebugCheck(order == sh.order); for(int i = 0; i < basisNum(); i++) { - m_elemArray[i] *= sh.m_elemArray[i]; + coef[i] *= sh.coef[i]; } } @@ -160,17 +155,17 @@ namespace nv void operator*= ( float f ) { for(int i = 0; i < basisNum(); i++) { - m_elemArray[i] *= f; + coef[i] *= f; } } /// Add scaled spherical harmonics. void addScaled( const Sh & sh, float f ) { - nvDebugCheck(order() == sh.order()); + nvDebugCheck(order == sh.order); for(int i = 0; i < basisNum(); i++) { - m_elemArray[i] += sh.m_elemArray[i] * f; + coef[i] += sh.coef[i] * f; } } @@ -188,7 +183,7 @@ namespace nv /// Evaluate void eval(const Vector3 & dir) { - for(int l = 0; l <= m_order; l++) { + for(int l = 0; l <= order; l++) { for(int m = -l; m <= l; m++) { elem(l, m) = shBasis(l, m, dir); } @@ -199,17 +194,15 @@ namespace nv /// Evaluate the spherical harmonic function. float sample(const Vector3 & dir) const { - Sh sh(order()); + Sh sh(order); sh.eval(dir); return dot(sh, *this); } - protected: - - const int m_order; - float * m_elemArray; + const int order; + float * coef; }; @@ -217,10 +210,10 @@ namespace nv /// Compute dot product of the spherical harmonics. inline float dot(const Sh & a, const Sh & b) { - nvDebugCheck(a.order() == b.order()); + nvDebugCheck(a.order == b.order); float sum = 0; - for( int i = 0; i < Sh::basisNum(a.order()); i++ ) { + for( int i = 0; i < Sh::basisNum(a.order); i++ ) { sum += a.elemAt(i) * b.elemAt(i); } @@ -239,9 +232,34 @@ namespace nv /// Copy constructor. Sh2(const Sh2 & sh) : Sh(sh) {} + // Fast evaluation from: PPS' Efficient Spherical Harmonic Evaluation http://jcgt.org/published/0002/02/06/ + void eval(const Vector3 & dir) { + float fZ2 = dir.z * dir.z; + coef[0] = 0.2820947917738781f; + coef[2] = 0.4886025119029199f * dir.z; + coef[6] = 0.9461746957575601f * fZ2 + -0.3153915652525201f; + + float fC0 = dir.x; + float fS0 = dir.y; + + float fTmpA = -0.48860251190292f; + coef[3] = fTmpA * fC0; + coef[1] = fTmpA * fS0; + + float fTmpB = -1.092548430592079f * dir.z; + coef[7] = fTmpB * fC0; + coef[5] = fTmpB * fS0; + + float fC1 = dir.x * fC0 - dir.y * fS0; + float fS1 = dir.x * fS0 + dir.y * fC0; + + float fTmpC = 0.5462742152960395f; + coef[8] = fTmpC * fC1; + coef[4] = fTmpC * fS1; + } + /// Spherical harmonic resulting from projecting the clamped cosine transfer function to the SH basis. - void cosineTransfer() - { + 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) @@ -256,17 +274,17 @@ namespace nv const float const4 = c4 * normalization * (1.0f / 4.0f); const float const5 = c5 * normalization * (1.0f / 4.0f); - m_elemArray[0] = const1; + coef[0] = const1; - m_elemArray[1] = -const2; - m_elemArray[2] = const2; - m_elemArray[3] = -const2; + coef[1] = -const2; + coef[2] = const2; + coef[3] = -const2; - m_elemArray[4] = const3; - m_elemArray[5] = -const3; - m_elemArray[6] = const4; - m_elemArray[7] = -const3; - m_elemArray[8] = const5; + coef[4] = const3; + coef[5] = -const3; + coef[6] = const4; + coef[7] = -const3; + coef[8] = const5; } }; @@ -352,8 +370,8 @@ namespace nv /// 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 ); + nvCheck( dest->order <= order ); + nvCheck( order <= source.order ); if (m_identity) { *dest = source; @@ -361,7 +379,7 @@ namespace nv } // Loop through each band. - for (int l = 0; l <= dest->m_order; l++) { + for (int l = 0; l <= dest->order; l++) { for (int mo = -l; mo <= l; mo++) { diff --git a/src/nvtt/CubeSurface.cpp b/src/nvtt/CubeSurface.cpp index f79aa14..2b42a89 100644 --- a/src/nvtt/CubeSurface.cpp +++ b/src/nvtt/CubeSurface.cpp @@ -529,6 +529,24 @@ void CubeSurface::clamp(int channel, float low/*= 0.0f*/, float high/*= 1.0f*/) CubeSurface CubeSurface::irradianceFilter(int size, EdgeFixup fixupMethod) const { + // Evaluate spherical harmonic for each output texel. + CubeSurface output; + output.m->allocate(size); + + Sh2 shr, shg, shb; + + computeIrradianceSH3(0, shr.coef); + computeIrradianceSH3(1, shg.coef); + computeIrradianceSH3(2, shb.coef); + + // @@ Sample spherical harmonic from every direction. + + return CubeSurface(); +} + + +void CubeSurface::computeLuminanceIrradianceSH3(float coef[9]) const{ + m->allocateTexelTable(); // Transform this cube to spherical harmonic basis @@ -537,6 +555,10 @@ CubeSurface CubeSurface::irradianceFilter(int size, EdgeFixup fixupMethod) const // For each texel of the input cube. const uint edgeLength = m->edgeLength; for (uint f = 0; f < 6; f++) { + + const Surface & inputFace = m->face[f]; + const FloatImage * inputImage = inputFace.m->image; + for (uint y = 0; y < edgeLength; y++) { for (uint x = 0; x < edgeLength; x++) { @@ -546,24 +568,56 @@ CubeSurface CubeSurface::irradianceFilter(int size, EdgeFixup fixupMethod) const Sh2 shDir; shDir.eval(dir); - sh.addScaled(sh, solidAngle); + float r = inputImage->pixel(0, x, y, 0); + float g = inputImage->pixel(1, x, y, 0); + float b = inputImage->pixel(2, x, y, 0); + float lum = 0.333f * (r + g + b); // @@ use the proper luminance formula. + + sh.addScaled(shDir, lum * solidAngle); } } } + for (int i = 0; i < 9; i++) { + coef[i] = sh.coef[i]; + } +} - // Evaluate spherical harmonic for each output texel. - CubeSurface output; - output.m->allocate(size); +void CubeSurface::computeIrradianceSH3(int channel, float coef[9]) const { + m->allocateTexelTable(); + // Transform this cube to spherical harmonic basis + Sh2 sh; - // @@ TODO - return CubeSurface(); -} + // For each texel of the input cube. + const uint edgeLength = m->edgeLength; + for (uint f = 0; f < 6; f++) { + const Surface & inputFace = m->face[f]; + const FloatImage * inputImage = inputFace.m->image; + + for (uint y = 0; y < edgeLength; y++) { + for (uint x = 0; x < edgeLength; x++) { + Vector3 dir = m->texelTable->direction(f, x, y); + float solidAngle = m->texelTable->solidAngle(f, x, y); + + Sh2 shDir; + shDir.eval(dir); + + float c = inputImage->pixel(channel, x, y, 0); + + sh.addScaled(shDir, c * solidAngle); + } + } + } + + for (int i = 0; i < 9; i++) { + coef[i] = sh.elemAt(i); + } +} // Convolve filter against this cube. @@ -832,7 +886,7 @@ CubeSurface CubeSurface::cosinePowerFilter(int size, float cosinePower, EdgeFixu CubeSurface filteredCube; filteredCube.m->allocate(size); - // Texel table is stored along with the surface so that it's compute only once. + // Texel table is stored along with the surface so that it's computed only once. m->allocateTexelTable(); const float threshold = 0.001f; diff --git a/src/nvtt/CubeSurface.h b/src/nvtt/CubeSurface.h index b5e3757..ed1b6c4 100644 --- a/src/nvtt/CubeSurface.h +++ b/src/nvtt/CubeSurface.h @@ -88,6 +88,9 @@ namespace nvtt void allocateTexelTable() { + if (edgeLength == 0) { + edgeLength = face[0].width(); + } if (texelTable == NULL) { texelTable = new TexelTable(edgeLength); } diff --git a/src/nvtt/nvtt.h b/src/nvtt/nvtt.h index 303a1b4..1ca48ae 100644 --- a/src/nvtt/nvtt.h +++ b/src/nvtt/nvtt.h @@ -665,6 +665,10 @@ namespace nvtt NVTT_API CubeSurface fastResample(int size, EdgeFixup fixupMethod) const; + // Spherical Harmonics: + NVTT_API void computeLuminanceIrradianceSH3(float sh[9]) const; + NVTT_API void computeIrradianceSH3(int channel, float sh[9]) const; + /* NVTT_API void resize(int w, int h, ResizeFilter filter); NVTT_API void resize(int w, int h, ResizeFilter filter, float filterWidth, const float * params = 0);