sunLoadingImage
whowedImag
decoration left 1
decoration left 2
transhome
transprojects
transgallery
transarticles
decoration rigth
English

Show/Hide search bar
black cat logo variable logo
[12 Sep 2013]

Compute coefficients of spherical harmonics from cube texture (lighting map)

This tutorial describes how to compute coefficients for first three bands of spherical harmonics from arbitrary cube texture. Coefficients of spherical harmonics are then used to reconstruct diffuse and ambient lighting in a shader. Process of calculation of coefficients for spherical harmonics is a lot quicker than calculation of blurred cube texture for image-based lighting. Due to this, spherical harmonics can be used for real-time dynamic updates to parameters of lighting after changes to a scene. To reconstruct diffuse lighting from first three bands of spherical harmonics you need only nine coefficients (for each basis function). As coefficients are calculated separately for each RGB channel, there are 27 floating point numbers to pass to a shader (9 vec3 values).
Function sphericalHarmonicsFromTexture() computes coefficients of spherical harmonics from cube texture. First it allocates and initializes required amount of memory to save coefficients and intermediate data. Then it binds input cube texture as current OpenGL texture. For each face of cube texture, data about dimensions and texels is acquired. Then, for each texel, function computes normalized direction from the imaginary center of the cube texture to the texel (on the unit cube) and scale factor for the texel. Direction from current texel is projected on each basis function of spherical harmonics. Then function retrieves the color of current texel and maps value from 8 bit representation to floating point representation in the range [0, 1]. Color (it is treated as intensity) is scaled depending on the distance from current texel to the center of cube map's face (to properly take into account greater ratio of texels to unit of space in places that are closer to edges of the cube texture. Calculated coefficients for each RGB channel are then saved to the appropriate buffer. In the end, accumulated coefficients of spherical harmonics are normalized to the range, appropriate for reconstruction of lighting.
Function colorFromSphericalHarmonics() implements reconstruction of diffuse lighting from coefficients of spherical harmonics. Function takes 9 coefficients of spherical harmonics, that are calculated by sphericalHarmonicsFromTexture() function.
Function to calculate coefficients of spherical harmonics from cube texture:
void sphericalHarmonicsFromTexture(TextureCube * cubeTexture,
                                    std::vector<glm::vec3> & output, const uint order)
{
   const uint sqOrder = order*order;

   // allocate memory for calculations
   output.resize(sqOrder);
   std::vector<float> resultR(sqOrder);
   std::vector<float> resultG(sqOrder);
   std::vector<float> resultB(sqOrder);

   // variables that describe current face of cube texture
   std::unique_ptr data;
   GLint width, height;
   GLint internalFormat;
   GLint numComponents;

   // initialize values
   float fWt = 0.0f;
   for(uint i=0; i < sqOrder; i++)
   {
      output[i].x = 0;
      output[i].y = 0;
      output[i].z = 0;
      resultR[i] = 0;
      resultG[i] = 0;
      resultB[i] = 0;
   }
   std::vector<float> shBuff(sqOrder);
   std::vector<float> shBuffB(sqOrder);

   // bind current texture
   glBindTexture(GL_TEXTURE_CUBE_MAP, cubeTexture->texture());
   // for each face of cube texture
   for(int face=0; face < 6; face++)
   {
      // get width and height
      glGetTexLevelParameteriv(cubeSides[face], 0, GL_TEXTURE_WIDTH, &width);
      glGetTexLevelParameteriv(cubeSides[face], 0, GL_TEXTURE_HEIGHT, &height);

      if(width != height)
      {
         return;
      }

      // get format of data in texture
      glGetTexLevelParameteriv(cubeSides[face], 0,
            GL_TEXTURE_INTERNAL_FORMAT, &internalFormat);

      // get data from texture
      if(internalFormat == GL_RGBA)
      {
         numComponents = 4;
         data = std::unique_ptr(new GLubyte[numComponents * width * width]);
      }
      else if(internalFormat == GL_RGB)
      {
         numComponents = 3;
         data = std::unique_ptr(new GLubyte[numComponents * width * width]);
      }
      else
      {
         return;
      }
      glGetTexImage(cubeSides[face], 0, internalFormat, GL_UNSIGNED_BYTE, data.get());

      // step between two texels for range [0, 1]
      float invWidth = 1.0f / float(width);
      // initial negative bound for range [-1, 1]
      float negativeBound = -1.0f + invWidth;
      // step between two texels for range [-1, 1]
      float invWidthBy2 = 2.0f / float(width);

      for(int y=0; y < width; y++)
      {
         // texture coordinate V in range [-1 to 1]
         const float fV = negativeBound + float(y) * invWidthBy2;

         for(int x=0; x < width; x++)
         {
            // texture coordinate U in range [-1 to 1]
            const float fU = negativeBound + float(x) * invWidthBy2;

            // determine direction from center of cube texture to current texel
            glm::vec3 dir;
            switch(cubeSides[face])
            {
               case GL_TEXTURE_CUBE_MAP_POSITIVE_X:
                  dir.x = 1.0f;
                  dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
                  dir.z = 1.0f - (invWidthBy2 * float(x) + invWidth);
                  dir = -dir;
                  break;
               case GL_TEXTURE_CUBE_MAP_NEGATIVE_X:
                  dir.x = -1.0f;
                  dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
                  dir.z = -1.0f + (invWidthBy2 * float(x) + invWidth);
                  dir = -dir;
                  break;
               case GL_TEXTURE_CUBE_MAP_POSITIVE_Y:
                  dir.x = - 1.0f + (invWidthBy2 * float(x) + invWidth);
                  dir.y = 1.0f;
                  dir.z = - 1.0f + (invWidthBy2 * float(y) + invWidth);
                  dir = -dir;
                  break;
               case GL_TEXTURE_CUBE_MAP_NEGATIVE_Y:
                  dir.x = - 1.0f + (invWidthBy2 * float(x) + invWidth);
                  dir.y = - 1.0f;
                  dir.z = 1.0f - (invWidthBy2 * float(y) + invWidth);
                  dir = -dir;
                  break;
               case GL_TEXTURE_CUBE_MAP_POSITIVE_Z:
                  dir.x = - 1.0f + (invWidthBy2 * float(x) + invWidth);
                  dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
                  dir.z = 1.0f;
                  break;
               case GL_TEXTURE_CUBE_MAP_NEGATIVE_Z:
                  dir.x = 1.0f - (invWidthBy2 * float(x) + invWidth);
                  dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
                  dir.z = - 1.0f;
                  break;
               default:
                  return;
            }

            // normalize direction
            dir = glm::normalize(dir);

            // scale factor depending on distance from center of the face
            const float fDiffSolid = 4.0f / ((1.0f + fU*fU + fV*fV) *
                                       sqrtf(1.0f + fU*fU + fV*fV));
            fWt += fDiffSolid;

            // calculate coefficients of spherical harmonics for current direction
            sphericalHarmonicsEvaluateDirection(shBuff.data(), order, dir);

            // index of texel in texture
            uint pixOffsetIndex = (x + y * width) * numComponents;
            // get color from texture and map to range [0, 1]
            glm::vec3 clr(
                  float(data[pixOffsetIndex]) / 255,
                  float(data[pixOffsetIndex+1]) / 255,
                  float(data[pixOffsetIndex+2]) / 255
               );

            // scale color and add to previously accumulated coefficients
            sphericalHarmonicsScale(shBuffB.data(), order,
                  shBuff.data(), clr.r * fDiffSolid);
            sphericalHarmonicsAdd(resultR.data(), order,
                  resultR.data(), shBuffB.data());
            sphericalHarmonicsScale(shBuffB.data(), order,
                  shBuff.data(), clr.g * fDiffSolid);
            sphericalHarmonicsAdd(resultG.data(), order,
                  resultG.data(), shBuffB.data());
            sphericalHarmonicsScale(shBuffB.data(), order,
                  shBuff.data(), clr.b * fDiffSolid);
            sphericalHarmonicsAdd(resultB.data(), order,
                  resultB.data(), shBuffB.data());
         }
      }
   }

   // final scale for coefficients
   const float fNormProj = (4.0f * M_PI) / fWt;
   sphericalHarmonicsScale(resultR.data(), order, resultR.data(), fNormProj);
   sphericalHarmonicsScale(resultG.data(), order, resultG.data(), fNormProj);
   sphericalHarmonicsScale(resultB.data(), order, resultB.data(), fNormProj);

   // save result
   for(uint i=0; i < sqOrder; i++)
   {
      output[i].r = resultR[i];
      output[i].g = resultG[i];
      output[i].b = resultB[i];
   }
}
Function that restores lighting value from coefficients of spherical harmonics:
glm::vec3 sphericalHarmonicsFromTexture(glm::vec3 & N, std::vector & coef)
{
   return

      // constant term, lowest frequency //////
      C4 * coef[0] +

      // axis aligned terms ///////////////////
      2.0 * C2 * coef[1] * N.y +
      2.0 * C2 * coef[2] * N.z +
      2.0 * C2 * coef[3] * N.x +

      // band 2 terms /////////////////////////
      2.0 * C1 * coef[4] * N.x * N.y +
      2.0 * C1 * coef[5] * N.y * N.z +
      C3 * coef[6] * N.z * N.z - C5 * coef[6] +
      2.0 * C1 * coef[7] * N.x * N.z +
      C1 * coef[8] * (N.x * N.x - N.y * N.y);
}
Function that calculates coefficients of spherical harmonics for direction:
void sphericalHarmonicsEvaluateDirection(float * result, int order,
                              const glm::vec3 & dir)
{
   // calculate coefficients for first 3 bands of spherical harmonics
   double p_0_0 = 0.282094791773878140;
   double p_1_0 = 0.488602511902919920 * dir.z;
   double p_1_1 = -0.488602511902919920;
   double p_2_0 = 0.946174695757560080 * dir.z * dir.z - 0.315391565252520050;
   double p_2_1 = -1.092548430592079200 * dir.z;
   double p_2_2 = 0.546274215296039590;
   result[0] = p_0_0;
   result[1] = p_1_1 * dir.y;
   result[2] = p_1_0;
   result[3] = p_1_1 * dir.x;
   result[4] = p_2_2 * (dir.x * dir.y + dir.y * dir.x);
   result[5] = p_2_1 * dir.y;
   result[6] = p_2_0;
   result[7] = p_2_1 * dir.x;
   result[8] = p_2_2 * (dir.x * dir.x - dir.y * dir.y);
}
Functions to add and scale coefficients of spherical harmonics:
void sphericalHarmonicsAdd(float * result, int order,
               const float * inputA, const float * inputB)
{
   const int numCoeff = order * order;
   for(int i=0; i < numCoeff; i++)
   {
      result[i] = inputA[i] + inputB[i];
   }
}

void sphericalHarmonicsScale(float * result, int order,
               const float * input, float scale)
{
   const int numCoeff = order * order;
   for(int i=0; i < numCoeff; i++)
   {
      result[i] = input[i] * scale;
   }
}



Sun and Black Cat- Igor Dykhta (igor dykhta email) 2007-2014