解析空中三角测量程序设计

已于 2024-11-28 15:45:50 修改
阅读量1k 收藏 12
点赞数 45
分类专栏： 某油地理科专业复习参考 文章标签： 算法 人工智能 c++
于 2024-11-22 20:42:54 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/weixin_74942631/article/details/143823483
版权
  • 目的要求
  1. 通过程序的设计深入理解解析空中三角测量的整个过程
  2. 提高应用程序设计语言解决问题的能力。
  • 资料及用具

一台微机，一套调试程序所用的数据。

  • 实习内容及作业步骤

用VB（或者C语言）程序设计语言编写单航带区域网概算的程序，具有包含一些三方面的内容：

  • （一）连续法像对的相对定向

  目标：求解各模型点在各模型像空间辅助坐标系中的坐标（Xi,Yi,Zi）

  步骤：

  1. 用连续法相对定向求解相对定向元素（bx，by，bz，φ，ω，κ）。
  2. 用前方交会公式计算各模型点坐标。
  • （二）航带网模型的建立

   目标：求出航带内各模型点在航带统一坐标系（第一个像片的像空间坐标系）中的坐标（Xi‘,Yi‘,Zi‘）。

   步骤：

  1. 根据模型间的三个公共点求解个模型的缩放系数ki
  2. 通过各模型缩放系数ki将模型内各模型点转换到航带统一坐标系。
  • （三）模型的绝对定向

  目标：求出航带内各模型点在大地坐标系中的坐标（Xti,Yti,Zti）。

  步骤：

  1. 根据两个已知控制点确定大地坐标系与地面摄测坐标系之间的转换参数（a、b）。
  2. 将所有控制点的大地坐标（Xt,Yt,Zt）通过转换参数（a、b）转换为地面摄测坐标系（Xtp,Ytp,Ztp）。
  3. 根据已知控制点的地面摄测坐标（Xtp,Ytp,Ztp）求解绝对定向元素（ΔX，ΔY，ΔZ，Φ，Ω，Κ，λ）。
  4. 根据绝对定向元素（ΔX，ΔY，ΔZ，Φ，Ω，Κ，λ）将各模型点在航带统一坐标系中的坐标（Xi‘,Yi‘,Zi‘）转换为地面摄测坐标系（Xtp,Ytp,Ztp）。
  5. 将各模型点的地面摄测坐标（Xtp,Ytp,Ztp）通过转换参数（a、b）转换为所有控制点的大地坐标（Xt,Yt,Zt）。
  • 上交成果

调试完成的程序一份及计算结果。

附：调试数据：

  1. 基本数据

摄影机主距：f＝153.033mm

                     bx＝200mm

  1. 像对坐标数据（单位：微米）

像对1：    701         702   

    1201      648790      735260     -230980      788550 

    1818      113860      595800     -771380      679190

    1202      241050      586260     -644640      661950

    1204      578030      223960     -327700      277590

    1203      256140      187820     -648360      260160  

    1205      606230     -104550     -315660      -51440

    1206      278340     -565020     -660020     -487200

    1052      179220     -757800     -765430     -670400

     410      478510     -637940     -466930     -569840

     399      975930     -700850       16690     -658940

     192      917380      -16160       -4600       18540

     370      803150      758570      -76440      802960

    1118       94870      709670     -785850      795830

     194       35030      -34550     -877140       50930

-398       19790     -843710     -925660     -745370

像对2：     703         702   

     400     -568240     -631500      426790     -634450

    1207     -601170      642970      323920      637060

     399     -980300     -630600       16690     -658940

     192     -963100       51030       -4600       18540

     370     -989450      836700      -76410      802880

    1209       -8790      641530      921460      679700------20

     401      -29060     -915900      981740     -884160

-190        1460        5490      965780       38900

像对3：    703         704   

    1826      931930      729560      -26020      680090

     402      907100     -785750       20330     -836630

    1055      753660     -838360     -134160     -896670

     411      397770     -725220     -500070     -791380

    1211       49010       50160     -875490      -17320

     188      918060       48330      -10150       11950

    1225      890340      544420      -58540      499270

    1210      624000      444520     -317940      392160

    1209       -8690      641540     -945780      560530---------31

     401      -29020     -915940     -930070    -1004070

    -190        1460        5460     -921980      -63410

3、控制点数据

点号        X坐标          Y坐标            Z坐标

==================================================

 1201      24204.689      46604.652         46.251

 1205      24343.683      45111.263         48.198

 1206      24965.047      44309.253         49.340

 1209      22167.904      46432.515         46.457

4、检查点数据

点号        X坐标          Y坐标          Z坐标

==================================================

 1118      25192.533      46608.059         48.936

 1202      24941.046      46375.998         46.539

 1203      24945.705      45662.638         49.052

 1204      24369.486      45700.421         49.020

 1207      23218.292      46347.142         48.993

注意：上述检查点数据、控制点数据、像对数据可能每次的不一样，本程序只适用上述数据的特定情况，并不能直接套用进行计算，但具体原理在代码中已注释，非全原创，但步骤详细清晰，均按照PPT公式进行。需要PPT的小伙伴可以在评论区留言，获取一整套摄影测量复习PPT。

代码程序（C语言）：

#define _CRT_SECURE_NO_WARNINGS 1
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

// 定义必要的常量和数据结构
#define f 0.153033 // 摄影机主距，单位：米
#define bx 0.200              // 基线，单位：米

typedef struct//输入 像点坐标
{
	double index;
	double x1;
	double y1;
	double x2;
	double y2;
}data1;

typedef struct      //输入 控制点坐标
{
	double index;
	double X;
	double Y;
	double Z;
}data3;

typedef struct  //相对定向求得的像空辅
{
	double index;
	double X1;
	double Y1;
	double Z1;
	double X2;
	double Y2;
	double Z2;
}data2;

typedef struct//像控点和像对的对应
{
	double index;
	double x;//左片的相同点
	double y;
	double z;
	int dex;
}data4;

typedef struct//相对定向后利用前方交会公式求得的模型点坐标
{
	double index;
	double X1;
	double Y1;
	double Z1;
}data5;


typedef struct //检查点数据
{
	double index;
	double X_;
	double Y_;
	double Z_;
} data6;

void matrix_multiply(double* A, int A_rows, int A_cols, double* B, int B_rows, int B_cols, double* result)
{ // A的列数必须等于B的行数
	if (A_cols != B_rows) {
		printf("Matrix dimensions do not match for multiplication!\n");
		return;
	}

	// 矩阵乘法计算
	for (int i = 0; i < A_rows; i++) {
		for (int j = 0; j < B_cols; j++) {
			result[i * B_cols + j] = 0; // 初始化结果矩阵元素
			for (int k = 0; k < A_cols; k++) {
				result[i * B_cols + j] += A[i * A_cols + k] * B[k * B_cols + j];
			}
		}
	}
}

#pragma region InverseMatrix

int InverseMatrix(double* B, int n, double* inverse) {
	// 创建增广矩阵，将 B 扩展成 [B | I]
	double* aug = (double*)malloc(n * 2 * n * sizeof(double));
	if (aug == NULL) {
		printf("内存分配失败\n");
		return 0;
	}

	// 构造增广矩阵
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			aug[i * 2 * n + j] = B[i * n + j];  // 左边是矩阵 B
			aug[i * 2 * n + j + n] = (i == j) ? 1.0 : 0.0;  // 右边是单位矩阵 I
		}
	}

	// 高斯-约当消元法
	for (int i = 0; i < n; i++) {
		// 找到当前列中绝对值最大的元素并交换行
		double maxElement = fabs(aug[i * 2 * n + i]);
		int maxRow = i;
		for (int k = i + 1; k < n; k++) {
			if (fabs(aug[k * 2 * n + i]) > maxElement) {
				maxElement = fabs(aug[k * 2 * n + i]);
				maxRow = k;
			}
		}

		// 如果无法求逆（矩阵不可逆），返回错误代码
		if (maxElement == 0.0) {
			printf("矩阵不可逆\n");
			free(aug);
			return 0;
		}

		// 交换当前行与最大元素所在行
		if (maxRow != i) {
			for (int k = 0; k < 2 * n; k++) {
				double temp = aug[i * 2 * n + k];
				aug[i * 2 * n + k] = aug[maxRow * 2 * n + k];
				aug[maxRow * 2 * n + k] = temp;
			}
		}

		// 归一化当前行，使对角线元素变为 1
		double diag = aug[i * 2 * n + i];
		for (int k = 0; k < 2 * n; k++) {
			aug[i * 2 * n + k] /= diag;
		}

		// 消去其他行的当前列
		for (int j = 0; j < n; j++) {
			if (j != i) {
				double factor = aug[j * 2 * n + i];
				for (int k = 0; k < 2 * n; k++) {
					aug[j * 2 * n + k] -= factor * aug[i * 2 * n + k];
				}
			}
		}
	}

	// 提取右侧的单位矩阵部分作为逆矩阵
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			inverse[i * n + j] = aug[i * 2 * n + j + n];
		}
	}

	// 释放增广矩阵内存
	free(aug);
	return 1;  // 成功返回 1
}

#pragma endregion

int main()
{
#pragma region import data
	data1 ImageData[34];//共有34对像点 15/8/11

	//导入像点坐标
	FILE* file = fopen("./像点数据.csv", "r"); // 打开文件
	if (!file) {
		printf("无法打开文件\n");
		return 1;
	}

	char buffer[100];
	fgets(buffer, sizeof(buffer), file);// 跳过第一行（标题行）

	int i = 0;
	while (fscanf(file, "%lf,%lf,%lf,%lf,%lf",
		&ImageData[i].index, &ImageData[i].x1, &ImageData[i].y1, &ImageData[i].x2, &ImageData[i].y2) == 5) {
		ImageData[i].x1 = ImageData[i].x1 / 1000000;
		ImageData[i].y1 = ImageData[i].y1 / 1000000;
		ImageData[i].x2 = ImageData[i].x2 / 1000000;
		ImageData[i].y2 = ImageData[i].y2 / 1000000;
		i++;
		if (i >= 34) break; // 确保不超过数组的大小
	}

	fclose(file);

	for (int i = 0; i < 34; i++) {	 //输出读取的数据验证结果
		printf("点号: %lf, 左片x: %lf, 左片y: %lf, 右片x: %lf, 右片y: %lf\n",
			ImageData[i].index, ImageData[i].x1, ImageData[i].y1, ImageData[i].x2, ImageData[i].y2);
	}

	//导入控制点坐标
	data3 ControlData[4];//共有4组控制点

	FILE* file3 = fopen("./控制点数据.csv", "r"); // 打开文件
	if (!file3) {
		printf("无法打开文件\n");
		return 1;
	}

	char buffer3[100];
	fgets(buffer3, sizeof(buffer3), file3);// 跳过第一行（标题行）

	int j = 0;
	while (fscanf(file3, "%lf,%lf,%lf,%lf",
		&ControlData[j].index, &ControlData[j].X, &ControlData[j].Y, &ControlData[j].Z) == 4) {
		j++;
		if (j >= 4) break; // 确保不超过数组的大小
	}

	fclose(file3);
	printf("\n");

	for (int j = 0; j < 4; j++) {	 //输出读取的数据验证结果
		printf("点号: %lf, 控制点X: %lf, 控制点Y: %lf, 控制点Z: %lf\n",
			ControlData[j].index, ControlData[j].X, ControlData[j].Y, ControlData[j].Z);
	}

	//导入检查点坐标
	data6 CheckData[5];//共有4组控制点

	FILE* file6 = fopen("./检查点数据.csv", "r"); // 打开文件
	if (!file6) {
		printf("无法打开文件\n");
		return 1;
	}

	char buffer6[100];
	fgets(buffer6, sizeof(buffer6), file6);// 跳过第一行（标题行）

	int k = 0;
	while (fscanf(file6, "%lf,%lf,%lf,%lf",
		&CheckData[k].index, &CheckData[k].X_, &CheckData[k].Y_, &CheckData[k].Z_) == 4) {
		k++;
		if (k >= 5) break; // 确保不超过数组的大小
	}

	fclose(file6);
	printf("\n");
	for (int k = 0; k < 5; k++) {	 //输出读取的数据验证结果
		printf("点号: %lf, 检查点X: %lf, 检查点Y: %lf, 检查点Z: %lf\n",
			CheckData[k].index, CheckData[k].X_, CheckData[k].Y_, CheckData[k].Z_);
	}
#pragma endregion

#pragma region 读取文件结束，开始相对定向----第一个像对

	double result1[5] = { 0 }, result2[5] = { 0 }, result3[5] = { 0 };
	int i1, i2, i3;
	//投影系数
	double N11[15], N21[15];
	double N12[8], N22[8];
	double N13[11], N23[11];

	//相对定向元素
	double μ1 = 0, γ1 = 0, φ21 = 0, ω21 = 0, κ21 = 0;
	double μ2 = 0, γ2 = 0, φ22 = 0, ω22 = 0, κ22 = 0;
	double μ3 = 0, γ3 = 0, φ23 = 0, ω23 = 0, κ23 = 0;
	double a4, a5, a6, b4, b5, b6, c4, c5, c6;

	data2 ImageData1[15]; data2 ImageData2[8]; data2 ImageData3[11];

	double by1, bz1;
	double by2, bz2;
	double by3, bz3;

	double A1[15][5];//误差方程——矩阵系数
	double B1[5][5];
	double R11[5][15];//误差方程——A的转置矩阵
	double A11[5][15];//误差方程——
	double Q1[15];
	double BB1[5][5];//误差方程——B的逆矩阵

	double A2[8][5];//误差方程——矩阵系数
	double B2[5][5];//误差方程——过渡矩阵
	double R12[5][8];//误差方程——A的转置矩阵
	double A12[5][8];//误差方程——
	double Q2[8];
	double BB2[5][5];//误差方程——B的逆矩阵

	double A3[11][5];
	double B3[5][5];
	double R13[5][11];
	double A13[5][11];
	double Q3[11];
	double BB3[5][5];

	for (i1 = 1;; i1++)//第一个文件相对定向元素
	{
		//初始化相对定向元素
		μ1 += result1[0];
		γ1 += result1[1];
		φ21 += result1[2];
		ω21 += result1[3];
		κ21 += result1[4];
		by1 = bx * μ1;
		bz1 = bx * γ1;

		//旋转矩阵的方向元素
		double
			a1 = cos(φ21) * cos(κ21) - sin(φ21) * sin(ω21) * sin(κ21),
			a2 = -cos(φ21) * sin(κ21) - sin(φ21) * sin(ω21) * cos(κ21),
			a3 = -sin(φ21) * cos(ω21),
			b1 = cos(ω21) * sin(κ21),
			b2 = cos(ω21) * cos(κ21),
			b3 = -sin(ω21),
			c1 = sin(φ21) * cos(κ21) + cos(ω21) * sin(ω21) * sin(κ21),
			c2 = -sin(φ21) * sin(κ21) + cos(ω21) * sin(ω21) * cos(κ21),
			c3 = cos(φ21) * cos(ω21);

		for (int i = 0; i < 15; i++)
		{
			//像点的像空间辅助坐标
			ImageData1[i].index = ImageData[i].index;
			ImageData1[i].X1 = ImageData[i].x1;//R1=E
			ImageData1[i].Y1 = ImageData[i].y1;
			ImageData1[i].Z1 = -f;
			ImageData1[i].X2 = a1 * ImageData[i].x2 + a2 * ImageData[i].y2 + a3 * (-f);//R2
			ImageData1[i].Y2 = b1 * ImageData[i].x2 + b2 * ImageData[i].y2 + b3 * (-f);
			ImageData1[i].Z2 = c1 * ImageData[i].x2 + c2 * ImageData[i].y2 + c3 * (-f);

			//计算投影系数
			N11[i] = (bx * ImageData1[i].Z2 - bz1 * ImageData1[i].X2) / (ImageData1[i].X1 * ImageData1[i].Z2 - ImageData1[i].X2 * ImageData1[i].Z1);
			N21[i] = (bx * ImageData1[i].Z1 - bz1 * ImageData1[i].X1) / (ImageData1[i].X1 * ImageData1[i].Z2 - ImageData1[i].X2 * ImageData1[i].Z1);

			//计算线性化矩阵A（系数）
			A1[i][0] = bx;
			A1[i][1] = -(ImageData1[i].Y2 / ImageData1[i].Z2) * bx;
			A1[i][2] = -(ImageData1[i].X2 * ImageData1[i].Y2 * N21[i]) / ImageData1[i].Z2;
			A1[i][3] = -(ImageData1[i].Z2 + (ImageData1[i].Y2 * ImageData1[i].Y2) / ImageData1[i].Z2) * N21[i];
			A1[i][4] = ImageData1[i].X2 * N21[i];

			//计算上下视差矩阵Q
			Q1[i] = N11[i] * ImageData1[i].Y1 - N21[i] * ImageData1[i].Y2 - by1;
		}

		//构建误差方程及法方程，最小二乘法求相对定向元素
		//计算转置矩阵
		for (int i = 0; i < 15; i++)
		{
			for (int j = 0; j < 5; j++)
			{
				R11[j][i] = A1[i][j];
			}
		}

		matrix_multiply(&R11[0][0],5,15,&A1[0][0],15,5, &B1[0][0]);//得到B
		InverseMatrix(&B1[0][0],5, &BB1[0][0]);

		//计算(ATA）-1AT
		for (int i = 0; i < 5; i++)
		{
			for (int m = 0; m < 15; m++)
			{
				A11[i][m] = 0;
				for (int j = 0; j < 5; j++)
				{
					A11[i][m] += BB1[i][j] * R11[j][m];
				}
			}
		}

		for (int i = 0; i < 5; i++)
		{
			result1[i] = 0;
			for (int j = 0; j < 15; j++)
			{
				result1[i] += A11[i][j] * Q1[j];
			}

		}//得到改正数

		if (fabs(result1[0]) < 0.000001 && fabs(result1[1]) < 0.000001 && fabs(result1[2]) < 0.000001 && fabs(result1[3]) < 0.000001 && fabs(result1[4]) < 0.000001)
		{
			φ21 += result1[2];
			ω21 += result1[3];
			κ21 += result1[4];
			a1 = cos(φ21) * cos(κ21) - sin(φ21) * sin(ω21) * sin(κ21),
			a2 = -cos(φ21) * sin(κ21) - sin(φ21) * sin(ω21) * cos(κ21),
			a3 = -sin(φ21) * cos(ω21),
			b1 = cos(ω21) * sin(κ21),
			b2 = cos(ω21) * cos(κ21),
			b3 = -sin(ω21),
			c1 = sin(φ21) * cos(κ21) + cos(ω21) * sin(ω21) * sin(κ21),
			c2 = -sin(φ21) * sin(κ21) + cos(ω21) * sin(ω21) * cos(κ21),
			c3 = cos(φ21) * cos(ω21);

			a4 = a1;
			a5 = a2;
			a6 = a3;
			b4 = b1;
			b5 = b2;
			b6 = b3;
			c4 = c1;
			c5 = c2;
			c6 = c3;
			break;
		}
	}

	printf("\n第一个像对的相对定向结果\n");
	printf("μ=%lf  γ=%1f  φ2=%1f  ω2=%1f  κ2=%1f  by=%1f  bz=%1f  循环次数=%d", μ1, γ1, φ21, ω21, κ21, by1, bz1, i1);
#pragma endregion

#pragma region 第二个像对的相对定向
	for (i2 = 1;; i2++)//第二个文件相对定向元素
	{
		//初始化相对定向元素
		μ2 += result2[0];
		γ2 += result2[1];
		φ22 += result2[2];
		ω22 += result2[3];
		κ22 += result2[4];
		by2 = bx * μ2;
		bz2 = bx * γ2;

		//旋转矩阵的方向元素
		double
			a1 = cos(φ22) * cos(κ22) - sin(φ22) * sin(ω22) * sin(κ22),
			a2 = -cos(φ22) * sin(κ22) - sin(φ22) * sin(ω22) * cos(κ22),
			a3 = -sin(φ22) * cos(ω22),
			b1 = cos(ω22) * sin(κ22),
			b2 = cos(ω22) * cos(κ22),
			b3 = -sin(ω22),
			c1 = sin(φ22) * cos(κ22) + cos(ω22) * sin(ω22) * sin(κ22),
			c2 = -sin(φ22) * sin(κ22) + cos(ω22) * sin(ω22) * cos(κ22),
			c3 = cos(φ22) * cos(ω22);

		for (int i = 0; i < 8; i++)
		{
			//像点的像空辅坐标系
			ImageData2[i].index = ImageData[i+15].index;
			ImageData2[i].X1 = a4 * ImageData[i + 15].x1 + a5 * ImageData[i + 15].y1 + a6 * (-f);
			ImageData2[i].Y1 = b4 * ImageData[i + 15].x1 + b5 * ImageData[i + 15].y1 + b6 * (-f);
			ImageData2[i].Z1 = c4 * ImageData[i + 15].x1 + c5 * ImageData[i + 15].y1 + c6 * (-f);
			ImageData2[i].X2 = a1 * ImageData[i + 15].x2 + a2 * ImageData[i + 15].y2 + a3 * (-f);//R2
			ImageData2[i].Y2 = b1 * ImageData[i + 15].x2 + b2 * ImageData[i + 15].y2 + b3 * (-f);
			ImageData2[i].Z2 = c1 * ImageData[i + 15].x2 + c2 * ImageData[i + 15].y2 + c3 * (-f);

			//计算投影系数
			N12[i] = (bx * ImageData2[i].Z2 - bz2 * ImageData2[i].X2) / (ImageData2[i].X1 * ImageData2[i].Z2 - ImageData2[i].X2 * ImageData2[i].Z1);
			N22[i] = (bx * ImageData2[i].Z1 - bz2 * ImageData2[i].X1) / (ImageData2[i].X1 * ImageData2[i].Z2 - ImageData2[i].X2 * ImageData2[i].Z1);

			//计算线性化矩阵A（系数）
			A2[i][0] = bx;
			A2[i][1] = -(ImageData2[i].Y2 / ImageData2[i].Z2) * bx;
			A2[i][2] = -(ImageData2[i].X2 * ImageData2[i].Y2 * N22[i]) / ImageData2[i].Z2;
			A2[i][3] = -(ImageData2[i].Z2 + (ImageData2[i].Y2 * ImageData2[i].Y2) / ImageData2[i].Z2) * N22[i];
			A2[i][4] = ImageData2[i].X2 * N22[i];

			//计算上下视差矩阵Q
			Q2[i] = N12[i] * ImageData2[i].Y1 - N22[i] * ImageData2[i].Y2 - by2;
		}

		//构建误差方程及法方程，最小二乘法求相对定向元素
		//计算转置矩阵
		for (int i = 0; i < 8; i++)
		{
			for (int j = 0; j < 5; j++)
			{
				R12[j][i] = A2[i][j];
			}
		}

		matrix_multiply(&R12[0][0], 5, 8, &A2[0][0], 8, 5, &B2[0][0]);//得到B
		InverseMatrix(&B2[0][0], 5, &BB2[0][0]);

		//计算(ATA）-1AT
		for (int i = 0; i < 5; i++)
		{
			for (int m = 0; m < 8; m++)
			{
				A12[i][m] = 0;
				for (int j = 0; j < 5; j++)
				{
					A12[i][m] += BB2[i][j] * R12[j][m];
				}
			}
		}

		for (int i = 0; i < 5; i++)
		{
			result2[i] = 0;
			for (int j = 0; j < 8; j++)
			{
				result2[i] += A12[i][j] * Q2[j];
			}

		}//得到改正数

		if (fabs(result2[0]) < 0.000001 && fabs(result2[1]) < 0.000001 && fabs(result2[2]) < 0.000001 && fabs(result2[3]) < 0.000001 && fabs(result2[4]) < 0.000001)
		{
			φ22 += result2[2];
			ω22 += result2[3];
			κ22 += result2[4];
			a1 = cos(φ22) * cos(κ22) - sin(φ22) * sin(ω22) * sin(κ22),
				a2 = -cos(φ22) * sin(κ22) - sin(φ22) * sin(ω22) * cos(κ22),
				a3 = -sin(φ22) * cos(ω22),
				b1 = cos(ω22) * sin(κ22),
				b2 = cos(ω22) * cos(κ22),
				b3 = -sin(ω22),
				c1 = sin(φ22) * cos(κ22) + cos(ω22) * sin(ω22) * sin(κ22),
				c2 = -sin(φ22) * sin(κ22) + cos(ω22) * sin(ω22) * cos(κ22),
				c3 = cos(φ22) * cos(ω22);

			a4 = a1;
			a5 = a2;
			a6 = a3;
			b4 = b1;
			b5 = b2;
			b6 = b3;
			c4 = c1;
			c5 = c2;
			c6 = c3;
			break;
		}
	}

	printf("\n第二个像对的相对定向结果\n");
	printf("μ=%lf  γ=%1f  φ2=%1f  ω2=%1f  κ2=%1f  by=%1f  bz=%1f  循环次数=%d", μ2, γ2, φ22, ω22, κ22, by2, bz2, i2);
#pragma endregion

#pragma region 第三个像对的相对定向
	for (i3 = 1;; i3++)//第三个文件相对定向元素
	{
		//初始化相对定向元素
		μ3 += result3[0];
		γ3 += result3[1];
		φ23 += result3[2];
		ω23 += result3[3];
		κ23 += result3[4];
		by3 = bx * μ3;
		bz3 = bx * γ3;

		//旋转矩阵的方向元素
		double
			a1 = cos(φ23) * cos(κ23) - sin(φ23) * sin(ω23) * sin(κ23),
			a2 = -cos(φ23) * sin(κ23) - sin(φ23) * sin(ω23) * cos(κ23),
			a3 = -sin(φ23) * cos(ω23),
			b1 = cos(ω23) * sin(κ23),
			b2 = cos(ω23) * cos(κ23),
			b3 = -sin(ω23),
			c1 = sin(φ23) * cos(κ23) + cos(ω23) * sin(ω23) * sin(κ23),
			c2 = -sin(φ23) * sin(κ23) + cos(ω23) * sin(ω23) * cos(κ23),
			c3 = cos(φ23) * cos(ω23);

		for (int i = 0; i < 11; i++)
		{
			//像点的像空间辅助坐标系
			ImageData3[i].index = ImageData[i + 23].index;
			ImageData3[i].X1 = a4 * ImageData[i + 23].x1 + a5 * ImageData[i + 23].y1 + a6 * (-f);
			ImageData3[i].Y1 = b4 * ImageData[i + 23].x1 + b5 * ImageData[i + 23].y1 + b6 * (-f);
			ImageData3[i].Z1 = c4 * ImageData[i + 23].x1 + c5 * ImageData[i + 23].y1 + c6 * (-f);
			ImageData3[i].X2 = a1 * ImageData[i + 23].x2 + a2 * ImageData[i + 23].y2 + a3 * (-f);//R2
			ImageData3[i].Y2 = b1 * ImageData[i + 23].x2 + b2 * ImageData[i + 23].y2 + b3 * (-f);
			ImageData3[i].Z2 = c1 * ImageData[i + 23].x2 + c2 * ImageData[i + 23].y2 + c3 * (-f);

			//计算投影系数
			N13[i] = (bx * ImageData3[i].Z2 - bz3 * ImageData3[i].X2) / (ImageData3[i].X1 * ImageData3[i].Z2 - ImageData3[i].X2 * ImageData3[i].Z1);
			N23[i] = (bx * ImageData3[i].Z1 - bz3 * ImageData3[i].X1) / (ImageData3[i].X1 * ImageData3[i].Z2 - ImageData3[i].X2 * ImageData3[i].Z1);

			//计算线性化矩阵A（系数）
			A3[i][0] = bx;
			A3[i][1] = -(ImageData3[i].Y2 / ImageData3[i].Z2) * bx;
			A3[i][2] = -(ImageData3[i].X2 * ImageData3[i].Y2 * N23[i]) / ImageData3[i].Z2;
			A3[i][3] = -(ImageData3[i].Z2 + (ImageData3[i].Y2 * ImageData3[i].Y2) / ImageData3[i].Z2) * N23[i];
			A3[i][4] = ImageData3[i].X2 * N23[i];

			//计算上下视差矩阵Q
			Q3[i] = N13[i] * ImageData3[i].Y1 - N23[i] * ImageData3[i].Y2 - by3;
		}

		//构建误差方程及法方程，最小二乘法求相对定向元素
		//计算转置矩阵
		for (int i = 0; i < 11; i++)
		{
			for (int j = 0; j < 5; j++)
			{
				R13[j][i] = A3[i][j];
			}
		}

		matrix_multiply(&R13[0][0], 5, 11, &A3[0][0], 11, 5, &B3[0][0]);//得到B
		InverseMatrix(&B3[0][0], 5, &BB3[0][0]);

		//计算(ATA）-1AT
		for (int i = 0; i < 5; i++)
		{
			for (int m = 0; m < 11; m++)
			{
				A13[i][m] = 0;
				for (int j = 0; j < 5; j++)
				{
					A13[i][m] += BB3[i][j] * R13[j][m];
				}
			}
		}

		for (int i = 0; i < 5; i++)
		{
			result3[i] = 0;
			for (int j = 0; j < 11; j++)
			{
				result3[i] += A13[i][j] * Q3[j];
			}

		}//得到改正数

		if (fabs(result3[0]) < 0.000001 && fabs(result3[1]) < 0.000001 && fabs(result3[2]) < 0.000001 && fabs(result3[3]) < 0.000001 && fabs(result3[4]) < 0.000001)
		{
			φ23 += result3[2];
			ω23 += result3[3];
			κ23 += result3[4];
			break;
		}
	}

	printf("\n第三个像对的相对定向结果\n");
	printf("μ=%lf  γ=%1f  φ2=%1f  ω2=%1f  κ2=%1f  by=%1f  bz=%1f  循环次数=%d", μ3, γ3, φ23, ω23, κ23, by3, bz3, i3);
#pragma endregion

#pragma region 求比例尺m
	double distance1[3], distance2[3];
	distance1[0] = sqrt((ImageData[0].x1 - ImageData[5].x1) * (ImageData[0].x1 - ImageData[5].x1) + (ImageData[0].y1 - ImageData[5].y1) * (ImageData[0].y1 - ImageData[5].y1));
	distance1[1] = sqrt((ImageData[0].x1 - ImageData[6].x1) * (ImageData[0].x1 - ImageData[6].x1) + (ImageData[0].y1 - ImageData[6].y1) * (ImageData[0].y1 - ImageData[6].y1));
	distance1[2] = sqrt((ImageData[5].x1 - ImageData[6].x1) * (ImageData[5].x1 - ImageData[6].x1) + (ImageData[5].y1 - ImageData[6].y1) * (ImageData[5].y1 - ImageData[6].y1));

	distance2[0] = sqrt((ControlData[0].X - ControlData[1].X) * (ControlData[0].X - ControlData[1].X) + (ControlData[0].Y - ControlData[1].Y) * (ControlData[0].Y - ControlData[1].Y));
	distance2[1] = sqrt((ControlData[0].X - ControlData[2].X) * (ControlData[0].X - ControlData[2].X) + (ControlData[0].Y - ControlData[2].Y) * (ControlData[0].Y - ControlData[2].Y));
	distance2[2] = sqrt((ControlData[2].X - ControlData[1].X) * (ControlData[2].X - ControlData[1].X) + (ControlData[2].Y - ControlData[1].Y) * (ControlData[2].Y - ControlData[1].Y));

	double M[3], m;
	M[0] = distance2[0] / distance1[1];
	M[1] = distance2[1] / distance1[1];
	M[2] = distance2[2] / distance1[2];

	m = (M[0] + M[1] + M[2]) / 3;

	printf("\n");
	printf("比例尺：m=%lf ", m);

#pragma endregion

#pragma region 航带网模型的建立
	//上述已求得像点的像空间辅助坐标系
	// 此时各模型的像空间辅助坐标系坐标轴向都保持平行
	//但模型比，例尺各不相同，坐标原点也不一致

	data5 ModelspaceAid1[15];//模型点的像空间辅助坐标系
	data5 ModelspaceAid2[8];
	data5 ModelspaceAid3[11];

	for (int i = 0; i < 15; i++)
	{
		ModelspaceAid1[i].index = ImageData1[i].index;
		ModelspaceAid1[i].X1 = ImageData1[i].X1 * N11[i];
		ModelspaceAid1[i].Y1 = (ImageData1[i].Y1 * N11[i] + ImageData1[i].Y2 * N21[i] + by1) / 2;
		ModelspaceAid1[i].Z1 = ImageData1[i].Z1 * N11[i];
	}

	for (int i = 0; i < 8; i++)
	{
		ModelspaceAid2[i].index = ImageData2[i].index;
		ModelspaceAid2[i].X1 = ImageData2[i].X1 * N12[i];
		ModelspaceAid2[i].Y1 = (ImageData2[i].Y1 * N12[i] + ImageData2[i].Y2 * N22[i] + by2) / 2;
		ModelspaceAid2[i].Z1 = ImageData2[i].Z1 * N12[i];
	}

	for (int i = 0; i < 11; i++)
	{
		ModelspaceAid3[i].index = ImageData3[i].index;
		ModelspaceAid3[i].X1 = ImageData3[i].X1 * N13[i];
		ModelspaceAid3[i].Y1 = (ImageData3[i].Y1 * N13[i] + ImageData3[i].Y2 * N23[i] + by3) / 2;
		ModelspaceAid3[i].Z1 = ImageData3[i].Z1 * N13[i];
	}

	//航带内各立体模型利用公共点进行连接，建立起统一的航带网模型。
	//若要将航带内所有模型连接构成航带网，必须要进行各模型之间比例尺的归算，即K1，K2
	double K1[3] = { 0 }, K2[3] = { 0 };
	int count1 = 0, count2 = 0;
	double k1 = 0, k2 = 0;
	for (int i = 0; i < 15; i++)
	{
		for (int j = 0; j < 8; j++)
		{
			if (ModelspaceAid1[i].index == ModelspaceAid2[j].index && count1 < 3)
			{
				K1[count1] = (ModelspaceAid1[i].Z1 - bz1) / (ModelspaceAid2[j].Z1);
				k1 += K1[count1];
				count1++;
			}
		}
	}
	k1 = k1 / 3;
	for (int i = 0; i < 8; i++)
	{
		for (int j = 0; j < 11; j++)
		{
			if (ModelspaceAid2[i].index == ModelspaceAid3[j].index && count2 < 3)
			{
				K2[count2] = (ModelspaceAid2[i].Z1 - bz2) / (ModelspaceAid3[j].Z1);
				k2 += K2[count2];
				count2++;
			}
		}
	}
	k2 = k2 / 3;

	printf("\n");
	printf("各模型之间的比例尺的归算：k1 = %1f k2 = %1f ", k1, k2);
	printf("\n");

	//航带内各单个模型建立以后，以相邻两模型重叠范围内三个连接点的高度应相等为条件，
	// 沿航带从左至右，逐个模型的归化比例尺，统一坐标原点，
	// 使全航带各模型连接成一个统一的自由航带网模型，将模型点的坐标换算到摄影测量坐标系中

	//此时需要使用K1、K2，继续计算模型点摄影测量坐标
	double Xps[4], Yps[4], Zps[4];
	Xps[0] = 0;
	Yps[0] = 0;
	Zps[0] = m * f;
	Xps[1] = Xps[0] + m * bx;
	Yps[1] = Yps[0] + m * by1;
	Zps[1] = Zps[0] + m * bz1;
	Xps[2] = Xps[1] + k1 * m * bx;
	Yps[2] = Yps[1] + k1 * m * by2;
	Zps[2] = Zps[1] + k1 * m * bz2;
	Xps[3] = Xps[2] + k1 * k2 * m * bx;
	Yps[3] = Yps[2] + k1 * k2 * m * by3;
	Zps[3] = Zps[2] + k1 * k2 * m * bz3;

	//计算最终的模型点坐标(摄影测量坐标系)
	data5 Photography1[15];
	data5 Photography2[8];
	data5 Photography3[11];

	for (int i = 0; i < 15; i++)
	{
		Photography1[i].index = ModelspaceAid1[i].index;
		Photography1[i].X1 = ModelspaceAid1[i].X1 * m;
		Photography1[i].Y1 = ModelspaceAid1[i].Y1 * m;
		Photography1[i].Z1 = Zps[0] + ModelspaceAid1[i].Z1 * m;
	}

	for (int i = 0; i < 8; i++)
	{
		Photography2[i].index = ModelspaceAid2[i].index;
		Photography2[i].X1 = Xps[1] + k1 * ModelspaceAid2[i].X1 * m;
		Photography2[i].Y1 = (Yps[1] + k1 * m * N12[i] * ImageData2[i].Y1 + Yps[2] + k1 * m * N22[i] * ImageData2[i].Y2) / 2;//公式有点问题？？？
		Photography2[i].Z1 = Zps[1] + k1 * ModelspaceAid2[i].Z1 * m;
	}

	for (int i = 0; i < 11; i++)
	{
		Photography3[i].index = ModelspaceAid3[i].index;
		Photography3[i].X1 = Xps[2] + k1 * k2 * ModelspaceAid3[i].X1 * m;
		Photography3[i].Y1 = (Yps[2] + k1 * k2 * m * N13[i] * ImageData3[i].Y1 + Yps[3] + k1 * k2 * m * N23[i] * ImageData3[i].Y2 ) / 2;
		Photography3[i].Z1 = Zps[2] + k1 * k2 * ModelspaceAid3[i].Z1 * m;
	}
#pragma endregion

#pragma region 绝对定向求解七参数
	//在航带网的两端，分别选定1和2两个控制点，根据这两点的地面测量坐标和摄影测量坐标，
	// 将测区内所有地面控制点的地面测量坐标和对应的摄影测量坐标都换算为以1点为坐标原点的坐标。
	// 同时，将自由航带网内所有加密点的摄影测量坐标也换算为以1点为原点的坐标。

	// 求dx dy
	double dXt, dYt, dXp, dYp;
	dXt = ControlData[0].X - ControlData[3].X;//控制点的大地坐标差
	dYt = ControlData[0].Y - ControlData[3].Y;
	dXp = Photography1[0].X1 - Photography3[8].X1;//摄影测量坐标差
	dYp = Photography1[0].Y1 - Photography3[8].Y1;

	//计算转换系数a b T
	double a, b, T;
	a = (dXp * dYt + dYp * dXt) / (dXt * dXt + dYt * dYt);
	b = (dXp * dXt - dYp * dYt) / (dXt * dXt + dYt * dYt);
	T = sqrt(a * a + b * b);

	printf("\n");
	printf("三个转换系数：%lf  %lf  %lf  ", a, b, T);

	//求得三参数后，将所有地面控制点的大地坐标系计算变换成的地面摄影测量坐标Xtp 
	data3 controlXt_to_Xtp[4];
	for (int i = 0; i < 4; i++)
	{
		controlXt_to_Xtp[i].X = b * (ControlData[i].X - ControlData[0].X) + (ControlData[i].Y - ControlData[0].Y) * a;
		controlXt_to_Xtp[i].Y = a * (ControlData[i].X - ControlData[0].X) + (ControlData[i].Y - ControlData[0].Y) * (-b);
		controlXt_to_Xtp[i].Z = T * (ControlData[i].Z - ControlData[0].Z);
	}

	//摄影测量坐标Xp
	data4 Photography1_controlXp[10] = { 0 };//与控制点对应的模型点摄影测量坐标
	int point4 = 0;
	for (int i = 0; i < 15; i++)
	{
		for (int j = 0; j < 4; j++)
		{
			if (Photography1[i].index == ControlData[j].index)
			{
				Photography1_controlXp[point4].index = Photography1[i].index;
				Photography1_controlXp[point4].x = Photography1[i].X1;
				Photography1_controlXp[point4].y = Photography1[i].Y1;
				Photography1_controlXp[point4].z = Photography1[i].Z1;
				point4++;
			}
		}
	}

	data4 Photography3_controlXp[10] = { 0 };
	int point5 = 0;
	for (int i = 0; i < 11; i++)
	{
		for (int j = 0; j < 4; j++)
		{
			if (Photography3[i].index == ControlData[j].index)
			{
				Photography3_controlXp[point5].index = Photography3[i].index;
				Photography3_controlXp[point5].x = Photography3[i].X1;
				Photography3_controlXp[point5].y = Photography3[i].Y1;
				Photography3_controlXp[point5].z = Photography3[i].Z1;
				point5++;
			}
		}
	}


	//计算地面摄影测量坐标和摄影测量坐标
	data3 tp[4], p[4];
	for (int i = 0; i < 4; i++)
	{
		tp[i].X = controlXt_to_Xtp[i].X;
		tp[i].Y = controlXt_to_Xtp[i].Y;
		tp[i].Z = controlXt_to_Xtp[i].Z;

	}

	for (int i = 0; i < 3; i++)
	{
		p[i].X = Photography1_controlXp[i].x;
		p[i].Y = Photography1_controlXp[i].y;
		p[i].Z = Photography1_controlXp[i].z;
	}
	p[3].X = Photography3_controlXp[0].x;
	p[3].Y = Photography3_controlXp[0].y;
	p[3].Z = Photography3_controlXp[0].z;

	//printf("\n------------------------\n");
	//for (int i = 0; i < 4; i++)
	//{
	//	printf("%lf  ", controlXt_to_Xtp[i].X);
	//}

	//计算绝对定向元素
	double T0 = 1, X0 = 0, Y0 = 0, Z0 = 0, φ = 0, ω = 0, κ = 0;//七参数
	double result4[7] = { 0 };//改正数
	double L[12] = { 0 };//4个控制点，Lx,Ly,Lz有三个常数项
	
	//误差方程与法方程
	double A[12][7] = { 0 };
	double AT[7][12];
	double ATA[7][7];
	double ATA1[7][7];
	double AA[7][12];

	for (int i5 = 1;; i5++)
	{
		//printf("%d ", i5);
		//初始化

		X0 += result4[0];
		Y0 += result4[1];
		Z0 += result4[2];
		T0 += result4[3];
		φ += result4[4];
		ω += result4[5];
		κ += result4[6];

		double
			a1 = cos(φ) * cos(κ) - sin(φ) * sin(ω) * sin(κ),
			a2 = -cos(φ) * sin(κ) - sin(φ) * sin(ω) * cos(κ),
			a3 = -sin(φ) * cos(ω),
			b1 = cos(ω) * sin(κ),
			b2 = cos(ω) * cos(κ),
			b3 = -sin(ω),
			c1 = sin(φ) * cos(κ) + cos(ω) * sin(ω) * sin(κ),
			c2 = -sin(φ) * sin(κ) + cos(ω) * sin(ω) * cos(κ),
			c3 = cos(φ) * cos(ω);

		for (int i = 0; i < 4; i++)
		{
			int j = i * 3;
			L[j] = tp[i].X - T0 * (a1 * p[i].X + a2 * p[i].Y + a3 * p[i].Z) - X0;
			L[j + 1] = tp[i].Y - T0 * (b1 * p[i].X + b2 * p[i].Y + b3 * p[i].Z) - Y0;
			L[j + 2] = tp[i].Z - T0 * (c1 * p[i].X + c2 * p[i].Y + c3 * p[i].Z) - Z0;
		}

		//计算A矩阵
		// 计算A矩阵
		for (int i = 0; i < 12; i++) {
			for (int j = 0; j < 7; j++) {
				int idx = i / 3;  // p数组的索引

				// 处理对角线元素
				if (i % 3 == j) {
					A[i][j] = 1;
				}
				// 处理特定列的赋值
				else if ((i == 0 || i == 3 || i == 6 || i == 9) && j == 3) {
					A[i][j] = p[idx].X;
				}
				else if ((i == 2 || i == 5 || i == 8 || i == 11) && j == 4) {
					A[i][j] = p[idx].X;
				}
				else if ((i == 1 || i == 4 || i == 7 || i == 10) && j == 6) {
					A[i][j] = p[idx].X;
				}
				else if ((i == 1 || i == 4 || i == 7 || i == 10) && j == 3) {
					A[i][j] = p[idx].Y;
				}
				else if ((i == 2 || i == 5 || i == 8 || i == 11) && j == 5) {
					A[i][j] = p[idx].Y;
				}
				else if ((i == 0 || i == 3 || i == 6 || i == 9) && j == 6) {
					A[i][j] = -p[idx].Y;
				}
				else if ((i == 2 || i == 5 || i == 8 || i == 11) && j == 3) {
					A[i][j] = p[idx].Z;
				}
				else if ((i == 1 || i == 4 || i == 7 || i == 10) && j == 5) {
					A[i][j] = -p[idx].Z;
				}
				else if ((i == 0 || i == 3 || i == 6 || i == 9) && j == 4) {
					A[i][j] = -p[idx].Z;
				}
			}
		}

		//计算A的转置
		for (int i = 0; i < 12; i++)
		{
			for (int j = 0; j < 7; j++)
			{
				AT[j][i] = A[i][j];
			}
		}

		for (int i = 0; i < 7; i++)
		{
			for (int m = 0; m < 7; m++)
			{
				ATA[i][m] = 0;
				for (int j = 0; j < 12; j++)
				{
					ATA[i][m] += AT[i][j] * A[j][m];
				}
			}
		}

		//求逆
		InverseMatrix(&ATA[0][0], 7, &ATA1[0][0]);

		//计算(ATA）-1AT
		for (int i = 0; i < 7; i++)
		{
			for (int m = 0; m < 12; m++)
			{
				AA[i][m] = 0;
				for (int j = 0; j < 7; j++)
				{
					AA[i][m] += ATA1[i][j] * AT[j][m];
				}
			}
		}

		//计算最终结果
		for (int i = 0; i < 7; i++)
		{
			result4[i] = 0;
			for (int j = 0; j < 12; j++)
			{
				result4[i] += AA[i][j] * L[j];
			}
		}

		if (fabs(result4[0]) < 0.0000001 && fabs(result4[1]) < 0.0000001 && fabs(result4[5]) < 0.0000001 && fabs(result4[2]) < 0.0000001 && fabs(result4[3]) < 0.0000001 && fabs(result4[4]) < 0.0000001 && fabs(result4[6]) < 0.0000001)
		{
			printf("\n");
			printf("绝对定向结束");
			break;
		}
	}
	printf("\n");
	printf("%lf  %lf  %lf  %lf  %lf  %lf  %lf  ", X0, Y0, Z0, T0, φ, ω, κ);

#pragma endregion


#pragma region 绝对定向摄影转地面摄影
	// 计算三个模型的摄影测量坐标——>地面摄影测量坐标
	data3 modeltp1[15], modeltp2[8], modeltp3[11];
	for (int i = 0; i < 15; i++)
	{
		modeltp1[i].index = Photography1[i].index;
		modeltp1[i].X = T0 * (Photography1[i].X1 - κ * Photography1[i].Y1 - φ * Photography1[i].Z1) + X0;
		modeltp1[i].Y = T0 * (κ * Photography1[i].X1 + Photography1[i].Y1 - ω * Photography1[i].Z1) + Y0;
		modeltp1[i].Z = T0 * (φ * Photography1[i].X1 + ω * Photography1[i].Y1 + Photography1[i].Z1) + Z0;
	}

	for (int i = 0; i < 8; i++)
	{
		modeltp2[i].index = Photography2[i].index;
		modeltp2[i].X = T0 * (Photography2[i].X1 - κ * Photography2[i].Y1 - φ * Photography2[i].Z1) + X0;
		modeltp2[i].Y = T0 * (κ * Photography2[i].X1 + Photography2[i].Y1 - ω * Photography2[i].Z1) + Y0;
		modeltp2[i].Z = T0 * (φ * Photography2[i].X1 + ω * Photography2[i].Y1 + Photography2[i].Z1) + Z0;
	}
	for (int i = 0; i < 11; i++)
	{
		modeltp3[i].index = Photography3[i].index;
		modeltp3[i].X = T0 * (Photography3[i].X1 - κ * Photography3[i].Y1 - φ * Photography3[i].Z1) + X0;
		modeltp3[i].Y = T0 * (κ * Photography3[i].X1 + Photography3[i].Y1 - ω * Photography3[i].Z1) + Y0;
		modeltp3[i].Z = T0 * (φ * Photography3[i].X1 + ω * Photography3[i].Y1 + Photography3[i].Z1) + Z0;
	}

	//计算最终的大地坐标
	data3 finish1[15], finish2[8], finish3[11];
	for (int i = 0; i < 15; i++)
	{
		finish1[i].index = modeltp1[i].index;
		finish1[i].X = (b * modeltp1[i].X + a * modeltp1[i].Y) / T / T + ControlData[0].X;
		finish1[i].Y = (a * modeltp1[i].X - b * modeltp1[i].Y) / T / T + ControlData[0].Y;
		finish1[i].Z = (modeltp1[i].Z) / T + ControlData[0].Z;
	}

	printf("\n---------------------------\n");
	for (int i = 0; i < 15; i++)
	{
		printf("%lf  ", finish1[i].X);
	}


	for (int i = 0; i < 8; i++)
	{
		finish2[i].index = modeltp2[i].index;
		finish2[i].X = (b * modeltp2[i].X + a * modeltp2[i].Y) / T / T + ControlData[0].X;
		finish2[i].Y = (a * modeltp2[i].X - b * modeltp2[i].Y) / T / T + ControlData[0].Y;
		finish2[i].Z = (modeltp2[i].Z) / T + ControlData[0].Z;
	}


	for (int i = 0; i < 11; i++)
	{
		finish3[i].index = modeltp3[i].index;
		finish3[i].X = (b * modeltp3[i].X + a * modeltp3[i].Y) / T / T + ControlData[0].X;
		finish3[i].Y = (a * modeltp3[i].X - b * modeltp3[i].Y) / T / T + ControlData[0].Y;
		finish3[i].Z = (modeltp3[i].Z) / T + ControlData[0].Z;
	}


	//检查点 计算差值	
	double  D[5][4] = { 0 };
	for (i = 0; i < 5; i++)
	{
		for (int j = 0; j < 15; j++)
		{
			if (CheckData[i].index == finish1[j].index)
			{
				D[i][0] = CheckData[i].index;
				D[i][1] = CheckData[i].X_ - finish1[j].X;
				D[i][2] = CheckData[i].Y_ - finish1[j].Y;
				D[i][3] = CheckData[i].Z_ - finish1[j].Z;
				break;
			}

		}
		for (int j = 0; j < 8; j++)
		{
			if (CheckData[i].index == finish2[j].index)
			{
				D[i][0] = CheckData[i].index;
				D[i][1] = CheckData[i].X_ - finish2[j].X;
				D[i][2] = CheckData[i].Y_ - finish2[j].Y;
				D[i][3] = CheckData[i].Z_ - finish2[j].Z;
				break;
			}

		}
		for (int j = 0; j < 11; j++)
		{
			if (CheckData[i].index == finish3[j].index)
			{
				D[i][0] = CheckData[i].index;
				D[i][1] = CheckData[i].X_ - finish3[j].X;
				D[i][2] = CheckData[i].Y_ - finish3[j].Y;
				D[i][3] = CheckData[i].Z_ - finish3[j].Z;
				break;
			}
		}
	}
#pragma endregion



#pragma region 写入文件
	//写入文件
	FILE* fp2;
	fp2 = fopen("result.txt", "w");
	fprintf(fp2, "****************相对定向元素*******************\n");
	fprintf(fp2, "\n");
	fprintf(fp2, " 三个像对最终的相对定向元素结果:\n ");
	fprintf(fp2, "\n");
	fprintf(fp2, "μ   γ   φ   ω   κ   by   bz  循环次数\n");
	fprintf(fp2, "第一个像对：\n");
	fprintf(fp2, " %10.8f  %f   %f  %f  %f  %f  %f   %d\n", μ1, γ1, φ21, ω21, κ21, by1, bz1, i1);
	fprintf(fp2, "\n");
	fprintf(fp2, "第一个像对：\n");
	fprintf(fp2, " %10.8f  %f   %f  %f  %f  %f  %f   %d\n", μ2, γ2, φ22, ω22, κ22, by2, bz2, i2);
	fprintf(fp2, "\n");
	fprintf(fp2, "第一个像对：\n");
	fprintf(fp2, " %10.8f  %f   %f  %f  %f  %f  %f   %d\n", μ3, γ3, φ23, ω23, κ23, by3, bz3, i3);
	fprintf(fp2, "\n");

	fprintf(fp2, "****************航带网的建立*******************\n");
	fprintf(fp2, "\n");

	fprintf(fp2, "比例尺m = %f    \n", m);
	fprintf(fp2, "\n");

	fprintf(fp2, "比例尺缩放系数： k1 = %f   k2 = %f  \n", k1, k2);
	fprintf(fp2, "\n");

	fprintf(fp2, "三个转换系数： a = %f  b = %f  T = %f\n", a, b, T);

	fprintf(fp2, "\n绝对定向元素  \n");
	fprintf(fp2, " T0 = %f   X0 = %f   Y0 = %f   Z0 = %f   φ = %f   ω = %f   κ = %f  \n", T0, X0, Y0, Z0, φ, ω, κ);

	fprintf(fp2, "\n");
	fprintf(fp2, "****************模型连接*******************\n");

	fprintf(fp2, "\n 最后的大地坐标 \n");
	for (int i = 0; i < 15; i++)
	{
		fprintf(fp2, " %f   %f    %f \n", finish1[i].X, finish1[i].Y, finish1[i].Z);
	}
	fprintf(fp2, "\n");
	for (int i = 0; i < 8; i++)
	{
		fprintf(fp2, " %f   %f    %f\n ", finish2[i].X, finish2[i].Y, finish2[i].Z);
	}
	fprintf(fp2, "\n");
	for (int i = 0; i < 11; i++)
	{
		fprintf(fp2, " %f   %f    %f \n", finish3[i].X, finish3[i].Y, finish3[i].Z);
	}
	fprintf(fp2, "\n 与待定点的差值 \n");
	for (int j = 0; j < 5; j++)
	{
		fprintf(fp2, " %f   %f    %f  \n", D[j][1], D[j][2], D[j][3]);
	}
	fprintf(fp2, "\n");
	fclose(fp2);
	//写入文件结束
#pragma endregion


	return 0;
}

注意：上述代码可以直接运行，只需要将数据导入CSV表格即可，格式如下：

最终结果：

附：其他版本的像点数据之间可能不一样，今年版本的数据稍有改动，但是不知道为什么检查点的误差在Z轴上这么大。

连续像对相对定向
01-02
利用连续像对相对定向的方法进行立体像对相对定向，得到像对方位元素
连续相对定向 matlab
04-08
连续法相对定向程序，使用matlab编写
Python 卖空算法教程（三）
龙哥盟
05-09 1072
原文：zh.annas-archive.org/md5/ceefdd89e585c59c20db6a7760dc11f1 译者：飞龙 协议：CC BY-NC-SA 4.0 第十章：优化投资范围 市场参与者通常觉得市场的广大相当令人生畏。因此，在我们开始将我们在前几章中涵盖的想法转化为结合的投资策略之前，我们将专门撰写一个简短的章节，将市场简化为可管理的投资范围。 在本章中，我们将从一些概念性的盲点开始，讨论长/短期业务，试图为你将要进行交易的现实世界提供一些有价值的背景信息。接下来，我们将跟随资金的脚步
摄影测量--相对定向-绝对定向（C++实现）
热门推荐
Jensen Lee的博客
09-06 2万+
一． 目的 理解并掌握空间相对定向-原理。 体会相对定向-绝对定向与前后方交会的异同，理解各参数含义。 熟悉计算流程，并通过编程运用到实践上。 提高编程计算能力，并将算式转换为程序，体会编程计算的特点。 二． 实验要求 提交实习报告：程序框图、程序源代码、计算结果、体会。 计算结果：相对定向五参数及精度评定、计算模型点坐标、绝对定向七参...
C#代码-摄影测量实习-解析空中三角测量
wuboyang0430的博客
11-24 2555
解析空中三角测量-连续法相对定向
解析空中三角测量程序设计.zip
06-03
本压缩包“解析空中三角测量程序设计.zip”包含了实现这一过程的相关代码和数据，旨在帮助用户理解并实践空中三角测量的算法。 首先，我们要理解空中三角测量的基本概念。它基于多张从不同角度拍摄的航片，通过匹配...
解析空中三角测量
09-24
空中三角测量是摄影测量学中的一个关键步骤，用于通过多张照片重建三维地理空间信息。在数字摄影测量中，它通常涉及使用计算机视觉算法来确定相机的位置和姿态，以及图像中特征点的空间坐标。本代码资源是用C#语言...
解析空中三角测量C#程序
11-18
在实现空中三角测量的C#程序时，我们首先需要理解以下几个核心概念： 1. **影像匹配**：这是整个过程的基础，通过特征检测（如SIFT、SURF等）找到不同照片之间的对应点，这些对应点提供了照片间重叠区域的信息。 2...
用c#编写了解析空中三角测量程序
08-14
5. **安装程序封装**：完成所有计算后，将程序打包成安装程序，通常会用到Visual Studio或其他类似工具，如Inno Setup，这样用户就可以在自己的计算机上轻松安装和运行空中三角测量程序。 6. **用户界面设计**：...
摄影测量解析三角测量程序
11-30
在本压缩包中，包含的是一个用于摄影测量实习的解析空中三角测量程序，它由C++和C语言编写，表明该程序具有高度的灵活性和效率。 解析三角测量的核心在于通过图像匹配找到同一场景在不同视点下的对应点，然后利用...
相对定向 摄影测量 C++完整能运行代码
11-30
相对定向 摄影测量 C++完整能运行代码，才写的，完整能跑，有问题问我
摄影测量连续法相对定向
10-27
摄影测量连 续法相对定向 测绘 遥感 摄影测量连 续法相对定向 测绘 遥感 摄影测量连 续法相对定向 测绘 遥感
用python实现相对定向
02-23
用python实现相对定向，里面有详细的注释和实验所用的数据，代码简洁易懂
相对定向代码
10-27
在摄影测量课程中一定要做的一节课，相对定向和绝对定向的代码实现matlab
连续法相对定向
06-14
摄影测量中的连续法相对定向，MATLAB实现，亲测没问题
测绘-空中三角测量程序设计
想上岸的鱼
06-05 7037
实习解析空中三角测量程序设计一、 目的要求1、 通过程序的设计深入理解解析空中三角测量的整个过程2、 提高应用程序设计语言解决问题的能力。二、 资料及用具一台微机，一套调试程序所用的数据。三、 实习内容及作业步骤用VB（或者C语言）程序设计语言编写单航带区域网概算的程序，具有包含一些三方面的内容：（一） 连续法像对的相对定向  目标：求解各模型点在各模型像空间辅助坐标系中的坐标（Xi,Yi,Zi）...
摄影测量简答题和作图题
09-08 5074
单后：获取已知数据；量测控制点的像点坐标；确定未知数的初始值；计算旋转矩阵R；逐点计算像点坐标的近似值；组成误差方程；组成法方程式；解外方位元素；检查计算是否收敛； 双前：获取已知数据；量测像点坐标；由像片的角元素计算各像片的旋转矩阵R1，R2；计算摄影基线分量；计算像点在左右像片上的像空间辅助坐标系；计算点投影系数N1，N2；计算待定点的地面摄影测量坐标； 相对定向步骤：获取已知数据X0，Y...
探索摄影测量连续法相对定向：测绘与遥感的利器
最新发布
gitblog_09726的博客
10-28 344
探索摄影测量连续法相对定向：测绘与遥感的利器 【下载地址】摄影测量连续法相对定向资源下载 本仓库提供了一个关于“摄影测量连续法相对定向”的资源文件下载。该资源文件详细介绍了摄影测量中的连续法相对定向技术，适用于测绘和遥感领域的专业人士和学生。通过本资源，您可以深入了解摄影测量连续法相对定向的基本原理、操作步骤以及实际应用...
空中三角测量加密原理及4D产品制作流程
u011322358的博客
02-13 1万+
空中三角测量加密原理及4D产品制作流程 1. 航空摄影测量技术 航空摄影测量技术作为空间信息技术体系的两大分支之一，是空间数据获取的重要工具之一。由于其运行成本低、执行任务灵活性高、安全性高、测量精度高等优点，在各国得到了广泛应用。我国以中科遥感院为代表的无人机在国内也取得了一系列重大的突破，研制出许多航空摄影测量设备（如武大生产的全数字摄影测量系统VirtuoZo等）。其航测主要以生产4D产品：...
C#实现空中三角测量解析与封装
在编程语言C#中编写解析空中三角测量程序，需要掌握以下几个方面的知识点： 1. **C# 基础知识**：了解C#语言的基本语法、数据类型、控制流语句、数组和集合、面向对象编程等核心概念。 2. **线性代数与数值方法**...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值