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

附录一、快速付里哀变换与反变换程序实例#include #include #define pi (double)3.14159265359/*复数定义*/typedef structdouble 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=1power;/*分配运算所需存储器*/W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);/*计算加权系数*/for(i=0;icount/2;i+)angle=-i*pi*2/count;Wi.re=cos(angle);Wi.im=sin(angle);/*将时域点写入存储器*/memcpy(X1,TD,sizeof(COMPLEX)*count);/*蝶形运算*/for(k=0;kpower;k+)for(j=0;j1k;j+)bfsize=1power-k;for(i=0;ibfsize/2;i+)p=j*bfsize;X2i+p=Add(X1i+p,X1i+p+bfsize/2);X2i+p+bfsize/2=Mul(Sub(X1i+p,X1i+p+bfsize/2),Wi*(1k);X=X1;X1=X2;X2=X;/*重新排序*/for(j=0;jcount;j+)p=0;for(i=0;ipower;i+)if (j&(1i) p+=1power-i-1;FDj=X1p;/*释放存储器*/free(W);free(X1);free(X2);/*快速付里哀反变换,利用快速付里哀变换FD为频域值,TD为时域值,power为2的幂数*/void IFFT(COMPLEX *FD, COMPLEX *TD, int power)int i,count;COMPLEX *x;/*计算付里哀反变换点数*/count=1power;/*分配运算所需存储器*/x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);/*将频域点写入存储器*/memcpy(x,FD,sizeof(COMPLEX)*count);/*求频域点的共轭*/for(i=0;icount;i+)xi.im=-xi.im;/*调用快速付里哀变换*/FFT(x,TD,power);/*求时域点的共轭*/for(i=0;icount;i+)TDi.re/=count;TDi.im=-TDi.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=1power;/*分配运算所需存储器*/X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);/*延拓*/memset(X,0,sizeof(COMPLEX)*count*2);/*将时域点写入存储器*/for(i=0;icount;i+)Xi.re=fi;/*调用快速付里哀变换*/FFT(X,X,power+1);/*调整系数*/s=1/sqrt(count);F0=X0.re*s;s*=sqrt(2);for(i=1;icount;i+)Fi=(Xi.re*cos(i*pi/(count*2)+Xi.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=1power;/*分配运算所需存储器*/X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);/*延拓*/memset(X,0,sizeof(COMPLEX)*count*2);/*调整频域点并写入存储器*/for(i=0;icount;i+)Xi.re=Fi*cos(i*pi/(count*2);Xi.im=Fi*sin(i*pi/(count*2);/*调用快速付里哀反变换*/IFFT(X,X,power+1);/*调整系数*/s=1/sqrt(count);for(i=0;icount;i+)fi=(1-sqrt(2)*s*F0+sqrt(2)*s*Xi.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=1power;/*分配运算所需存储器*/X1=(double *)malloc(sizeof(double)*count);X2=(double *)malloc(sizeof(double)*count);/*将时域点写入存储器*/memcpy(X1,f,sizeof(double)*count);/*蝶形运算*/for(k=0;kpower;k+)for(j=0;j1k;j+)bfsize=1power-k;for(i=0;ibfsize/2;i+)p=j*bfsize;X2i+p=X1i+p+X1i+p+bfsize/2;X2i+p+bfsize/2=X1i+p-X1i+p+bfsize/2;X=X1;X1=X2;X2=X;/*调整系数*/for(i=0;icount;i+)Wi=X1i/count;/*释放存储器*/free(X1);free(X2);/*快速沃尔什-哈达玛反变换,利用快速沃尔什-哈达玛变换F为变换域值,f为时域值,power为2的幂数*/void IWALh(double *W, double *f, int power)int i,count;count=1power;/*调用快速沃尔什-哈达玛变换*/WALh(W,f,power);/*调整系数*/for(i=0;icount;i+)fi*=count;附录四、二维小波分解与重构程序实例/ *（基于Windows的二维小波分解与重构BC+4。5）*/#include #include #include #include #include #include #include #include #include #include #include #include #include stdio.h#include stdlib.h#include bmpview.h#define PI 3.1415926#define SIZE 256#define DD 13float hDD=-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 gDD=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 hiDD,giDD;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 *imgSIZE,*imgxSIZE,*imgySIZE;/*for evolution agents*/int AgentNum,ActiveAgent;int SolutionX1000,SolutionY1000;int SearchNum,AgentLife;int Agent10002;int isAgentAlive1000;int SearchTime;unsigned short Image3232;#define MAXAPPNAME 20static 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. SubWinPtrType is set to 0 to indicate that the window/ has be closed./TSubWindow:TSubWindow() SubWinPtr = 0;voidTSubWindow:Paint(TDC& dc, bool, TRect&)class TBmpViewWindow : virtual public TWindow, public TClipboardViewer public: char FileNameMAXPATH; 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 *f256,*g256; 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 ss5; int i,j,k,x,y; AgentNum=5;ActiveAgent=5;SearchNum=0; AgentLife=2; SearchTime=0; for(i=0;iAgentNum;i+)Agenti0=5*i+1;Agenti1=5*i+1;isAgentAlivei=1;/MessageBox(ss,ss,MB_OK);if(pMemDC)for(x=0;x32;x+) for(y=0;yGetPixel(x,y);/ imgij=fij; while(1)MessageBox(ss,ss,MB_OK); if(ActiveAgent=0)break; for( i=0;iSetPixel(x,y, (unsigned char)Imagexy); AdjustScroller();SearchNum+=1;/*judge direction*/ if(x+1)=0 & Imagex-1y=0) x=x-1;Agenti0=x; elseif(y-1)=0&Imagexy-1=0)y=y-1;Agenti1=y; else if(y+1)32&Imagexy+1=0)y=y+1;Agenti1=y;else isAgentAlivei=0; ActiveAgent-=1;else /*judge direction*/ if(x+1)=0 & Imagex-1y=0) x=x-1;Agenti0=x; elseif(y-1)=0&Imagexy-1=0)y=y-1;Agenti1=y; else if(y+1)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;iGetName(), MB_OK);CloseWindow(0);if(imgyi=(float far*)farmalloc(SIZE*sizeof(float)=NULL) MessageBox(Cannot malloc memory for displaying bitmap file, GetApplication()-GetName(), MB_OK);CloseWindow(0);if(imgi=(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;iSIZE;i+) fread(buf,1,SIZE,inputfile); for(j=0;jSIZE;j+) imgij=(float)(bufj)-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;iSIZE;i+) for(j=0;jSIZE;j+) bufj=(unsigned char)(s(imgij)+128); imgij=(float)bufj-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;iSIZE;i+) for(j=0;jSIZE;j+) bufj=(unsigned char)(s(imgij)+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(imgyi); farfree(imgxi); farfree(imgi); if(wavelet_direction=1)wavelet_direction=0; else wavelet_direction=1;/Edge Process/int a(int x,int xsize) if(x=xsize)x=xsize*2-x-2; return(x); /Threshold/int s(float x) if(x127)return(127); if(x-128)return(-128); return(x); /Set Inverse Filter Coefficients/void coef() int i; for(i=0;iDD;i+) hii=hDD-1-i;gii=gDD-1-i; /*二维小波分解*/void wt(int xs,int ys,long xsize,long ysize) int i,j,k,n; float temp1,temp2; float bufferxSIZE,bufferySIZE; for(n=0;nys;n+) for(i=0;ixs;i+=2) temp1=0; temp2=0; for(j=-(DD-1)/2;j=(DD-1)/2;j+) temp1=temp1+hj+(DD-1)/2*imgna(i+j,xs); for(j=-(DD-1)/2+1;j=(DD-1)/2+1;j+) temp2=temp2+gj+(DD-1)/2-1*imgna(i+j,xs); bufferxi/2=temp1; bufferxi/2+xs/2=temp2; for(k=0;kxs;k+)imgnk=bufferxk; /*Filter in the vertical direction*/ for(n=0;nxs;n+) for(i=0;iys;i+=2) temp1=0; temp2=0; for(j=-(DD-1)/2;j=(DD-1)/2;j+) temp1=temp1+hj+(DD-1)/2*imga(i+j,ys)n; for(j=-(DD-1)/2+1;j=(DD-1)/2+1;j+) temp2=temp2+gj+(DD-1)/2-1*imga(i+j,ys)n; bufferyi/2=temp1; bufferyi/2+ys/2=temp2; for(k=0;kys;k+)imgkn=bufferyk; /*二维小波重构*/void iwt(int xs,int ys,long xsi