快速中值滤波算法

快速中值滤波算法
————————————————————————————————作者：————————————————————————————————日期：
南昌大学实验报告
学生姓名：洪僡婕学号：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_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<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) { //如果做傅立叶逆变换，则必须对数列除以N
for(l=0;l<N;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<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 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<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 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<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 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<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 here
int 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 here
int 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;
}
六、实验数据
七、思考及体会
在设计算法编制程序的时候，我们充分考虑到程序运行的时间复杂度和空间复杂度问题，在解决问题的前提下，使算法尽量简单，使程序运行占有的空间尽量的小，这样来减少不必要的时问浪费和空间浪费，从而太大的提高程序执行的效率。

采用迭代，逐次逼近的方法得到本次的中值，在一行处理完毕后转人下一行也采用走S型的方法．这样除第一个窗口采用了一伏排序得到中值外，其它的窗口都利用上伏的窗口的象素删除无用的3个象素后再加人新的3个象素，利用迭代的方法得到本次窗口的中值．这样太大地提高了程序执行的效率。

本实验充分考虑到程序运行的时间复杂度和空间复杂度问题。

解决了图象大而内存不足的问题。

对中值滤波的普通算法采用一般的常规算法，而对其快速算法采用走S型并且迭代的方法。

采用对程序运行计时的方法．对中值滤波的快速算法和普通算法进行精确的比较，结论可靠。