dwtobj.cpp

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

// DwtObj.cpp: implementation of the CDwtObj class.
//////////////////////////////////////////////////////////////////////
/*
 * DWT-GPU License Version 1.0
 *
 * Copyright (c) 2004-2006  The Chinese University of Hong Kong
 *
 * Developed by: 
 *    Jianqing Wang, Tien-Tsin Wong, Pheng-Ann Heng, Chi-Sing Leung,
 *    and Liang Wan
 *
 * All rights reserved.
 *
 * References:
 *    Tien-Tsin Wong, Chi-Sing Leung, Pheng-Ann Heng, and Jianqing Wang, 
 *    "Discrete Wavelet Transform on Consumer-Level Graphics Hardware,"
 *    IEEE Transaction on Multimedia, Vol. 9, No. 3, April 2007, pp. 668-673.
 *
 * Modification History:
 * From 0.9a to 0.9b:  (by Tien-Tsin Wong)
 * - Make "dwtsource" and "dwtcolsource" consistent with published paper
 * - Make "invdwtsource" and "invdwtcolsource" consistent with published paper
 * - Simplify the creation of indirect address table for inverse DWT (createIDATexture)
 * - Make the code more neat and tidy
 *
 * From 0.9c to 1.0:  (by Liang Wan)
 * - Provide Frame Buffer Objects (FBO) as an alternative for RenderTexture
 * - Use the latest fragment profile available on the GPU
 *
*/

/*
 * Permission is hereby granted, free of charge, to any person (the
 * "User") obtaining a copy of this software and associated documentation
 * files (the "Software"), to deal in the Software without restriction,
 * including without limitation the rights to use, copy, modify, merge,
 * publish, distribute, and/or sell copies of the Software, and to permit
 * persons to whom the Software is furnished to do so, subject to the
 * following conditions:
 *
 * 1.  If User publishes work based on the Software or its derivative, User
 * agrees to cite the following reference in the publication:
 *
 *     Tien-Tsin Wong, Chi-Sing Leung, Pheng-Ann Heng, and Jianqing Wang, 
 *     "Discrete Wavelet Transform on Consumer-Level Graphics Hardware,"
 *     IEEE Transaction on Multimedia, Vol. 9, No. 3, April 2007, pp. 668-673.
 *     (http://www.cse.cuhk.edu.hk/~ttwong/papers/dwtgpu/dwtgpu.html)
 *
 * 2.  The above copyright notices and this permission notice (which
 * includes the disclaimer below) shall be included in all copies or
 * substantial portions of the Software.
 *
 * 3.  The name of a copyright holder shall not be used to endorse or
 * promote products derived from the Software without specific prior
 * written permission.
 *
 * 4.  The Software is free for both non-commercial and commercial usages.
 * Only the commercial usage requires online registration through the 
 * following webpage:
 *
 * http://www.cse.cuhk.edu.hk/~ttwong/software/dwtgpu/dwtgpu.html 
 * 
 *
 * THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL PART OF THIS
 * LICENSE.  NO USE OF THE SOFTWARE IS AUTHORIZED HEREUNDER EXCEPT UNDER
 * THIS DISCLAIMER.  THE SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS
 * "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
 * BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 * PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.  IN NO
 * EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL
 * INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING
 * FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
 * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.  NO ASSURANCES ARE
 * PROVIDED BY THE COPYRIGHT HOLDERS THAT THE SOFTWARE DOES NOT INFRINGE
 * THE PATENT OR OTHER INTELLECTUAL PROPERTY RIGHTS OF ANY OTHER ENTITY.
 * EACH COPYRIGHT HOLDER DISCLAIMS ANY LIABILITY TO THE USER FOR CLAIMS
 * BROUGHT BY ANY OTHER ENTITY BASED ON INFRINGEMENT OF INTELLECTUAL
 * PROPERTY RIGHTS OR OTHERWISE.  AS A CONDITION TO EXERCISING THE RIGHTS
 * GRANTED HEREUNDER, EACH USER HEREBY ASSUMES SOLE RESPONSIBILITY TO SECURE
 * ANY OTHER INTELLECTUAL PROPERTY RIGHTS NEEDED, IF ANY.  THE SOFTWARE
 * IS NOT FAULT-TOLERANT AND IS NOT INTENDED FOR USE IN MISSION-CRITICAL
 * SYSTEMS, SUCH AS THOSE USED IN THE OPERATION OF NUCLEAR FACILITIES,
 * AIRCRAFT NAVIGATION OR COMMUNICATION SYSTEMS, AIR TRAFFIC CONTROL
 * SYSTEMS, DIRECT LIFE SUPPORT MACHINES, OR WEAPONS SYSTEMS, IN WHICH
 * THE FAILURE OF THE SOFTWARE OR SYSTEM COULD LEAD DIRECTLY TO DEATH,
 * PERSONAL INJURY, OR SEVERE PHYSICAL OR ENVIRONMENTAL DAMAGE ("HIGH
 * RISK ACTIVITIES").  THE COPYRIGHT HOLDERS SPECIFICALLY DISCLAIM ANY
 * EXPRESS OR IMPLIED WARRANTY OF FITNESS FOR HIGH RISK ACTIVITIES.
 * 
 * Credit:
 * This C++ class uses GLEW for OpenGL Extensions.
 */


#include <math.h>
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <time.h>
#include <GL/glew.h>
#ifdef _WIN32
#include <GL/wglew.h>
#else
#include <GL/glxew.h>
#endif

#define  CHECKSUCC(X)  if(!(X)){return false;};
#define  DELTEXTURES(n,t)  if(glIsTexture((*t))){glDeleteTextures(n,t);};




#include "DwtObj.h"

// Wavelet filters kernel


const float decomp_filter[18] = {	
  0.026749, -0.016864, -0.078223, 0.266864, 0.602949, 0.266864, -0.078223, -0.016864, 0.026749,
  0,  0.045636, -0.028772, -0.29564,  0.55754, -0.29564, -0.028772,  0.045636,        0
};


// The Uninterleaved filter
/*
//const float scal=2.0;
float recon_filter[18]={  0, -0.091272f, -0.057544f, 0.591272f, 1.115087f,
                    0.591272f, -0.057544f, -0.091272f, 0 ,
                    0.026749f*scal, 0.016864f*scal, -0.078223f*scal, -0.266864f*scal,
                    0.602949f*scal, -0.266864f*scal, -0.078223f*scal, 0.016864f*scal, 0.026749f*scal
                    };
*/


// Interleaved reconstruction filter
const float recon_filter[18] = {
  0,  0.033728f, -0.057544f, -0.533728f, 1.115087f,-0.533728f, -0.057544f,  0.033728f, 0, 
  0.053498f, -0.091272f, -0.156446f,  0.591272f, 1.205898f, 0.591272f, -0.156446f, -0.091272f, 0.053498f
};




/* Global cgcontext  */
CGcontext   cgContext;


/**********************************************************************/
/*******     Cg program string					   			 **********/
/**********************************************************************/


/*************  Simple Pass Through  **************************/

char fillsource[]="	fragout_float main(vf30 IN,"
"             uniform samplerRECT texture) \n"
"{"
"  fragout_float OUT;\n"

"  OUT.col =f4texRECT(texture, IN.TEX0.xy); \n"
"  OUT.col.a=1.0; \n"
"  return OUT;	\n"
"} \n";


/********************** dwt horizontal decomposition ************************/
char dwtsource[] = 
"fragout_float main(vf30 IN,"
"                   uniform samplerRECT dwt,	  "
"                   uniform samplerRECT filter,   "
"                   uniform samplerRECT lut,      "
"                   uniform float       level) \n"
"{ \n"
"  float  base, offset, level_center=level+0.5;	\n"
"  float2 neighbor;                             \n"
"  float3 lookup, sum=float3(0,0,0);            \n"

"  // Look up alpha (lut.g) and beta (lut.b)                        \n"
"  lookup     = f3texRECT(lut,float2(IN.TEX0.x+4.0, level_center)); \n"
"  offset     = lookup.g*9.0;  // 0 - low-pass, 1 - high-pass       \n"	
"  base       = lookup.b;   \n"
"  neighbor.y = IN.TEX0.y;  \n" 
	
"  for (int i=0; i<9; i++)	\n"
"  {"
"    // Boundary extension is handled by the indirect addressing with lut.r         \n"
"    neighbor.x = f3texRECT(lut,float2(base+i,level_center)).r;	                    \n"
"    // Convolve corresponding filter values with data from level j                 \n"
"    sum += f3texRECT(filter,float2(i+offset+0.5,0.5)).x * f3texRECT(dwt,neighbor); \n"
"  } \n"

"  fragout_float OUT;          \n"
"  OUT.col = float4(sum, 1.0); \n"
"  return OUT;                 \n"
"}" ;


/********************** dwt vertical decomposition ************************/
char dwtcolsource[] =
"fragout_float main(vf30 IN,"
"                   uniform samplerRECT dwt,     "  
"                   uniform samplerRECT filter,  "
"                   uniform samplerRECT lut,     "
"                   uniform float       level) \n"
"{"
"  float  base, offset, level_center=level+0.5; \n"
"  float2 neighbor;                             \n"
"  float3 lookup, sum=float3(0,0,0);            \n"

"  // Look up alpha (lut.g) and beta (lut.b)    \n"
"  lookup     = f3texRECT(lut,float2(IN.TEX0.y+4.0, level_center)); \n"
"  offset     = lookup.g*9.0;  // 0 - low-pass, 1 - high-pass       \n"
"  base       = lookup.b;     \n"
"  neighbor.x = IN.TEX0.x;    \n"

"  for (int i=0; i<9; i++)                   \n"
"  {"
"    // Boundary extension is handled by the indirect addressing with lut.r         \n"
"    neighbor.y = f3texRECT(lut,float2(base+i,level_center)).r;                     \n"
"    // Convolve corresponding filter values with data from level j                 \n"
"    sum += f3texRECT(filter,float2(i+offset+0.5,0.5)).x * f3texRECT(dwt,neighbor); \n"
"  }"

"  fragout_float OUT;          \n"
"  OUT.col = float4(sum, 1.0); \n"
"  return OUT;                 \n"
"} ";



/********************** idwt horizontal reconstruction ************************/
char invdwtsource[] =
"fragout_float main(vf30 IN,"
"                   uniform samplerRECT dwt,     "
"                   uniform samplerRECT filter,  "
"                   uniform samplerRECT lut,     "
"                   uniform float       level) \n"
"{"
"  float  lookup, level_center=level+0.5; \n"
"  float2 st, filter_st;                  \n"
"  float3 sum=float3(0,0,0);              \n"

"  // Look up filter selector (lut.g)                           \n"
"  lookup    = f4texRECT(lut,float2(IN.TEX0.x,level_center)).g; \n"
"  filter_st = float2(lookup*9.0+0.5,0.5);                      \n"
"  st        = float2(0.0,IN.TEX0.y);                           \n"

"  // Lookup indirect address (lut.r) & filter values (filter.x), then convolve \n"
"  for (int i=0 ; i<9 ; i++)                                           \n"
"  {                                                                   \n"
"    st.x         = f4texRECT(lut,float2(IN.TEX0.x+i,level_center)).r; \n"
"    sum         += f3texRECT(filter,filter_st).x * f3texRECT(dwt,st); \n"
"    filter_st.x += 1.0;	                                           \n"
"  }                                                                   \n"

"  fragout_float OUT;	        \n"
"  OUT.col = float4(sum, 1.0);	\n"
"  return OUT;	                \n"
"}";



/********************** idwt vertical reconstruction ************************/
char invdwtcolsource[]=
"fragout_float main(vf30 IN,"
"                   uniform samplerRECT dwt,      "
"                   uniform samplerRECT filter,   "
"                   uniform samplerRECT lut,      "
"                   uniform float       texwidth, "
"                   uniform float       texheight,"
"                   uniform float       level)  \n"
"{"
"  float  lookup, level_center=level+0.5; \n"
"  float2 st, filter_st;                  \n"
"  float3 sum=float3(0,0,0);              \n"

"  // Look up filter selector (lut.g)                                               \n"
"  lookup    = f4texRECT(lut,float2(texwidth-IN.TEX0.y, texheight-level_center)).g; \n"
"  filter_st = float2(0.0+lookup*9.0+0.5,0.5);                                      \n"
"  st        = float2(IN.TEX0.x,0.0);                                               \n"

"  // Lookup indirect address (lut.r) & filter values (filter.x), then convolve            \n"
"  for (int i=0 ; i<9 ; i++) {                                                             \n"
"    st.y         = f4texRECT(lut,float2(texwidth-IN.TEX0.y-i, texheight-level_center)).r; \n"
"    sum         += f3texRECT(filter,filter_st).x * f3texRECT(dwt,st);                     \n"
"    filter_st.x += 1.0;	                                                           \n"
"  }                                                                                       \n"

"  fragout_float OUT;	        \n"
"  OUT.col = float4(sum, 1.0);	\n"
"  return OUT;	                \n"
"}";




//------------------------------------------------------------------------------
// Function			: CDwtObj::ext()
// Description	    : boundary extension function
//------------------------------------------------------------------------------
// [in] index, datalen -- the item index and the whole data length for boundary extension
// [in] mode --boundary extension mode (symper, per)
//  return value --  the index after perform boundary extension


int CDwtObj::ext(int index, int datalen, extmode mode)
{
  int id = index;
  int eo = datalen % 2;

  // symmetric padding
  switch (mode)
  {
    case symper:

      if (id < 0)
        id = -id;

      if (id >= datalen)
        id = datalen - 1 - (id - datalen + 1);

      return id;
      break;


    //periodic padding

    case per:

      id = (id + datalen + eo) % (datalen + eo);
      //id=(id+datalen)%(datalen);

      if (eo == 0)
        return id;

      else
      {
        if (id != datalen)
          return id;
        else
          return datalen - 1;
      }
      break;

    default:
      id = (id + datalen) % datalen;
      return id;
      break;
  }
}



//------------------------------------------------------------------------------
// Function         : CDwtObj::calLevels()
// Description      : calculate maximum levels
//------------------------------------------------------------------------------
// [in] startind, endind -- start index and end index in JPEG2000 tiling
//  return value --  total levels

int CDwtObj::calLevels(int startind, int endind)
{
  int level, lstart, lend;

  lstart = startind;
  lend   = endind;
  level  = 1;

  while (lend - lstart > 1)