dwtobj.cpp

This code demonstrate Fast Wavelet Transform. Executable and can be run on Visual C++ platform

  {
    lstart = ceil(lstart / 2.0f);
    lend = ceil(lend / 2.0f);
    level++;
  }
  return level;

}



//------------------------------------------------------------------------------
// Function         : CDwtObj::calLength()
// Description      : calculate the downsampling length and sampling offset for lowpassed part
//------------------------------------------------------------------------------
// [in] startind, endind, level -- start index and end index, and the desired level in JPEG2000 tiling
// [out]  llenghth , loffset -- downsampling length and sampling offset for lowpassed part

void CDwtObj::calLength(int startind, int endind, int level, int *llength, int *loffset)
{

  int hstart, lstart, hend, lend, i;

  lstart = startind;
  lend   = endind;
  hstart = startind;
  hend   = endind;

  for (i = 0; i < level; i++)
  {
    hstart = floor(lstart / 2.0f);
    hend   = floor(lend   / 2.0f);
    lstart = ceil (lstart / 2.0f);
    lend   = ceil (lend   / 2.0f);
  }

  // low-passed signal length
  *llength = lend - lstart;

  // controls the downsampling, i.e, the even/oddness of the filtered samples to be taken for lowpassed-part
  if (hstart < lstart)
    *loffset = 1;
  else
    *loffset = 0;
}




//------------------------------------------------------------------------------
// Function			: CDwtObj::createIDATexture()
// Description	    : create the indirect addressing table for inverse DWT
//------------------------------------------------------------------------------
// [out]  itex1, -- the texture for holding the indirect addressing table
//		 -- memory allocation is done inside the function, user is responsible to free the memory using delete[]
// [in]  width, height   -- the width and height of the image (not for the indirect addressing table)
// [in]  startx, endx, starty, endy -- the indices of x and y direction
// [out] texwidth, texheight --  the dimensions of the indirect addressing table texture
// 
// Modification history:
//
// 1 June 2005  (By Tien-Tsin Wong) 
// -  Remove the 9D indirect address tables (tex1, tex2, tex3) and combine them into 
//    1 indirect address table. Because they are redundant.
// -  Due to the removal of 9D table, the computation of indirect addresses can be reduced by 9 folds
//
bool CDwtObj::createIDATexture(extmode mode, float **itex1, int width, int height,
                               int &texwidth, int &texheight, int startx, int starty, int endx, int endy)
{
  if (width < 1 || height < 1)
    return false;

  int i, j, k;
  int levellength, halflength, loffset;
  int L = width > height ? width : height;
  int xlevels = calLevels(startx, endx) + 1;
  int ylevels = calLevels(starty, endy) + 1;
  int maxlevels = xlevels > ylevels ? xlevels : ylevels;    //maximum number of decomposition levels for this image

  texwidth = L + 1 + 8 + 8;     // the (8 + 8) is due to boundary extension
  texheight = maxlevels;

  // lookup table holding indices for indirect addressing
  float   (*tex1)[4] = new float[texwidth * texheight][4];

  // create the indirect addressing texture
  for (i = 0; i < xlevels; i++)
  {
    calLength(startx, endx, i, &levellength, &loffset);
    calLength(startx, endx, i + 1, &halflength, &loffset);

    for (j = 0; j < levellength; j++)
    {
      int   eo = j % 2;
      float neighbours[9];

      // put filter selector to tex3
      if (loffset == 1)
        tex1[i * texwidth + j][1] = 1 - eo;
      else
        tex1[i * texwidth + j][1] = eo;
      tex1[i * texwidth + j][2] = 0;


      // Calculate neighbor indirect address (offset = -4), and put it to tex1
      if (loffset == 1)
        neighbours[0] = 0.5 + halflength * (1.0 - eo) + (ext(j - 4, levellength, mode) - 1.0 * eo) / 2.0;
      else
        neighbours[0] = 0.5 + halflength * eo + (ext(j - 4, levellength, mode) - 1.0 * eo) / 2.0;
      tex1[i * texwidth + j][0] = neighbours[0];

      // Compute 9 neighbors' addresses only for the last element (boundary extension)
      if (j == levellength - 1)
      {
        //  these 5 neighbours share same parity (even/oddness) 
        for (k = -4; k <= 4; k += 2)
        {
          if (loffset == 1)
            neighbours[k + 4] = 0.5 + halflength * (1.0 - eo) + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
          else
            neighbours[k + 4] = 0.5 + halflength * eo + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
        }

        eo = 1 - eo;

        //  these 4 neighbours share same parity (even/oddness) 
        for (k = -3; k <= 3; k += 2)
        {
          if (loffset == 1)
            neighbours[k + 4] = 0.5 + halflength * (1.0 - eo) + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
          else
            neighbours[k + 4] = 0.5 + halflength * eo + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
        }

        for (k = 0; k < 9; k++) // compute 9 neighbor's indirect address
          tex1[i * texwidth + j + k][0] = neighbours[k];
      }
    }
  }


  for (i = 0; i < ylevels; i++)
  {
    calLength(starty, endy, i, &levellength, &loffset);
    calLength(starty, endy, i + 1, &halflength, &loffset);

    for (j = 0; j < levellength; j++)
    {
      int   eo = j % 2;
      float neighbours[9];

      // put filter selector to tex1
      if (loffset == 1)
        tex1[(maxlevels - i - 1) * texwidth + (texwidth - 1) - j][1] = 1 - eo;
      else
        tex1[(maxlevels - i - 1) * texwidth + (texwidth - 1) - j][1] = eo;
      tex1[(maxlevels - i - 1) * texwidth + (texwidth - 1) - j][2] = 0;


      // Calculate neighbor indirect address (offset = -4), and put it to tex1            
      if (loffset == 1)
        neighbours[0] = 0.5 + halflength * (1.0 - eo) + (ext(j - 4, levellength, mode) - 1.0 * eo) / 2.0;
      else
        neighbours[0] = 0.5 + halflength * eo + (ext(j - 4, levellength, mode) - 1.0 * eo) / 2.0;
      tex1[(maxlevels - i - 1) * texwidth + (texwidth - 1) - j][0] = neighbours[0];


      // Compute 9 neighbors' addresses only for the last element (boundary extension)
      if (j == levellength - 1)
      {
        //  these 5 neighbour share same parity (even/oddness) 
        for (k = -4; k <= 4; k += 2)
        {
          if (loffset == 1)
            neighbours[k + 4] = 0.5 + halflength * (1.0 - eo) + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
          else
            neighbours[k + 4] = 0.5 + halflength * eo + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
        }

        eo = 1 - eo;

        //  these 4 neighbour share same parity (even/oddness) 
        for (k = -3; k <= 3; k += 2)
        {
          if (loffset == 1)
            neighbours[k + 4] = 0.5 + halflength * (1.0 - eo) + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
          else
            neighbours[k + 4] = 0.5 + halflength * eo + (ext(j + k, levellength, mode) - 1.0 * eo) / 2.0;
        }

        for (k = 0; k < 9; k++)
          tex1[(maxlevels - i - 1) * texwidth + (texwidth - 1) - j - k][0] = neighbours[k];
      }
    }
  }

  *itex1 = (float *) tex1;
  return true;
}








//------------------------------------------------------------------------------
// Function         : CDwtObj::CDwtObj()
// Description      : Constructor
//------------------------------------------------------------------------------

CDwtObj::CDwtObj()
{
  m_imageWidth  = 0;
  m_imageHeight = 0;
  m_iReadBuffer = 0;
  m_rt[0]       = NULL;
  m_rt[1]       = NULL;	
}




//------------------------------------------------------------------------------
// Function         : CDwtObj::~CDwtObj()
// Description      : Destructor
//------------------------------------------------------------------------------

CDwtObj::~CDwtObj()
{
  finalcleanup();
}




//------------------------------------------------------------------------------
// Function         : cgErrorCallback()
// Description      : Cg Error Callback
//------------------------------------------------------------------------------
void cgErrorCallback()
{
  CGerror lastError = cgGetError();
  
  if(lastError)
  {
    printf("%s\n\n", cgGetErrorString(lastError));
    printf("%s\n", cgGetLastListing(cgContext));
    printf("Cg error, exiting...\n");
    
    exit(0);
  }
} 





//------------------------------------------------------------------------------
// Function         : PrintGLerror()
// Description      : Print OpenGL Error Msg
//------------------------------------------------------------------------------

void PrintGLerror( char *msg )
{
 GLenum errCode;
 const GLubyte *errStr;

 if ((errCode = glGetError()) != GL_NO_ERROR) 
 {
    errStr = gluErrorString(errCode);
    fprintf(stderr,"OpenGL ERROR: %s: %s\n", errStr, msg);
 }
}




//------------------------------------------------------------------------------
// Function         : CDwtObj::createLookupTex
// Description      : Creat the  Ext,HL,BasePos Lookup Texture 
//------------------------------------------------------------------------------
// [in] mode --boundary extension mode (symper, per)
// [in] startx, endx, starty, endy -- the indices of x and y direction

bool CDwtObj::createLookupTex(extmode mode, int startx, int starty, int endx, int endy)
{
  //uninitialized imagedimension
  if (m_imageWidth <= 0 || m_imageHeight <= 0)
    return false;

  int   i;
  int   maxrowlevels = calLevels(startx, endx) + 3;
  int   maxcollevels = calLevels(starty, endy) + 3;
  float (*rowlookup)[3];
  float (*collookup)[3];

  // Designed for filter kernel with maximum length 9, so adding 4 at each side for boundary extension, total 8 more entries 
  rowlookup = new float[(m_imageWidth  + 8) * maxrowlevels][3];
  collookup = new float[(m_imageHeight + 8) * maxcollevels][3];

  // wrapping
  //***********************************************************************************
  //*    -4 -3 -2 -1 0 1 2.........m_imagewidth-2, m_imagewidth-1, m_imagewidth-2.... *
  //***********************************************************************************


  // fill up the  lookup table for horizontal decomposition
  for (i = 0 ; i < (m_imageWidth + 8) * maxrowlevels ; i++)
  {
    // determine current row--current level
    int fullrow, halfrow, loffset;
    int level = i / (m_imageWidth + 8);
    int pos   = i % (m_imageWidth + 8);

    calLength(startx, endx, level, &fullrow, &loffset);
    calLength(startx, endx, level + 1, &halfrow, &loffset);

    // The 3 components: boundary extension / High|low / Base Pos
    
    // boundary extension
    rowlookup[i][0] = pos - 4;

    switch (mode)
    {
      case symper:
        if (rowlookup[i][0] < 0)
          rowlookup[i][0] = -rowlookup[i][0];
        if (rowlookup[i][0] > fullrow - 1)
          rowlookup[i][0] = fullrow - 1 - (rowlookup[i][0] - fullrow + 1);
        break;

      case per:
        rowlookup[i][0] = fmod((rowlookup[i][0] + fullrow), fullrow);
        break;
    }

    // High pass or low pass?
    if (rowlookup[i][0] >= halfrow)
      rowlookup[i][1] = 1;      // for Highpass
    else
      rowlookup[i][1] = 0;      // for Lowpass

    // Base Position
    // loffset determine even or odd samples to choose for dowsampling
    if (loffset == 1)
      rowlookup[i][2] = (rowlookup[i][0] - halfrow * rowlookup[i][1]) * 2.0 + 0.5 + 1.0 - rowlookup[i][1];
    else
      rowlookup[i][2] = (rowlookup[i][0] - halfrow * rowlookup[i][1]) * 2.0 + 0.5 + rowlookup[i][1];

    rowlookup[i][0] += 0.5;

  }




  // fill up the lookup table for vertical decomposition
  for (i = 0; i < (m_imageHeight + 8) * maxcollevels; i++)
  {
    //determine current row--current level
    int fullcol, halfcol, loffset;
    int level = i / (m_imageHeight + 8);
    int pos   = i % (m_imageHeight + 8);

    calLength(starty, endy, level, &fullcol, &loffset);
    calLength(starty, endy, level + 1, &halfcol, &loffset);

    // The 3 component:  boundary extension/ High|low / Base Pos