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