摄影测量前方交会编程实现

在最近用摄影测量前方交会的知识过程中，需要将前方交会用编程来实现，因为我在设计本次代码的过程中用到了OpenCV这个库来进行矩阵运算，相比之前的用纯C++代码来实现矩阵来说不知道方便了多少，提供给大家一种新的思想，供大家参考，我的OpenCV库是3.4.0版本的。摄影测量前方交会的理论在相关的专业书里都有叙述，在这里我就不多说了，下面是代码：

//利用像片外方位元素进行前方交会
#include<iostream>
#include<fstream>
#include <iomanip>
#include "core/core.hpp"

using namespace std;
using namespace cv;

int main()
{
	//从文件中读取同名像点坐标到数组
	double T_Z_B[2][2] = {0.0 };
	ifstream infile1;//定义文件流对象
	infile1.open("C:\\Users\\ASUS\\Desktop\\6\\15T.txt");//打开文档
	double *ptr1 = &T_Z_B[0][0];
	while (!infile1.eof())
	{
		infile1 >> *ptr1;//这是把文档里面的数对应放在ptr位置上的数值上
		ptr1++;
	}
	infile1.close();
	double x1 = T_Z_B[0][0];
	double y1 = T_Z_B[0][1];
	double x2 = T_Z_B[1][0];
	double y2 = T_Z_B[1][1];
	//定义内方位元素
	double x0 = -0.00000;//mm
	double y0 = 0.00000;//mm
	double f = 70.50;//mm
	//从文件中读取外方位元素到数组
	/*int i, j;*/
	double data[2][6] = { 0.0 };
	ifstream infile2;//定义文件流对象
	infile2.open("C:\\Users\\ASUS\\Desktop\\6\\2037_2038.txt");//打开文档
	double *ptr2 = &data[0][0];
	while (!infile2.eof())
	{
		infile2 >> *ptr2;//这是把文档里面的数对应放在ptr位置上的数值上
		ptr2++;
	}
	infile2.close();
	/*for (i = 0; i < 2; i++)
	{
		for (j = 0; j < 6; j++)
		{
			cout.precision(12);
		  cout << data[i][j] << "  ";
		}
		cout << '\n';
	}*/
	//左右像片的外方位元素
	double Xs1 = data[0][0];
	double Ys1 = data[0][1];
	double Zs1 = data[0][2];
	double phi1 = data[0][3];
	double omig1 = data[0][4];
	double kappa1 = data[0][5];
	double Xs2 = data[1][0];
	double Ys2 = data[1][1];
	double Zs2 = data[1][2];
	double phi2 = data[1][3];
	double omig2 = data[1][4];
	double kappa2 = data[1][5];
	//cout.precision(12);//控制输出的小数点位数
//计算摄影基线的三个分量
double Bx = Xs2 - Xs1;
double By = Ys2 - Ys1;
double Bz = Zs2 - Zs1;
//利用外方位角元素计算左右像片的旋转矩阵R1和R2，用OpenCV矩阵，方便
double a1 = cos(phi1)*cos(kappa1) - sin(phi1)*sin(omig1)*sin(kappa1);
double a2 = -cos(phi1)*sin(kappa1) - sin(phi1)*sin(omig1)*cos(kappa1);
double a3 = -sin(phi1)*cos(omig1);
double b1 = cos(omig1)*sin(kappa1);
double b2 = cos(omig1)*cos(kappa1);
double b3 = -sin(omig1);
double c1 = sin(phi1)*cos(kappa1) + cos(phi1)*sin(omig1)*sin(kappa1);
double c2 = -sin(phi1)*sin(kappa1) + cos(phi1)*sin(omig1)*cos(kappa1);
double c3 = cos(phi1)*cos(omig1);
Mat R1 = Mat::ones(3, 3, CV_64F);
R1.at<double>(0, 0) = a1;
R1.at<double>(0, 1) = a2;
R1.at<double>(0, 2) = a3;
R1.at<double>(1, 0) = b1;
R1.at<double>(1, 1) = b2;
R1.at<double>(1, 2) = b3;
R1.at<double>(2, 0) = c1;
R1.at<double>(2, 1) = c2;
R1.at<double>(2, 2) = c3;
double a11 = cos(phi2)*cos(kappa2) - sin(phi2)*sin(omig2)*sin(kappa2);
double a22 = -cos(phi2)*sin(kappa2) - sin(phi2)*sin(omig2)*cos(kappa2);
double a33 = -sin(phi2)*cos(omig2);
double b11 = cos(omig2)*sin(kappa2);
double b22 = cos(omig2)*cos(kappa2);
double b33 = -sin(omig2);
double c11 = sin(phi2)*cos(kappa2) + cos(phi2)*sin(omig2)*sin(kappa2);
double c22 = -sin(phi2)*sin(kappa2) + cos(phi2)*sin(omig2)*cos(kappa2);
double c33 = cos(phi2)*cos(omig2);
Mat R2 = Mat::ones(3, 3, CV_64F);
R2.at<double>(0, 0) = a11;
R2.at<double>(0, 1) = a22;
R2.at<double>(0, 2) = a33;
R2.at<double>(1, 0) = b11;
R2.at<double>(1, 1) = b22;
R2.at<double>(1, 2) = b33;
R2.at<double>(2, 0) = c11;
R2.at<double>(2, 1) = c22;
R2.at<double>(2, 2) = c33;
//计算同名像点的像空间辅助坐标系(X1,Y1,Z1)与(X2,Y2,Z2)
Mat RR1, RR2;
Mat R3 = Mat::ones(3, 1, CV_64F);
R3.at<double>(0, 0) = x1;
R3.at<double>(1, 0) = y1;
R3.at<double>(2, 0) = -f;
Mat R33 = Mat::ones(3, 1, CV_64F);
R33.at<double>(0, 0) = x2;
R33.at<double>(1, 0) = y2;
R33.at<double>(2, 0) = -f;
RR1 = R1 * R3;
double X1 = RR1.at<double>(0, 0);
double Y1 = RR1.at<double>(1, 0);
double Z1 = RR1.at<double>(2, 0);
RR2 = R2 * R33;
double X2 = RR2.at<double>(0, 0);
double Y2 = RR2.at<double>(1, 0);
double Z2 = RR2.at<double>(2, 0);
//计算投影系数N1,N2
double N1 = (Bx*Z2 - Bz * X2) / (X1*Z2 - X2 * Z1);
double N2 = (Bx*Z1 - Bz * X1) / (X1*Z2 - X2 * Z1);
计算投影系数N1,N2
//double N1 = (Bx*Z2 - Bz * X2) / (X1*Z2 - X2 * Z1);
//double N2 = (Bx*Z1 - Bz * X1) / (X1*Z2 - X2 * Z1);
//计算地面点的左像辅系坐标(deteX,deteY,deteZ)
double deteX = N1 * X1;
double deteY = 0.5*(N1*Y1 + N2 * Y2 + By);
double deteZ = N1 * Z1;
//计算地面点的地面坐标(X,Y,Z)
double X = Xs1 + deteX;
double Y = Ys1 + deteY;
double Z = Zs1 + deteZ;
cout.precision(11);
cout << "计算得出地面点的地面坐标分别是：" << endl
	<< "X=" << X<< endl
	<< "Y=" << Y << endl
	<< "Z=" << Z<< endl;
cout << "You have finished the work  !" << endl;
system("pause");
return 0;

}