www.gusucode.com > 基于dwt的视频数字水印源代码 > 基于dwt的视频数字水印源代码\code\HaYDWT\dwt.cpp

    //Download by http://www.NewXing.com
/*
	Daubechies 9/7 wavelet transform [default]

    Daubechies 9/7 wavelet transform with Lifting

		Written by C. Bloom.
		Modified by Chasan Chouse, 2004-2006.
*/
//---------------------------------------------------------------------------


#pragma hdrstop

#include "dwt.h"
#include "assert.h"

//---------------------------------------------------------------------------
#pragma package(smart_init)


#define WT_ROUND_IN	0
#define WT_ROUND_OUT	0.5

#define F_EXTPAD 4
#define D_EXTPAD 2

// For L97 routines
#define MULT(x,y)	((((x)*(y))+16383)>>15)
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))

#define	C_A	51974
#define C_B	1736
#define C_C	28931
#define C_D	14533


static double *x_alloc = NULL;	//** temp work

void forwardDWT(double *x_in, short N)
{
short i, n, half, j;
double *x,*r,*d;

	x = x_alloc + F_EXTPAD;
	
	// copy array x_in into array x
	for(j=0;j<N;j++)
		*x++ = *x_in++;
	x     -=N; // back to start adress
	x_in  -=N; // back to start adress

	for(i=1;i<=F_EXTPAD;i++) 
	{
		x[-i] = x[i];
		x[(N-1) + i] = x[(N-1) - i];
	}

	half = N>>1;  

	r = x_in; 
	d = x_in + half;
	for (n=half;n--;) 
	{

		*r++ = 0.0378285 * (x[ 4] + x[-4]) - 0.0238495 * (x[ 3] + x[-3]) - 0.110624 * (x[ 2] + x[-2]) + 
					0.377403 * (x[ 1] + x[-1]) + 0.852699 * x[0];

		x++;

		*d++ =  - 0.064539 * (x[ 3] + x[-3]) + 0.040689 * (x[ 2] + x[-2]) + 
				0.418092 * (x[ 1] + x[-1]) - 0.788486 * x[ 0] ;

		x++;
    	}
}

void inverseDWT(double *x, short N)
{
short i, j, n, half;
double *r, *d;

	half = N>>1;

	r = x_alloc + D_EXTPAD;  
	d = x_alloc + D_EXTPAD + half+D_EXTPAD+D_EXTPAD;  

	// copy data to new
	for(j=0; j < half; j++)
		*r++ = *x++;
	r  -=half; // back to start adress

	for(j=0; j < half; j++)
		*d++ = *x++; // above x was in the middle of array. we continue to use remaining elments
	d  -=half; // back to start adress
	x  -=N; // back to start adress


	for(i=1;i<=D_EXTPAD;i++) 
	{
		r[-i] = r[i];
		r[(half-1)+i] = r[half - i];
		d[-i] = d[i-1];
		d[(half-1)+i] = d[(half-1) - i];
	}
 
	for (n = half; n--;) 
	{

		*x++ = 0.788486 * r[0] - 0.0406894 * ( r[1] + r[-1] )
		     - 0.023849 * (d[1] + d[-2]) + 0.377403 * (d[0] + d[-1]);

		*x++ = 0.418092 * (r[1] + r[0]) - 0.0645389 * (r[2] + r[-1])
		     - 0.037829 * (d[2] + d[-2]) + 0.110624 * (d[1] + d[-1]) - 0.852699 * d[0];

		d++; 
		r++;
	}
}

void waveletTransform2D(UCHAR **raw, double **wave, long width, long height, int levels, bool forward)
{
	long i, j, x, y, w, h, L, c;
	double *buffer;

    /* Check the dimensions for compatability. */

    if (width%(1 << levels) || height%(1 << levels))
	{
              	MessageDlg("width and height must be divisible by 2^levels", mtError, TMsgDlgButtons() << mbOK, 0);
                return;
	}

	if ( (x_alloc = (double *)malloc(sizeof(double)*(width+height+F_EXTPAD+F_EXTPAD))) == NULL )  
	{
              	MessageDlg("Malloc failed", mtError, TMsgDlgButtons() << mbOK, 0);
                return;
	}

    /* Allocate a work array (for transposing columns) */

   	if ( (buffer = (double *)calloc(height,sizeof(double))) == NULL )
	{
              	MessageDlg("Malloc failed", mtError, TMsgDlgButtons() << mbOK, 0);
                return;
	}

    /* Compute the wavelet transform. */
  
	if ( forward ) /* forward transform. */
	{ 

		/** raw -> wave **/
		//  Copy raw ---> wave matrix
		for(i=0;i<height;i++)
		{
			for(j=0;j<width;j++)
			{
				wave[i][j] = raw[i][j] + WT_ROUND_IN;
			}
		}

		for (L = 0; L < levels; L++) 
		{
			w = width >> L;
			h = height >> L;
      
			/* Rows. */
	
			for (y = 0; y < h; y++)
				forwardDWT(wave[y], (short)w);
    
			/* Columns. */


			for (x = 0; x < w; x++) 
			{
				// place one column from wave matrix to "buffer" array
				for(y=0; y < h; y++)
					buffer[y] = wave[y][x];

				forwardDWT(buffer, (short)h);

				// write back the buffer array to wave matrix
				for(y=0; y < h; y++)
					wave[y][x] = buffer[y];
			} // for
		} // for L
    } // if forward
	else  // inverse
	{
		/* backward transform. */
		for (L = levels-1; L >= 0; L--) {
			w = width >> L;
			h = height >> L;

			/* Columns. */
	
			for (x = 0; x < w; x++) 
			{
				for(y=0; y < h; y++)
					buffer[y] = wave[y][x];

				inverseDWT(buffer, (short)h);

				for(y=0; y < h; y++)
					wave[y][x] = buffer[y];
			}

			/* Rows. */
	
			for (y = 0; y < h; y++)
				inverseDWT(wave[y], (short)w);
    	}

		//  Copy wave ---> raw matrix
		for(i=0;i<height;i++)
		{
			for(j=0;j<width;j++)
			{
				c = (long)(wave[i][j] + WT_ROUND_OUT);
				if ( c > 255)
					raw[i][j] = 255;
				else if ( c < 0)
					raw[i][j] = 0;
				else
					raw[i][j] = (UCHAR)c;
			}
		}

	}

	free(x_alloc);
	free(buffer);
}

// This routine for YUV image DWT transform
void waveletTransform2DDouble(double **raw, double **wave, long width, long height, int levels, bool forward)
{
	long i, j, x, y, w, h, L, c;
	double *buffer;

    /* Check the dimensions for compatability. */

    if (width%(1 << levels) || height%(1 << levels))
	{
              	MessageDlg("width and height must be divisible by 2^levels", mtError, TMsgDlgButtons() << mbOK, 0);
                return;
	}

	if ( (x_alloc = (double *)malloc(sizeof(double)*(width+height+F_EXTPAD+F_EXTPAD))) == NULL )
	{
              	MessageDlg("Malloc failed", mtError, TMsgDlgButtons() << mbOK, 0);
                return;
	}

    /* Allocate a work array (for transposing columns) */

   	if ( (buffer = (double *)calloc(height,sizeof(double))) == NULL )
	{
              	MessageDlg("Malloc failed", mtError, TMsgDlgButtons() << mbOK, 0);
                return;
	}

    /* Compute the wavelet transform. */
  
	if ( forward ) /* forward transform. */
	{ 

		/** raw -> wave **/
		//  Copy raw ---> wave matrix
		for(i=0;i<height;i++)
		{
			for(j=0;j<width;j++)
			{
				wave[i][j] = raw[i][j] + WT_ROUND_IN;
			}
		}

		for (L = 0; L < levels; L++) 
		{
			w = width >> L;
			h = height >> L;
      
			/* Rows. */
	
			for (y = 0; y < h; y++)
				forwardDWT(wave[y], (short)w);

			/* Columns. */


			for (x = 0; x < w; x++) 
			{
				// place one column from wave matrix to "buffer" array
				for(y=0; y < h; y++)
					buffer[y] = wave[y][x];

				forwardDWT(buffer, (short)h);

				// write back the buffer array to wave matrix
				for(y=0; y < h; y++)
					wave[y][x] = buffer[y];
			} // for
		} // for L
    } // if forward
	else  // inverse
	{
		/* backward transform. */
		for (L = levels-1; L >= 0; L--) {
			w = width >> L;
			h = height >> L;

			/* Columns. */
	
			for (x = 0; x < w; x++) 
			{
				for(y=0; y < h; y++)
					buffer[y] = wave[y][x];

				inverseDWT(buffer, (short)h);

				for(y=0; y < h; y++)
					wave[y][x] = buffer[y];
			}

			/* Rows. */
	
			for (y = 0; y < h; y++)
				inverseDWT(wave[y], (short)w);
    	}

		//  Copy wave ---> raw matrix
		for(i=0;i<height;i++)
		{
			for(j=0;j<width;j++)
			{
				raw[i][j] = (long)(wave[i][j] + WT_ROUND_OUT);
			}
		}

	}

	free(x_alloc);
	free(buffer);
}


/***********************************************

  Daubechies 9/7 wavelet transform with Lifting

************************************************/


static void do_tdec_line(int * to,int *from,int len)
{
int x,*ptr,*low,*high,half;

	half = len>>1;

//	assert(len & 1 == 0);
//	assert(half >= 2);

	low  = to; high = to + half; ptr = from;

		*high = ptr[1] - MULT(C_A,( ptr[0] + ptr[2] ));
		*low  = ptr[0] - MULT(C_B,( high[0] + high[0]));
		low++; high++; ptr += 2;

	for(x=half-2;x--;) {
		*high = ptr[1] - MULT(C_A,( ptr[0] + ptr[2] ));
		*low  = ptr[0] - MULT(C_B,( high[0] + high[-1]));
		low++; high++; ptr += 2;
	}

		*high = ptr[1] - MULT(C_A,( ptr[0] + ptr[0] ));
		*low  = ptr[0] - MULT(C_B,( high[0] + high[-1]));

	low  = to; high = to + half;

		*high += MULT(C_C,( low[0] + low[1] ));
		*low  += MULT(C_D,( high[0] + high[0]));
		low++; high++;

	for(x=half-2;x--;) {
		*high += MULT(C_C,( low[0] + low[1] ));
		*low  += MULT(C_D,( high[0] + high[-1]));
		low++; high++;
	}

		*high += MULT(C_C,( low[0] + low[0] ));
		*low  += MULT(C_D,( high[0] + high[-1]));

}

static void un_tdec_line(int *to,int *from,int len)
{
int x,*ptr,*low,*high,half;

	half = len>>1;

	low  = from; high = from + half;
	for(x=0;x<half;x++) {
		if ( x == 0 )		*low  -= MULT(C_D,( high[0] + high[0]));
		else				*low  -= MULT(C_D,( high[0] + high[-1]));
		low++; high++;
	}

	low  = from; high = from + half;
	for(x=0;x<half;x++) {
		if ( x == half-1 )	*high -= MULT(C_C,( low[0] + low[0] ));
		else				*high -= MULT(C_C,( low[0] + low[1] ));
		low++; high++;
	}

	low  = from; high = from + half; ptr = to;
	for(x=0;x<half;x++) {
		if ( x == 0 )		ptr[0] = *low  + MULT(C_B,( high[0] + high[0]));
		else				ptr[0] = *low  + MULT(C_B,( high[0] + high[-1]));
		low++; high++; ptr += 2;
	}

	low  = from; high = from + half; ptr = to;
	for(x=0;x<half;x++) {
		if ( x == half-1 )	ptr[1] = *high + MULT(C_A,( ptr[0] + ptr[0] ));
		else				ptr[1] = *high + MULT(C_A,( ptr[0] + ptr[2] ));
		low++; high++; ptr += 2;
	}

}

void l97_2D(int **rows, int width, int height, int levels,bool inverse)
{
int x, y, w, h, l;
int *buffer,*tempbuf,*temprow;


    /* Check the dimensions for compatability. */

    if (width%(1 << (levels+1)) || height%(1 << (levels+1)))
	{
       	MessageDlg("width and height must be divisible by 2^levels", mtError, TMsgDlgButtons() << mbOK, 0);
        return;
	}

   	if ( (buffer = (int *)malloc( sizeof(int) * (height+max(width,height)+height))) == NULL )
    {
       	MessageDlg("Malloc failed", mtError, TMsgDlgButtons() << mbOK, 0);
        return;
    }

	temprow = buffer+height;
	tempbuf = buffer+height+height;

	if ( !inverse ) {  // forward transform

		for (l = 0; l < levels; l++) {
			w = width >> l;
			h = height >> l;

			/* Rows */
	
			do_tdec_line(temprow,rows[h-1],w);
			for (y = h-2; y >=0; y--) {
				do_tdec_line(rows[y+1],rows[y],w);
			}
	
			/* Columns */
	
			for (x = 0; x < w; x++) {
					for (y = 1; y < h; y++) buffer[y-1] = rows[y][x];
					buffer[h-1] = temprow[x];
				do_tdec_line(tempbuf,buffer,h);
					for (y = 0; y < h; y++) rows[y][x] = tempbuf[y];
			}
		}

    } else {

		for (l = levels-1; l >= 0; l--) { /** backwards in scale **/
			w = width >> l;
			h = height >> l;

			/* Columns */
	
			for (x = 0; x < w; x++) {
					for (y = 0; y < h; y++) buffer[y] = rows[y][x];
				un_tdec_line(tempbuf,buffer,h);
					for (y = 0; y < h-1; y++) rows[y+1][x] = tempbuf[y];
					temprow[x] = tempbuf[h-1];
			}

			/* Rows */
			for (y = 0; y < h-1; y++) {
				un_tdec_line(rows[y],rows[y+1],w);
			} 
			un_tdec_line(rows[h-1],temprow,w);
		}
	}

	free(buffer);
}


void l97Quad(int *band,int w,int h,int fullw,bool inverse)
{
int x, y;
int *buffer,*tempbuf,*bptr,*temprow;


   	if ( (buffer = (int *)malloc( sizeof(int) * h+h+max(w,h))) == NULL )
    {
       	MessageDlg("Malloc failed", mtError, TMsgDlgButtons() << mbOK, 0);
        return;
    }

	temprow = buffer+h;
	tempbuf = buffer+h+h;
  
	if ( !inverse ) { /* forward transform. */
	
		bptr = band + (h-1)*fullw;
		do_tdec_line(temprow,bptr,w);
		for (y = (h-1); y--;) {
			bptr -= fullw;
			do_tdec_line(bptr+fullw,bptr,w);
		}
    
		for (x = 0; x < w; x++) {
			bptr = band + x + fullw;
			for (y = 0; y < (h-1); y++) { buffer[y] = *bptr; bptr += fullw; }
			buffer[h-1] = temprow[x];
			do_tdec_line(tempbuf,buffer,h);
			bptr = band + x;
			for (y = 0; y < h; y++) { *bptr = tempbuf[y]; bptr += fullw; }
		}

    } else {

		for (x = 0; x < w; x++) {
			bptr = band + x;
			for (y = 0; y < h; y++) { buffer[y] = *bptr; bptr += fullw; }
			un_tdec_line(tempbuf,buffer,h);
			bptr = band + x + fullw;
			for (y = 0; y < h-1; y++) { *bptr = tempbuf[y]; bptr += fullw; }
			temprow[x] = tempbuf[h-1];
		}

		bptr = band;
		for (y = (h-1); y--; ){
			un_tdec_line(bptr,bptr+fullw,w);
			bptr += fullw;
		}
		un_tdec_line(bptr,temprow,w);

	}

	free(buffer);
}