Compute spherical harmonics from cube maps. Work in progress.

src/nvmath/SphericalHarmonic.h

@@ -4,6 +4,7 @@
 #define NV_MATH_SPHERICALHARMONIC_H
 
 #include "nvmath.h"
+#include "Vector.h"
 
 #include <string.h> // 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++) {
 

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;

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);
             }

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);