//********************************************************************************
 //    開頭            程序名:中值濾波 
                                                
 //    作者            王統金                                                                                                            
 //    目的            對位場數據進行中值濾波                                                               
 //    用法            填寫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 ;
}