文章出处：http://outdoorsports.blog.hexun.com/18282012_d.html
int WaveletFilter(unsigned char *pData, int heightOld, int widthOld){int i,j; // 循环变量int height,width; // 扩展后图像的高和宽int scale; // 分解尺度int m; // 整除数, 2^scalescale = 3;//表示分解3层m = (int)pow(2,scale);// 计算扩展后图像的尺寸，使得图像的高和宽能被m整除if(0 == heightOld%m){height = heightOld;}else{height = (m - heightOld%m) + heightOld;}if(0 == widthOld%m){width = widthOld;}else{width = (m - widthOld%m) + widthOld;}double **tempData; // 临时变量，存放图像数据tempData = new double *[height];for(i=0; i<height; i++){tempData[i] = new double[width];}// 初始化,扩展部分采用周期延拓for(i=0; i<height; i++){for(j=0; j<width; j++){if(i >= heightOld){if(j >= widthOld)tempData[i][j] = (double)pData[(i-heightOld)*widthOld + j-widthOld];elsetempData[i][j] = (double)pData[(i-heightOld)*widthOld + j];}else{if(j >= widthOld)tempData[i][j] = (double)pData[i*widthOld + j-widthOld];elsetempData[i][j] = (double)pData[i*widthOld + j];} }}// 取得滤波器系数char* wavelet = "db4"; // 取何种小波int filterlength; // 滤波器长度//形成小波系数if(1 != Wfilters(wavelet, &filterlength)){return 0;}// 分解int heightNew = height;int widthNew = width;for(i=1; i<=scale; i++){DWT2(tempData,heightNew,widthNew,filterlength);// 低频子图像的高和宽heightNew /= 2;widthNew /= 2;}// int starttime = highGetTime();// 滤波Denoise(tempData,height,width,scale);// int endtime = highGetTime();// m_filtertime = endtime - starttime;// 重构for(i=scale; i>=1; i--){//重构图像的高和宽heightNew += heightNew;widthNew += widthNew;
IDWT2(tempData,heightNew,widthNew,filterlength);}// int endtime = highGetTime();// int totaltime = endtime - starttime;// 输出滤波结果for(i=0; i<heightOld; i++){for(j=0; j<widthOld; j++){if(tempData[i][j] > 255){pData[i*widthOld + j] = (byte)255;}else if(tempData[i][j] < 0){pData[i*widthOld + j] = (byte)0;}else{pData[i*widthOld + j] = (byte)tempData[i][j];}}}// 删除分配的内存for(i=0; i<height; i++){delete []tempData[i];}delete []tempData;return 1;}
// int Wfilters(char* wavelet, int* length) 函数
//
//
int Wfilters(char* wavelet, int* length){char* waveletName; // 小波名称，按照Matlab中waveletName = "db4"; // db4小波，消失矩4if(0 == strcmp(wavelet,waveletName)){double db[8] = {0.23037781330886, 0.71484657055254, 0.63088076792959, -0.02798376941698,-0.18703481171888, 0.03084138183599, 0.03288301166698, -0.01059740178500};*length = 8; // 滤波器长度SetFilters(db,*length);}return 1;}
// int DWT2 函数
int DWT2(double **data, int height, int width, int filterlength){int i,j; // 循环变量double *temp = new double[height]; // 临时存放每一列的数据// 行分解for(i=0; i<height; i++)DWT(data[i],width,filterlength);// 列分解for(j=0; j<width; j++){// 列转置for(i=0; i<height; i++){temp[i] = data[i][j];}DWT(temp,height,filterlength);// 反转置for(i=0; i<height; i++){data[i][j] = temp[i];}}delete []temp;return 1;}
//
//
// int DWT(double *data, int length, int filterlength)
/
int DWT(double *data, int length, int filterlength){int i,k; // 循环变量int n = filterlength/2; // 滤波器中心位置double sum1,sum2; // 临时变量double *temp = new double[length]; // 临时变量，存放变换结果for(i=0; i<length; i+=2) //下采样{// 初始化sum1 = 0.0;sum2 = 0.0;// 计算滤波器卷积结果，滤波器中心位置和i重合for(k=0; k<filterlength; k++){if(i-n+k < 0) // 边界处理，采用周期延拓{sum1 += data[i-n+k + length]*m_lo_D[k]; // 低通sum2 += data[i-n+k + length]*m_hi_D[k]; // 高通}else if(i-n+k >=length) // 边界处理，采用周期延拓{sum1 += data[i-n+k - length]*m_lo_D[k]; // 低通sum2 += data[i-n+k - length]*m_hi_D[k]; // 高通}else{sum1 += data[i-n+k]*m_lo_D[k]; // 低通sum2 += data[i-n+k]*m_hi_D[k]; // 高通}}// 结果temp[i/2] = sum1; // 概貌信息temp[length/2 + i/2] = sum2; // 细节信息}// 更新for(i=0; i<length; i++)data[i] = temp[i];delete []temp;return 1;}
//
// int IDWT2 函数
int IDWT2(double **coeff, int height, int width, int filterlength){int i,j; // 循环变量double *temp = new double[height]; // 临时存放每一列的数据// 列重构for(j=0; j<width; j++){// 列转置for(i=0; i<height; i++){temp[i] = coeff[i][j];}IDWT(temp,height,filterlength);// 反转置for(i=0; i<height; i++){coeff[i][j] = temp[i];}}// 行重构for(i=0; i<height; i++)IDWT(coeff[i],width,filterlength);return 1;}
///
//int IDWT(double *coeff, int length, int filterlength) 函数
//
/// int IDWT(double *coeff, int length, int filterlength){int i,k; // 循环变量int n = filterlength/2; // 滤波器中心位置double sum1,sum2; // 临时变量double *temp = new double[length]; // 临时变量，存放变换结果for(i=0; i<length; i++){// 初始化sum1 = 0.0;sum2 = 0.0;// 计算滤波器卷积结果for(k=0; k<filterlength; k++){if(i-n+k+1 < 0) // 边界处理，采用周期延拓{if(0 == (i-n+k+1 + length)%2) // 上采样{sum1 += coeff[(i-n+k+1 + length)/2]*m_lo_R[k]; // 低通sum2 += coeff[length/2 + (i-n+k+1 + length)/2]*m_hi_R[k]; // 高通}}else if(i-n+k+1 >=length) // 边界处理，采用周期延拓{if(0 == (i-n+k+1 - length)%2) // 上采样{sum1 += coeff[(i-n+k+1 - length)/2]*m_lo_R[k]; // 低通sum2 += coeff[length/2 + (i-n+k+1 - length)/2]*m_hi_R[k]; // 高通}}else{if(0 == (i-n+k+1)%2) // 上采样{sum1 += coeff[(i-n+k+1)/2]*m_lo_R[k]; // 低通sum2 += coeff[length/2 + (i-n+k+1)/2]*m_hi_R[k]; // 高通}}}// 结果temp[i] = sum1 + sum2;}// 更新for(i=0; i<length; i++)coeff[i] = temp[i];delete []temp;return 1;}
///
// int Denoise(double **coeff, int height, int width, int scale) 函数
//
///
int Denoise(double **coeff, int height, int width, int scale){int i,j,k; // 循环变量int r,c; // 行和列double delta; // 噪声系数的标准差for(k=1; k<=scale; k++){if(1 == k){// 利用对角子带（HH）求噪声标准差r = (int)(height/pow(2,k)); // 子带的行c = (int)(width/pow(2,k)); // 子带的列double *temp = new double[r * c]; // 临时存放子带数据，求取中值for(i=0; i<r; i++)for(j=0; j<c; j++)temp[i*c + j] = ABS(coeff[r+i][c+j]);delta = prhap(temp,r*c) / 0.6745;delete []temp;}// 域值处理r = (int)(height/pow(2,k-1)); // LL,LH,HL,HH四个字图像组成的图像的行c = (int)(width/pow(2,k-1)); // LL,LH,HL,HH四个字图像组成的图像的列Thresholding(coeff, r/2, c/2, r, c, delta); //HHThresholding(coeff, r/2, 0, r, c/2, delta); //HLThresholding(coeff, 0, c/2, r/2, c, delta); //LH}return 1;}
//
//
//int Thresholding(double **coeff, int rTop, int cLeft, int rBottom,int cRight, double delta)
/
int Thresholding(double **coeff, int rTop, int cLeft, int rBottom,int cRight, double delta){int i,j,m,n;int nr = 3; // 窗口半径double sum1,sum2;int num;double mean,var,T;for(i=rTop; i<rBottom; i++){for(j=cLeft; j<cRight; j++){sum1 = 0.0; sum2 = 0.0; num = 0;for(m=-nr; m<=nr; m++){for(n=-nr; n<=nr; n++){if((i+m >= rTop) && (i+m < rBottom) && (j+n >= cLeft) && (j+n < cRight)){num += 1;sum1 += coeff[i+m][j+n];sum2 += coeff[i+m][j+n]*coeff[i+m][j+n];}}}mean = sum1/num;var = sum2/num - mean*mean;if(var <= delta*delta)T = 1000; // 取很大的值，使处理后的系数为0elseT = sqrt(2) * delta*delta/sqrt(var-delta*delta);// 软域值处理mean = ABS(coeff[i][j]) - T;if(mean <= 0)coeff[i][j] = 0;elsecoeff[i][j] = coeff[i][j]*mean / (mean+T);}}return 1;}
// //
// double prhap(double *p,int n) 堆排序,比冒泡快得多
double prhap(double *p,int n){int i,mm;double t;mm=n/2;for(i=mm-1; i>=0; i--)rsift(p,i,n-1);for(i=n-1; i>=1; i--){t=p[0];p[0]=p[i];p[i]=t;rsift(p,0,i-1);}return p[n/2];}
