单像空间后方交会程序报告

单像空间后方交会程序报告指导老师：刘老师班级：测绘101姓名：尚锋学号： 19号1、应用程序的主入口部分的代码：using System;using System.Collections.Generic;using System.Linq;using System.Windows.Forms;namespace单像空间后方交会{static class Program{///<summary>///应用程序的主入口点。

///</summary>[STAThread]static void Main(){Application.EnableVisualStyles();Application.SetCompatibleTextRenderingDefault(false);Application.Run(new Form1());}}}2、方法解算类（通用）部分的代码：using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace单像空间后方交会{class Tongyong{struct image_point//一个像点结构，包含像点坐标和地面点坐标 {public double x;public double y;public double X;public double Y;public double Z;}private double f; //主距private double u; //u为外方位元素，下面5个相同private double w;private double k;private double Xs;private double Ys;private double Zs;private image_point[] p = new image_point[4]; //四个控制点private double[] R = new double[9]; //旋转矩阵private double[] a = new double[8]; //像点坐标近似值private double[,] A = new double[8, 6]; //误差方程式系数private double[] L = new double[8]; //误差方程式常数项private int count = 0; //统计代次数public Tongyong(double g, double[] q) //构造函数，初始化各变量,单位m{f = g;for (int i = 0; i < 4; i++){int j = i * 5;p[i].x = q[j];p[i].y = q[j + 1];p[i].X = q[j + 2];p[i].Y = q[j + 3];p[i].Z = q[j + 4];}double ave = 0, sum = 0; //求比例尺分母for (int i = 0; i < 3; i++){for (int j = i + 1; j < 4; j++){sum += Math.Sqrt(Math.Pow(p[i].X - p[j].X, 2) + Math.Pow(p[i].Y - p[j].Y, 2)) / Math.Sqrt(Math.Pow(p[i].x - p[j].x, 2) + Math.Pow(p[i].y - p[j].y, 2));}}ave = sum / 6;u = 0; //给定外方位元素的初始值,角度均设置为0w = 0;k = 0;Xs = (p[0].X + p[1].X + p[2].X + p[3].X) / 4; //Xs为四个控制点X的平均值，Ys类似Ys = (p[0].Y + p[1].Y + p[2].Y + p[3].Y) / 4;Zs = (p[0].Z + p[1].Z + p[2].Z + p[3].Z) / 4 + ave * f; }private double sin(double m) //正弦，为简化而写,下同{return Math.Sin(m);}private double cos(double m){return Math.Cos(m);}private void calcos() //计算旋转矩阵{R[0] = cos(u) * cos(k) - sin(u) * sin(w) * sin(k);R[1] = -cos(u) * sin(k) - sin(u) * sin(w) * cos(k);R[2] = -sin(u) * cos(w);R[3] = cos(w) * sin(k);R[4] = cos(w) * cos(k);R[5] = -sin(w);R[6] = sin(u) * cos(k) + cos(u) * sin(w) * sin(k);R[7] = cos(u) * sin(w) * cos(k) - sin(u) * sin(k);R[8] = cos(u) * cos(w);}private void calabout() //像点坐标的近似值{int i;for (i = 0; i < 4; i++){a[2 * i] = -f * (R[0] * (p[i].X - Xs) + R[3] * (p[i].Y - Ys) + R[6] * (p[i].Z - Zs)) / (R[2] * (p[i].X - Xs) + R[5] * (p[i].Y - Ys) + R[8] * (p[i].Z - Zs));a[2 * i + 1] = -f * (R[1] * (p[i].X - Xs) + R[4] * (p[i].Y - Ys) + R[7] * (p[i].Z - Zs)) / (R[2] * (p[i].X - Xs) + R[5] * (p[i].Y - Ys) + R[8] * (p[i].Z - Zs));}}private void calxx() //误差方程式的系数和常数项{int i;for (i = 0; i < 4; i++) //系数{double z = R[2] * (p[i].X - Xs) + R[5] * (p[i].Y - Ys) + R[8] * (p[i].Z - Zs);int n = i * 2;A[n, 0] = (R[0] * f + R[2] * p[i].x) / z;A[n, 1] = (R[3] * f + R[5] * p[i].x) / z;A[n, 2] = (R[6] * f + R[8] * p[i].x) / z;A[n, 3] = p[i].y * sin(w) - f * cos(w) * cos(k) - p[i].x / f * (p[i].x * cos(w) * cos(k) - p[i].y * cos(w) * sin(k));A[n, 4] = -f * sin(k) - p[i].x / f * (p[i].x * sin(k) + p[i].y * cos(k));A[n, 5] = p[i].y;A[n + 1, 0] = (R[1] * f + R[2] * p[i].y) / z;A[n + 1, 1] = (R[4] * f + R[5] * p[i].y) / z;A[n + 1, 2] = (R[7] * f + R[8] * p[i].y) / z;A[n + 1, 3] = -p[i].x * sin(w) + f * cos(w) * sin(k) - p[i].x / f * (p[i].x * cos(w) * cos(k) - p[i].y * sin(k) * cos(w));A[n + 1, 4] = -f * cos(k) - p[i].y / f * (p[i].x * sin(k) + p[i].y * cos(k));A[n + 1, 5] = -p[i].x;}for (i = 0; i < 4; i++) //常数项{L[2 * i] = p[i].x - a[2 * i];L[2 * i + 1] = p[i].y - a[2 * i + 1];}}private double calAdd(){double[,] temp = new double[6, 6]; //A的转置与A相乘的积double[,] ANew = new double[6, 8]; //A的转置double[] t = new double[6]; //A的转置与L相乘的积double[] X = new double[6]; //改正数int i, j, n;for (i = 0; i < 8; i++) //求A的转置ANew{for (j = 0; j < 6; j++){ANew[j, i] = A[i, j];}}for (i = 0; i < 6; i++) //求A的转置与A相乘的积temp{for (j = 0; j < 6; j++){temp[i, j] = 0;for (n = 0; n < 8; n++){temp[i, j] += ANew[i, n] * A[n, j];}}}MATINV(temp); //temp的逆，保存在本身矩阵中for (i = 0; i < 6; i++) //求A的转置与L的乘积t{t[i] = 0;for (j = 0; j < 8; j++){t[i] += ANew[i, j] * L[j];}}for (i = 0; i < 6; i++) //求改正数X{X[i] = 0;for (j = 0; j < 6; j++){X[i] += temp[i, j] * t[j];}}Xs += X[0]; //外方位元素初始值加上改正数Ys += X[1];Zs += X[2];u += X[3];w += X[4];k += X[5];return maxone(X); //返回判断条件,最大的改正数的值}public void makeSure() //计算流程控制函数{calcos();calabout();calxx();double VALUE = calAdd();count++;while (VALUE > 0.00001) //迭代至最大改正数<=0.00001为止 {calcos();calabout();calxx();VALUE = calAdd();count++;}}private void MATINV(double[,] c) //求6阶矩阵的逆{int i, j, h, m;const int n = 6;double l;double[,] q = new double[n, 12];for (i = 0; i < n; i++) //构造高斯矩阵{for (j = 0; j < n; j++){q[i, j] = c[i, j];}}for (i = 0; i < n; i++) //单位矩阵{for (j = n; j < 12; j++){if (i + 6 == j)q[i, j] = 1;elseq[i, j] = 0;}}for (h = 0, m = 0; m < n - 1; m++, h++) //消去对角线以下的数据{for (i = m + 1; i < n; i++){if (q[i, h] == 0d)continue;l = q[m, h] / q[i, h];for (j = 0; j < 12; j++){q[i, j] *= l;q[i, j] -= q[m, j];}}}for (h = n - 1, m = n - 1; m > 0; m--, h--) // 消去对角线以上的数据{for (i = m - 1; i >= 0; i--){if (q[i, h] == 0d)continue;l = q[m, h] / q[i, h];for (j = 0; j < 12; j++){q[i, j] *= l;q[i, j] -= q[m, j];}}}for (i = 0; i < n; i++)//将对角线上数据化为1{l = 1.0 / q[i, i];for (j = 0; j < 12; j++){q[i, j] *= l;}}for (i = 0; i < n; i++) //提取逆矩阵{for (j = 0; j < n; j++){c[i, j] = q[i, j + 6];}}}private double maxone(double[] Arr) //返回六个元素中的最大值 {double[] ARR = new double[6];for (int i = 0; i < 6; i++){ARR[i] = Arr[i];if (ARR[i] < 0d) //取正ARR[i] = -ARR[i];}Array.Sort(ARR, 0, 6);return ARR[5];}public int COUNT{get{return count;}}public double U{get{return u; }}public double W{get{return w; }}public double K{get{return k; }}public double XS {get{return Xs; }}public double YS {get{return Ys; }}public double ZS {get{return Zs; }}}3、窗体一部分的代码：using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Windows.Forms;using System.Diagnostics;using System.Data.OleDb;using System.IO;namespace单像空间后方交会{public partial class Form1 : Form{private double[] data = new double[20]; //保存表中的数据private double f = 0; //主距public Form1(){InitializeComponent();}private void Form1_Load(object sender, EventArgs e){// TODO: 这行代码将数据加载到表“database1DataSet1.data”中。