快速中值滤波及c语言实现

快速中值滤波及c语言实现学生姓名：刘勇学号：6100410218 专业班级：数媒101【摘要】本文讨论了用c语言在微机上实现中值滤波及快速算法，在程序设计的过程中充分考虑到程序运行的时间复杂度和空间复杂度的问题．解决了由于图像太大而内存不够的问题，运用对程序运行时的方法，得出在PENTIUM-S100MHz上中值滤渡的一般算法运行4．23秒．而快速算法运行2 58秒。

【关键词】c语言；中值滤波；快速算法1 引言中值滤波是涂基发明的一种非线性信号处理技术，对抑制图像的噪声非常有效，在二维形式下，中值滤渡器是一个古有奇数个像素的滑动窗口，窗口正中的象素的灰度值用窗口内各个象素的中值代替窗口的中值为窗口中象素按大小顺序排列后处于中间位置的象素；本文讨论中值滤的一般算法并比较其运算速度。

2 用C语言实现算法的若干问题在设计算法编制程序的时候，我们充分考虑到程序运行的时间复杂度和空间复杂度问题，在解决问题的前提下，使算法尽量简单，使程序运行占有的空间尽量的小，这样来减少不必要的时问浪费和空间浪费，从而太大的提高程序执行的效率。

首先考虑到的内存问题。

由于在本文算法中用的图像是512+512 8bit，这就存在一个内存不够大一整幅图像不能一次性调入的问题。

为了解受此问题，可以只开辟一个3"512的缓冲区n，将原图像采用分批调入缓冲区，使内存不够的问题得到了圆满的解决。

另外为了对中值滤波的快速算法和普通算法进行精确的比较，采用对程序运行计时的方法，并精确计算每个算法运行的时间，使得出的结论更可靠。

3 中值滤波算法的C语言程序实现本算法采用对开辟的3*512的缓冲区从左到右依次形成一个3*3的窗口．然后将此3*3的窗口放人一个一维数组中，调用求中值子函数．通过排序得出中值，当此中值不等于窗口中间位置的象素时．用此中值来代替窗VI中间位置的象素灰度值．若此缓冲区处理完毕后，将缓冲区的第一行存入新建的文件中，将第二、第三行分别向上移动一行，若存人新建的文件中的行数小于或等于511(即这样处理的行数小于或等于511)，则从原文件中调入一行作为缓冲区第三行，按上述方法进行直到处理的总行数等于511为止，最后，将缓冲区的第二、三行存人新建的文件，程序流程框图如图14 中值滤波快速算法的C语言程序实现本算法充分利用了上一次处理的结果．采用迭代，逐次逼近的方法得到本次的中值，在一行处理完毕后转人下一行也采用走S型的方法．这样除第一个窗口采用了一伏排序得到中值外，其它的窗口都利用上伏的窗口的象素删除无用的3个象素后再加人新的3个象素，利用迭代的方法得到本次窗口的中值．这样太大地提高了程序执行的效率。

4．1算法的解释首先是开辟一个3*512的缓冲区a，在初始化缓冲区时考虑到时间复杂度的问题，所以只初始化了第二、三行，而对第一行只初始化了前三个象素，这样便在缓冲区中可以得到一个3*3的窗口，对此窗口进行排序求中值后得出第一个窗口的中间象素的值．将文件指针定位在2*512处。

然后开始循环，当处理的行数小于或等于511时，将缓冲区a中的第二、三行分别向上移动一行变为第一、二行，从文件中读人512个字节作为缓冲区的第三行，并用行数模2的方法设置方向标志k．当k为0时，从左向右移动窗口．当k为1时．从右向左移动窗口．而每一窗口都利用上次的窗口的像素删除无用的3个象素后再加入新的3个象素．利用迭代的方法从上次的中值得到本次的中值。

当处理完一行后将缓冲区的第一行存入新建的文件中，最后将缓冲区的第二、三行存入文件中。

4．2 算法代码// ImageProcessingDoc.cpp : implementation of the CImageProcessingDoc class//#include "stdafx.h"#include "ImageProcessing.h"#include "ImageProcessingDoc.h"#include "GreyRatio.h"#include <math.h>#define PI (acos(0.0) * 2)#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif///////////////////////////////////////////////////////////////////// ////////// CImageProcessingDocIMPLEMENT_DYNCREATE(CImageProcessingDoc, CDocument)BEGIN_MESSAGE_MAP(CImageProcessingDoc, CDocument)//{{AFX_MSG_MAP(CImageProcessingDoc)ON_COMMAND(ID_HISTOGRAM_ADJUSTIFCATION, OnHistogramAdjustifcation) ON_COMMAND(ID_FFT, OnFft)ON_COMMAND(ID_SALT_PEPPER_NOICE, OnSaltPepperNoice)ON_COMMAND(ID_RANDOM_NOISE, OnRandomNoise)ON_COMMAND(ID_MEDIAN_FILTERING, OnMedianFiltering)ON_COMMAND(ID_DCT, OnDct)ON_COMMAND(ID_FWT, OnFwt)ON_COMMAND(ID_DHT, OnDht)ON_COMMAND(ID_WAVELET_TRANSFORM, OnWaveletTransform)ON_COMMAND(ID_GREY_ADJUSTIFCATION, OnGreyAdjustifcation)ON_COMMAND(ID_GREY_LINEAR_ADJUSTIFCATION, OnGreyLinearAdjustifcation)ON_COMMAND(ID_GREY_SEGLINEAR_ADJUSTIFCATION, OnGreySeglinearAdjustifcation)ON_COMMAND(ID_2DGRAD, On2dgrad)ON_COMMAND(ID_ROBERT, OnRobert)//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////// ////////// CImageProcessingDoc construction/destructionCImageProcessingDoc::CImageProcessingDoc(){// TODO: add one-time construction code heremImageFile = NULL;bFileIsLoad = FALSE;nRows = 256;nCols = 256;mSourceData = NULL;pSourceData = NULL;bDataIsProcessed = FALSE;mResultData = FALSE;pResultData = FALSE;FourierDataR = NULL;FourierDataI = NULL;}CImageProcessingDoc::~CImageProcessingDoc(){}BOOL CImageProcessingDoc::OnNewDocument(){if (!CDocument::OnNewDocument())return FALSE;// TODO: add reinitialization code here// (SDI documents will reuse this document)return TRUE;}///////////////////////////////////////////////////////////////////// ////////// CImageProcessingDoc serializationvoidCImageProcessingDoc::Serialize(CArchive&ar){if (ar.IsStoring()){// TODO: add storing code here}else{// TODO: add loading code here}}///////////////////////////////////////////////////////////////////// ////////// CImageProcessingDoc diagnostics#ifdef _DEBUGvoidCImageProcessingDoc::AssertValid() const{CDocument::AssertValid();}voidCImageProcessingDoc::Dump(CDumpContext& dc) const{CDocument::Dump(dc);}#endif //_DEBUG///////////////////////////////////////////////////////////////////// ////////// CImageProcessingDoc commandsBOOL CImageProcessingDoc::OnOpenDocument(LPCTSTR lpszPathName) {int x;int y;if (!CDocument::OnOpenDocument(lpszPathName))return FALSE;// TODO: Add your specialized creation code hereif(mSourceData){free(mSourceData);mSourceData = NULL;}if (!(mSourceData = (unsigned char*)malloc(nRows*nCols*sizeof(unsigned char))))return FALSE;if (pSourceData){free(pSourceData);pSourceData = NULL;}if (!(pSourceData = (unsigned char*)malloc(3*nRows*nCols*sizeof(unsigned char))))return FALSE;if (mResultData){free(mResultData);mResultData = NULL;}if (!(mResultData = (unsigned char*)malloc(nRows*nCols*sizeof(unsigned char))))return FALSE;if (pResultData){free(pResultData);pResultData = NULL;}if (!(pResultData = (unsigned char*)malloc(3*nRows*nCols*sizeof(unsigned char))))return FALSE;if (mImageFile){fclose(mImageFile);mImageFile = NULL;}if (!(mImageFile = fopen(lpszPathName,"rb"))){free(mSourceData);return FALSE;}if (fread(mSourceData,sizeof(unsignedchar),nRows*nCols,mImageFile) != (unsigned)nCols*nRows) {free(mSourceData);fclose(mImageFile);mImageFile = NULL;bFileIsLoad = false;return FALSE;}for(y = 0; y <nRows; y++)for(x = 0; x <nCols; x++){pSourceData[3*y*nCols+3*x] = mSourceData[y*nCols+x]; pSourceData[3*y*nCols+3*x+1] = mSourceData[y*nCols+x]; pSourceData[3*y*nCols+3*x+2] = mSourceData[y*nCols+x];}bFileIsLoad = TRUE;return TRUE;}voidCImageProcessingDoc::OnHistogramAdjustifcation(){// TODO: Add your command handler code hereintx,y;double *mR;double *mS;mR = new double[256];mS = new double[256];for(x=0;x<256;x++){mR[x] = mS[x] = 0.0;}//统计直方图for(y = 0; y <nRows; y++)for(x = 0; x <nCols; x++){mR[mSourceData[y*nCols+x]] ++;}for(x=0;x<256;x++){for(y=0;y<x;y++)mS[x] += mR[y];mS[x] /= nRows*nCols;}//直方图变换for(y = 0; y <nRows; y++)for(x = 0; x <nCols; x++)mResultData[y*nRows+x] = (char) (255* mS[mSourceData[y*nRows+x]]);//灰度计算for(y = 0; y <nRows; y++)for(x = 0; x <nCols; x++){pResultData[3*y*nCols+3*x] = mResultData[y*nCols+x]; pResultData[3*y*nCols+3*x+1] = mResultData[y*nCols+x]; pResultData[3*y*nCols+3*x+2] = mResultData[y*nCols+x];}//更新显示UpdateAllViews(NULL);}// FFTandIFFT 一维傅立叶变换与你变换函数// 输入时域数据实部Tr,虚部Ti// 输出频域数据实部Tr,虚部Ti// 序列长度N，N等于2的r次幂// FFTorIFFT,逻辑变量，非零做正变换，零做反变换voidCImageProcessingDoc::FFTandIFFT(float *Tr, float *Ti, int N, boolFFTorIFFT){int r; //迭代次数intl,j,k;//循环变量int p; //用于蝶形计算加权系数的指数int B; //对偶结点距离float X,Y,XX,YY;float w;float cosw,sinw;if (!FFTorIFFT){ //如果做傅立叶你变换，则必须对数列除以N for(l=0;l<N;l++){Tr[l] /= N;Ti[l] /= N;}}//计算循环次数rr = 0; l = N;while(l /= 2) r++;//倒序int LH = N/2;inti;float temp;j = 0;for (i=1;i<N-1;i++){k = LH;while(j>=k){j = j-k;k = k/2;}j = j + k;if (i<=j){temp = Tr[i]; Tr[i] = Tr[j]; Tr[j] = temp;temp = Ti[i]; Ti[i] = Ti[j]; Ti[j] = temp;}}for(l=0; l <= r; l++) //共r级{B = 1<<(l-1); // 第l层对偶结点距离为2^(l-1)for(j=0; j <B;j++){p = j*(1<<(r-l));w = 2*PI*p/N;for(k=j;k<N-1;k+=(1<<l)){if (FFTorIFFT){ //若做傅立叶正变换cosw =cos(-w);sinw =sin(-w);}else{ //傅立叶反变换cosw =cos(w);sinw =sin(w);}X = Tr[k] + Tr[k+B]*cosw - Ti[k+B] * sinw;Y = Ti[k] +Tr[k+B]*sinw + Ti[k+B] * cosw; XX = Tr[k] - Tr[k+B]*cosw + Ti[k+B] * sinw;YY = Ti[k] - Tr[k+B]*sinw - Ti[k+B] * cosw;Tr[k] = X;Ti[k] = Y;Tr[k+B] = XX;Ti[k+B] = YY;}}}}voidCImageProcessingDoc::OnFft(){// TODO: Add your command handler code hereinti,j;intii,jj;float temp;float *Tr;float *Ti;Tr = new float[nCols];Ti = new float[nCols];if ( FourierDataR){deleteFourierDataR;FourierDataR = NULL;}if ( FourierDataI){deleteFourierDataI;FourierDataR = NULL;}FourierDataR = new float[nRows*nCols];FourierDataI = new float[nRows*nCols];for(i=0;i<nRows;i++)for(j=0;j<nCols;j++){ //图像数据先给傅立叶变换数组FourierDataR[i*nCols+j] = (float) mSourceData[i*nCols+j];FourierDataI[i*nCols+j] = 0.0;}for (i=0;i<nRows;i++){ //每行进行傅立叶变换for (j=0;j<nCols;j++){Tr[j] = FourierDataR[i*nCols + j];Ti[j] = FourierDataI[i*nCols + j];}FFTandIFFT(Tr,Ti,nCols,1);for (j=0;j<nCols;j++){FourierDataR[i*nCols + j] = Tr[j];FourierDataI[i*nCols + j] = Ti[j];}}deleteTr;delete Ti;Tr = new float[nRows];Ti = new float[nRows];for(j=0;j<nCols;j++){ //每列进行傅立叶变换for (i=0;i<nRows;i++){Tr[i] = FourierDataR[i*nCols + j];Ti[i] = FourierDataI[i*nCols + j];}FFTandIFFT(Tr,Ti,nRows,1);for (i=0;i<nRows;i++){FourierDataR[i*nCols + j] = Tr[i];FourierDataI[i*nCols + j] = Ti[i];}}for (i=0;i<nRows;i++)for (j=0;j<nCols;j++){temp = sqrt(FourierDataR [i*nCols+j]*FourierDataR [i*nCols+j] +FourierDataI [i*nCols+j]*FourierDataI[i*nCols+j] );temp /= 100;if(temp > 255.0)temp = 255.0;ii = nRows - 1 - (i<nRows/2?i+nRows/2:i-nRows/2);jj = (j<nCols/2)?(j+nCols/2):(j-nCols/2);//将变换后现实的原点调整在中心位置pResultData[3*ii*nCols+3*jj] = (int) temp;pResultData[3*ii*nCols+3*jj+1] = (int) temp;pResultData[3*ii*nCols+3*jj+2] = (int) temp;}//更新显示UpdateAllViews(NULL);deleteFourierDataR;deleteFourierDataI;FourierDataI = NULL;FourierDataR = NULL;return;}voidCImageProcessingDoc::OnSaltPepperNoice(){// TODO: Add your command handler code here // TODO: Add your command handler code here int x;int y;Salt_Pepper_Noise(mSourceData,nCols,nRows);for(y = 0; y <nRows; y++)for(x = 0; x <nCols; x++){pSourceData[3*y*nCols+3*x] = (unsigned char) mSourceData[y*nCols+x];pSourceData[3*y*nCols+3*x+1] = (unsigned char) mSourceData[y*nCols+x];pSourceData[3*y*nCols+3*x+2] = (unsigned char) mSourceData[y*nCols+x];}UpdateAllViews(NULL);}voidCImageProcessingDoc::OnRandomNoise(){// TODO: Add your command handler code hereint x;int y;Random_Noise(mSourceData,nRows,nCols);for(y = 0; y <nRows; y++)for(x = 0; x <nCols; x++){pSourceData[3*y*nCols+3*x] = (unsigned char) mSourceData[y*nCols+x];pSourceData[3*y*nCols+3*x+1] = (unsigned char) mSourceData[y*nCols+x];pSourceData[3*y*nCols+3*x+2] = (unsigned char) mSourceData[y*nCols+x];}UpdateAllViews(NULL);}voidCImageProcessingDoc::Salt_Pepper_Noise(unsigned char *mdata, intnHeight, intnWidth){unsigned char* lpSrc;//循环变量longi;long j;//生成伪随机种子srand((unsigned)time(NULL));//在图像中加噪for (j = 0;j <nHeight ;j++){for(i = 0;i <nWidth ;i++){if(rand()>31500){// 指向源图像倒数第j行，第i个象素的指针lpSrc = (unsigned char *)&mdata[j*nWidth + i];//图像中当前点置为黑*lpSrc = 0;}}}// 返回return ;}voidCImageProcessingDoc::Random_Noise(unsigned char *mdata, intnHeight, intnWidth){// 指向源图像的指针unsigned char* lpSrc;//循环变量longi;long j;//像素值unsigned char pixel;//噪声BYTE NoisePoint;//生成伪随机种子srand((unsigned)time(NULL));//在图像中加噪for (j = 0;j <nHeight ;j++){for(i = 0;i <nWidth ;i++){NoisePoint=rand()/1024;// 指向源图像倒数第j行，第i个象素的指针lpSrc = (unsigned char *)&mdata[nWidth * j + i];//取得像素值pixel = (unsigned char)*lpSrc;*lpSrc = (unsigned char)(pixel*224/256 + NoisePoint);}}// 返回return ;}voidCImageProcessingDoc::MedianFiltering(unsigned char *sourcedata, unsigned char *resultdata,intnHeight, intnWidth, intnR){inti,j,m,n,r;unsignedtmp;unsigned char* mdata = new unsigned char[(2*nR+1)*(2*nR+1)];for (i=0;i<nRows;i++)for (j=0;j<nCols;j++){if((i<nR) || (i>nHeight-nR-1) || (j<nR) || (j>nWidth-nR-1)) resultdata[i*nWidth+j] = 0;else{for(m=-nR;m<=nR;m++)for(n=-nR;n<=nR;n++)mdata[(m+nR)*(2*nR+1)+n+nR]=sourcedata[(i+m)*nWidth+(j+n)];//排序for(m=0;m<(2*nR+1)*(2*nR+1)-2;m++){r = 1;for (n=m+1;n<(2*nR+1)*(2*nR+1);n++){if (mdata[n]<mdata[n+1]){tmp = mdata[n];mdata[n]=mdata[n+1];mdata[n+1] =tmp;r=0;}}if (r)break;}mResultData[i*nWidth+j] = mdata[nR*(2*nR+1)+nR];}}}voidCImageProcessingDoc::OnMedianFiltering(){// TODO: Add your command handler code hereint x;int y;MedianFiltering(mSourceData,mResultData,nRows,nCols,1);for(y = 0; y <nRows; y++)for(x = 0; x <nCols; x++){pResultData[3*y*nCols+3*x] = (unsigned char) mResultData[y*nCols+x];pResultData[3*y*nCols+3*x+1] = (unsigned char) mResultData[y*nCols+x];pResultData[3*y*nCols+3*x+2] = (unsigned char) mResultData[y*nCols+x];}UpdateAllViews(NULL);}voidCImageProcessingDoc::OnDct(){// TODO: Add your command handler code here}voidCImageProcessingDoc::OnFwt(){// TODO: Add your command handler code here}voidCImageProcessingDoc::OnDht(){// TODO: Add your command handler code here}voidCImageProcessingDoc::OnWaveletTransform(){// TODO: Add your command handler code here}voidCImageProcessingDoc::OnGreyAdjustifcation(){// TODO: Add your command handler code here}voidCImageProcessingDoc::OnGreyLinearAdjustifcation(){// TODO: Add your command handler code hereint x;int y;inttmp;CGreyRatiomdlg;mdlg.DoModal();for(y=0;y<nRows;y++)for(x=0;x<nCols;x++){tmp =(int)(mdlg.m_GreyRatio*mSourceData[y*nCols+x]);tmp = tmp>255?255:tmp;pResultData[3*y*nCols+3*x] = tmp;pResultData[3*y*nCols+3*x+1] = tmp;pResultData[3*y*nCols+3*x+2] = tmp;}UpdateAllViews(NULL);}voidCImageProcessingDoc::OnGreySeglinearAdjustifcation() {// TODO: Add your command handler code here}voidCImageProcessingDoc::On2dgrad(){// TODO: Add your command handler code hereint x;int y;int dx;intdy;inttmp;for(y=0;y<nRows-1;y++)for(x=0;x<nCols-1;x++){dx = mSourceData[y*nCols+x] - mSourceData[y*nCols+x+1];dy = mSourceData[y*nCols+x] - mSourceData[(y+1)*nCols+x];tmp = (int) sqrt(dx*dx+dy*dy);tmp = tmp>255?255:tmp;pResultData[3*y*nCols+3*x] = tmp;pResultData[3*y*nCols+3*x+1] = tmp;pResultData[3*y*nCols+3*x+2] = tmp;}UpdateAllViews(NULL);}voidCImageProcessingDoc::OnRobert(){// TODO: Add your command handler code hereint x;int y;int dx;intdy;inttmp;for(y=0;y<nRows-1;y++)for(x=0;x<nCols-1;x++){dx = mSourceData[y*nCols+x] - mSourceData[(y+1)*nCols+x+1];dy = mSourceData[y*nCols+x+1] - mSourceData[(y+1)*nCols+x];tmp = (int) sqrt(dx*dx+dy*dy);tmp = tmp>255?255:tmp;pResultData[3*y*nCols+3*x] = tmp;pResultData[3*y*nCols+3*x+1] = tmp;pResultData[3*y*nCols+3*x+2] = tmp;}UpdateAllViews(NULL);}voidCImageProcessingDoc::DCTandIDCT(float *Ff, int N, boolbDctIDct) {float *mR;float *mI;inti;float Ff0 = 0;mR = new float[N*2];mI = new float[N*2];if(bDctIDct){for(i=0;i<2*N;i++){if(i<N)mR[i] = Ff[i];elsemR[i] = 0;mI[i] = 0;}for(i=0;i<N;i++)Ff0 += Ff[i];Ff0 = Ff0/sqrt(N);FFTandIFFT(mR,mI,2*N,true);Ff[0] = Ff0;for(i=0;i<N;i++)Ff[i] = (mR[i]*cos(i*PI/(2*N)) + mI[i]*sin(i*PI/(2*N))) *sqrt(2.0/N);}else{for(i=0;i<2*N;i++){if(i<N){mR[i] = Ff[i]*cos(i*PI/(2*N));mI[i] = Ff[i]*sin(i*PI/(2*N));}else{mR[i] = 0;mI[i] = 0;}}for(i=0;i<N;i++)Ff0 += Ff[i];Ff0 = Ff0/sqrt(N);FFTandIFFT(mR,mI,2*N,false);for(i=0;i<N;i++)Ff[i] = 1/sqrt(N) - sqrt(2.0/N) + sqrt(2.0/N)*mR[i];}return;}结果截图;5 结论本文充分考虑到程序运行的时间复杂度和空间复杂度问题。