快速中值滤波及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为止，最后，将缓冲区的第二、三行存人新建的文件，程序流程框图如图1
4 中值滤波快速算法的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
#define PI (acos(0.0) * 2)
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
///////////////////////////////////////////////////////////////////// ////////
// CImageProcessingDoc
IMPLEMENT_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_MAP
END_MESSAGE_MAP()
///////////////////////////////////////////////////////////////////// ////////
// CImageProcessingDoc construction/destruction
CImageProcessingDoc::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 serialization
void CImageProcessingDoc::Serialize(CArchive& ar)
{
if (ar.IsStoring())
{
// TODO: add storing code here
}
else
{
// TODO: add loading code here
}
}
///////////////////////////////////////////////////////////////////// ////////
// CImageProcessingDoc diagnostics
#ifdef _DEBUG
void CImageProcessingDoc::AssertValid() const
{
CDocument::AssertValid();
}
void CImageProcessingDoc::Dump(CDumpContext& dc) const
{
CDocument::Dump(dc);
}
#endif //_DEBUG
///////////////////////////////////////////////////////////////////// ////////
// CImageProcessingDoc commands
BOOL 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++)
{
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;
}
void CImageProcessingDoc::OnHistogramAdjustifcation() {
// TODO: Add your command handler code here
int x,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 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,逻辑变量，非零做正变换，零做反变换 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) { //如果做傅立叶你变换，则必须对数列除以N for(l=0;l { Tr[l] /= N; Ti[l] /= N; } } //计算循环次数r r = 0; l = N; while(l /= 2) r++; //倒序 int LH = N/2; int i; float temp; j = 0; for (i=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 { 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; } } } } void CImageProcessingDoc::OnFft() { // TODO: Add your command handler code here int i,j; int ii,jj; float temp; float *Tr; float *Ti; Tr = new float[nCols]; Ti = new float[nCols]; if ( FourierDataR) { delete FourierDataR; FourierDataR = NULL; } if ( FourierDataI) { delete FourierDataI; FourierDataR = NULL; } FourierDataR = new float[nRows*nCols]; FourierDataI = new float[nRows*nCols]; for(i=0;i for(j=0;j { //图像数据先给傅立叶变换数组 FourierDataR[i*nCols+j] = (float) mSourceData[i*nCols+j]; FourierDataI[i*nCols+j] = 0.0; } for (i=0;i { //每行进行傅立叶变换 for (j=0;j { Tr[j] = FourierDataR[i*nCols + j]; Ti[j] = FourierDataI[i*nCols + j]; } FFTandIFFT(Tr,Ti,nCols,1); for (j=0;j { FourierDataR[i*nCols + j] = Tr[j]; FourierDataI[i*nCols + j] = Ti[j]; } } delete Tr; delete Ti; Tr = new float[nRows]; Ti = new float[nRows]; for(j=0;j { //每列进行傅立叶变换 for (i=0;i { Tr[i] = FourierDataR[i*nCols + j]; Ti[i] = FourierDataI[i*nCols + j]; } FFTandIFFT(Tr,Ti,nRows,1); for (i=0;i { FourierDataR[i*nCols + j] = Tr[i]; FourierDataI[i*nCols + j] = Ti[i]; } } for (i=0;i for (j=0;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 jj = (j //将变换后现实的原点调整在中心位置 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); 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++) { 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); } 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++) { 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); } 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 < nWidth ;i++) { if(rand()>31500) { // 指向源图像倒数第j行，第i个象素的指针 lpSrc = (unsigned char *)&mdata[j*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 *)&mdata[nWidth * 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;i for (j=0;j { if((i 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] { 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]; } } } 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++) { 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); } void CImageProcessingDoc::OnDct() { // TODO: Add your command handler code here } void CImageProcessingDoc::OnFwt() { // TODO: Add your command handler code here } void CImageProcessingDoc::OnDht() { // TODO: Add your command handler code here } void CImageProcessingDoc::OnWaveletTransform() { // TODO: Add your command handler code here } void CImageProcessingDoc::OnGreyAdjustifcation() { // TODO: Add your command handler code here } void CImageProcessingDoc::OnGreyLinearAdjustifcation() { // TODO: Add your command handler code here int x; int y; int tmp; CGreyRatio mdlg; mdlg.DoModal(); for(y=0;y for(x=0;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); } void CImageProcessingDoc::OnGreySeglinearAdjustifcation() { // TODO: Add your command handler code here } void CImageProcessingDoc::On2dgrad() { // TODO: Add your command handler code here int x; int y; int dx; int dy; int tmp; for(y=0;y for(x=0;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); } void CImageProcessingDoc::OnRobert() { // TODO: Add your command handler code here int x; int y; int dx; int dy; int tmp; for(y=0;y for(x=0;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); } void CImageProcessingDoc::DCTandIDCT(float *Ff, int N, bool bDctIDct) { float *mR; float *mI; int i; float Ff0 = 0; mR = new float[N*2]; mI = new float[N*2]; if(bDctIDct) { for(i=0;i<2*N;i++) { if(i mR[i] = Ff[i]; else mR[i] = 0; mI[i] = 0; } for(i=0;i Ff0 += Ff[i]; Ff0 = Ff0/sqrt(N); FFTandIFFT(mR,mI,2*N,true); Ff[0] = Ff0; for(i=0;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 { 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 Ff0 += Ff[i]; Ff0 = Ff0/sqrt(N); FFTandIFFT(mR,mI,2*N,false); for(i=0;i Ff[i] = 1/sqrt(N) - sqrt(2.0/N) + sqrt(2.0/N)*mR[i]; } return; } 结果截图; 5 结论 本文充分考虑到程序运行的时间复杂度和空间复杂度问题。解决了图象大而内存不足的问题。对中 值滤波的普通算法采用一般的常规算法，而对其快速算法采用走S型并且迭代的方法。采用对程序运行 计时的方法．对中值滤波的快速算法和普通算法进行精确的比较，结论可靠。【参考文献】 [1]王积分，张新荣．计算机图象识别IM]．北京：中国铁道出版社．1988 [2]徐金梧，杨德斌．徐科．Turbo C 实用大全[M]．北京：机械工业出版社，1997