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

快速中值滤波及c语言实现学生姓名： 刘 勇 学 号： 6100410218 专业班级： 数媒101 【摘要】本文讨论了用c语言在微机上实现中值滤波及快速算法，在程序设计的过程中充分考虑到程序运行的时间复杂度和空间复杂度的问题解决了由于图像太大而内存不够的问题，运用对程序运行时的方法，得出在PENTIUM-S100MHz上中值滤渡的一般算法运行423秒而快速算法运行2 58秒。【关键词】c语言；中值滤波；快速算法1 引言中值滤波是涂基发明的一种非线性信号处理技术，对抑制图像的噪声非常有效，在二维形式下，中值滤渡器是一个古有奇数个像素的滑动窗口，窗口正中的象素的灰度值用窗口内各个象素的中值代替窗口的中值为窗口中象素按大小顺序排列后处于中间位置的象素；本文讨论中值滤的一般算法并比较其运算速度。2 用C语言实现算法的若干问题在设计算法编制程序的时候，我们充分考虑到程序运行的时间复杂度和空间复杂度问题，在解决问题的前提下，使算法尽量简单，使程序运行占有的空间尽量的小，这样来减少不必要的时问浪费和空间浪费，从而太大的提高程序执行的效率。首先考虑到的内存问题。由于在本文算法中用的图像是512+512 8bit，这就存在一个内存不够大一整幅图像不能一次性调入的问题。为了解受此问题，可以只开辟一个3512的缓冲区n，将原图像采用分批调入缓冲区，使内存不够的问题得到了圆满的解决。另外为了对中值滤波的快速算法和普通算法进行精确的比较，采用对程序运行计时的方法，并精确计算每个算法运行的时间，使得出的结论更可靠。3 中值滤波算法的C语言程序实现本算法采用对开辟的3*512的缓冲区从左到右依次形成一个3*3的窗口然后将此3*3的窗口放人一个一维数组中，调用求中值子函数通过排序得出中值，当此中值不等于窗口中间位置的象素时用此中值来代替窗VI中间位置的象素灰度值若此缓冲区处理完毕后，将缓冲区的第一行存入新建的文件中，将第二、第三行分别向上移动一行，若存人新建的文件中的行数小于或等于511(即这样处理的行数小于或等于511)，则从原文件中调入一行作为缓冲区第三行，按上述方法进行直到处理的总行数等于511为止，最后，将缓冲区的第二、三行存人新建的文件，程序流程框图如图14 中值滤波快速算法的C语言程序实现本算法充分利用了上一次处理的结果采用迭代，逐次逼近的方法得到本次的中值，在一行处理完毕后转人下一行也采用走S型的方法这样除第一个窗口采用了一伏排序得到中值外，其它的窗口都利用上伏的窗口的象素删除无用的3个象素后再加人新的3个象素，利用迭代的方法得到本次窗口的中值这样太大地提高了程序执行的效率。41算法的解释首先是开辟一个3*512的缓冲区a，在初始化缓冲区时考虑到时间复杂度的问题，所以只初始化了第二、三行，而对第一行只初始化了前三个象素，这样便在缓冲区中可以得到一个3*3的窗口，对此窗口进行排序求中值后得出第一个窗口的中间象素的值将文件指针定位在2*512处。然后开始循环，当处理的行数小于或等于511时，将缓冲区a中的第二、三行分别向上移动一行变为第一、二行，从文件中读人512个字节作为缓冲区的第三行，并用行数模2的方法设置方向标志k当k为0时，从左向右移动窗口当k为1时从右向左移动窗口而每一窗口都利用上次的窗口的像素删除无用的3个象素后再加入新的3个象素利用迭代的方法从上次的中值得到本次的中值。当处理完一行后将缓冲区的第一行存入新建的文件中，最后将缓冲区的第二、三行存入文件中。42 算法代码/ ImageProcessingDoc.cpp : implementation of the CImageProcessingDoc class/#include stdafx.h#include ImageProcessing.h#include ImageProcessingDoc.h#include GreyRatio.h#include #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 here mImageFile = 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 serializationvoid CImageProcessingDoc:Serialize(CArchive& ar)if (ar.IsStoring()/ TODO: add storing code hereelse/ TODO: add loading code here/ CImageProcessingDoc diagnostics#ifdef _DEBUGvoid CImageProcessingDoc:AssertValid() constCDocument:AssertValid();void CImageProcessingDoc:Dump(CDumpContext& dc) constCDocument: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 here if(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(unsigned char),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+) pSourceData3*y*nCols+3*x = mSourceDatay*nCols+x; pSourceData3*y*nCols+3*x+1 = mSourceDatay*nCols+x; pSourceData3*y*nCols+3*x+2 = mSourceDatay*nCols+x; bFileIsLoad = TRUE; return TRUE;void CImageProcessingDoc:OnHistogramAdjustifcation() / TODO: Add your command handler code hereint x,y;double *mR;double *mS;mR = new double256;mS = new double256; for(x=0;x256;x+)mRx = mSx = 0.0;/统计直方图 for(y = 0; y nRows; y+) for(x = 0; x nCols; x+) mRmSourceDatay*nCols+x +; for(x=0;x256;x+) for(y=0;yx;y+) mSx += mRy; mSx /= nRows*nCols; /直方图变换 for(y = 0; y nRows; y+) for(x = 0; x nCols; x+) mResultDatay*nRows+x = (char) (255* mSmSourceDatay*nRows+x); /灰度计算 for(y = 0; y nRows; y+) for(x = 0; x nCols; x+) pResultData3*y*nCols+3*x = mResultDatay*nCols+x; pResultData3*y*nCols+3*x+1 = mResultDatay*nCols+x; pResultData3*y*nCols+3*x+2 = mResultDatay*nCols+x; /更新显示 UpdateAllViews(NULL);/ FFTandIFFT 一维傅立叶变换与你变换函数 / 输入时域数据实部Tr,虚部Ti/ 输出频域数据实部Tr,虚部Ti/ 序列长度N，N等于2的r次幂/ FFTorIFFT,逻辑变量，非零做正变换，零做反变换void CImageProcessingDoc:FFTandIFFT(float *Tr, float *Ti, int N, bool FFTorIFFT)int r; /迭代次数int l,j,k;/循环变量int p; /用于蝶形计算加权系数的指数int B; /对偶结点距离 float X,Y,XX,YY;float w;float cosw,sinw; if (!FFTorIFFT) /如果做傅立叶你变换，则必须对数列除以Nfor(l=0;lN;l+)Trl /= N;Til /= N; /计算循环次数rr = 0;l = N;while(l /= 2) r+; /倒序int LH = N/2; int i;float temp; j = 0;for (i=1;i=k)j = j-k;k = k/2;j = j + k;if (i=j) temp = Tri; Tri = Trj; Trj = temp; temp = Tii; Tii = Tij; Tij = 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;kN-1;k+=(1l) if (FFTorIFFT) /若做傅立叶正变换 cosw =cos(-w); sinw =sin(-w); else /傅立叶反变换 cosw =cos(w); sinw =sin(w); X = Trk + Trk+B*cosw - Tik+B * sinw; Y = Tik + Trk+B*sinw + Tik+B * cosw; XX = Trk - Trk+B*cosw + Tik+B * sinw; YY = Tik - Trk+B*sinw - Tik+B * cosw; Trk = X; Tik = Y; Trk+B = XX; Tik+B = YY; void CImageProcessingDoc:OnFft() / TODO: Add your command handler code here int i,j;int ii,jj; float temp;float *Tr;float *Ti;Tr = new floatnCols;Ti = new floatnCols; if ( FourierDataR) delete FourierDataR;FourierDataR = NULL; if ( FourierDataI) delete FourierDataI;FourierDataR = NULL; FourierDataR = new floatnRows*nCols;FourierDataI = new floatnRows*nCols; for(i=0;inRows;i+)for(j=0;jnCols;j+) /图像数据先给傅立叶变换数组 FourierDataRi*nCols+j = (float) mSourceDatai*nCols+j; FourierDataIi*nCols+j = 0.0; for (i=0;inRows;i+) /每行进行傅立叶变换 for (j=0;jnCols;j+)Trj = FourierDataRi*nCols + j; Tij = FourierDataIi*nCols + j;FFTandIFFT(Tr,Ti,nCols,1); for (j=0;jnCols;j+) FourierDataRi*nCols + j = Trj; FourierDataIi*nCols + j = Tij; delete Tr;delete Ti;Tr = new floatnRows;Ti = new floatnRows;for(j=0;jnCols;j+) /每列进行傅立叶变换 for (i=0;inRows;i+)Tri = FourierDataRi*nCols + j; Tii = FourierDataIi*nCols + j;FFTandIFFT(Tr,Ti,nRows,1); for (i=0;inRows;i+) FourierDataRi*nCols + j = Tri; FourierDataIi*nCols + j = Tii; for (i=0;inRows;i+)for (j=0;j 255.0)temp = 255.0; ii = nRows - 1 - (inRows/2?i+nRows/2:i-nRows/2);jj = (jnCols/2)?(j+nCols/2):(j-nCols/2); /将变换后现实的原点调整在中心位置 pResultData3*ii*nCols+3*jj = (int) temp; pResultData3*ii*nCols+3*jj+1 = (int) temp; pResultData3*ii*nCols+3*jj+2 = (int) temp; /更新显示 UpdateAllViews(NULL); delete FourierDataR;delete FourierDataI;FourierDataI = NULL;FourierDataR = NULL; return;void CImageProcessingDoc: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+) pSourceData3*y*nCols+3*x = (unsigned char) mSourceDatay*nCols+x; pSourceData3*y*nCols+3*x+1 = (unsigned char) mSourceDatay*nCols+x; pSourceData3*y*nCols+3*x+2 = (unsigned char) mSourceDatay*nCols+x; UpdateAllViews(NULL);void CImageProcessingDoc:OnRandomNoise() / TODO: Add your command handler code here int x;int y;Random_Noise(mSourceData,nRows,nCols);for(y = 0; y nRows; y+) for(x = 0; x nCols; x+) pSourceData3*y*nCols+3*x = (unsigned char) mSourceDatay*nCols+x; pSourceData3*y*nCols+3*x+1 = (unsigned char) mSourceDatay*nCols+x; pSourceData3*y*nCols+3*x+2 = (unsigned char) mSourceDatay*nCols+x; UpdateAllViews(NULL);void CImageProcessingDoc:Salt_Pepper_Noise(unsigned char *mdata, int nHeight, int nWidth)unsigned char*lpSrc;/循环变量long i;long j;/生成伪随机种子srand(unsigned)time(NULL);/在图像中加噪for (j = 0;j nHeight ;j+)for(i = 0;i 31500)/ 指向源图像倒数第j行，第i个象素的指针lpSrc = (unsigned char *)&mdataj*nWidth + i;/图像中当前点置为黑*lpSrc = 0;/ 返回return ;void CImageProcessingDoc:Random_Noise(unsigned char *mdata, int nHeight, int nWidth)/ 指向源图像的指针unsigned char*lpSrc;/循环变量long i;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 *)&mdatanWidth * j + i;/取得像素值pixel = (unsigned char)*lpSrc;*lpSrc = (unsigned char)(pixel*224/256 + NoisePoint);/ 返回return ;void CImageProcessingDoc:MedianFiltering(unsigned char *sourcedata, unsigned char *resultdata,int nHeight, int nWidth, int nR) int i,j,m,n,r; unsigned tmp; unsigned char* mdata = new unsigned char(2*nR+1)*(2*nR+1); for (i=0;inRows;i+) for (j=0;jnCols;j+) if(inHeight-nR-1) | (jnWidth-nR-1) resultdatai*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 (mdatanmdatan+1) tmp = mdatan;mdatan=mdatan+1;mdatan+1 =tmp;r=0; if (r) break; mResultDatai*nWidth+j = mdatanR*(2*nR+1)+nR; void CImageProcessingDoc:OnMedianFiltering() / TODO: Add your command handler code here int x;int y;MedianFiltering(mSourceData,mResultData,nRows,nCols,1);for(y = 0; y nRows; y+) for(x = 0; x nCols; x+) pResultData3*y*nCols+3*x = (unsigned char) mResultDatay*nCols+x; pResultData3*y*nCols+3*x+1 = (unsigned char) mResultDatay*nCols+x; pResultData3*y*nCols+3*x+2 = (unsigned char) mResultDatay*nCols+x; UpdateAllViews(NULL); void CImageProcessingDoc:OnDct() / TODO: Add your command handler code herevoid CImageProcessingDoc:OnFwt() / TODO: Add your command handler code herevoid CImageProcessingDoc:OnDht() / TODO: Add your command handler code herevoid CImageProcessingDoc:OnWaveletTransform() / TODO: Add your command handler code herevoid CImageProcessingDoc:OnGreyAdjustifcation() / TODO: Add your command handler code herevoid CImageProcessingDoc:OnGreyLinearAdjustifcation() / TODO: Add your command handler code hereint x;int y;int tmp; CGreyRatio mdlg;mdlg.DoModal(); for(y=0;ynRows;y+)for(x=0;x255?255:tmp;pResultData3*y*nCols+3*x = tmp;pResultData3*y*nCols+3*x+1 = tmp;pResultData3*y*nCols+3*x+2 = tmp; UpdateAllViews(NULL);void CImageProcessingDoc:OnGreySeglinearAdjustifcation() / TODO: Add your command handler code herevoid CImageProcessingDoc:On2dgrad() / TODO: Add your command handler code hereint x;int y;int dx;int dy;int tmp; f