图像处理之基础---不同小波基的小波变换（卷积）
#include <stdio.h>
#include <stdlib.h>
#include <cv.h>
#include <highgui.h>
#include <math.h>
#define WIN_WIDTH 1//1~10,选择小波基
double *Ld = new double[2*WIN_WIDTH]; //分解尺度函数
double *Hd = new double[2*WIN_WIDTH]; //分解母函数
double *Lr = new double[2*WIN_WIDTH]; //重建尺度函数
double *Hr = new double[2*WIN_WIDTH]; //重建母函数
int m_preoffset;
int m_aftoffset;
double m_GrayMax=0;
double m_GrayMin=255;
BOOL Wavelet(double *image,int FWidth,int nLevel, int width, int hight);
BOOL DisWavelet(double *image,int FWidth,int nLevel, int width, int hight);
void main()
{
int FWidth =2*WIN_WIDTH;
// 小波变换层数
int nLayer = 1;
// 输入彩色图像
IplImage *pSrc = cvLoadImage("1.bmp", CV_LOAD_IMAGE_GRAYSCALE);
if(!pSrc)
exit(1);
//cvSaveImage("222.bmp",pSrc);
// 计算小波图象大小
int height = ((pSrc->height)>>nLayer)<<nLayer;
int width = ((pSrc->width)>>nLayer)<<nLayer;
// 创建小波图象
cvNamedWindow("1",1);
cvShowImage("1",pSrc);
//cvWaitKey(0);
IplImage *pWavelet = cvCreateImage(cvGetSize(pSrc), IPL_DEPTH_32F, 1);
IplImage *image = cvCreateImage(cvGetSize(pSrc), IPL_DEPTH_8U, 1);
double *pMatrix = new double[height*width];
int temp = 0;
int temp1 = 0;
for(int i = 0; i < height; i++)
{
for(int j = 0; j < width; j++)
{
pMatrix[i*width + j] = 0.0;
temp1 = (unsigned char)(*(pSrc->imageData+i*pSrc->widthStep+j));
pMatrix[i*width + j] = (double)(temp1);
}
}
{
Wavelet(pMatrix, FWidth,1, width, height);
for(int i = 0; i < height; i++)
{
for(int j = 0; j < width; j++)
{
temp = (int)(pMatrix[i*width + j]);
if(temp < 0)
temp = 0;//abs(temp);
if(temp > 255)
temp = 255;
*(image->imageData+i*image->widthStep+j) = (unsigned char)temp;
}
}
cvNamedWindow("12",1);
cvShowImage("12",image);
cvWaitKey(0);
DisWavelet(pMatrix, FWidth,1, width, height);
}
for(int i = 0; i < height; i++)
{
for(int j = 0; j < width; j++)
{
temp = (int)(pMatrix[i*width + j]);
if(temp < 0)
temp = 0;//abs(temp);
if(temp > 255)
temp = 255;
*(image->imageData+i*image->widthStep+j) = (unsigned char)temp;
}
}
delete pMatrix;
cvNamedWindow("123",1);
cvShowImage("123",image);
cvWaitKey(0);
cvReleaseImage(&pWavelet);
return ;
}
/*************************************************************************
*
* 函数名称：
* Convolution()
*
* 参数:
* double * LF HF - 指向小波的指针，是常量
* FR - 小波窗的宽度
* double * f - 指向时域值的指针和返回的小波变换频域的指针
* fr －原图象宽
* 说明:
* 该函数用来实现卷积运算。
*
************************************************************************/
void Convolution(double *LF,double *HF,int FR, double *f, int fr)
{
int i,j,m;// 循环变量
double *X;
X = new double[fr]; // 分配运算所需的数组
// 卷积运算
for(i=0;i<fr/2;i++)
{
X[i]=0;
X[i+fr/2]=0;
for(j=0;j<FR;j++)
{
m=(2*i+j+fr-m_preoffset)%fr;
X[i]+=f[m]*LF[j];
X[i+fr/2]+=f[m]*HF[j];
}
}
//运算结果反传给f。
for(i= 0; i <fr; i++)
{
f[i]=X[i];
}
delete X;// 释放内存
}
/*************************************************************************
*
* 函数名称：
* DisConvolution()
*
* 参数:
* double * F - 指向小波的指针，是常量
* FR - 小波窗的宽度
* double * f - 指向时域值的指针和返回的小波变换频域的指针
* fr －原图象宽
* 说明:
* 该函数用来实现解卷积运算。
*
************************************************************************/
void DisConvolution(double *LF,double *HF,int FR, double *f0,double *f1, int fr)
{
int i,j;// 循环变量
double *X,*Y;
// 分配运算所需的数组
X = new double[fr];
Y = new double[fr];
// 解卷积运算
for(i=0;i<fr;i++)
{
X[i]=0;
Y[i]=0;
for(j=0;j<FR;j++)
{
X[i]+=f0[(i+j)%fr]*LF[j];
Y[i]+=f1[(i+j)%fr]*HF[j];
}
}
//运算结果反传给f0。
for(i= 0; i <fr; i++)
{
j=(i+fr-m_aftoffset)%fr;//循环移位
f0[i]=X[j]+Y[j];
}
delete X,Y; // 释放内存
}
/*************************************************************************
*
* 函数名称：
* DIBWavelet()
* 参数:
* double *LF - 使用的小波尺度函数，是常量
* double *HF - 使用的小波母函数，是常量
* int FWidth - 小波窗的宽度
* int nLevel - 小波分解的层数
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用来对图像进行小波变换分解。
************************************************************************/
BOOL Wavelet( double *f,int FWidth,int nLevel, int width, int hight)
{
//生成小波分解和重建的尺度函数和母函数
switch(FWidth)
{
case 2:
Ld[0]=1/sqrt(2.0f); Ld[1]=1/sqrt(2.0f);
m_preoffset=0;
m_aftoffset=1;
break;
case 4:
Ld[0]=0.4829629131445341; Ld[1]=0.8365163037378077; Ld[2]=0.2241438680420134;
Ld[3]=-0.1294095225512603;
m_preoffset=1;
m_aftoffset=2;
break;
case 6:
Ld[0]=0.3326705529500825; Ld[1]=0.8068915093110924; Ld[2]=0.4598775021184914;
Ld[3]=-0.1350110200102546; Ld[4]=-0.0854412738820267; Ld[5]=0.035226218857095;
m_preoffset=2;
m_aftoffset=3;
break;
case 8:
Ld[0]=-0.107148901418/sqrt(2.0); Ld[1]=-0.041910965126/sqrt(2.0); Ld[2]=0.703739068656/sqrt(2.0);
Ld[3]=1.136658243409/sqrt(2.0); Ld[4]=0.421234534204/sqrt(2.0); Ld[5]=-0.140317624179/sqrt(2.0);
Ld[6]=-0.017824701442/sqrt(2.0); Ld[7]=0.045570345896/sqrt(2.0);
m_preoffset=3;
m_aftoffset=4;
break;
case 10:
Ld[0]=0.038654795955/sqrt(2.0); Ld[1]=0.041746864422/sqrt(2.0); Ld[2]=-0.055344186117/sqrt(2.0);
Ld[3]=0.281990696854/sqrt(2.0); Ld[4]=1.023052966894/sqrt(2.0); Ld[5]=0.896581648380/sqrt(2.0);
Ld[6]=0.023478923136/sqrt(2.0); Ld[7]=-0.247951362613/sqrt(2.0); Ld[8]=-0.029842499869/sqrt(2.0);
Ld[9]=0.027632152958/sqrt(2.0);
m_preoffset=4;
m_aftoffset=5;
break;
case 12:
Ld[0]=0.021784700327/sqrt(2.0); Ld[1]=0.004936612372/sqrt(2.0); Ld[2]=-0.166863215412/sqrt(2.0);
Ld[3]=-0.068323121587/sqrt(2.0); Ld[4]=0.0694457972958/sqrt(2.0); Ld[5]=1.113892783926/sqrt(2.0);
Ld[6]=0.477904371333/sqrt(2.0); Ld[7]=-0.102724969862/sqrt(2.0); Ld[8]=-0.029783751299/sqrt(2.0);
Ld[9]=0.063250562660/sqrt(2.0); Ld[10]=0.002499922093/sqrt(2.0); Ld[11]=-0.011031867509/sqrt(2.0);
m_preoffset=5;
m_aftoffset=6;
break;
case 14:
Ld[0]=0.003792658534/sqrt(2.0); Ld[1]=-0.001481225915/sqrt(2.0); Ld[2]=-0.017870431651/sqrt(2.0);
Ld[3]=0.043155452582/sqrt(2.0); Ld[4]=0.096014767936/sqrt(2.0); Ld[5]=-0.070078291222/sqrt(2.0);
Ld[6]=0.024665659489/sqrt(2.0); Ld[7]=0.758162601964/sqrt(2.0); Ld[8]=1.085782709814/sqrt(2.0);
Ld[9]=0.408183939725/sqrt(2.0); Ld[10]=-0.198056706807/sqrt(2.0); Ld[11]=-0.152463872896/sqrt(2.0);
Ld[12]=0.005671342686/sqrt(2.0); Ld[13]=0.014521394762/sqrt(2.0);
m_preoffset=7;
m_aftoffset=6;
break;
case 16:
Ld[0]=0.002672793393/sqrt(2.0); Ld[1]=-0.000428394300/sqrt(2.0); Ld[2]=-0.021145686528/sqrt(2.0);
Ld[3]=0.005386388754/sqrt(2.0); Ld[4]=0.069490465911/sqrt(2.0); Ld[5]=-0.038493521263/sqrt(2.0);
Ld[6]=-0.073462508761/sqrt(2.0); Ld[7]=0.515398670374/sqrt(2.0); Ld[8]=1.099106630537/sqrt(2.0);
Ld[9]=0.680745347190/sqrt(2.0); Ld[10]=-0.086653615406/sqrt(2.0); Ld[11]=-0.202648655286/sqrt(2.0);
Ld[12]=0.010758611751/sqrt(2.0); Ld[13]=0.044823623042/sqrt(2.0); Ld[14]=-0.000766690896/sqrt(2.0);
Ld[15]=-0.004783458512/sqrt(2.0);
m_preoffset=8;
m_aftoffset=7;
break;
case 18:
Ld[0]=0.001512487309/sqrt(2.0); Ld[1]=-0.000669141509/sqrt(2.0); Ld[2]=-0.014515578553/sqrt(2.0);
Ld[3]=0.012528896242/sqrt(2.0); Ld[4]=0.087791251554/sqrt(2.0); Ld[5]=-0.025786445930/sqrt(2.0);
Ld[6]=-0.270893783503/sqrt(2.0); Ld[7]=0.049882830959/sqrt(2.0); Ld[8]=0.873048407349/sqrt(2.0);
Ld[9]=1.015259790832/sqrt(2.0); Ld[10]=0.337658923602/sqrt(2.0); Ld[11]=-0.077172161097/sqrt(2.0);
Ld[12]=0.000825140929/sqrt(2.0); Ld[13]=0.042744433602/sqrt(2.0); Ld[14]=-0.016303351226/sqrt(2.0);
Ld[15]=-0.018769396836/sqrt(2.0); Ld[16]=0.000876502539/sqrt(2.0); Ld[17]=0.001981193736/sqrt(2.0);
m_preoffset=8;
m_aftoffset=9;
break;
case 20:
Ld[0]=0.001089170447/sqrt(2.0); Ld[1]=0.000135245020/sqrt(2.0); Ld[2]=-0.012220642630/sqrt(2.0);
Ld[3]=-0.002072363923/sqrt(2.0); Ld[4]=0.064950924579/sqrt(2.0); Ld[5]=0.016418869426/sqrt(2.0);
Ld[6]=-0.225558972234/sqrt(2.0); Ld[7]=-0.100240215031/sqrt(2.0); Ld[8]=0.667071338154/sqrt(2.0);
Ld[9]=1.088251530500/sqrt(2.0); Ld[10]=0.542813011213/sqrt(2.0); Ld[11]=-0.050256540092/sqrt(2.0);
Ld[12]=-0.045240772218/sqrt(2.0); Ld[13]=0.070703567550/sqrt(2.0); Ld[14]=0.008152816799/sqrt(2.0);
Ld[15]=-0.028786231926/sqrt(2.0); Ld[16]=-0.001137535314/sqrt(2.0); Ld[17]=0.006495728375/sqrt(2.0);
Ld[18]=0.000080661204/sqrt(2.0); Ld[19]=-0.000649589896/sqrt(2.0);
m_preoffset=9;
m_aftoffset=10;
break;
default:
printf("错误的窗口尺寸\n");
//break;
}
for(int i=0;i<FWidth;i++)
{
Hd[i]=pow(-1.0,i+1)*Ld[-i-1+FWidth];
}
for(int i=0;i<FWidth;i++)
{
Lr[i]=Ld[-i-1+FWidth];
Hr[i]=Hd[-i-1+FWidth];
}
//小波分解
LONG lWidth, lHeight;
lWidth=width;
lHeight=hight;
LONG i,j;//循环变量
{
int n;//层数循环变量
//小波变换分解过程循环
for(n=0;n<nLevel;n++)
{
LONG Height,Width;//第n层图象的高度和宽度
Height=long(lHeight>>n);
Width=long(lWidth>>n);
double *LH=new double[Width]; //存放每一行元素
for(i = 0; i < Height; i++)
{
for(j=0;j<Width;j++)
{
LH[j]=f[i*lWidth+j];
}
Convolution( Ld,Hd, FWidth,LH, Width);// 对x方向进行卷积运算
for(j=0;j<Width;j++)
{
f[i*lWidth+j]=LH[j];
}
}
delete LH;
LH=new double[Height]; //存放每一列元素
for(i = 0; i < Width; i++)
{
for(j=0;j<Height;j++)
{
LH[j]=f[i+j*lWidth];
}
Convolution( Ld,Hd, FWidth,LH, Height);// 对y方向进行卷积运算
for(j=0;j<Height;j++)
{
f[i+j*lWidth]=LH[j];
}
}
delete LH;//释放内存
}
{//显示小波变换后的图像需要将数值归一化
//将分解后的值规划处理
for(i=0;i<lHeight;i++)
{
for(j=0;j<lWidth;j++)
{
m_GrayMax=m_GrayMax>f[i * lWidth + j]?m_GrayMax:f[i * lWidth + j];
m_GrayMin=m_GrayMin<f[i * lWidth + j]?m_GrayMin:f[i * lWidth + j];
}
}
// 更新源图像
double Temp;// 中间变量
Temp=255.0/(m_GrayMax-m_GrayMin);
for(i = 0; i < lHeight; i++)// 每列
{
for(j = 0; j < lWidth; j++)// 每行
{
double dTemp;// 中间变量
dTemp = f[i * lWidth + j]; // 计算频谱
dTemp=Temp*(dTemp-m_GrayMin);
f[i * lWidth + j] = dTemp;// 更新源图像
}
}
}
}
return TRUE;// 返回
}
/*************************************************************************
*
* 函数名称：
* DIBDisWavelet()
* 参数:
* double *LF - 使用的小波尺度函数，是常量
* double *HF - 使用的小波母函数，是常量
* int FWidth - 小波窗的宽度
* int nLevel - 小波分解的层数
* 说明:
* 该函数用来对图像进行小波变换重建。
************************************************************************/
BOOL DisWavelet(double *f,int FWidth,int nLevel, int width, int hight)
{
double dTemp;// 中间变量
LONG lWidth,lHeight;
lWidth=width;
lHeight=hight;
LONG i,j;//循环变量
{
// 从源图像中读取数据。
for(i = 0; i < lHeight; i++)// 每列
{
for(j = 0; j < lWidth; j++)// 每行
{
//f[i*lWidth+j] = image[i*lWidth+j];// 给时域赋值
//将规划处理后的值变回原样
f[i*lWidth+j]=(m_GrayMax-m_GrayMin)/255*f[i*lWidth+j]+m_GrayMin;
}
}
int n = 0;//层数循环变量
{
LONG Height,Width;
Height=long(lHeight>>n);
Width=long(lWidth>>n);
double *H00=new double[Height]; //按列存放低低元素
double *H01=new double[Height]; //按列存放低高元素
double *H10=new double[Height]; //按列存放高低元素
double *H11=new double[Height]; //按列存放高高元素
for(i = 0; i < Width/2; i++)
{
for(j=0;j<Height/2;j++)
{
H00[2*j]=f[i+j*lWidth];
H00[2*j+1]=0;
}
for(j=Height/2;j<Height;j++)
{
H01[2*j-Height]=f[i+j*lWidth];
H01[2*j-Height+1]=0;
}
DisConvolution( Lr,Hr, FWidth,H00,H01 ,Height);// 对y方向进行解内积运算
for(j=0;j<Height;j++)
{
f[i+j*lWidth]=H00[j];
}
}
for(i =Width/2 ; i < Width; i++)
{
for(j=0;j<Height/2;j++)
{
H10[2*j]=f[i+j*lWidth];
H10[2*j+1]=0;
}
for(j=Height/2;j<Height;j++)
{
H11[2*j-Height]=f[i+j*lWidth];
H11[2*j-Height+1]=0;
}
DisConvolution( Lr,Hr, FWidth,H10,H11, Height); // 对y方向进行解内积运算
for(j=0;j<Height;j++)
{
f[i+j*lWidth]=H10[j];
}
}
delete H00,H01,H10,H11;//释放内存
double *H0=new double[Width]; //按行存放低元素
double *H1=new double[Width]; //按行存放高元素
for(i = 0; i < Height; i++)
{
for(j=0;j<Width/2;j++)
{
H0[2*j]=f[i*lWidth+j];
H0[2*j+1]=0;
}
for(j=Width/2;j<Width;j++)
{
H1[2*j-Width]=f[i*lWidth+j];
H1[2*j-Width+1]=0;
}
DisConvolution( Lr,Hr, FWidth,H0,H1, Width);// 对x方向进行解卷积运算
for(j=0;j<Width;j++)
{
f[i*lWidth+j]=H0[j];
}
}
delete H0,H1;
}
}
return TRUE;// 返回
}
