摄影测量后方交会算法C#实现

最新推荐文章于 2023-12-07 16:56:22 发布
最新推荐文章于 2023-12-07 16:56:22 发布
阅读量3.4k 收藏 75
点赞数 11
文章标签： C# 摄影测量后方交会
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/wudiyuehao/article/details/100598940
版权
  1. 界面设计
    在这里插入图片描述
  2. 控制点
PID,x(mm),y(mm),X(m),Y(m),Z(m)
1,-86.15,-68.99,36589.41,25273.32,2195.17
2,-53.4,82.21,37631.08,31324.51,728.69
3,-14.78,-76.63,39100.97,24934.98,2386.5
4,10.46,64.43,40426.54,30319.81,757.31
  1. 导入数据
OpenFileDialog dialog = new OpenFileDialog();
            dialog.Filter = "文本文件.txt|*.txt";
            if (dialog.ShowDialog() == DialogResult.OK)
            {
                path = Path.GetFullPath(dialog.FileName);
                //从路径中获取文件名
            }
          
            string data;
           
            DataTable dt = new DataTable();
            StreamReader reader = new StreamReader(path);
            while (true)
            {
                data = reader.ReadLine();
                if (data == null)
                    break;
                string[] v = data.Split(new string[] { "," }, StringSplitOptions.RemoveEmptyEntries);
                if (data.Contains("PID"))
                {
                    dt.Columns.Add(v[0]);
                    dt.Columns.Add(v[1]);
                    dt.Columns.Add(v[2]);
                    dt.Columns.Add(v[3]);
                    dt.Columns.Add(v[4]);
                    dt.Columns.Add(v[5]);
                    continue;
                }
                DataRow dr = dt.NewRow();
                MyPoint p = new MyPoint();
                p.id = Convert.ToInt32(v[0]);
                p.x = (Convert.ToDouble(v[1])-x0) / 1000;//减去内方位元素
                p.y = (Convert.ToDouble(v[2])-y0) / 1000;
                p.X = Convert.ToDouble(v[3]);
                p.Y = Convert.ToDouble(v[4]);
                p.Z = Convert.ToDouble(v[5]);
                dr[0] = p.id;
                dr[1] = p.x;
                dr[2] = p.y;
                dr[3] = p.X;
                dr[4] = p.Y;
                dr[5] = p.Z;
                dt.Rows.Add(dr);
                ConPoint.Add(p);
            }
            n = ConPoint.Count;
            
            dataGridView1.DataSource = dt;

4.迭代计算

double error = Convert.ToDouble(txterror.Text);
            //外方位元素赋初值
            double sumx = 0, sumy = 0;
            for (int i = 0; i < ConPoint.Count; i++)
            {
                sumx += ConPoint[i].X;
                sumy += ConPoint[i].Y;
            }
            Xs = sumx / n;
            Ys = sumy / n;
            Zs = m * f;
            //旋转矩阵
            double[,] R = new double[3, 3];
            //近似（x）（y）
            double[] xAppro = new double[n];
            double[] yAppro = new double[n];
            double[] XX = new double[n];
            double[] YY = new double[n];
            double[] ZZ = new double[n];
            double[,] V = new double[2 * n, 1];
            double[,] A = new double[2 * n, 6];
            double[,] X = new double[6, 1];
            double[,] L = new double[2 * n, 1];
            int num = 0;
            do
            {
                num++;
                R[0, 0] = Math.Cos(phi) * Math.Cos(kappa) - Math.Sin(phi) * Math.Sin(omega) * Math.Sin(kappa);//a1
                R[0, 1] = -Math.Cos(phi) * Math.Sin(kappa) - Math.Sin(phi) * Math.Sin(omega) * Math.Cos(kappa);//a2
                R[0, 2] = -Math.Sin(phi) * Math.Cos(omega);//a3
                R[1, 0] = Math.Cos(omega) * Math.Sin(kappa);//b1
                R[1, 1] = Math.Cos(omega) * Math.Cos(kappa);//b2
                R[1, 2] = -Math.Sin(omega);//b3
                R[2, 0] = Math.Sin(phi) * Math.Cos(kappa) + Math.Cos(phi) * Math.Sin(omega) * Math.Sin(kappa);//c1
                R[2, 1] = -Math.Sin(phi) * Math.Sin(kappa) + Math.Cos(phi) * Math.Sin(omega) * Math.Cos(kappa);//c2
                R[2, 2] = Math.Cos(phi) * Math.Cos(omega);//c3
                //计算(x)(y),n个控制点的
                for (int i = 0; i < n; i++)
                {
                    XX[i] = R[0, 0] * (ConPoint[i].X - Xs) + R[1, 0] * (ConPoint[i].Y - Ys) + R[2, 0] * (ConPoint[i].Z - Zs);
                    YY[i] = R[0, 1] * (ConPoint[i].X - Xs) + R[1, 1] * (ConPoint[i].Y - Ys) + R[2, 1] * (ConPoint[i].Z - Zs);
                    ZZ[i] = R[0, 2] * (ConPoint[i].X - Xs) + R[1, 2] * (ConPoint[i].Y - Ys) + R[2, 2] * (ConPoint[i].Z - Zs);
                    xAppro[i] = -f * XX[i] / ZZ[i];
                    yAppro[i] = -f * YY[i] / ZZ[i];
                }
                //计算A L矩阵
                for (int i = 0; i < n; i++)
                {
                    A[2 * i, 0] = (R[0, 0] * f + R[0, 2] * ConPoint[i].x) / ZZ[i];
                    A[2 * i, 1] = (R[1, 0] * f + R[1, 2] * ConPoint[i].x) / ZZ[i];
                    A[2 * i, 2] = (R[2, 0] * f + R[2, 2] * ConPoint[i].x) / ZZ[i];
                    A[2 * i, 3] = ConPoint[i].y * Math.Sin(omega) - (ConPoint[i].x * (ConPoint[i].x * Math.Cos(kappa) - ConPoint[i].y * Math.Sin(kappa)) / f + f * Math.Cos(kappa)) * Math.Cos(omega);
                    A[2 * i, 4] = -f * Math.Sin(kappa) - ConPoint[i].x * (ConPoint[i].x * Math.Sin(kappa) + ConPoint[i].y * Math.Cos(kappa)) / f;
                    A[2 * i, 5] = ConPoint[i].y;
                    A[2 * i + 1, 0] = (R[0, 1] * f + R[0, 2] * ConPoint[i].y) / ZZ[i];
                    A[2 * i + 1, 1] = (R[1, 1] * f + R[1, 2] * ConPoint[i].y) / ZZ[i];
                    A[2 * i + 1, 2] = (R[2, 1] * f + R[2, 2] * ConPoint[i].y) / ZZ[i];
                    A[2 * i + 1, 3] = -ConPoint[i].x * Math.Sin(omega) - (ConPoint[i].y * (ConPoint[i].x * Math.Cos(kappa) - ConPoint[i].y * Math.Sin(kappa)) / f - f * Math.Sin(kappa)) * Math.Cos(omega);
                    A[2 * i + 1, 4] = -f * Math.Cos(kappa) - ConPoint[i].y * (ConPoint[i].x * Math.Sin(kappa) + ConPoint[i].y * Math.Cos(kappa)) / f;
                    A[2 * i + 1, 5] = -ConPoint[i].x;
                    L[2 * i, 0] = ConPoint[i].x - xAppro[i];
                    L[2 * i + 1, 0] = ConPoint[i].y - yAppro[i];
                }
//矩阵运算求解
MatrixO matrixw = new MatrixO();
                double[,] AT = matrixw.MatrixTrans(A);
                double[,] ATA = matrixw.Matrixmulti(AT, A);
                double[,] Inv = matrixw.MatrixInvByCom(ATA);
                double[,] InvAT = matrixw.Matrixmulti(Inv, AT);
                X = matrixw.Matrixmulti(InvAT, L);
                
                Xs += X[0, 0];
                Ys += X[1, 0];
                Zs += X[2, 0];
                phi += X[3, 0];
                omega += X[4, 0];
                kappa += X[5, 0];
            }
            while (Math.Abs(X[0, 0]) >= error || Math.Abs(X[1, 0]) >= error || Math.Abs(X[2, 0]) >= error || Math.Abs(X[3, 0]) >= error || Math.Abs(X[4, 0]) >= error || Math.Abs(X[5, 0]) >= error);

//整理输出内容
MatrixO matrix1 = new MatrixO();
            V = matrix1.MatrixSub(matrix1.Matrixmulti(A, X), L);
            double[,] at = matrix1.MatrixTrans(A);
            double[,] ata = matrix1.Matrixmulti(at, A);
            double[,] inv = matrix1.MatrixInvByCom(ata);

double m0 = Math.Sqrt((matrix1.Matrixmulti(matrix1.MatrixTrans(V), V)[0,0]) / (2 * Convert.ToDouble(n) - 6));//中误差
            double[] mi = new double[6];
            for (int i = 0; i < 6; i++)
            {
                mi[i] = m0 * Math.Sqrt(inv[i, i]);
            }
            txtXs.Text = Math.Round(Xs, 4).ToString();
            txtYs.Text = Math.Round(Ys, 4).ToString();
            txtZs.Text = Math.Round(Zs, 4).ToString();
            RadianAngle angle = new RadianAngle();
            txtPhi.Text = angle.TransArcToDMS(phi);
            txtOmega.Text = angle.TransArcToDMS(omega);
            txtKappa.Text = angle.TransArcToDMS(kappa);
            labnum.Text = "迭代次数"+num.ToString();
//写出到文件
StreamWriter writer = new StreamWriter(@"C: \Userdata\结果.txt");
            writer.WriteLine("像片外方位元素为：");
            writer.WriteLine("Xs={0}  m={1}", Xs, mi[0]);
            writer.WriteLine("Ys={0}  m={1}", Ys, mi[1]);
            writer.WriteLine("Zs={0}  m={1}", Zs, mi[2]);
            writer.WriteLine("φ={0}  m={1}", phi, mi[3]);
            writer.WriteLine("ω={0}  m={1}", omega, mi[4]);
            writer.WriteLine("κ={0}  m={1}", kappa, mi[5]);
            writer.WriteLine();
            writer.WriteLine("旋转矩阵为：");
            writer.WriteLine("{0}    {1}    {2}", R[0, 0], R[0, 1], R[0, 2]);
            writer.WriteLine("{0}    {1}    {2}", R[1, 0], R[1, 1], R[1, 2]);
            writer.WriteLine("{0}    {1}    {2}", R[2, 0], R[2, 1], R[2, 2]);
            writer.Flush();
            writer.Close();

5、矩阵运算类

//矩阵乘法
        public double [,] Matrixmulti(double[,] A, double[,] B)
        {
            int m1 = A.GetLength(0);
            int n1 = A.GetLength(1);
            int m2 = B.GetLength(0);
            int n2 = B.GetLength(1);
            if (n1 != m2)
            {
                Exception myException = new Exception("数组维数不匹配");
                throw myException;
            }
            double[,] c = new double[m1, n2];
            double[,] a = A;
            double[,] b = B;
            for (int i = 0; i < m1; i++)
                for (int j = 0; j < n2; j++)
                {
                    c[i, j] = 0;
                    for (int k = 0; k < n1; k++)
                        c[i, j] += a[i, k] * b[k, j];
                }
            return c;
        }
        //矩阵乘法
        public static double[,] MatrixMulti(double[,] A, double[,] B)
        {
            int m1 = A.GetLength(0);
            int n1 = A.GetLength(1);
            int m2 = B.GetLength(0);
            int n2 = B.GetLength(1);
            if (n1 != m2)
            {
                Exception myException = new Exception("数组维数不匹配");
                throw myException;
            }
            double[,] c = new double[m1, n2];
            double[,] a = A;
            double[,] b = B;
            for (int i = 0; i < m1; i++)
                for (int j = 0; j < n2; j++)
                {
                    c[i, j] = 0;
                    for (int k = 0; k < n1; k++)
                        c[i, j] += a[i, k] * b[k, j];
                }
            return c;
        }
        //矩阵转置
        public double[,] MatrixTrans(double[,]A)
        {
            int m = A.GetLength(0);
            int n = A.GetLength(1);
            double[,] c = new double[n, m];
            double[,] a = A;
            for (int i = 0; i < n; i++)
                for (int j = 0; j < m; j++)
                    c[i, j] = a[j, i];
            return c;
        }
        //矩阵求逆
        public double[,] MatrixInvByCom( double[,]A)
        {
            double d = MatrixO.MatrixDet(A);
            if (d == 0)
            {
                Exception myException = new Exception("没有逆矩阵");
                throw myException;
            }
            double[,] Ax = MatrixO.MatrixCom(A);
            double[,] An = MatrixO.MatrixSimpleMulti((1.0 / d), Ax);
            return An;
        }
        //对应行列式的代数余子式矩阵
        public static double[,] MatrixSpa(double[,] A, int ai, int aj)
        {
            int m = A.GetLength(0);
            int n = A.GetLength(1);
            if (m != n)
            {
                Exception myException = new Exception("矩阵不是方阵");
                throw myException;
            }
            int n2 = n - 1;
            double[,] a = A;
            double[,] b = new double[n2,n2];
            //左上
            for (int i = 0; i < ai; i++)
                for (int j = 0; j < aj; j++)
                {
                    b[i, j] = a[i, j];
                }
            //右下
            for (int i = ai; i < n2; i++)
                for (int j = aj; j < n2; j++)
                {
                    b[i, j] = a[i + 1, j + 1];
                }
            //右上
            for (int i = 0; i < ai; i++)
                for (int j = aj; j < n2; j++)
                {
                    b[i, j] = a[i, j + 1];
                }
            //左下
            for (int i = ai; i < n2; i++)
                for (int j = 0; j < aj; j++)
                {
                    b[i, j] = a[i + 1, j];
                }
            //符号位
            if ((ai + aj) % 2 != 0)
            {
                for (int i = 0; i < n2; i++)
                    b[i, 0] = -b[i, 0];
            }
            return b;
        }

        //矩阵的行列式,矩阵必须是方阵
        public static double MatrixDet(double[,]A)
        {
            int m = A.GetLength(0);
            int n = A.GetLength(1);
            if (m != n)
            {
                Exception myException = new Exception("数组维数不匹配");
                throw myException;
            }
            double[,] a = A;
            if (n == 1) return a[0, 0];
            double D = 0;
            for (int i = 0; i < n; i++)
            {
                D += a[1, i] * MatrixDet(MatrixSpa(A, 1, i));
            }
            return D;
        }
        //矩阵的伴随矩阵
        public static double[,] MatrixCom(double[,]A)
        {
            int m = A.GetLength(0);
            int n = A.GetLength(1);
            double[,] c = new double[m, n];
            double[,] a = A;
            for (int i = 0; i < m; i++)
                for (int j = 0; j < n; j++)
                    c[i, j] = MatrixDet(MatrixSpa(A, j, i));
            return c;
        }

  1. 结果
    在这里插入图片描述
    在这里插入图片描述

安平桥
关注
C#版后方交会程序
12-08
自己当时写的一个空间后方交会程序，已调试，能通过，当年过关的程序
C#简单实现摄影测量后方交会
Lynn_whu的博客
09-16 1万+
C#简单实现摄影测量后方交会                                                                                              ——题目来自武大版《摄影测量学》习题2.27                   一、    实验目的和要求 通过对题目已知的4对影像坐标和地面坐标，通过编程实现计算近似垂直摄...
摄影测量后方交会-前方交会(C#)
热门推荐
Geo-xudong的博客
09-05 1万+
双像解析摄影测量中的后交-前交解法，常用于已知像片的外方位元素、需确定少量待定点坐标的情况。解算过程分为两个阶段：后方交会、前方交会。 思路为利用后方交会计算外方位元素，利用前方交会解算待定地面点坐标。 后方交会的解算过程为： （1）获取已知数据：比例尺分母m，内方位元素x0,y0,f，以及控制点的地面摄影测量坐标X、Y、Z。 （2）量测控制点的像点坐标：得到像点坐标x、y。 （3）确定未
后方交会
Hanford的专栏
12-15 7536
目 录 第1章 计算公式    1 1.1 计算公式一    1 1.1.1 角度取法    1 1.1.2 代码实现    3 1.2 计算公式二    5 1.2.1 公式推导    5 1.2.2 代码实现    6 1.3 危险圆    7 第2章 精度评定    8   第1章
测量-后方交会C#代码.zip_C#_C#后方_后方交会_测绘_测量程序
07-13
从测量工程原理出发，利用C#编程工具，编写后方交会测量程序，实现简单类型的后方交会测量计算。
单片空间后方交会[C#]
Giser.Lu
10-15 5344
  using System;using System.Collections.Generic;using System.Text;namespace _33...{    class Program    ...{        public class Data        ...{            //内方位元素            double f = 0.1
摄影测量_单张像片的空间后方交会(C# windows form)
11-29
总的来说，"摄影测量_单张像片的空间后方交会(C# windows form)"项目结合了摄影测量学的理论知识、C#编程技巧和图像处理算法，为图像的地理定位提供了有效工具。开发者需要对摄影测量的基本原理有深入理解，并具备...
摄影测量空间后方交会.zip
09-14
空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。 算法：由于每...
摄影测量后方交会_摄影测量后方交会_源码
10-04
"后方交会"是摄影测量中的一种关键算法，主要用于从已知控制点恢复未知像点的地面坐标。在此，我们详细探讨一下后方交会的基本原理、实现方法以及可能遇到的问题。 后方交会的理论基础在于几何成像模型，即摄影机...
摄影测量空间后方交会
03-28
根据提供的信息，我们可以深入探讨“摄影测量空间后方交会”这一主题。这涉及到摄影测量学中的一个关键领域：外方位元素解算。 ### 摄影测量空间后方交会简介 摄影测量空间后方交会是一种在摄影测量学中确定摄影机...
单片空间后方交会C#源代码
11-26
（1）已知航摄仪内方位元素f＝153.24mm，Xo＝Yo＝0。限差0.1秒 （2）已知4对点的影像坐标和地面坐标：
C#前后方交会定点计算程序
10-28
利用C#编写的交会定点计算程序，可以实现计算、结果输出、图形绘制与输出等功能。
C#代码实现空间后方交会
11-06
C#代码实现空间后方交会，可以直接使用，无需配置
距离交会程序源代码C#
07-11
测绘工程专业，距离交会计算程序，C#+VS2017编写。测绘工程专业，距离交会计算程序，C#+VS2017编写。测绘工程专业，距离交会计算程序，C#+VS2017编写。
相机内定向参数计算程序设计（C#）
10-31
该文档介绍了数字摄影测量中相机内定向参数平差计算的基本方法和程序设计，附程序代码和运行截图。
【c#】vs2017摄影测量影像后方交会
05-06
用vs2017调试，c#编写的后方交会程序，读取txt数据来计算（内附测试数据），读取左片或者右片坐标数据一键计算外方位元素（我解释下怎么用这个程序，打开一个左片文件只能计算左片外方位，要计算右片那就还要再打开...
C#单像空间后方交会计算
最新发布
m0_68160088的博客
12-07 769
摄影测量学
用C#编写摄影测量后方交会求外方位元素
u011713916的博客
09-14 2681
实验原理 **详细题目： ** 编程思想： 具体代码 我使用了MathNet.Numerics.dll矩阵库，功能强大，还可以求解病态矩阵 MathNet.Numerics.dll 完整源码： System; using System.Net.Http.Headers; using System.Xml.Linq; using MathNet.Numerics.LinearAlgebra...
单像空间后方交会
qq_64775230的博客
04-22 897
定义：像片覆盖范围内的地面三个以上的控制点坐标 + 对应点的像点坐标 -----------------》像片的6个外方位元素。7.如果角元素的改正值均小于0.1’，不用再次计算，如果不满足条件，就要将新值重新带回第一步重新计算，直到改正值满足条件。但是，如果想继续研究，那就要转成线性的（按泰勒级数展开的对共线方程线性化）将这个真实值带入到上述已经线性化的共线方程中，那么就能得到关于改正数的方程。如果我们将（x，y）设计为观测值，vx，vy设计为改正数。6.利用逆矩阵求解改正数的值。2.根据比例尺求航高。
评论 9
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值