hardy_lineardlg.cpp

利用多面函数拟合法来拟合似大地水准面

// Hardy_LinearDlg.cpp : implementation file
//

#include "stdafx.h"
#include "GeoinModel.h"
#include "Hardy_LinearDlg.h"
//#include "MathBase.h"
#include "matrix.h" 

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif


////////////////////////////////
#ifndef _NO_NAMESPACE
using namespace std;
using namespace math;
#define STD std
#else
#define STD
#endif

#ifndef _NO_TEMPLATE
typedef matrix<double> Matrix;
#else
typedef matrix Matrix;
#endif

#ifndef _NO_EXCEPTION
#  define TRYBEGIN()	try {
#  define CATCHERROR()	} catch (const STD::exception& e) { \
						cerr << "Error: " << e.what() << endl; }
#else
#  define TRYBEGIN()
#  define CATCHERROR()
#endif
///////////////////////////////

extern double CenterL(double dl);
extern  double CenterB(double db);
extern  double ReturnCenterL(double dl);
extern  double ReturnCenterB(double db);
extern  double Arc2DMS(double dd);
extern  double int_radian(double intdeg);
extern  double deg_int(double gms);
extern double myu(double a, double b);

/////////////////////////////////////////////////////////////////////////////
// CHardy_LinearDlg dialog


CHardy_LinearDlg::CHardy_LinearDlg(CWnd* pParent /*=NULL*/)
	: CDialog(CHardy_LinearDlg::IDD, pParent)
{
	//{{AFX_DATA_INIT(CHardy_LinearDlg)
	m_strKnownFilePath = _T("");
	m_strOutputFilePath = _T("");
	m_strCoreFunction = _T("不顾及高程项");
	m_strB = _T("0.5");
	m_strP = _T("不定权");
	m_bOverwrite = true;
	m_dFlat = 0.01;
	//}}AFX_DATA_INIT
}


void CHardy_LinearDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CHardy_LinearDlg)
	DDX_Text(pDX, IDC_EDTKnownFilePath, m_strKnownFilePath);
	DDX_Text(pDX, IDC_EDTOutputFilePath, m_strOutputFilePath);
	DDX_CBString(pDX, IDC_COMBOFunction, m_strCoreFunction);
	DDX_CBString(pDX, IDC_COMBOB, m_strB);
	DDX_CBString(pDX, IDC_COMBOP, m_strP);
	DDX_Check(pDX, IDC_CHKOverWrite, m_bOverwrite);
	DDX_Text(pDX, IDC_EDTFlat, m_dFlat);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CHardy_LinearDlg, CDialog)
	//{{AFX_MSG_MAP(CHardy_LinearDlg)
	ON_BN_CLICKED(IDCalc, OnCalc)
	ON_BN_CLICKED(IDC_BTNOpenKnownFile, OnBTNOpenKnownFile)
	ON_BN_CLICKED(IDC_BTNOpenOutputPath, OnBTNOpenOutputPath)
	ON_EN_CHANGE(IDC_EDTFlat, OnChangeEDTFlat)
	ON_CBN_SELCHANGE(IDC_COMBOFunction, OnSelchangeCOMBOFunction)
	ON_CBN_SELCHANGE(IDC_COMBOP, OnSelchangeCombop)
	ON_CBN_SELCHANGE(IDC_COMBOB, OnSelchangeCombob)
	ON_BN_CLICKED(IDC_CHKOverWrite, OnCHKOverWrite)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CHardy_LinearDlg message handlers

//计算
void CHardy_LinearDlg::OnCalc() 
{
	int i,j;
	int iPtNum,iCodeNum;  //拟合点数，节点数
	double dB[300],dL[300],dH[300],dE[300];  //拟合点
	double dCodeB[200],dCodeL[200],dCodeH[200],dCodeE[200];  //节点
	
	double dTemp;	
	CString strPtname,strB,strL,strH,strh,strType,strEnd;


	if(m_strKnownFilePath=="")
	{
		MessageBox("请先指定已知数据文件。");
		return;
	}
	if(m_strOutputFilePath=="")
	{
		MessageBox("请先指定结果输出文件。");
		return;
	}

	//◎◎读入已知坐标文件
	char buffer[256];
	FILE *File=fopen(m_strKnownFilePath,"r");

	//读取已知坐标文件,点名,B,L,H,h,i
	iPtNum=0;iCodeNum=0;
	while(fgets(buffer,256,File)!=NULL)		
	{
		strPtname=strtok(buffer,",， ");			
		strB=strtok('\0',",， ");
		strL=strtok(  '\0',",， ");
		strH=strtok('\0',",， ");
		strh=strtok('\0',",， ");
		strType=strtok('\0',",， ");
		strEnd=strtok('\0',",， \r\n");
			
		if(atoi(strType)<2)  //0-拟和点，1-节点，2-不用
		{
			iPtNum++;	

			dTemp =atof(strB);       //B格式：DD.MMSSS
			dB[iPtNum]=int_radian(deg_int(dTemp));
			dB[iPtNum]=CenterB(dB[iPtNum]);    //中心化
			
			dTemp=atof(strL);          //L格式：DD.MMSSS
			dL[iPtNum]=int_radian(deg_int(dTemp));    //弧度
			dL[iPtNum]=CenterL(dL[iPtNum]);    //中心化

			dH[iPtNum]=atof(strH);  //大地高
			dE[iPtNum]=dH[iPtNum]-atof(strh);    //高程异常
		
			if (atoi(strType)==1)  //节点数据
			{
				iCodeNum++;
				dCodeB[iCodeNum]=dB[iPtNum];
				dCodeL[iCodeNum]=dL[iPtNum];
				dCodeH[iCodeNum]=dH[iPtNum];
				dCodeE[iCodeNum]=dE[iPtNum];
			}
		}
	}
	fclose(File);

	if(iPtNum<=0)
	{
		MessageBox("已知数据文件格式有误，请查看。");
		return;
	}
	if(iCodeNum<=0)
	{
		MessageBox("已知数据文件中没有节点（点类型＝1），无法进行多面函数拟合，请查看。");
		return;
	}

	//◎◎
	Matrix mQ,mA,mE,mP;
	mE.SetSize(iPtNum,1);     //高程异常
	mQ.SetSize(iPtNum,iCodeNum);
	mA.SetSize(iCodeNum,1);
	mP.SetSize(iPtNum,iPtNum);

	for(i=0;i<iPtNum;i++) mE(i,0)=dE[i+1];
    //Q矩阵，与核函数相关
	for(i=0;i<iPtNum;i++)
	{
		for(j=0;j<iCodeNum;j++)
		{
			if(m_strCoreFunction=="不顾及高程项")  //不顾及高程项
			{
				mQ(i,j)=sqrt(pow(dL[i+1]-dCodeL[j+1],2)+
					         pow(dB[i+1]-dCodeB[j+1],2)+
							 m_dFlat);
			}
			else  //顾及高程项
			{
				mQ(i,j)=sqrt(pow(dL[i+1]-dCodeL[j+1],2)+
					         pow(dB[i+1]-dCodeB[j+1],2)+
							 pow(dH[i+1]-dCodeH[j+1],2)+
							 m_dFlat);
			}

			//正双曲，或者倒双曲
			if(m_strB=="-0.5")	mQ(i,j)=1.0/mQ(i,j);
		}
	}
	//权矩阵
	double dSum;
	if(m_strP=="距离平方的倒数")
	{		
		for(i=0;i<iPtNum;i++)
		{
			dSum=0;
			for(j=0;j<iPtNum;j++)
			{
				mP(i,j)=pow(dL[i+1]-dL[j+1],2)+pow(dB[i+1]-dB[j+1],2);
				dSum+=mP(i,j);
			}
			dSum=dSum/iPtNum;
			for(j=0;j<iPtNum;j++)
			{
				mP(i,j)=mP(i,j)/dSum;
				if(mP(i,j)>0) mP(i,j)=1.0/mP(i,j);
			}
		}
		//
	}
	else if( m_strP=="距离的倒数")
	{
		for(i=0;i<iPtNum;i++)
		{
			for(j=0;j<iPtNum;j++)
			{
				mP(i,j)=sqrt((pow(dL[i+1]-dL[j+1],2)+pow(dB[i+1]-dB[j+1],2)));
				if(mP(i,j)>0) mP(i,j)=1.0/mP(i,j);
			}
		}
	}
	else if(m_strP== "距离的指数形式")
	{
		for(i=0;i<iPtNum;i++)
		{
			for(j=0;j<iPtNum;j++)
			{
				mP(i,j)=exp(-1.0*(pow(dL[i+1]-dL[j+1],2)+pow(dB[i+1]-dB[j+1],2))/0.02);  
				//还要除以数据点平均间距的两倍
			}
		}
	}

	//◎◎解算
	if (m_strP== "不定权")
	{	
		mA=!((~mQ)*mQ)*(~mQ)*mE;
	}
	else
	{
		mA=!((~mQ)*mP*mQ)*(~mQ)*mP*mE;
	}

	
	


	//◎◎输出结果
	FILE *fp;

/*
“rt”　　　　　　只读打开一个文本文件，只允许读数据 
“wt”　　　　　　只写打开或建立一个文本文件，只允许写数据
“at”　　　　　　追加打开一个文本文件，并在文件末尾写数据
“rb”　　　　　　只读打开一个二进制文件，只允许读数据
“wb”　　　　 　 只写打开或建立一个二进制文件，只允许写数据
“ab” 　　　　 　追加打开一个二进制文件，并在文件末尾写数据
“rt+”　　　　　 读写打开一个文本文件，允许读和写
“wt+”　　　　　 读写打开或建立一个文本文件，允许读写
“at+”　　　　　 读写打开一个文本文件，允许读，或在文件末追加数 据
“rb+”　　　　　 读写打开一个二进制文件，允许读和写 
“wb+”　　　　　 读写打开或建立一个二进制文件，允许读和写
“ab+” 　　　　　读写打开一个二进制文件，允许读，或在文件末追加数据
*/
	CString str;
	if(m_bOverwrite==0)
	{
		str="at";
	}
	else
	{
		str="wt";
	}
	if((fp=fopen(m_strOutputFilePath,str))!=NULL)
	{
		//起始标志（-9999.9999），节点个数
		//fprintf(fp,"-9999.9999,%s,%s\n",m_strCoreFunction,m_strP);
		fprintf(fp,"-9999.9999,%d\n",iCodeNum);
		//转换参数a，平滑参数q，曲面参数b（0.5/-0.5），大地经度，大地纬度，高程异常
		for(j=0;j<iCodeNum;j++) fprintf(fp,"%.15f,%.3f,%.1f,%.9f,%.9f,%.4f\n",
			                   mA(j,0),m_dFlat,atof(m_strB),
							   dCodeB[j+1],dCodeL[j+1],dCodeE[j+1]);
		fclose(fp);
		MessageBox("计算完成，请查看结果文件Hardy_LMS");
		//打开文本文件
		CString Command;
		Command = "notepad ";
		Command += m_strOutputFilePath;
		WinExec(Command,SW_SHOWMAXIMIZED);

	}
	else
	{
		MessageBox("写入结果文件出错，请检查。。。");
	}

}

void CHardy_LinearDlg::OnCancel() 
{	
	CDialog::OnCancel();
}

//打开已知文件
void CHardy_LinearDlg::OnBTNOpenKnownFile() 
{
	CString FileName;
	CFileDialog file(TRUE,NULL,NULL,OFN_HIDEREADONLY|OFN_OVERWRITEPROMPT,
		     TEXT("Chart Files (*.txt)|*.txt"),NULL);
	//打开文件
	if(file.DoModal()==IDOK)
	{
		m_strKnownFilePath=file.m_ofn.lpstrFile;
		UpdateData(false);
	}
	
}

//打开结果输出文件
void CHardy_LinearDlg::OnBTNOpenOutputPath() 
{	
	CString FileName;
	CFileDialog file(TRUE,NULL,NULL,OFN_HIDEREADONLY|OFN_OVERWRITEPROMPT,
		     TEXT("Chart Files (*.txt)|*.txt"),NULL);
	//打开文件
	if(file.DoModal()==IDOK)
	{
		m_strOutputFilePath=file.m_ofn.lpstrFile;
		UpdateData(false);
	}
	
}

void CHardy_LinearDlg::OnChangeEDTFlat() 
{
	UpdateData(true);	
}
void CHardy_LinearDlg::OnSelchangeCOMBOFunction() 
{
	UpdateData(true);	
}
void CHardy_LinearDlg::OnSelchangeCombop() 
{
	UpdateData(true);	
}
void CHardy_LinearDlg::OnSelchangeCombob() 
{
	UpdateData(true);
}
void CHardy_LinearDlg::OnCHKOverWrite() 
{
	UpdateData(true);
	
}