project.cpp

基于小波的SAR斑点处理

//#include <windows.h>
#include "iostream.h"
#include "stdafx.h"
#include "project.h"
#include "math.h"

CProject::CProject()
{
	//PRJPARAM  PrjParam;

	//PrjParam.m_bRadius_a=TRUE; 
	//PrjParam.m_bRadius_b=TRUE;
	//PrjParam.m_bLatitude1=FALSE; 
	//PrjParam.m_bLatitude2=FALSE;  
	//PrjParam.m_bCentralLatitude=FALSE;  
	//PrjParam.m_bCentralLongitude=FALSE; 
	//PrjParam.m_bCentralAngle=FALSE;  
	//PrjParam.m_bSouthLatitude=FALSE;  
	//PrjParam.m_bZonewidth=FALSE;
	//PrjParam.m_bMapScale=FALSE;
}

// 获取e、e2、e4、e6、e8的值
E_value CProject::Get_e()
{                                         
	E_value   E;
	
	E.e2 = (a*a-b*b) / a / a;
	E.e  = sqrt(E.e2);
	E.e4 = E.e2*E.e2;
	E.e6 = E.e4*E.e2;
	E.e8 = E.e6*E.e2;
	
	return E;
}

// 获取纬圈半径
double CProject::r_latitude(double latitude)
{
	double r,B;
	E_value E;

	E = Get_e();	
	B = Dms_to_arc(latitude);
	r = a*cos(B) / sqrt(1.0-E.e2*sin(B)*sin(B));
	
	return(r);
}      

// 获取子午圈曲率半径
double CProject::M_longitude(double latitude)
{
	double B,M;
	E_value  E;

	E = Get_e();
	B = Dms_to_arc(latitude);
	M = a*(1.0-E.e2) / pow((1-E.e2*sin(B)*sin(B)),1.5);
	
	return(M);
}	               

// 获取卯酉圈曲率半径
double CProject::N_latitude(double latitude)
{
	double N,B;
    E_value  E;
	
	E = Get_e();	
	B = Dms_to_arc(latitude);
	N = a / sqrt(1-E.e2*sin(B)*sin(B));
	
	return(N);
}	               

//-------------------------------------------------------------------------------------
double CProject::U_latitude(double latitude)
{
	double U,l,A,B,f;

    E_value  E;
	E = Get_e();	
	l = Dms_to_arc(latitude);
	f = asin(E.e*sin(l))*180.0/PI;	//angle in degree.
	
	A = (45.0+Dms_to_degree(latitude)/2.0)*PI/180.0;
	B = (45.0+f/2.0)*PI/180.0;
	U = tan(A)/pow(tan(B),E.e); 

	return(U);
}	


//-----------------------------------------------------------------------------------------------	
Radius CProject::Get_R(double latitude)
{
    Radius P;
    double M,N;

	E_value  E;
	E = Get_e();	
    
	M = M_longitude(latitude);
	N = N_latitude(latitude);
		    
	P.R  = sqrt(M*N);
	P.Rm = (a+a+b)/3.0;
	P.Rv = pow(pow(a,2.0)*b,1.0/3.0);
	P.Rf = sqrt(pow(a,2.0)/2.0+pow(b,2.0)/4.0/E.e*log((1.0+E.e)/(1.0-E.e)));
	
	return P;
}

//=================================================
// 等距
double CProject::Isodistance_S(double latitude)
{
	double A0,B0,C0,D0,E0;
	double  e,e2,e4,e6,e8;
	double  S,B;
	
    E_value  E;
	E  = Get_e();	
	e  = E.e;
	e2 = E.e2;
	e4 = E.e4;
	e6 = E.e6;
	e8 = E.e8;
	
/*	
	A0=1.0050517783495;
	B0=5.0623821196504e-003;
	C0=1.0624539772572e-005;
	D0=2.0808250963498e-008;
	E0=3.8590896397417e-011;
*/	
	
	A0 = 1.0+3.0/4.0*e2+45.0/64.0*e4+175.0/256.0*e6+11025.0/16384.0*e8;
	B0 = 3.0/4.0*e2+15.0/16.0*e4+525.0/512.0*e6+2205.0/2048.0*e8;
	C0 = 15.0/64.0*e4+105.0/256.0*e6+2205.0/4096.0*e8;
	D0 = 35.0/512.0*e6+315.0/2048.0*e8;
	E0 = 315.0/16384.0*e8;
//	cout<<"Please intput a number of the latitude:"<<endl<<"B=";
//	cin>>l1;	                                 
	B = Dms_to_arc(latitude);	                              
 	S = a*(1.0-e2)*(A0*B-B0/2.0*sin(2.0*B)+C0/4.0*sin(4.0*B)-D0/6.0*sin(6.0*B)+E0/8.0*sin(8.0*B));

// 	cout<<"S="<<S<<"(m)"<<endl;
	return(S);
}

// 等角
double CProject::Isoangle_Q(double latitude)	
{
	double  Q,B;
	double e,e2;

    E_value  E;
	E=Get_e();	
	e2=E.e2;
 	e=E.e;
/*	e4=e2*e2;
	e6=e4*e2;
	e8=e6*e2;
	
	A0=1.0/2.0*e2+5.0/24.0*e4+3.0/32.0*e6+1399.0/53760.0*e8;
	B0=5.0/48.0*e4+7.0/80.0*e6+689.0/17920.0*e8;
	C0=13.0/480.0*e6+1363.0/17920.0*e8;
	D0=677.0/107520.0*e8;*/
	
//	cout<<"Please intput a number of the latitude:"<<endl<<"B=";
//	cin>>l1;	                                 

	B=Dms_to_arc(latitude);	                              
// 	Q=B-A0*sin(2.0*B)+B0*sin(4.0*B)-C0*sin(6.0*B)+D0*sin(8.0*B);
	Q=log(tan(PI/4.0+B/2.0)*pow((1.0-e*sin(B))/(1.0+e*sin(B)),e/2.0));

// 	cout<<"Q="<<Q<<endl;
	return(Q);
}

// 等积
double CProject::Isosqare_F(double latitude)
{
	double  f,f2;
	double  F,B;
    E_value  E;
	E=Get_e();	

	B=Dms_to_arc(latitude);
	f=E.e*sin(B);
	f2=f*f;
	
 	F=a*a*(1.0-E.e2)*(1.0/2.0*sin(B)/(1.0-f2)+1.0/4.0/E.e*log((1.0+f)/(1.0-f)))*1.0e-06;

	return(F);
}

//================================================
// 等距
double CProject::Isodistance_B(double S)
{	 	 		
	const double A1=2.5184647783e-003,B1=3.6998873e-006;
 	const double C1=7.4449e-008,D1=1.828e-011;
 	double l,B_latitude,B;
 	
// 	cout<<"Please input the length of a longitude:"<<endl<<"S=";
// 	cin>>s;
 	                                       				
 	l=S*1.570460641219e-007;
 	B_latitude=l+A1*sin(2.0*l)+B1*sin(4.0*l)+C1*sin(6.0*l)+D1*sin(8.0*l);
	B=Arc_to_dms(B_latitude);
	 	
//	cout<<"B="<<B<<"(dms)"<<endl;
	return(B);
}

// 等角
double CProject::Isoangle_B(double Q)
{
	const double PI=3.1415926535897932384;
	const double A1=3.3560695588e-003,B1=6.5700353e-006;
 	const double C1=1.76221e-007,D1=6.08e-011;
 	double l,B_latitude,B;
 	
// 	cout<<"Please input the length of a longitude:"<<endl<<"Q=";
// 	cin>>q;
 	                                       				
 	l=2.0*atan(exp(Q))-PI/2.0;
 	B_latitude=l+A1*sin(2.0*l)+B1*sin(4.0*l)+C1*sin(6.0*l)+D1*sin(8.0*l);
	B=Arc_to_dms(B_latitude);
	 	
// 	cout<<"B="<<B<<"(dms)"<<endl;
 	return(B);
}                    

// 等积
double CProject::Isosqare_B(double F)
{
	const double PI = 3.1415926535897932385;
	const double A1 = 2.2388876665e-003,B1=2.8823802e-006;
 	const double C1 = 5.07852e-009,D1=1.198e-011;
 	const double R2 = 40591120.14123;
//	double A1,B1,C1,D1;
 	double l,B_latitude,B;
 	
// 	A1=1.0/3.0*e2+31.0/180.0*e4+59.0/560.0*e6+126853.0/518400.0*e8;
//	B1=17.0/360.0*e4+61.0/1260.0*e6+3622447.0/94089600.0*e8;
//	C1=383.0/45360.0*e6+6688039.0/658627200.0*e8;
//	D1=27787.0/23522400.0*e8;

// 	cout<<"Please input the length of a longitude:"<<endl<<"Q=";
 //	cin>>f;
 	l = asin(F/R2);
 	B_latitude = l+A1*sin(2.0*l)+B1*sin(4.0*l)+C1*sin(6.0*l)+D1*sin(8.0*l);
	B = Arc_to_dms(B_latitude);
	 	
// 	cout<<"B="<<B<<"(dms)"<<endl;
  	return(B);
}

//++++++++++++++++++++++++++++++++++++++++++++++++

// 将度转化为度分秒
double Degree_to_dms(double angle)
{
	double degree,minute,second,fraction;
	double Dms;
	
	fraction=modf(angle,&degree);
	second=modf(fraction*60.0,&minute)*60.0;
	Dms=degree+minute/100.0+second/10000.0;
	
	return(Dms);
}

// 将度分秒转化为度
double Dms_to_degree(double angle)
{
	double degree,minute,fraction,second;
	double Deg;
	
	fraction=modf(angle+1.0e-08,&degree);
	second=modf(fraction*100.0,&minute)*100.0;	
	Deg=degree+minute/60.0+second/3600.0;
	
	return(Deg);
}

// 将度分秒转化为弧度
double Dms_to_arc(double angle)
{
	double degree,arc;

//	cout<<"Please input the number of an angle:"<<endl;
//	cin>>angle;
	
/*	fraction=modf(angle,&degree);                             
	second=modf(fraction*100,&minute)*100;*/
		                                                                                                       
	degree=Dms_to_degree(angle);
	arc=degree*PI/180.0;
	
//	cout<<"The arc number is :"<<arc<<endl; 	
	return(arc);
}
                                                                            
// 将弧度转化为度分秒
double  Arc_to_dms(double angle)
{
	double degree,Dms;
	
//	cout<<"Please input athe numer of an angle"<<endl;
//	cin>>angle;
	
	degree=angle*180.0/PI;
	Dms=Degree_to_dms(degree);
/*	fraction=modf(angle*180.0/PI,&degree);
	second=modf(fraction*60.0,&minute)*60.0;
		
	Dms=degree+minute/100.0+second/10000.0;
*/	
//	cout<<"The number of the angle in degree is :"<<Deg<<endl;
	return(Dms);
}

//================================================
/*  void main (void)
 	{
 	CProject  P(6378245.0,6356863.0);

 	E_value    E;
	E=P.Get_e();

	cout.precision(12);
 	cout<<"e2="<<E.e2<<endl<<endl;
 	cout<<"M="<<P.M_longitude(30.0)<<endl;
 	cout<<"N="<<P.N_latitude(30.0)<<endl;
 	cout<<"r="<<P.r_latitude(30.0)<<endl<<endl;

	Radius    radius;
	radius=P.Get_R(31.0);	
	cout<<"R="<<radius.R<<endl;
	cout<<"Rm="<<radius.Rm<<endl;
	cout<<"Rv="<<radius.Rv<<endl;           
	cout<<"Rf="<<radius.Rf<<endl;

/*
	cout.precision(12);
 	cout<<"U="<<P.U_latitude(30.0)<<endl;
 	cout<<"S="<<P.Isodistance_S(30.0)<<endl;
 	cout<<"Q="<<P.Isoangle_Q(30.0)<<endl;
 	cout<<"F="<<P.Isosqare_F(30.0)<<endl<<endl;
*/ 	
/*
	cout.precision(12);
 	cout<<"Bs="<<P.Isodistance_B(552895.344)<<endl;
 	cout<<"Bq="<<P.Isoangle_B(0.0867940561)<<endl;
 	cout<<"Bf="<<P.Isosqare_B(3522057.4488)<<endl<<endl;
	}  	*/