中值滤波.cpp

中值滤波源程序可以滤除波形干扰有利于数据处理

//********************************************************************************
 //    开头            程序名:中值滤波 
                                                
 //    作者            王统金                                                                                                            
 //    目的            对位场数据进行中值滤波                                                               
 //    用法            填写par.txt文件，按照用户要求填写滤波参数，                                                                        
 //    文件格式        par.txt文件要读取数据文件sp.grd和输出结果数据spout.grd                                           
 //    限制            输入和输出数据格式均为*.grd格式。
                                        
 //    错误处理        程序中断，不再执行。                                                                                                        
 //    注释            见论文中对变量和标识符的注释                                                       
 //    注释            name：原始数据文件名    outname：输出数据文件名     //******************************************************************************
#include "fstream.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
void Midfiltr(float Data[], int M,int N, int Numx, int Numy);        //定义中值滤波子程序
void paiqi(float Data[],int m_TraceN,int m_PointN,float Arg);        //定义排齐子程序
void main()                                                          //main()主程序
{

	float *data,*Data,*Data2;
	int i,j,M,N,Numx,Numy;
	float Arg;
	char name[80],outname[80];
	char skipchar[100];
	ifstream fpar("midfiltr.par");
	fpar.seekg(40,ios::cur);
    fpar.getline(name,sizeof(name));
	fpar.seekg(40,ios::cur);
    fpar.getline(outname,sizeof(outname));
    fpar.seekg(40,ios::cur);
    fpar>>Numx;
	fpar.seekg(40,ios::cur);
    fpar>>Numy;	
	fpar.seekg(40,ios::cur);
    fpar>>Arg;
	fpar.close();
	ifstream fin(name);
	ofstream fout(outname);	
	for(i=0;i<9;i++)
	{
		fin>>skipchar;
	    fout<<skipchar;
		if(i%2==0)fout<<"\n";
		else fout<<"  ";
		if(i==1)N=atoi(skipchar);
		if(i==2)M=atoi(skipchar);
	}
    Arg=Arg/45*atan(1);  
	int Ar=fabs(N*tan(Arg));
	int m1=M+Ar;
	int m2=M+2*Ar;
	data=new float[N*M];
	Data=new float[N*m2];
	cout<<name<<"\n"<<outname<<"\n"<<M<<"\n"<<N<<"\n"<<Numx<<"\n"<<Numy<<"\n";
	for(i=0;i<M;i++)
	{
	  for(j=0;j<N;j++)
	  {
		 fin>>data[i*N+j];
		 
	  }
	}  

	for(i=0;i<m2;i++)
	{
	  for(j=0;j<N;j++)
	  {
		if(i>=Ar && i<m1)Data[i*N+j]=data[(i-Ar)*N+j];
		else Data[i*N+j]=0;		 
	  }
	}
	paiqi(Data,M,N,Arg);                                        //调用排齐子程序
	Midfiltr(Data,m2,N,Numx,Numy);                              //调用中值滤波子程序
	paiqi(Data,M,N,-Arg);	                                    //再次调用排齐子程序，进行反排齐（注意角度此时加了负号）
	for(i=Ar;i<m1;i++)
	{
	  for(j=0;j<N;j++)
	  {		 
		 fout<<Data[i*N+j]<<"  ";		 
	  }
	  fout<<"\n\n";
	}
	fin.close();
	fout.close();

	
}
void paiqi(float Data[],int M,int N,float Arg)                   //排齐子程序
{
	int i,j,temp;
	float step;	
	int Ar=abs(int(N*tan(Arg)));
	int m1=M+Ar;
	int m2=M+2*Ar;
	int start=0,T=1;
	float *A;
	A=new float[m2*N];
	if(tan(Arg)<0)
	{
		start=Ar;
		T=-1;
	}
	for(i=0;i<m2;i++)
	  for(j=0;j<N;j++)
	  {
		  A[i*N+j]=0;
	  }
	for(i=start;i<start+m1;i++)
	{		
		for(j=0;j<N;j++)
		{
			step=(float)tan(Arg)*j;			
	        temp=int(step)+T;			
			A[i*N+j]=Data[(i+temp)*N+j]-fabs(temp-step)*(Data[(i+temp)*N+j]-Data[(i+temp-T)*N+j]);		   

		}
		
	}
	for(i=0;i<m2;i++)
	  for(j=0;j<N;j++)
	  {
		  Data[i*N+j]=A[i*N+j];
	  }

	return;

}
void Midfiltr(float Data[], int m_TraceN,int m_PointN, int Numx, int Numy)             //中值滤波子程序
{
	int i,j,k,ki,kj,kx,ky,startx,starty;
	float Mid;
	float *A,*temp;
	A=new float[Numx*Numy];
	temp=new float[m_TraceN*m_PointN];
	int len;
    for(i=0;i<m_TraceN;i++)
	{
		if(i<(Numx-1)/2)
		{
		   startx=0;
		   kx=2*i+1;
		}
		else if((m_TraceN-i)<=(Numx-1)/2)
		{
		   kx=2*(m_TraceN-i-1)+1;
		   startx=m_TraceN-kx;
		   
		}
		else
		{
		   startx=i-(Numx-1)/2;
		   kx=Numx;
		}
			
		for(j=0;j<m_PointN;j++)
		{			
			if(j<(Numy-1)/2)
			{
			  starty=0;
			  ky=2*j+1;
			}
			else if((m_PointN-j)<=(Numy-1)/2)
			{
		      ky=2*(m_PointN-j-1)+1;
		      starty=m_PointN-ky;
		   
			}
		    else
			{
			   starty=j-(Numy-1)/2;
		       ky=Numy;
			}
			k=0;
			for(ki=startx;ki<startx+kx;ki++)
			{
			  for(kj=starty;kj<starty+ky;kj++)
			  {
		         A[k]=Data[ki*m_PointN+kj];	
				 k++;
			  }
			}
			len=kx*ky;
		    for(ki=0;ki<(len+1)/2+1;ki++)
			{
				for(kj=ki+1;kj<len;kj++)
				{
					if(A[ki]<A[kj])
					{
						Mid=A[kj];
						A[kj]=A[ki];
						A[ki]=Mid;
					}
				}
			}
		    temp[i*m_PointN+j]=(A[len/2]+A[len-1-(len/2)])/2;		
			
		}
	}
	for(i=0;i<m_TraceN*m_PointN;i++)
	{
		Data[i]=temp[i];
	}
	delete []temp;
	delete []A;
	return ;
}