数字图像处理C语言程序实例.doc

附录一、快速付里哀变换与反变换程序实例

#include <math.h>
#include <malloc.h>

#define pi (double)3.14159265359

/*复数定义*/
typedef struct
{
double re;
double im;
}COMPLEX;

/*复数加运算*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
COMPLEX c;
c.re=c1.re+c2.re;
c.im=c1.im+c2.im;
return c;
}

/*复数减运算*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
COMPLEX c;
c.re=c1.re-c2.re;
c.im=c1.im-c2.im;
return c;
}

/*复数乘运算*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
COMPLEX c;
c.re=c1.re*c2.re-c1.im*c2.im;
c.im=c1.re*c2.im+c2.re*c1.im;
return c;
}

/*快速付里哀变换
TD为时域值,FD为频域值,power为2的幂数*/
void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{
int count;
int i,j,k,bfsize,p;
double angle;
COMPLEX *W,*X1,*X2,*X;

/*计算付里哀变换点数*/
count=1<<power;
/*分配运算所需存储器*/
W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);
X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
/*计算加权系数*/
for(i=0;i<count/2;i++)
{
angle=-i*pi*2/count;
W[i].re=cos(angle);
W[i].im=sin(angle);
}
/*将时域点写入存储器*/
memcpy(X1,TD,sizeof(COMPLEX)*count);
/*蝶形运算*/
for(k=0;k<power;k++)
{
for(j=0;j<1<<k;j++)
{
bfsize=1<<power-k;
for(i=0;i<bfsize/2;i++)
{
p=j*bfsize;
X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
}
}
X=X1;
X1=X2;
X2=X;
}
/*重新排序*/
for(j=0;j<count;j++)
{
p=0;
for(i=0;i<power;i++)
{
if (j&(1<<i)) p+=1<<power-i-1;
}
FD[j]=X1[p];
}
/*释放存储器*/
free(W);
free(X1);
free(X2);
}

/*快速付里哀反变换,利用快速付里哀变换
FD为频域值,TD为时域值,power为2的幂数*/
void IFFT(COMPLEX *FD, COMPLEX *TD, int power)
{
int i,count;
COMPLEX *x;

/*计算付里哀反变换点数*/
count=1<<power;
/*分配运算所需存储器*/
x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
/*将频域点写入存储器*/
memcpy(x,FD,sizeof(COMPLEX)*count);
/*求频域点的共轭*/
for(i=0;i<count;i++)
{
x[i].im=-x[i].im;
}
/*调用快速付里哀变换*/
FFT(x,TD,power);
/*求时域点的共轭*/
for(i=0;i<count;i++)
{
TD[i].re/=count;
TD[i].im=-TD[i].im/count;
}
/*释放存储器*/
free(x);
}
附录二、快速余弦变换与反变换程序实例
（利用2N点付里哀变换实现快速余弦变换）
/*快速离散余弦变换,利用快速付里哀变换
f为时域值,F为变换域值,power为2的幂数*/
void DCT(double *f, double *F, int power)
{
int i,count;
COMPLEX *X;
double s;

/*计算离散余弦变换点数*/
count=1<<power;
/*分配运算所需存储器*/
X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);
/*延拓*/
memset(X,0,sizeof(COMPLEX)*count*2);
/*将时域点写入存储器*/
for(i=0;i<count;i++)
{
X[i].re=f[i];
}
/*调用快速付里哀变换*/
FFT(X,X,power+1);
/*调整系数*/
s=1/sqrt(count);
F[0]=X[0].re*s;
s*=sqrt(2);
for(i=1;i<count;i++)
{
F[i]=(X[i].re*cos(i*pi/(count*2))+X[i].im*sin(i*pi/(count*2)))*s;
}
/*释放存储器*/
free(X);
}

/*快速离散余弦反变换,利用快速付里哀反变换
F为变换域值,f为时域值,power为2的幂数*/
void IDCT(double *F, double *f, int power)
{
int i,count;
COMPLEX *X;
double s;

/*计算离散余弦反变换点数*/
count=1<<power;
/*分配运算所需存储器*/
X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);
/*延拓*/
memset(X,0,sizeof(COMPLEX)*count*2);
/*调整频域点并写入存储器*/
for(i=0;i<count;i++)
{
X[i].re=F[i]*cos(i*pi/(count*2));
X[i].im=F[i]*sin(i*pi/(count*2));
}
/*调用快速付里哀反变换*/
IFFT(X,X,power+1);
/*调整系数*/
s=1/sqrt(count);
for(i=0;i<count;i++)
{
f[i]=(1-sqrt(2))*s*F[0]+sqrt(2)*s*X[i].re*count*2;
}
/*释放存储器*/
free(X);
}

附录三、快速Walsh-Hadamard变换与反变换程序实例

/*快速沃尔什-哈达玛变换
f为时域值,F为变换域值,power为2的幂数*/
void WALh(double *f, double *W, int power)
{
int count;
int i,j,k,bfsize,p;
double *X1,*X2,*X;

/*计算快速沃尔什变换点数*/
count=1<<power;
/*分配运算所需存储器*/
X1=(double *)malloc(sizeof(double)*count);
X2=(double *)malloc(sizeof(double)*count);
/*将时域点写入存储器*/
memcpy(X1,f,sizeof(double)*count);
/*蝶形运算*/
for(k=0;k<power;k++)
{
for(j=0;j<1<<k;j++)
{
bfsize=1<<power-k;
for(i=0;i<bfsize/2;i++)
{
p=j*bfsize;
X2[i+p]=X1[i+p]+X1[i+p+bfsize/2];
X2[i+p+bfsize/2]=X1[i+p]-X1[i+p+bfsize/2];
}
}
X=X1;
X1=X2;
X2=X;
}
/*调整系数*/
for(i=0;i<count;i++)
{
W[i]=X1[i]/count;
}
/*释放存储器*/
free(X1);
free(X2);
}

/*快速沃尔什-哈达玛反变换,利用快速沃尔什-哈达玛变换
F为变换域值,f为时域值,power为2的幂数*/
void IWALh(double *W, double *f, int power)
{
int i,count;
count=1<<power;

/*调用快速沃尔什-哈达玛变换*/
WALh(W,f,power);
/*调整系数*/
for(i=0;i<count;i++)
{
f[i]*=count;
}
}

附录四、二维小波分解与重构程序实例

// *********（基于Windows的二维小波分解与重构BC++4。5）**********//
#include <owl/owlpch.h>
#include <owl/applicat.h>
#include <owl/dc.h>
#include <owl/menu.h>
#include <owl/framewin.h>
#include <owl/scroller.h>
#include <owl/opensave.h>
#include <owl/clipview.h>
#include <string.h>
#include <alloc.h>
#include <dir.h>
#include <math.h>
#include "stdio.h"
#include "stdlib.h"
#include "bmpview.h"

#define PI 3.1415926
#define SIZE 256

#define DD 13
float h[DD]={-0.00332761,0.00569794,0.0196637,-0.0482603,-0.0485391,
0.292562,0.564406,0.292562,-0.0485391,-0.0482602,-0.0196637,
0.00569794,-0.0033276};
float g[DD]={0.00332761,0.00569794,-0.0196637,-0.0482603,0.0485391,
0.292562,-0.564406,0.292562,0.0485391,-0.0482602,0.0196637,
0.00569794,0.0033276};
float hi[DD],gi[DD];

int wavelet_direction=1;

int a(int x,int xsize);
//Threshold//
int s(float x);
//Set Inverse Filter Coefficients//
void coef();
//********Wavelet Transform************************//

void wt(int xs,int ys,long xsize,long ysize);
//*****Inverse Wavelet Transform*********************//
void iwt(int xs,int ys,long xsize,long ysize);

float *img[SIZE],*imgx[SIZE],*imgy[SIZE];

//****for evolution agents*****//
int AgentNum,ActiveAgent;
int SolutionX[1000],SolutionY[1000];
int SearchNum,AgentLife;
int Agent[1000][2];
int isAgentAlive[1000];
int SearchTime;
unsigned short Image[32][32];

#define MAXAPPNAME 20
static const char AppName[] = "Image Processing";

//
// TBmpViewWindow, a Bitmap displaying window derived from TClipboardViewer to
// facilitate receiving of clipboard change notifications. Could mix it in if
// an additional base was desired.
//
class TSubWindow : public TFrameWindow {
public:

TSubWindow(TWindow* parent);
~TSubWindow();

protected:
void EvSize(UINT sizeType, TSize& size)
{Invalidate(); TFrameWindow::EvSize(sizeType, size);}
void Paint(TDC& dc, bool, TRect&);

DECLARE_RESPONSE_TABLE(TSubWindow);
};

DEFINE_RESPONSE_TABLE1(TSubWindow, TFrameWindow)
EV_WM_SIZE,
END_RESPONSE_TABLE;

// pointers to different child windows.
//
TWindow* SubWinPtr = 0;
TSubWindow::TSubWindow(TWindow* parent)
: TFrameWindow(parent)
{
Attr.Style |= WS_VISIBLE | WS_POPUP | WS_OVERLAPPEDWINDOW,
Attr.X = 100;
Attr.Y = 100;
Attr.W = 300;
Attr.H = 400;
}

//
// Destroy window. SubWinPtr[Type] is set to 0 to indicate that the window
// has be closed.
//
TSubWindow::~TSubWindow()
{
SubWinPtr = 0;
}
void
TSubWindow::Paint(TDC& dc, bool, TRect&)
{
}

class TBmpViewWindow : virtual public TWindow, public TClipboardViewer {
public:
char FileName[MAXPATH];
TDib* Dib;
TBitmap* Bitmap;
TMemoryDC * pMemDC;
TPalette* Palette;
TBrush* BkgndBrush;
long Rop;
int PixelWidth;
int PixelHeight;
WORD Colors;
bool Fit;
bool AutoClipView;
// unsigned char far *f[256],*g[256];

TBmpViewWindow();
~TBmpViewWindow();
void ShowSubWindow(TWindow* parent);
protected:
void CmFileOpen();
void CmCopy();
void CmPaste();
void CmFit();
void CmAutoClipView();
void CeCopy(TCommandEnabler& ce);
void CePaste(TCommandEnabler& ce);
void CeFit(TCommandEnabler& ce);
void CeAutoClipView(TCommandEnabler& ce);
//*** Image Processing**//
void CmWavelet();
void CmAgent();
//**********************//
void Paint(TDC&, bool erase, TRect&);
void EvSize(UINT sizeType, TSize&);

void EvPaletteChanged(HWND hWndPalChg);
bool EvQueryNewPalette();
void EvSetFocus(HWND); // should use above when working
void EvDrawClipboard();
void EvDestroy();

bool UpdatePalette(bool alwaysRepaint);
void AdjustScroller();
void SetCaption(const char*);
void SetupFromDib(TDib* dib);
bool LoadBitmapFile(const char*);
bool LoadBitmapResource(WORD ResId);

DECLARE_RESPONSE_TABLE(TBmpViewWindow);
};

DEFINE_RESPONSE_TABLE2(TBmpViewWindow, TClipboardViewer, TWindow)
EV_COMMAND(CM_FILEOPEN, CmFileOpen),
EV_COMMAND(CM_EDITCOPY, CmCopy),
EV_COMMAND(CM_EDITPASTE, CmPaste),
EV_COMMAND(CM_FIT, CmFit),
EV_COMMAND(CM_AUTOCLIPVIEW, CmAutoClipView),
EV_COMMAND_ENABLE(CM_EDITCOPY, CeCopy),
EV_COMMAND_ENABLE(CM_EDITPASTE, CePaste),
EV_COMMAND_ENABLE(CM_FIT, CeFit),
EV_COMMAND_ENABLE(CM_AUTOCLIPVIEW, CeAutoClipView),

//*****Image Processing***************//
EV_COMMAND(CM_WAVELET,CmWavelet),
EV_COMMAND(CM_AGENT,CmAgent),

//***************************************//
EV_WM_SIZE,
EV_WM_PALETTECHANGED,
EV_WM_QUERYNEWPALETTE,
EV_WM_SETFOCUS,
EV_WM_DRAWCLIPBOARD,
EV_WM_DESTROY,
END_RESPONSE_TABLE;
//****agent search******//
//****小波变换菜单Wavelet*******//
//****原图象存在lena.bmp***//
//****运行后,分解图象存在lenawt.bmp******//
//****重建图象存在lenaiwt.bmp***********//
void TBmpViewWindow::CmAgent()
{
char ss[5];
int i,j,k,x,y;
AgentNum=5;ActiveAgent=5;SearchNum=0;
AgentLife=2;
SearchTime=0;
for(i=0;i<AgentNum;i++){Agent[i][0]=5*i+1;Agent[i][1]=5*i+1;isAgentAlive[i]=1;}

// MessageBox("ss","ss",MB_OK);
if(pMemDC)
{
for(x=0;x<32;x++)
for(y=0;y<32;y++)
Image[x][y]=(unsigned short)pMemDC->GetPixel(x,y);
// img[i][j]=f[i][j];
while(1){
MessageBox("ss","ss",MB_OK);
if(ActiveAgent<=0)break;
for( i=0;i<AgentNum;i++)
{
if(isAgentAlive[i]==1)
{
x=Agent[i][0];y=Agent[i][1];
if(Image[x][y]== 0)
{
SolutionX[SearchNum]=x;SolutionY[SearchNum]=y;
Image[x][y]=255;//***marking this point
pMemDC->SetPixel(x,y,
(unsigned char)Image[x][y]);
AdjustScroller();
SearchNum+=1;
//**judge direction**//
if((x+1)<32&&Image[x+1][y]==0){x=x+1;Agent[i][0]=x;}
else if((x-1)>=0 && Image[x-1][y]==0) {x=x-1;Agent[i][0]=x;}
else if((y-1)>=0&&Image[x][y-1]==0){y=y-1;Agent[i][1]=y;}
else if((y+1)<32&&Image[x][y+1]==0){y=y+1;Agent[i][1]=y;}
else {isAgentAlive[i]=0; ActiveAgent-=1;}
}
else
{ //**judge direction**//
if((x+1)<32&&Image[x+1][y]==0){x=x+1;Agent[i][0]=x;}
else if((x-1)>=0 && Image[x-1][y]==0) {x=x-1;Agent[i][0]=x;}
else if((y-1)>=0&&Image[x][y-1]==0){y=y-1;Agent[i][1]=y;}
else if((y+1)<32&&Image[x][y+1]==0){y=y+1;Agent[i][1]=y;}
else {isAgentAlive[i]=0; ActiveAgent-=1;}
}
SearchTime++;
}//**end if isAgentAlive==1
}//**end for
}//**end while
sprintf(ss,"%d",SearchTime);
MessageBox(ss,"Total Search times",MB_OK);
}

}

void
TBmpViewWindow::CmWavelet()
{
unsigned i,j,k;int xs,ys;
unsigned char buf[SIZE+1],buf1[1078];
FILE *inputfile,*mfile,*outputfile;
inputfile=fopen("lena.bmp","rb");
if(!inputfile)
{ MessageBox("Cannot open file lena.raw", GetApplication()->GetName(), MB_OK);
CloseWindow(0);
}
mfile=fopen("lenawt.bmp","wb");
if(!mfile)
{ MessageBox("Cannot open file lenawt.raw", GetApplication()->GetName(), MB_OK);
CloseWindow(0);
}
outputfile=fopen("lenaiwt.bmp","wb");
if(!outputfile)
{ MessageBox("Cannot open file lenawt.raw", GetApplication()->GetName(), MB_OK);
CloseWindow(0);
}
for(i=0;i<=SIZE-1;i++)
{
if((imgx[i]=(float far*)farmalloc(SIZE*sizeof(float)))==NULL)
{ MessageBox("Cannot malloc memory for displaying bitmap file", GetApplication()->GetName(), MB_OK);
CloseWindow(0);
}
if((imgy[i]=(float far*)farmalloc(SIZE*sizeof(float)))==NULL)
{ MessageBox("Cannot malloc memory for displaying bitmap file", GetApplication()->GetName(), MB_OK);
CloseWindow(0);
}
if((img[i]=(float far*)farmalloc(SIZE*sizeof(float)))==NULL)
{ MessageBox("Cannot malloc memory for displaying bitmap file", GetApplication()->GetName(), MB_OK);
CloseWindow(0);
}

}

//read input file
fread(buf1,1,1078,inputfile);
for(i=0;i<SIZE;i++)
{
fread(buf,1,SIZE,inputfile);
for(j=0;j<SIZE;j++)
img[i][j]=(float)(buf[j])-128.0;
}
fclose(inputfile);

//start wt....
coef();

xs=SIZE;
ys=SIZE;
wt(xs,ys,SIZE,SIZE);
xs=xs/2;ys=ys/2;
wt(xs,ys,SIZE,SIZE);
xs=xs/2;ys=ys/2;
wt(xs,ys,SIZE,SIZE);

//*****write the wt result file*****//
fwrite(buf1,1,1078,mfile);
for(i=0;i<SIZE;i++)
{
for(j=0;j<SIZE;j++)
{
buf[j]=(unsigned char)(s(img[i][j])+128);
img[i][j]=(float)buf[j]-128.0;
}
fwrite(buf,1,SIZE,mfile);
}
fclose(mfile);

//****start inverse wt....***//
iwt(xs,ys,SIZE,SIZE);
xs=2*xs;ys=2*ys;
iwt(xs,ys,SIZE,SIZE);
xs=2*xs;ys=2*ys;
iwt(xs,ys,SIZE,SIZE);
fwrite(buf1,1,1078,outputfile);
for(i=0;i<SIZE;i++)
{
for(j=0;j<SIZE;j++)
{
buf[j]=(unsigned char)(s(img[i][j])+128);
}
fwrite(buf,1,SIZE,outputfile);
}
fclose(outputfile);
/* LoadBitmapFile("lenawt.bmp");
SetCaption(strlwr("lenawt.bmp"));
LoadBitmapFile("lenaiwt.bmp");
SetCaption(strlwr("lenaiwt.bmp"));
*/
for(i=0;i<=SIZE-1;i++)
{
farfree(imgy[i]);
farfree(imgx[i]);
farfree(img[i]);
}
if(wavelet_direction==1)wavelet_direction=0;
else wavelet_direction=1;
}
//Edge Process//
int a(int x,int xsize)
{
if(x<0)x=-x;
if(x>=xsize)x=xsize*2-x-2;
return(x);
}

//Threshold//
int s(float x)
{
if(x>127)return(127);
if(x<-128)return(-128);
return(x);
}
//Set Inverse Filter Coefficients//
void coef()
{
int i;
for(i=0;i<DD;i++)
{
hi[i]=h[DD-1-i];
gi[i]=g[DD-1-i];
}
}
//********二维小波分解************************//

void wt(int xs,int ys,long xsize,long ysize)
{
int i,j,k,n;
float temp1,temp2;
float bufferx[SIZE],buffery[SIZE];
for(n=0;n<ys;n++)
{

for(i=0;i<xs;i+=2)
{
temp1=0;
temp2=0;
for(j=-(DD-1)/2;j<=(DD-1)/2;j++)
temp1=temp1+h[j+(DD-1)/2]*img[n][a(i+j,xs)];
for(j=-(DD-1)/2+1;j<=(DD-1)/2+1;j++)
temp2=temp2+g[j+(DD-1)/2-1]*img[n][a(i+j,xs)];
bufferx[i/2]=temp1;
bufferx[i/2+xs/2]=temp2;
}
for(k=0;k<xs;k++)
img[n][k]=bufferx[k];
}
//*****Filter in the vertical direction*****//
for(n=0;n<xs;n++)
{

for(i=0;i<ys;i+=2)
{
temp1=0;
temp2=0;
for(j=-(DD-1)/2;j<=(DD-1)/2;j++)
temp1=temp1+h[j+(DD-1)/2]*img[a(i+j,ys)][n];
for(j=-(DD-1)/2+1;j<=(DD-1)/2+1;j++)
temp2=temp2+g[j+(DD-1)/2-1]*img[a(i+j,ys)][n];
buffery[i/2]=temp1;
buffery[i/2+ys/2]=temp2;
}
for(k=0;k<ys;k++)
img[k][n]=buffery[k];
}
}
//*****二维小波重构*********************//
void iwt(int xs,int ys,long xsi