控制网平差（C++实现）

最新推荐文章于 2024-04-24 21:12:34 发布
最新推荐文章于 2024-04-24 21:12:34 发布
阅读量3.5k 收藏 81
点赞数 6
分类专栏： 测绘程序设计 MFC 文章标签： c++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/shaoxian_Li/article/details/126436826
版权

目录

一、需求分析

二、系统设计

三、原理

1、坐标概算

2、最小二乘平差

 四、代码

类设计：

 CAngle:

CCTRLPOINT：

CMatrix:

CErrorElipseDrow：

ControlNet：

CPP文件：

CAngle.cpp

CMatrix.cpp

CErrorElipseDrow.cpp

ControlNet.cpp

五、界面设计

六、总结

一、需求分析

导线网平差计算步骤十分繁琐且计算量大，人工计算不仅耗时长而且还容易运算错误，并且对于一些较为复杂的控制网，随着未知数的个数增多，误差方程的系数阵和常数阵也随之增大，此时人工计算就很难进行下去了。

随着计算机技术的普及和发展，人们把计算机的强大算力运用到我们的测绘学科里边，推动了测绘行业的飞速发展。为解决人工平差计算耗时长，计算量大，计算准确度不高的问题，平差程序也就随之孕育而生。

麻雀虽小，五脏俱全，作为一个简易的平差小程序，它仍应该具有以下功能：

1、平差计算：即通过输入边角观测数据列立误差方程，计算待定点坐标；

2、精度评定：评定观测精度以及待定点的点位精度；

3、控制网和误差椭圆绘制：根据输入数据的边角关系，绘制相应的控制网，让使用者明确导线网的布设形状及点位误差的分布规律。

二、系统设计

三、原理

1、坐标概算

​​​​​x_{k}=x_{a}+d\times \cos t=x_{a}+d\times \cos {T_{ak}}

y_{k}=y_{a}+d\times \sin t=x_{a}+d\times \sin {T_{ak}}

其中k为待定点，a为与k有距离和角度观测关系的已知点，b为与a有角度观测的另一已知点。

T_{ak}的计算为：

 逐点解算法：选择已知点或已计算出坐标的点作为起算点，根据具体的图形，逐个推算出各点的坐标。

计算步骤：

2、最小二乘平差

 四、代码

类设计：

 CAngle:

角度类，用于角度的DEG\RAD\DMS转换

#pragma once

//枚举数据类型，用于代表角度形式
enum AngleStyle
{
	DEG,
	DMS,
	RAD
};

//角度类，可用于角度不同类型的转换

class CAngle
{
public:
	CAngle(double value = 0, AngleStyle style = DMS);
	~CAngle(void);
private:
	double dValue;//角度值
	AngleStyle  nCurStyle;//当前角度值类型
private:
	//设置常成员函数的作用：1.类成员不会被改变
	//2.可以被常类变量调用
	double Deg(double dDms) const;
	double Dms(double dDeg) const;

public:
	//获取指定的类型获取角度值，
	//由于返回的是dValue的引用，所以该值大小可以改变，即可以进行赋值
	double& operator() (AngleStyle style);

	//重载，获取指定的类型获取角度值，该值不可改变，const CAngle类型变量调用
	double operator() (AngleStyle style) const;
	//重载运算符+/-
	friend CAngle operator + (const CAngle& m1, const CAngle& m2);
	friend CAngle operator - (const CAngle& m1, const CAngle& m2);
};

CCTRLPOINT：

观测边、观测角、测站、控制点的数据结构设计类

数据类的设计关系到后续程序编写的难易程度，各位读者可以在此类的基础上加以优化，好的数据结构可以达到事半功倍的的效果。推荐使用指针进行数据操作，避免后续对坐标进行修改时出现部分修改的情形。

#pragma once
#include "Angle.h"
#include"CErrorElipseDrow.h"
//控制点类，主要包括观测角度类、观测边长类及控制点类
class CCTRLPOINT
{
public:
	double X, Y;//控制点坐标
	CString strID;//控制点编号
	bool type=FALSE;//点的类型，若为已知点则为true,若为未知点则为FALSE
	bool tag = FALSE;//标识符，如果已经经过坐标概算得到了近似坐标，则tag值为真
	double dE, dF;//误差椭圆的长短半轴
	double dAlfa;//误差椭圆长半轴方位角
	double dP;//点位误差
	CErrorElipseDrow EDRAW;//误差椭圆绘制
	//构造函数
	CCTRLPOINT()
	{
		X = Y = 0;
		strID = _T("");
		dE = dF = dAlfa = 0;
		dP = 0;
	}
	~CCTRLPOINT() {};//析构函数
};

//方向观测类
class CObsDirection
{
public:
	CObsDirection() { ObsAngle = 0; cpST = NULL; cpOB = NULL; m = 0; };
	~CObsDirection()
	{
		if (cpST != NULL)
			cpST = NULL;
		if (cpOB != NULL)
			cpOB = NULL;
	};
public:
	CCTRLPOINT* cpST;//测站点
	CCTRLPOINT* cpOB;//照准点
	CAngle ObsAngle;//方向观测角度
	double m;

};

//距离观测类
class CObsDist
{
public:
	CObsDist() { Dist = 0; cpBG = NULL; cpEND = NULL; m = 0; };
	~CObsDist()
	{
		if (cpBG != NULL)
			cpBG = NULL;
		if (cpEND != NULL)
			cpEND = NULL;
	};
public:
	CCTRLPOINT* cpBG;//起始点
	CCTRLPOINT* cpEND;//终点
	double Dist;//观测距离
	double m;//中误差
};

//测站类（储存站点信息，以及初始方向角）
class STATION
{
public:
	CCTRLPOINT* station;//观测测站
	CCTRLPOINT* start;//定向边所在的照准站
	CAngle st_;
	STATION()
	{
		station = NULL;
		start = NULL;
	};
	~STATION()
	{
		if (station != NULL)
			station = NULL;
		if (start != NULL)
			start = NULL;
	};
};

CMatrix:

矩阵类，包括矩阵的转置、乘法、求逆计算

#pragma once
//矩阵类，用于矩阵运算（加、减、乘、求逆）
class CMatrix
{
public:
	CMatrix(int row = 3, int col = 3);
	// copy constructor

	CMatrix(const CMatrix& m);

	~CMatrix(void);
private:
	double** dMatData;//保存矩阵元素数据的二维数组
	int iRow;//矩阵的行
	int iCol;//矩阵的列

public:
	int Row() const { return iRow; }//返回行
	int Col() const { return iCol; }//返回列
	void SetSize(int row, int col);//调整数组的大小，原有数据不变（未测试）
	void ReSize(int row, int col);//调整数组大小，清空数据
	double& operator () (int row, int col);//获取矩阵元素
	double  operator () (int row, int col) const;//重载获取矩阵元素函数，只有const对象能访问
	CMatrix& operator = (const CMatrix& m);

	//注意：友元函数并不是类自己的成员函数
	friend CMatrix operator + (const CMatrix& m1, const CMatrix& m2);
	friend CMatrix operator - (const CMatrix& m1, const CMatrix& m2);
	friend CMatrix operator * (const CMatrix& m1, const CMatrix& m2);
	friend CMatrix operator * (const double& num, const CMatrix& m1);
	friend CMatrix operator * (const CMatrix& m1, const double& num);

	friend CMatrix operator ~ (const CMatrix& m);//矩阵转置
	CMatrix Inv();//矩阵求逆
	void Unit();//生成单位矩阵
	void Zero();//生成零矩阵
	CString Look();//矩阵元素查看
};

CErrorElipseDrow：

误差椭圆绘制类

#pragma once
//控制网和误差椭圆绘制
class CErrorElipseDrow
{
private:
  
    double dOrgX, dOrgY;//误差椭圆中心位置，即控制点平面坐标
    double dScale;//绘图比例
public:
    void SETORG(double X, double Y);
    void SETScale(double scale);
    void Draw(CDC* pDC, CRect& rect,double dAlfa,double dE,double dF);//dAlfa为长半轴方位角、dE为短半径，dF为长半径
};

ControlNet：

程序数据处理类，核心

#pragma once
#include "CTRLPOINT.h"
#include"Matrix.h"
//程序功能主要实现类
class ControlNet
{
public:
	ControlNet();
	~ControlNet();
private:
	int kpt_num;//已知点个数
	CCTRLPOINT* KPT;//已知点数组
	int ukpt_num;//已知点个数
	CCTRLPOINT* UKPT;//已知点数组
	int obsDirect_count;//方向观测数
	CObsDirection* OBSDIRECT;//方向观测数组
	int obsDist_count;//距离观测数
	CObsDist* OBSDIST;//距离观测数组
	int station_num;//测站数
	STATION *STA;//观测站
	CMatrix B;//系数矩阵
	CMatrix X;//未知数矩阵
	CMatrix P;//权矩阵
	CMatrix L;//常数项矩阵
	CMatrix Nbb;//法方程矩阵
	CMatrix Qxx;//协因数矩阵
	CMatrix V;//观测值改正数矩阵
	CMatrix QF;//观测值协因数矩阵
	double m0;//单位权中误差
	

public:int times;//迭代次数
public:
	
	bool LOADFILE(CString FileName,CString&strIn);//文件导入函数
	void Adjust();//平差计算函数
	void Estimate();//精度评定
	int SPILTSTRINF(const CString strLine, char split, CStringArray& strArray);//字符串分割函数
	void MainDraw(CDC* pDC, CRect& rect);//主绘图函数
	void OutPut(CString& result);//结果输出函数
	
private:
	bool SETKPT(int size);//给已知点数赋值并开辟存储空间
	bool SETUNKPT(int size);//给未知点数赋值并开辟存储空间
	bool SETDIRECT(int size);//给方向观测数赋值并开辟存储空间
	bool SETDIST(int size);//给距离观测赋值并开辟存储空间
	CCTRLPOINT* SearchPointUsingStrID(CString&strID);//按点编号查找点
	int SearchUKPosUsingID(CString ID);//按点编号查找未知数点号
	int SearchSTPosUsingID(CString ID);//按点编号查找测站点序号
	double AzimuthCal(const double x1, const double y1, const double x2, const double y2);//坐标方位角计算
	double DistCal(const double x1, const double y1, const double x2, const double y2);//距离计算函数
    void MAX_(double& minx, double& miny, double& maxx, double& maxy);//计算绘图区范围（经测试可用）
	void DrawFrame(CDC* pDC, CRect& rect,double &orgx,double &orgy,double &scalex,double&scaley);//框架绘制函数
	void checkspace();//当程序再次运行时进行空间检查，避免出现内存未释放的情况
	void RoughCal();//坐标概算函数
};

CPP文件：

CAngle.cpp

#include "pch.h"
#include "CAngle.h"
#include "math.h"
const double EPSILON = 1.0E-12;
const double PI = 4.0 * atan(1.0);

//重载构造函数，有缺省值
CAngle::CAngle(double value, AngleStyle style)
{
	dValue = value;
	nCurStyle = style;
}

CAngle::~CAngle(void)
{
}
//重载（）函数
double& CAngle::operator() (AngleStyle style) //指定的类型获取角度值
{
	//double dAngleValue;
	if (style == DMS)
	{
		if (nCurStyle == DEG)
		{
			dValue = Dms(dValue);
		}
		else if (nCurStyle == RAD)
		{
			dValue = Dms(dValue * 180.0 / PI);
		}
		nCurStyle = DMS;

	}
	else if (style == DEG)
	{
		if (nCurStyle == DMS)
		{
			dValue = Deg(dValue);
		}
		else if (nCurStyle == RAD)
		{
			dValue = dValue * 180.0 / PI;
		}
		nCurStyle = DEG;
	}
	else
	{
		if (nCurStyle == DMS)
		{
			dValue = Deg(dValue) * PI / 180;
		}
		else if (nCurStyle == DEG)
		{
			dValue = dValue * PI / 180;
		}
		nCurStyle = RAD;
	}
	return dValue;
}
//重载（）函数，该函数是常函数，只能被常CAngle对象使用
double CAngle::operator() (AngleStyle style) const //指定的类型获取角度值
{
	double dAngleValue;
	if (style == DMS)
	{
		if (nCurStyle == DEG)
		{
			dAngleValue = Dms(dValue);
		}
		else if (nCurStyle == RAD)
		{
			dAngleValue = Dms(dValue * 180.0 / PI);
		}
		else
		{
			dAngleValue = dValue;
		}

	}
	else if (style == DEG)
	{
		if (nCurStyle == DMS)
		{
			dAngleValue = Deg(dValue);
		}
		else if (nCurStyle == RAD)
		{
			dAngleValue = dValue * 180.0 / PI;
		}
		else
		{
			dAngleValue = dValue;
		}
	}
	else
	{
		if (nCurStyle == DMS)
		{
			dAngleValue = Deg(dValue) * PI / 180;
		}
		else if (nCurStyle == DEG)
		{
			dAngleValue = dValue * PI / 180;
		}
		else
		{
			dAngleValue = dValue;
		}
	}
	return dAngleValue;
}


//私有成员，度分秒向十进制度转换
double CAngle::Deg(double dDms) const
{
	int iDeg, iMin;
	double dSec;

	iDeg = int(dDms + EPSILON);//度//加一个很小的数，以防止取整时的出错
	iMin = int((dDms - iDeg) * 100 + EPSILON);//分
	dSec = ((dDms - iDeg) * 100 - iMin) * 100;//秒

	return iDeg + (double)iMin / 60 + dSec / 3600;
}

//私有成员，十进制度向度分秒转换
double CAngle::Dms(double dDeg) const
{
	int iDeg, iMin;
	double dSec;
	double dTmp;

	iDeg = int(dDeg + EPSILON);//整数部分度
	dTmp = (dDeg - iDeg) * 60;//小数部分转换成分
	iMin = int(dTmp + EPSILON);//取分的整数部分
	dSec = (dTmp - iMin) * 60;//截取秒

	return iDeg + (double)iMin / 100 + dSec / 10000;
}

//友元重载+函数
CAngle operator + (const CAngle& m1, const CAngle& m2)
{
	CAngle addAngle(0, RAD);
	addAngle(RAD) = m1(RAD) + m2(RAD);
	return addAngle;
}
//友元重载-函数
CAngle operator - (const CAngle& m1, const CAngle& m2)
{
	CAngle subAngle(0, RAD);
	subAngle(RAD) = m1(RAD) - m2(RAD);
	return subAngle;
}

CMatrix.cpp

#include "pch.h"
#include "Matrix.h"
#include "math.h"

CMatrix::CMatrix(int row, int col)
{
	iRow = row;
	iCol = col;
	dMatData = new double* [row];

	for (int i = 0; i < row; i++)
	{
		dMatData[i] = new double[col];
		for (int j = 0; j < col; j++)
		{
			dMatData[i][j] = 0;
		}
	}
}

// copy constructor，
//拷贝构造函数的作用：
//（1）以类对象作为函数参数传值调用时;
//（2）函数返回值为类对象；
//（3）用一个已定义的对象去初始化一个新对象时；

CMatrix::CMatrix(const CMatrix& m)
{
	iRow = m.Row();
	iCol = m.Col();
	dMatData = new double* [iRow];

	for (int i = 0; i < iRow; i++)
	{
		dMatData[i] = new double[iCol];
		//	for(int j=0;j<iCol;j++)
		{
			memcpy(dMatData[i], m.dMatData[i], sizeof(double) * iCol);
		}
	}

}

CMatrix::~CMatrix(void)
{
	for (int i = 0; i < iRow; i++)
	{
		delete[] dMatData[i];
	}
	delete[] dMatData;
}

//返回数组元素（引用返回）
double& CMatrix::operator () (int row, int col)
{
	if (row >= iRow || col >= iCol)
	{
		throw("CMatrix::operator(): Index out of range!");
	}

	return dMatData[row][col];
}

返回数组元素（重载）
double CMatrix::operator () (int row, int col) const
{
	if (row >= iRow || col >= iCol)
	{
		throw("CMatrix::operator(): Index out of range!");
	}

	return dMatData[row][col];
}


//重载预算符+
CMatrix operator + (const CMatrix& m1, const CMatrix& m2)
{

	if ((m1.Col() != m2.Col()) || (m1.Row() != m2.Row()))
	{
		throw("CMatrix::operator+: The two matrix have different size!");
	}

	CMatrix matTmp(m1.Row(), m1.Col());
	for (int i = 0; i < m1.Row(); i++)
	{
		for (int j = 0; j < m1.Col(); j++)
		{
			matTmp(i, j) = m1(i, j) + m2(i, j);
		}
	}
	return matTmp;
}

//重载赋值运算符=，当左右两边矩阵的大小不相等时，
//以右边的大小为基准,调整左边矩阵的大小

CMatrix& CMatrix::operator = (const CMatrix& m)
{
	//revised in 2011-4-1, by Daiwujiao
 //   if(iRow!=m.Row()||iCol!=m.Col())
	//{
 //       throw( "CMatrix::operator=: The two matrix have different size!");
	//}

	if (iRow != m.Row() || iCol != m.Col())
	{
		SetSize(m.Row(), m.Col());
	}
	for (int i = 0; i < iRow; i++)
	{

		for (int j = 0; j < iCol; j++)
		{
			dMatData[i][j] = m(i, j);
		}
	}
	return *this;
}


//调整矩阵大小，原有值不变
void CMatrix::SetSize(int row, int col)
{
	if (row == iRow && col == iCol)
	{
		return;
	}

	double** rsData = new double* [row];
	for (int i = 0; i < row; i++)
	{
		rsData[i] = new double[col];
		for (int j = 0; j < col; j++)
		{
			rsData[i][j] = 0;
		}
	}

	int minRow = (iRow > row) ? row : iRow;
	int minCol = (iCol > col) ? col : iCol;
	int  colSize = minCol * sizeof(double);


	for (int i = 0; i < minRow; i++)
	{
		memcpy(rsData[i], dMatData[i], colSize);
	}

	for (int i = 0; i < minRow; i++)
	{
		delete[] dMatData[i];
	}
	delete[] dMatData;
	dMatData = rsData;
	iRow = row;
	iCol = col;
	return;
}

void CMatrix::ReSize(int row, int col)
{
	for (int i = 0; i < iRow; i++)
	{
		delete[] dMatData[i];
	}
	delete[] dMatData;
	double** rsData = new double* [row];
	for (int i = 0; i < row; i++)
	{
		rsData[i] = new double[col];
		for (int j = 0; j < col; j++)
		{
			rsData[i][j] = 0;
		}
	}
	dMatData = rsData;
	iRow = row;
	iCol = col;
	return;
}



//重载预算符-
CMatrix operator - (const CMatrix& m1, const CMatrix& m2)
{

	if ((m1.Col() != m2.Col()) || (m1.Row() != m2.Row()))
	{
		throw("CMatrix::operator-: The two matrix have different size!");
	}

	CMatrix matTmp(m1.Row(), m1.Col());
	for (int i = 0; i < m1.Row(); i++)
	{
		for (int j = 0; j < m1.Col(); j++)
		{
			matTmp(i, j) = m1(i, j) - m2(i, j);
		}
	}
	return matTmp;
}

//重载预算符*,两个矩阵相乘，m1的列要等于m2的行
CMatrix operator * (const CMatrix& m1, const CMatrix& m2)
{

	if ((m1.Col() != m2.Row()))
	{
		throw("CMatrix::operator*: The col of matrix m1 doesn't equ to row of m2 !");
	}

	CMatrix matTmp(m1.Row(), m2.Col());
	for (int i = 0; i < m1.Row(); i++)
	{
		for (int j = 0; j < m2.Col(); j++)
		{
			for (int k = 0; k < m2.Row(); k++)
			{
				matTmp(i, j) += m1(i, k) * m2(k, j);
			}
		}
	}
	return matTmp;
}

//重载预算符*,矩阵右乘一个数
CMatrix operator * (const CMatrix& m1, const double& num)
{
	CMatrix matTmp(m1.Row(), m1.Col());
	for (int i = 0; i < m1.Row(); i++)
	{
		for (int j = 0; j < m1.Col(); j++)
		{
			matTmp(i, j) = m1(i, j) * num;

		}
	}
	return matTmp;
}

//重载预算符*,矩阵左乘一个数
CMatrix operator * (const double& num, const CMatrix& m1)
{
	CMatrix matTmp(m1.Row(), m1.Col());
	for (int i = 0; i < m1.Row(); i++)
	{
		for (int j = 0; j < m1.Col(); j++)
		{
			matTmp(i, j) = m1(i, j) * num;
		}
	}
	return matTmp;
}

//矩阵转置
CMatrix operator ~ (const CMatrix& m)
{
	CMatrix matTmp(m.Col(), m.Row());

	for (int i = 0; i < m.Row(); i++)
		for (int j = 0; j < m.Col(); j++)
		{
			matTmp(j, i) = m(i, j);
		}
	return matTmp;
}


//矩阵求逆
//采用选全主元法
CMatrix CMatrix::Inv()
{
	if (iRow != iCol)
	{
		throw("待求逆的矩阵行列不相等!");
	}

	int i, j, k, vv;

	CMatrix InvMat(iRow, iRow);

	//复制矩阵
	InvMat = *this;

	int* MainRow = new int[iRow];
	int* MainCol = new int[iRow];//用于记录主元素的行和列

	double dMainCell;//主元元素的值
	double dTemp;//临时变量

	for (k = 0; k < iRow; k++)
	{
		dMainCell = 0;
		//选全主元
		for (i = k; i < iRow; i++)
		{
			for (j = k; j < iRow; j++)
			{
				dTemp = fabs(InvMat(i, j));
				if (dTemp > dMainCell)
				{
					dMainCell = dTemp;
					MainRow[k] = i;
					MainCol[k] = j;
				}
			}
		}

		if (fabs(dMainCell) < 0.0000000000001)//矩阵秩亏，不能求逆
		{
			throw("矩阵秩亏");
		}

		if (MainRow[k] != k)//交换行
		{
			for (j = 0; j < iRow; j++)
			{
				vv = MainRow[k];
				dTemp = InvMat(k, j);
				InvMat(k, j) = InvMat(vv, j);
				InvMat(vv, j) = dTemp;
			}
		}

		if (MainCol[k] != k)//交换列
		{
			for (i = 0; i < iRow; i++)
			{
				vv = MainCol[k];
				dTemp = InvMat(i, k);
				InvMat(i, k) = InvMat(i, vv);
				InvMat(i, vv) = dTemp;
			}
		}
		InvMat(k, k) = 1.0 / InvMat(k, k);//计算乘数

		for (j = 0; j < iRow; j++) //计算主行
		{
			if (j != k)
			{
				InvMat(k, j) = InvMat(k, j) * InvMat(k, k);
			}
		}
		for (i = 0; i < iRow; i++)//消元
		{
			if (i != k)
			{
				for (j = 0; j < iRow; j++)
				{
					if (j != k)
					{
						InvMat(i, j) -= InvMat(i, k) * InvMat(k, j);
					}
				}
			}
		}
		for (i = 0; i < iRow; i++)//计算主列
		{
			if (i != k)
			{
				InvMat(i, k) = -InvMat(i, k) * InvMat(k, k);
			}
		}
	}

	for (k = iRow - 1; k >= 0; k--)
	{
		if (MainCol[k] != k)// 交换行
		{
			for (j = 0; j < iRow; j++)
			{
				vv = MainCol[k];
				dTemp = InvMat(k, j);
				InvMat(k, j) = InvMat(vv, j);
				InvMat(vv, j) = dTemp;
			}
		}

		if (MainRow[k] != k)//交换列
		{
			for (i = 0; i < iRow; i++)
			{
				vv = MainRow[k];
				dTemp = InvMat(i, k);
				InvMat(i, k) = InvMat(i, vv);
				InvMat(i, vv) = dTemp;
			}
		}
	}
	delete[] MainRow;
	delete[] MainCol;
	return InvMat;
}
//单位化矩阵
void CMatrix::Unit()
{
	for (int i = 0; i < iRow; i++)
	{
		for (int j = 0; j < iCol; j++)
		{
			dMatData[i][j] = (i == j) ? 1 : 0;
		}
	}
}

void CMatrix::Zero()
{
	for (int i = 0; i < iRow; i++)
	{
		for (int j = 0; j < iCol; j++)
		{
			dMatData[i][j] = 0;
		}
	}
}

CString CMatrix::Look()
{
	CString result, temp;
	result.Empty();
	for (int i = 0; i < iRow; i++)
	{
		for (int j = 0; j < iCol; j++)
		{
			temp.Empty();
			temp.Format(_T("%.4f\t"), dMatData[i][j]);
			result += temp;
			if (j == iCol - 1)
				result += "\r\n";
		}
	}
	return result;
}

CErrorElipseDrow.cpp

#include "pch.h"
#include "CErrorElipseDrow.h"
const double  PI = 3.1415926535;

//设置椭圆中心坐标
void CErrorElipseDrow::SETORG(double X, double Y)
{
	dOrgX = X;
	dOrgY = Y;
}
//设置椭圆绘图比例
void CErrorElipseDrow::SETScale(double scale)
{
	dScale = scale;
}

void CErrorElipseDrow::Draw(CDC* pDC, CRect& rect, double dAlfa, double dE, double dF)
{
	double dStartX, dStartY, dEndX, dEndY;
	//绘制短半轴
	dAlfa += PI / 2;
	dStartX = (dF * sin(dAlfa)) * dScale + dOrgX;
	dStartY = (-dF * cos(dAlfa)) * dScale + dOrgY;
	dEndX = (-dF * sin(dAlfa)) * dScale + dOrgX;
	dEndY = (dF * cos(dAlfa)) * dScale + dOrgY;
	CPen pen(PS_SOLID, 0, RGB(215, 0, 64));
	CPen* pOldPen = pDC->SelectObject(&pen);
	pDC->MoveTo(dStartX, dStartY);
	pDC->LineTo(dEndX, dEndY);
	//绘制长半轴
	dStartX = (-dE * cos(dAlfa)) * dScale + dOrgX;
	dStartY = (-dE * sin(dAlfa)) * dScale + dOrgY;
	dEndX = (dE * cos(dAlfa)) * dScale + dOrgX;
	dEndY = (dE * sin(dAlfa)) * dScale + dOrgY;
	pDC->MoveTo(dStartX, dStartY);
	pDC->LineTo(dEndX, dEndY);

	double ex, fy;
	ex = dE;
	fy = 0;
	//转换到长半轴方向上
	dStartX = (ex * cos(dAlfa)- fy * sin(dAlfa))* dScale + dOrgX;
	dStartY = (fy * cos(dAlfa)+ ex * sin(dAlfa)) * dScale + dOrgY;
	pDC->MoveTo(dStartX, dStartY);

	for (int i = 30; i <= 360; i += 30)
	{
		//在坐标轴方向的坐标
		ex = dE * cos((i / 180.0) * PI);
		fy = dF * sin((i / 180.0) * PI);

		//转换到长半轴方向上
		dEndX = (ex * cos(dAlfa)- fy * sin(dAlfa)) * dScale + dOrgX;
		dEndY = (fy * cos(dAlfa)+ ex * sin(dAlfa)) * dScale + dOrgY;
		pDC->LineTo(dEndX, dEndY);

	}
	pDC->SelectObject(pOldPen);
	pen.DeleteObject();
}

ControlNet.cpp

控制网平差主要实现类，代码很多，结合头文件来看，理解每个函数的作用

#include "pch.h"
#include "ControlNet.h"
#include"Angle.h"
const double Ro = 206264.806247;
const double EPSION = 1E-5;
const double PI = acos(-1);

ControlNet::ControlNet()
{
	kpt_num = ukpt_num = obsDirect_count = obsDist_count = station_num = times=0; 
	m0 = 0;
	KPT = NULL;
	UKPT = NULL;
	OBSDIRECT = NULL;
	OBSDIST = NULL;
	STA = NULL;
}

ControlNet::~ControlNet()
{
	if (KPT != NULL)
	{
		delete[]KPT;
		KPT = NULL;
	}
	if (UKPT != NULL)
	{
		delete[]UKPT;
		UKPT = NULL;
	}
	if (OBSDIRECT != NULL)
	{
		delete[]OBSDIRECT;
		OBSDIRECT = NULL;
	}
	if (OBSDIST != NULL)
	{
		delete[]OBSDIST;
		OBSDIST = NULL;
	}

}

bool ControlNet::SETKPT(int size)
{
	
		kpt_num = size;
		KPT = new CCTRLPOINT[kpt_num];
		return true;
	
}

bool ControlNet::SETUNKPT(int size)
{
	
	
		ukpt_num = size;
		UKPT = new CCTRLPOINT[ukpt_num];
		return true;
	
}

bool ControlNet::SETDIRECT(int size)
{
	
		obsDirect_count = size;
		OBSDIRECT = new CObsDirection[obsDirect_count];
		return true;
	
}

bool ControlNet::SETDIST(int size)
{
	
	
		obsDist_count = size;
		OBSDIST = new CObsDist[obsDist_count];
		return true;
	
}

//文件读取
bool ControlNet::LOADFILE(CString FileName, CString& strIn)
{
	CStdioFile sf;
	if (!sf.Open(FileName, CFile::modeRead))
		return false;
	checkspace();//再次运行程序时由于析构函数还未调用，所以需要释放内存空间，读入文件后重新开辟
	
	//读取已知点数据
	CString strLine;
	bool beof = sf.ReadString(strLine);
	strIn += strLine;
	strIn += "\r\n";
	//
	CString abc = _T("ABC2");
	int a = _ttoi(abc);//对于非数字，利用_ttoi会变为0
	//
	//数据格式检验
	if (_ttoi(strLine) == 0)
	{
		return false;
	}
	else

	{
		SETKPT(_ttoi(strLine));
		for (int i = 0; i < kpt_num; i++)
		{
			beof = sf.ReadString(strLine);
			strIn += strLine;
			strIn += "\r\n";
			CStringArray TEMP;
			int n = SPILTSTRINF(strLine, ',', TEMP);
			KPT[i].strID = TEMP[0];
			KPT[i].X = _ttof(TEMP[1]);
			KPT[i].Y = _ttof(TEMP[2]);
			KPT[i].tag = TRUE;
			KPT[i].type = TRUE;
		}
	}
		//读取未知点数据
	beof = sf.ReadString(strLine);
	strIn += strLine;
	strIn += "\r\n";
	if (_ttoi(strLine) == 0)
	{
		return false;
	}
	else
	{
		
		SETUNKPT(_ttoi(strLine));
		beof = sf.ReadString(strLine);
		strIn += strLine;
		strIn += "\r\n";
		CStringArray TEMP;
		int n = SPILTSTRINF(strLine, ',', TEMP);
		for (int i = 0; i < ukpt_num; i++)
		{
			UKPT[i].strID = TEMP[i];
			UKPT[i].tag = FALSE;
			UKPT[i].type = FALSE;
		}
	}
	//读取距离观测数据
	beof = sf.ReadString(strLine);
	if (_ttoi(strLine) == 0)
	{
		return false;
	}
	else
	{
		SETDIST(_ttoi(strLine));
		for (int i = 0; i < obsDist_count; i++)
		{
			beof = sf.ReadString(strLine);
			strIn += strLine;
			strIn += "\r\n";
			CStringArray TEMP;
			int n = SPILTSTRINF(strLine, ',', TEMP);
			OBSDIST[i].cpBG = SearchPointUsingStrID(TEMP[0]);
			OBSDIST[i].cpEND = SearchPointUsingStrID(TEMP[1]);
			OBSDIST[i].Dist = _ttof(TEMP[2]);
		}
	}
	//读取方位角观测数据
	beof = sf.ReadString(strLine);
	strIn += strLine;
	strIn += "\r\n";
	if (_ttoi(strLine) == 0)
	{
		return false;
	}
	else
	{
		
		SETDIRECT(_ttoi(strLine));
		station_num = 0;
		for (int i = 0; i < obsDirect_count; i++)
		{
			beof = sf.ReadString(strLine);
			strIn += strLine;
			strIn += "\r\n";
			CStringArray TEMP;
			int n = SPILTSTRINF(strLine, ',', TEMP);
			OBSDIRECT[i].cpST = SearchPointUsingStrID(TEMP[0]);
			OBSDIRECT[i].cpOB = SearchPointUsingStrID(TEMP[1]);
			OBSDIRECT[i].ObsAngle(DMS) = _ttof(TEMP[2]);
			if (OBSDIRECT[i].ObsAngle(RAD) == 0)
				station_num++;
		}
	}
	sf.Close();
	return true;
}

//字符串分割
int ControlNet::SPILTSTRINF(const CString strLine, char split, CStringArray& strArray)
{
	strArray.RemoveAll();//自带清空属性
	CString temp = strLine;
	int tag = 0;
	while (1)
	{
		tag = temp.Find(split);
		if (tag >= 0)
		{
			strArray.Add(temp.Left(tag));
			temp = temp.Right(temp.GetLength() - tag - 1);
		}
		else { break; }
	}
	strArray.Add(temp);
	return strArray.GetSize();
}

//按点编号查找点
CCTRLPOINT* ControlNet::SearchPointUsingStrID(CString& strID)
{
	for (int i = 0; i < kpt_num; i++)
	{
		if (strID == KPT[i].strID)
		{
			return &KPT[i];
		}
	}
	for (int i = 0; i < ukpt_num; i++)
	{
		if (strID == UKPT[i].strID)
		{
			return &UKPT[i];
		}
	}
	return NULL;
}
//按编号查找点的序号
int ControlNet::SearchUKPosUsingID(CString ID)
{
	for (int i = 0; i < ukpt_num; i++)
	{
		if (UKPT[i].strID == ID)
			return i;
	}
	return -1;
}

int ControlNet::SearchSTPosUsingID(CString ID)
{
	for (int i = 0; i < station_num; i++)
	{
		if (STA[i].station->strID == ID)
			return i;
	}
	return -1;
}

//方位角计算函数
double ControlNet::AzimuthCal(const double x1, const double y1, const double x2, const double y2)
{
	double dx, dy;
	dx = x2 - x1;
	dy = y2 - y1 + EPSION;
	if (dy >= 0)
		return PI - PI / 2 - atan(dx / dy);
	else
		return PI + PI / 2 - atan(dx / dy);
}
//距离计算函数
double ControlNet::DistCal(const double x1, const double y1, const double x2, const double y2)
{
	return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}
//控制网区域范围
void ControlNet::MAX_(double& minx, double& miny, double& maxx, double& maxy)
{
	minx = maxx = KPT[0].X;
	miny = miny = KPT[0].Y;

	for (int i = 0; i < ukpt_num + kpt_num; i++)
	{
		if (i < ukpt_num)
		{
			if (minx > UKPT[i].X)
				minx = UKPT[i].X;
			if (miny > UKPT[i].Y)
				miny = UKPT[i].Y;
			if (maxx < UKPT[i].X)
				maxx = UKPT[i].X;
			if (maxy < UKPT[i].Y)
				maxy = UKPT[i].Y;
		}
		else
		{
			if (minx > KPT[i-ukpt_num].X)
				minx = KPT[i - ukpt_num].X;
			if (miny > KPT[i - ukpt_num].Y)
				miny = KPT[i - ukpt_num].Y;
			if (maxx < KPT[i - ukpt_num].X)
				maxx = KPT[i - ukpt_num].X;
			if (maxy < KPT[i - ukpt_num].Y)
				maxy = KPT[i - ukpt_num].Y;
		}
	}
}

//控制网框架绘制
void ControlNet::DrawFrame(CDC*pDC,CRect&rect, double& orgx, double& orgy, double& scalex,double&scaley)
{
	
	double lowx, lowy, highx, highy;
	MAX_(lowx, lowy, highx, highy);//求得测图区的范围
	rect.DeflateRect(rect.Width()*0.1, rect.Height()*0.1);
	int num = 5;
	int dx = rect.Width() / (num);
	int dy = rect.Height() / (num);
	int tx = (highx - lowx) / num;
	int ty = (highy - lowy) / num;
	
	scalex = double(dy) / tx;//x方向的比例系数
	scaley = double(dx) / ty;//y方向的比例系数
	CRect grideRect(rect.left, rect.top, rect.left + dx * num, rect.top + dy * num);//调整矩形大小，以便于格网的绘制
	orgx = grideRect.bottom;
	orgy = grideRect.left;
	CPen gridpen(PS_DASH, 0, RGB(0, 100, 200));//创建格网画笔
	CPen* pOldPen = pDC->SelectObject(&gridpen);
	//绘制垂直线
	for (int i = 0; i <= num; i++)
	{
		pDC->MoveTo(grideRect.left + i * dx, grideRect.bottom);
		pDC->LineTo(grideRect.left + i * dx, grideRect.top);
		CString str;
		str.Format(_T("%d"), int(lowy)+i*ty);
		pDC->TextOutW(grideRect.left + i * dx-10, grideRect.bottom + 10, str);
		if (i == 0)
		{
			pDC->MoveTo(grideRect.left, grideRect.top);
			pDC->LineTo(grideRect.left, grideRect.top - 30);
			pDC->MoveTo(grideRect.left-10, grideRect.top - 20);
			pDC->LineTo(grideRect.left, grideRect.top - 30);
			pDC->LineTo(grideRect.left + 10, grideRect.top - 20);
		}
	}
	//绘制水平线
	for (int i = 0; i <= num; i++)
	{
		pDC->MoveTo(grideRect.left, grideRect.top+i*dy);
		pDC->LineTo(grideRect.right , grideRect.top + i * dy);
		CString str;
		str.Format(_T("%d"), int(lowx) + i * tx);
		pDC->TextOutW(grideRect.left-60, grideRect.bottom-i*dy-5, str);
		if (i == 0)
		{
			pDC->MoveTo(grideRect.right, grideRect.bottom);
			pDC->LineTo(grideRect.right+30, grideRect.bottom);
			pDC->MoveTo(grideRect.right+20, grideRect.bottom- 10);
			pDC->LineTo(grideRect.right+30, grideRect.bottom);
			pDC->LineTo(grideRect.right+20, grideRect.bottom+10);
		}
	}
	//绘制比例尺
	pDC->SelectObject(pOldPen);
	pDC->MoveTo(grideRect.left + grideRect.Width() / 2, grideRect.bottom + 40);
	pDC->LineTo(grideRect.left + grideRect.Width() / 2 + 30, grideRect.bottom + 40);
	CString sca;
	sca.Format(_T("1:%d"), tx / dx);
	pDC->TextOutW(grideRect.left + grideRect.Width() / 2 + 40, grideRect.bottom + 30, sca);
	pDC->TextOutW(grideRect.right+45, grideRect.bottom-5, _T("Y"));
	pDC->TextOutW(grideRect.left-5, grideRect.top-45, _T("X"));
	
	
	
}

void ControlNet::checkspace()
{
	if (UKPT != NULL)
	{
		delete[]UKPT;
		UKPT = NULL;
	}
	if (KPT != NULL)
	{
		delete[]KPT;
		KPT = NULL;
	}
	if (OBSDIRECT != NULL)
	{
		delete[]OBSDIRECT;
		OBSDIRECT = NULL;
	}
	if (OBSDIST != NULL)
	{
		delete[]OBSDIST;
		OBSDIST = NULL;
	}
}

//绘图主函数，绘制控制网以及误差椭圆
void ControlNet::MainDraw(CDC* pDC, CRect& rect)
{
	double scalex, scaley, orgx, orgy;
	double minx, miny, maxx, maxy;
	MAX_(minx, miny, maxx, maxy);//求得测图区的范围
	DrawFrame(pDC, rect,orgx,orgy,scalex,scaley);
	/*scalex = scalex * DScale;
	scaley= scaley * DScale;*/
	CBrush BTri;
	BTri.CreateSolidBrush(RGB(0, 100, 100));
	//绘制控制网
	for (int i = 0; i < obsDirect_count; i++)
	{
		int  drx1, dry1;
		drx1 = int((OBSDIRECT[i].cpST->Y - miny) * scaley + orgy);
		dry1 =int(- (OBSDIRECT[i].cpST->X - minx) * scalex + orgx);
		int drx2, dry2;
		drx2 =int( (OBSDIRECT[i].cpOB->Y - miny) * scaley + orgy);
		dry2 =int( - (OBSDIRECT[i].cpOB->X - minx) * scalex + orgx);
		if (OBSDIRECT[i].cpST->type == true && OBSDIRECT[i].cpOB->type == true)//两点均为已知点
		{
			
			pDC->MoveTo(drx1,dry1);
			pDC->LineTo(drx2, dry2);

			pDC->MoveTo(drx1+5, dry1+5);
			pDC->LineTo(drx2+5, dry2+5);
			//已知边绘制双线

			//已知点绘制三角形
			pDC->SelectObject(BTri);
			CPoint P1[3],P2[3];
			P1[0].SetPoint(drx1, dry1 - 10);
			P1[1].SetPoint(drx1 + 10, dry1 + 10);
			P1[2].SetPoint(drx1-10, dry1 + 10);

			P2[0].SetPoint(drx2, dry2 - 10);
			P2[1].SetPoint(drx2 + 10, dry2 + 10);
			P2[2].SetPoint(drx2-10, dry2 + 10);
			pDC->Polygon(P1,3);
			pDC->Polygon(P2, 3);
			BTri.DeleteObject();
			
		}
		else
		{
			pDC->MoveTo(drx1, dry1);
			pDC->LineTo(drx2, dry2);
		}
	}
	for (int i = 0; i < kpt_num; i++)
	{
		pDC->TextOutW((KPT[i].Y - miny) * scaley + orgy+15, -(KPT[i].X - minx) * scalex + orgx -15, KPT[i].strID);
	}
	//绘制误差椭圆
	for (int i = 0; i < ukpt_num; i++)
	{
		UKPT[i].EDRAW.SETORG((UKPT[i].Y - miny) * scaley + orgy, -(UKPT[i].X - minx) * scalex + orgx);
		UKPT[i].EDRAW.SETScale(30);
		UKPT[i].EDRAW.Draw(pDC,rect,UKPT[i].dAlfa,UKPT[i].dE,UKPT[i].dF);
		pDC->TextOutW((UKPT[i].Y - miny) * scaley + orgy+15, -(UKPT[i].X - minx) * scalex + orgx-15, UKPT[i].strID);
	}
	

}

//结果可视化输出
void ControlNet::OutPut(CString& result)
{
	CString AResult,DResult,PResult,temp;
	AResult.Format(_T("------------------------方向观测成果------------------------\r\n"));
	temp.Format(_T("测站\t照准\t方向值(dms)\t改正数(s)\t平差后值(dms)\t观测精度(″)\r\n"));
	AResult += temp;
	for (int i = 0; i < obsDirect_count; i++)
	{
		temp.Empty();
		CAngle pg;
		pg(RAD) = OBSDIRECT[i].ObsAngle(RAD) + V(i, 0) / Ro;
		temp.Format(_T("%s\t%s\t%.6f\t%.2f\t%.6f\t%.2f\r\n"), OBSDIRECT[i].cpST->strID, OBSDIRECT[i].cpOB->strID, OBSDIRECT[i].ObsAngle(DMS), V(i, 0), pg(DMS),OBSDIRECT[i].m);
		AResult += temp;
	}
	DResult.Format(_T("------------------------距离观测成果------------------------\r\n"));
	temp.Format(_T("测站\t照准\t距离(m)\t改正数(m)\t平差后值(m)\t观测精度(cm)\r\n"));
	DResult += temp;
	for (int i = 0; i < obsDist_count; i++)
	{
		temp.Empty();
		temp.Format(_T("%s\t%s\t%.4f\t%.4f\t%.4f\t%.2f\r\n"), OBSDIST[i].cpBG->strID, OBSDIST[i].cpEND->strID, OBSDIST[i].Dist, V(i + obsDirect_count, 0) / 100, OBSDIST[i].Dist + V(i + obsDirect_count, 0) / 100,OBSDIST[i].m);
		DResult += temp;
	}
	PResult.Format(_T("------------------------控制点平差成果------------------------\r\n"));
	temp.Format(_T("点名\tX(m)\tY(m)\t长轴(m)\t短轴(m)\t长轴方位(dms)\t点位中误差(m)\t备注\r\n"));
	PResult += temp;
	for (int i = 0; i < kpt_num; i++)
	{
		temp.Empty();
		CAngle pq;
		pq(RAD) = KPT[i].dAlfa;
		temp.Format(_T("%s\t%.4f\t%.4f\t\t\t\t\t已知点\r\n"), KPT[i].strID, KPT[i].X, KPT[i].Y, pq(DMS));
		PResult += temp;
	}
	for (int i = 0; i < ukpt_num; i++)
	{
		temp.Empty();
		CAngle pq;
		pq(RAD) = UKPT[i].dAlfa;
		temp.Format(_T("%s\t%.4f\t%.4f\t%.4f\t%.4f\t%.6f\t%.4f\r\n"), UKPT[i].strID, UKPT[i].X, UKPT[i].Y, UKPT[i].dE / 100, UKPT[i].dF / 100, pq(DMS),UKPT[i].dP);
		PResult += temp;
	}
	result.Format(_T("--------------------------平差成果--------------------------\r\n"));
	temp.Empty();
	temp.Format(_T("迭代次数:%d\r\n单位权中误差(″):%.4f\r\n"), times,m0);
	result += temp + AResult + DResult + PResult;
}



//概算函数，逐点计算法，需要递归计算
void ControlNet::RoughCal()
{

	for (int i = 0; i < ukpt_num; i++)
	{
		if (UKPT[i].tag == FALSE)//首先判断该点是否已经完成概算
		{
			for (int j = 0; j < obsDist_count; j++)
			{
				if (OBSDIST[j].cpEND->strID == UKPT[i].strID && OBSDIST[j].cpBG->tag == TRUE && UKPT[i].tag == FALSE)//寻找和该点有距离观测关系的已知点
				{
					for (int k = 0; k < obsDirect_count; k++)
					{
						if (OBSDIRECT[k].cpST == OBSDIST[j].cpBG && OBSDIRECT[k].cpOB->tag == TRUE)//寻找另一已知点
						{
							double T0 = AzimuthCal(OBSDIRECT[k].cpOB->X, OBSDIRECT[k].cpOB->Y, OBSDIRECT[k].cpST->X, OBSDIRECT[k].cpST->Y);//已知边的方位角
							for (int l = 0; l < obsDirect_count; l++)
							{
								if (OBSDIRECT[l].cpST == OBSDIRECT[k].cpST && OBSDIRECT[l].cpOB->strID == UKPT[i].strID)//寻找观测角度
								{
									double  A;
									double T1 = OBSDIRECT[l].ObsAngle(RAD) - OBSDIRECT[k].ObsAngle(RAD);//两方向的夹角
									if (T1 >= 0)
										A = T1 + T0 - PI;
									else
										A = PI + T1 + T0;
									UKPT[i].X = OBSDIST[j].cpBG->X + OBSDIST[j].Dist * cos(A);
									UKPT[i].Y = OBSDIST[j].cpBG->Y + OBSDIST[j].Dist * sin(A);
									UKPT[i].tag = TRUE;
									break;
								}
							}
							break;
						}

					}


				}
			}

		}
	}
	for (int n = 0; n < ukpt_num; n++)
	{
		if (UKPT[n].tag == FALSE)
		{

			RoughCal();
		}
	}
}
//平差函数
void ControlNet::Adjust()
{
	RoughCal();
	//计算定向边的近似方位角
	//概算检查
	CString checkg,tempg;
	for (int i = 0; i < kpt_num+ukpt_num; i++)
	{
		if (i < kpt_num)
		{
			tempg.Format(_T("%.6f\t%.6f\r\n"), KPT[i].X, KPT[i].Y);
			checkg += tempg;
		}
		else
		{
			tempg.Format(_T("%.6f\t%.6f\r\n"), UKPT[i-kpt_num].X, UKPT[i-kpt_num].Y);
			checkg += tempg;
		}
	}

	STA = new STATION[station_num];//为测站开辟空间，未释放

	B.ReSize(obsDirect_count + obsDist_count, 2 * ukpt_num + station_num);
	X.ReSize(2 * ukpt_num + station_num, 1);
	L.ReSize(obsDirect_count + obsDist_count, 1);
	P.ReSize(obsDirect_count + obsDist_count, obsDirect_count + obsDist_count);
	Nbb.ReSize(2 * ukpt_num + station_num, 2 * ukpt_num + station_num);
	//定权
	for (int i = 0; i < obsDirect_count + obsDist_count; i++)
	{
		if (i < obsDirect_count)
			P(i, i) = 1;
		else
			P(i, i) = 1E5 / (OBSDIST[i - obsDirect_count].Dist);//!!!如何定权存疑，单位cm
	}
	
	//迭代次数
	 times = 0;
	 double max_ = 0.0;//限差判断变量
	do
	{ 
		//构造B、L矩阵
		max_ = 0;
		int k = 0;
		//计算测站定向边的方位角
		for (int i = 0; i < obsDirect_count; i++)
		{
			if (OBSDIRECT[i].ObsAngle(RAD) == 0)
			{
				STA[k].station = OBSDIRECT[i].cpST;
				STA[k].start = OBSDIRECT[i].cpOB;
				STA[k].st_(RAD) = AzimuthCal(OBSDIRECT[i].cpST->X, OBSDIRECT[i].cpST->Y, OBSDIRECT[i].cpOB->X, OBSDIRECT[i].cpOB->Y);
				k++;
			}
		}

		//方向角观测误差方程

		for (int i = 0; i < obsDirect_count; i++)
		{
			//需要修改，按照测站点计算方位角

			double Tij = AzimuthCal(OBSDIRECT[i].cpST->X, OBSDIRECT[i].cpST->Y, OBSDIRECT[i].cpOB->X, OBSDIRECT[i].cpOB->Y);

			double  Dij = (DistCal(OBSDIRECT[i].cpOB->X, OBSDIRECT[i].cpOB->Y, OBSDIRECT[i].cpST->X, OBSDIRECT[i].cpST->Y));
			int s = SearchSTPosUsingID(OBSDIRECT[i].cpST->strID);

			if (OBSDIRECT[i].cpOB->type == FALSE && OBSDIRECT[i].cpST->type == FALSE)//如果两点均是未知点
			{
				int n = SearchUKPosUsingID(OBSDIRECT[i].cpOB->strID);
				int m = SearchUKPosUsingID(OBSDIRECT[i].cpST->strID);
				B(i, 2 * m) = Ro * sin(Tij) / (Dij * 100);//单位要统一，定权的单位分别是秒和厘米
				B(i, 2 * m + 1) = -Ro * cos(Tij) / (Dij * 100);
				B(i, 2 * n) = -Ro * sin(Tij) / (Dij * 100);
				B(i, 2 * n + 1) = Ro * cos(Tij) / (Dij * 100);
				B(i, 2 * ukpt_num + s) = -1;


			}
			else if (OBSDIRECT[i].cpOB->type == FALSE && OBSDIRECT[i].cpST->type == TRUE)//如果测站点为已知点
			{
				int n = SearchUKPosUsingID(OBSDIRECT[i].cpOB->strID);
				B(i, 2 * n) = -Ro * sin(Tij) / (Dij * 100);
				B(i, 2 * n + 1) = Ro * cos(Tij) / (Dij * 100);
				B(i, 2 * ukpt_num + s) = -1;
			}
			else if (OBSDIRECT[i].cpOB->type == TRUE && OBSDIRECT[i].cpST->type == FALSE)
			{
				int m = SearchUKPosUsingID(OBSDIRECT[i].cpST->strID);
				B(i, 2 * m) = Ro * sin(Tij) / (Dij * 100);
				B(i, 2 * m + 1) = -Ro * cos(Tij) / (Dij * 100);
				B(i, 2 * ukpt_num + s) = -1;

			}
			else
			{
				B(i, 2 * ukpt_num + s) = -1;
			}

			//L矩阵
			double dA = -STA[s].st_(RAD) + Tij;
			if (dA < 0)
			{
				L(i, 0) = -OBSDIRECT[i].ObsAngle(RAD) + 2 * PI + dA;
			}
			else if (dA > 0)
			{
				L(i, 0) = -OBSDIRECT[i].ObsAngle(RAD) + dA;
			}
			else
			{
				L(i, 0) = 0;
			}
		}


		for (int i = 0; i < obsDirect_count; i++)
		{
			L(i, 0) = Ro * L(i, 0);//弧度转为秒
		}

		//边长观测误差方程

		for (int i = 0; i < obsDist_count; i++)
		{
			double TT = AzimuthCal(OBSDIST[i].cpBG->X, OBSDIST[i].cpBG->Y, OBSDIST[i].cpEND->X, OBSDIST[i].cpEND->Y);
			L(obsDirect_count + i, 0) = (DistCal(OBSDIST[i].cpBG->X, OBSDIST[i].cpBG->Y, OBSDIST[i].cpEND->X, OBSDIST[i].cpEND->Y) - OBSDIST[i].Dist) * 100;
			if (OBSDIST[i].cpBG->type == FALSE && OBSDIST[i].cpEND->type == FALSE)//如果起始点和终点均为未知点
			{
				int n = SearchUKPosUsingID(OBSDIST[i].cpEND->strID);
				int m = SearchUKPosUsingID(OBSDIST[i].cpBG->strID);
				B(obsDirect_count + i, 2 * m) = -cos(TT);
				B(obsDirect_count + i, 2 * m + 1) = -sin(TT);
				B(obsDirect_count + i, 2 * n) = cos(TT);
				B(obsDirect_count + i, 2 * n + 1) = sin(TT);
			}
			else if (OBSDIST[i].cpBG->type == TRUE && OBSDIST[i].cpEND->type == FALSE)
			{
				int n = SearchUKPosUsingID(OBSDIST[i].cpEND->strID);
				B(obsDirect_count + i, 2 * n) = cos(TT);
				B(obsDirect_count + i, 2 * n + 1) = sin(TT);
			}
			else if (OBSDIST[i].cpBG->type == FALSE && OBSDIST[i].cpEND->type == TRUE)
			{
				int m = SearchUKPosUsingID(OBSDIST[i].cpBG->strID);
				B(obsDirect_count + i, 2 * m) = -cos(TT);
				B(obsDirect_count + i, 2 * m + 1) = -sin(TT);
			}

		}

		

		//解算法方程
		Nbb = (~B * P * B);
		X = Nbb.Inv() * ~B * P * (-1 * L);
		
		for (int i = 0; i < ukpt_num; i++)
		{
			UKPT[i].X += X(2 * i, 0) / 100;//改正未知数
			UKPT[i].Y += X(2 * i + 1, 0) / 100;
		}

		for (int i = 0; i < 2 * ukpt_num; i++)
		{
			if (max_ < fabs(X(i, 0)))
				max_ = fabs(X(i, 0));
		}

		times++;
	} while (max_ > 0.001);//当dX满足限差时结束迭代

	delete[]STA;
	STA = NULL;

}
//精度评定函数
void ControlNet::Estimate()
{
	V.ReSize(obsDirect_count + obsDist_count, 1);
	Qxx.ReSize(ukpt_num * 2 + station_num, ukpt_num * 2 + station_num);
	V = B * X + L;
	m0 = sqrt((~V * P * V)(0, 0) / (obsDirect_count + obsDist_count - ukpt_num * 2 - station_num));//单位权中误差
	Qxx = Nbb.Inv();//协因数矩阵
	//计算点位精度以及误差椭圆参数
	for (int i = 0; i < ukpt_num; i++)
	{
		double mx = m0 * sqrt(Qxx(i * 2, i * 2));
		double my = m0 * sqrt(Qxx(i * 2 + 1, i * 2 + 1));
		UKPT[i].dP = m0 * sqrt(mx * mx + my * my);//点位误差
		UKPT[i].dAlfa = 0.5 * atan(2 * Qxx(i * 2, i * 2 + 1) / (Qxx(i * 2, i * 2) - Qxx(i * 2 + 1, i * 2 + 1)));//误差椭圆元素赋值
		UKPT[i].dE = m0 * sqrt(Qxx(i * 2, i * 2) + Qxx(i * 2, i * 2 + 1) * tan(UKPT[i].dAlfa));//长轴
		UKPT[i].dF = m0 * sqrt(Qxx(i * 2, i * 2) + Qxx(i * 2, i * 2 + 1) * tan(UKPT[i].dAlfa + 0.5 * PI));//短轴
		double temp;//所求的E,F均为极值，需要判断大小
		if (UKPT[i].dF > UKPT[i].dE)
		{
			temp = UKPT[i].dE;
			UKPT[i].dE = UKPT[i].dF;
			UKPT[i].dF = temp;
			UKPT[i].dAlfa += 0.5 * PI;
		}
	}
	//计算观测精度
	CMatrix fT;
	fT.ReSize(obsDirect_count + obsDist_count, 2 * ukpt_num);//截取B矩阵
	for (int i = 0; i < obsDirect_count + obsDist_count; i++)
	{
		for (int j = 0; j < 2 * ukpt_num; j++)
			fT(i, j) = B(i, j);
	}

	CMatrix Qx;//矩阵维数不一致
	Qx.ReSize(2 * ukpt_num, 2 * ukpt_num);
	for (int i = 0; i < 2 * ukpt_num; i++)
	{
		for (int j = 0; j < 2 * ukpt_num; j++)
			Qx(i, j) = Qxx(i, j);
	}
	QF.ReSize(obsDirect_count + obsDist_count, obsDirect_count + obsDist_count);
	QF = fT * Qx * ~fT;
	for (int i = 0; i < obsDirect_count + obsDist_count; i++)
	{
		if (i < obsDirect_count)
			OBSDIRECT[i].m = m0 * sqrt(QF(i, i));//角度观测精度
		else
			OBSDIST[i - obsDirect_count].m = m0 * sqrt(QF(i, i));//距离观测精度
	}
	
}

至此，导线网平差的主要代码就已经编写完毕了，然后就是界面设计了，由于对MFC较为熟练，所以我采用了基于MFC的模板，各位读者也可以使用基于单个文档来写。这一部分代码我就不展示了（这部分写的不太行(* /ω＼*) ）。

五、界面设计

界面设计依个人的喜好而设计，尽量要大气美观，对使用者友好。 

以下是我设计的界面，比较简单，仅供参考(* /ω＼*)

 

 

 数据格式说明：

 六、总结

概算和平差是本次实习的两大难点，概算的关键在于找到与其匹配的计算元（V型结构），一次算不完就进行递归，直至所有待定点解算出来为止；平差部分是整个程序的核心，这一部分需要理解平差原理，列立误差方程，从而得到B,P,L三个矩阵，再利用平差公式解算。

以下是我编写程序时所遇见的一些bug及解决办法：

1、坐标概算无法跳出循环或得不到正确的计算结果

检查方位角的角度类型，计算过程需要使用DMS形式；

是否在循环嵌套里面忘记了加上break，跳不出循环，导致重复计算；

完成概算后是否修改状态，设为已知。

2、平差部分迭代次数过大甚至跳不出循环

检查是否改化了未知数；

检查判断条件；

检查坐标概算；

检查单位是否统一，P矩阵是有量纲的，角度和距离单位需要全程保持一致（统一化成″和cm）；

利用Look函数检查B、L矩阵的结构。

3、图像刷新问题

这里推荐http://t.csdn.cn/iWiPN，里面介绍了几种解决办法。

4、内存泄漏问题

内存泄漏问题主要就是通过new开辟的空间在程序运行完毕之后没有delete掉，解决办法就是去匹配new和delete。如果觉得利用new和delete动态开辟空间太复杂，可以用vector来给变量分配空间。vector的使用方法可以参考http://t.csdn.cn/EzBGV

 因本人水平有限，所设计的程序还不够完善，如有不正之处，还请各位看官批评指正，多多指教。

ShockeyRolli
关注
  • 6
    点赞
  • 踩
  • 81
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
专栏目录
控制网平差C++代码实现
12-24
C++实现的控制网平差实例 基于VS2008版本的 高版本的VS可以用
水平网平差c++实现
01-08
实先水平网平差，内有输入文件实例，结果文件实例，直接修改储存路径，原始数据，直接运行，可输入结果。设计了二维矩阵计算类，利用矩阵计算求取平差值。
间接平差代码实现（C++）（测量学）
最新发布
weixin_71272801的博客
04-24 1114
用C++语言，使用Eigen库实现，矩阵的运算。
测量平差程序源代码（已经调试通过）
04-14
这是军测的老师写的平差的程序的基础函数，包括求逆，矩阵传播，矩阵输出等功能的实现，比较有用。
导线网平差程序系统设计—CSU测绘实习
SEREND1P17Y的博客
06-14 2879
导线网平差程序系统设计，控制网平差程序对野外控制网观测数据进行平差数据处理，目的是根据最小二乘原理，消除网中的各种几何矛盾，求出全网各待定点元素如未知点的平面坐标。❤来自CSU测绘实习
GPS控制网网平差
gnsscome的博客
03-10 2096
GPS控制网平差概要
导线网平差程序系统设计
weixin_57056661的博客
08-08 1734
原始数据 界面设计 点击画图按钮后弹出的对话框 类的设计 Data类：存放点号、坐标、未知点误差椭圆元素 Gangle类：存储起测站、照准点、角度观测值 Length类：存储起测站、照准点、边长观测值 support类：存放各种函数和矩阵 Eigen类：用于矩阵运算 源代码 <Data.h> #pragma once #include <afxext.h> #include <math.h> #include <local....
如何理解条件平差
欧阳罢笔的博客
02-12 8194
为什么要平差，我们知道是因为我们测量的数据（在有多余观测的条件下）产生了矛盾。比如我们测一个三角形，如果只测了两个角，我们知道第三个角一定是180度减前两个角的和。但是如果我们测了三个角，但三个角加起来不等于180度，那就说明测得有误差，我们就需要用一些手段进行平差，消除这些误差 平差最主要的两种是间接平差和条件平差，它们的区别可以参考这里，间接平差使用的比条件平差多得多。简单的说条件平差就是以公...
测量平差之附有参数的条件平差
冯玮的博客
03-26 6348
附有参数的条件平差也可看作是一种间接平差和条件平差的混合模型。增加条件方程是“降维”，它可以把复杂的问题变得简单；添加参数是“升维”，它可以把问题变得理论上可分的概率更大。两者并不矛盾，而且相辅相成。
经典平差C++版（测绘专业课设）
04-25
该平差课设包括条件平差，间接平差，附有参数的条件平差和简介平差，这些将是测绘本科生课设必经之路，之后会有专门的一遍博客讲述它的用法
基于c++语言的附和导线测量程序
03-02
导线测量（traversesurvey）指的是测量导线长度、转角和高程，以及推算坐标等的作业。导线测量布设灵活，推进迅速，受地形限制小，边长精度分布均匀。如在平坦隐蔽、交通不便、气候恶劣地区，采用导线测量法布设大地控制网是有利的。但导线测量控制面积小、检核条件少，方位传算误差大。
条件平差的 C++代码
04-10
条件平差的
导线网测量平差C++
09-13
运用C++语言，MFC对话框结构写导线网平差
水准网平差C++语言实现
12-01
为解决水准网平差，利用间接平差模型，采用最小二乘进行精度评定，可直接利用VS2019以及其他版本运行运行。程序能较好的计算出平差结果
ControlNetAdjustment_导线网平差_C++_测绘程序_
10-01
"ControlNetAdjustment_导线网平差_C++_测绘程序_"这个项目是利用C++编程语言实现的一种针对导线网进行平差计算的软件工具。导线网平差通常应用于大地测量、工程测量和地籍测量等，它涉及到多个控制点间的观测数据...
基于MATLAB的高程控制网平差系统的设计与应用
05-30
文中利用MATLAB开发的计算机程序代码短小高效,能够快速的实现图形绘制和精度分析,与利用C++等高级编程语言开发的平差处理程序相比,可以在短时间内完成对平差的设计和验证。
《测绘程序设计实习》实验报告(MFC,C++)
qq_54353206的博客
07-08 2806
此文章为测绘程序设计实习中所设计的一个平差绘图项目。
测量平差MATLAB方法
yitian_00的博客
03-19 2327
测量平差MATLAB方法
python语言间接平差程序
07-15
### 回答1： Python语言可以编写一个间接平差程序。间接平差是一种通过观测数据进行测量结果精度分析和优化的方法，可以用于平差控制点的坐标和误差分析。 首先，我们可以使用Python中的NumPy库来进行矩阵运算。通过构建系数矩阵和常数向量，可以将间接平差问题转化为线性方程组求解。系数矩阵中的每一行对应一个观测方程，常数向量中的每个元素对应观测数据的观测值和已知值的差。 然后，我们可以使用线性代数库中的函数来求解线性方程组。例如，可以使用NumPy库中的`numpy.linalg.solve`函数来求解最小二乘解。这个函数可以输入系数矩阵和常数向量，返回最小二乘解的结果。 最后，我们可以使用Python中的Matplotlib库来可视化间接平差的结果。可以绘制观测数据的残差图，观察平差结果的拟合程度和误差情况。可以绘制平差后的控制点坐标图，以直观地展示平差结果的准确性。 通过以上步骤，我们可以编写一个间接平差程序，用Python语言实现。这个程序可以输入观测数据和已知值，进行间接平差计算，并输出平差结果和误差分析。使用Python编写间接平差程序可以提供灵活性和扩展性，方便进行参数调整和结果可视化。同时，Python语言的简洁性和易学性也使得编写和使用平差程序更加方便。 ### 回答2： Python语言可以用来编写间接平差程序。间接平差是一种常用的测量数据处理方法，用于通过多次测量来提高测量结果的精度。 在编写间接平差程序时，可以使用Python的数学计算库（如numpy）来进行矩阵运算和数值计算。通过矩阵运算，可以将观测数据转化为平差条件方程，并构建系数矩阵和常数向量。然后，可以使用线性代数方法（如最小二乘法）求解平差方程，得到未知参数的估计值。 Python的强大的数据处理和分析库（如pandas）也可以用于数据的输入、输出和分析。可以将观测数据从外部文件读取到Python程序中，进行数据预处理、筛选和清洗。然后，可以将结果输出到文件或数据库中，以供其他程序或人工使用。 此外，Python还有丰富的绘图库（如matplotlib）可以用于绘制观测数据的图表，帮助分析数据的分布和趋势，并验证平差结果的合理性。 编写间接平差程序时，还可以利用Python的模块化特性，将代码分为多个函数或类，提高代码的可读性和复用性。可以编写不同功能模块的单元测试，确保代码的正确性。 总之，Python语言具有丰富的数学计算、数据处理、绘图和模块化等功能，非常适合用于编写间接平差程序。通过Python的灵活性和易用性，可以快速开发出高效、稳定和可靠的平差程序，提高测量数据处理的效率和准确性。 ### 回答3： Python语言是一种通用的编程语言，可以用于编写各种类型的程序，包括间接平差程序。所谓间接平差，是一种通过测量数据之间的关系来估计未知量的技术。以下是一个简单的用Python编写的间接平差程序的示例： 首先，我们需要导入一些必要的库，如numpy和scipy。然后，定义一个函数来执行间接平差的计算。该函数接受测量数据和已知量作为输入，并返回估计的未知量。 在函数内部，我们可以使用最小二乘法来进行间接平差计算。通过建立一个误差模型和目标函数，我们可以最小化测量数据和模型之间的残差，并得到最优的未知量估计。 接下来，我们可以编写一个主程序来使用这个函数。在主程序中，我们可以定义一些测量数据和已知量，并调用间接平差函数来获得估计的未知量。 最后，我们可以在屏幕上打印出估计的未知量，并进行进一步的分析和处理。 这只是一个简单的例子，实际的间接平差程序可能会更加复杂，涉及更多的数学和统计方法。但是，使用Python编写间接平差程序可以大大简化程序的开发和调试过程，并且具有很好的可读性和可扩展性。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值