快速中值滤波算法

快速中值滤波算法————————————————————————————————作者：————————————————————————————————日期：南昌大学实验报告学生姓名：洪僡婕学号：6100411159 专业班级：数媒111班实验类型：■验证□综合□设计□创新实验日期： 4.29 实验成绩：一、实验项目名称数字图像处理二、实验目的实现快速中值滤波算法三、实验内容用VC++实现中值滤波的快速算法四、主要仪器设备及耗材PC机一台五、实验步骤// 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 serializationvoid CImageProcessingDoc::Serialize(CArchive& ar) {if (ar.IsStoring()) {// TODO: add storing code here}else{// TODO: add loading code here}}///////////////////////////////////////////////////////////////////////////// // CImageProcessingDoc diagnostics#ifdef _DEBUGvoid CImageProcessingDoc::AssertValid() const{CDocument::AssertValid();}void CImageProcessingDoc::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(unsignedchar))))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 hereint 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<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,逻辑变量，非零做正变换，零做反变换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;l<N;l++){Tr[l] /= N;Ti[l] /= N;}}//计算循环次数rr = 0; l = N;while(l /= 2) r++;//倒序int LH = N/2;int i;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;}}}}void CImageProcessingDoc::OnFft(){// TODO: Add your command handler code hereint 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<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];}}delete Tr;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);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 hereint 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 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);}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<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];}}}void CImageProcessingDoc::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);}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 hereint x;int y;int tmp;CGreyRatio mdlg;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);}void CImageProcessingDoc::OnGreySeglinearAdjustifcation() {// TODO: Add your command handler code here}void CImageProcessingDoc::On2dgrad() {// TODO: Add your command handler code hereint x;int y;int dx;int dy;int tmp;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);}void CImageProcessingDoc::OnRobert() {// TODO: Add your command handler code hereint x;int y;int dx;int dy;int tmp;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);}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<N)mR[i] = Ff[i];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,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;}六、实验数据七、思考及体会在设计算法编制程序的时候，我们充分考虑到程序运行的时间复杂度和空间复杂度问题，在解决问题的前提下，使算法尽量简单，使程序运行占有的空间尽量的小，这样来减少不必要的时问浪费和空间浪费，从而太大的提高程序执行的效率。