精华内容
下载资源
问答
  • 后方交会前方交会

    2018-02-19 16:25:26
    双像解析摄影测量中的后交-前交解法,常用于已知像片的外方位元素、需...解算过程分为两个阶段:后方交会前方交会。 思路为利用后方交会计算外方位元素,利用前方交会解算待定地面点坐标。 后方交会的解算过程为:
  • 后方交会_前方交会

    2015-09-05 20:51:45
    用C#实现的摄影测量后方交会前方交会的整体代码,利用网络的样例数据已验证代码没问题,请大家免费下载
  • 交会定点(前方交会、测边交会、后方交会)程序C++代码

    一、设计思路

    1.新建一个Function.h文件,在其中声明二维坐标结构体、角度转换函数、方位角计算函数、两点间距离计算函数和三个定点交会函数,并在新建一个Function.cpp文件对函数进行定义(Function.cpp文件中要加入#include "Function.h");

    2.在CExercise_4Dlg.cpp文件头处添加 #include "Function.h"

    3.计算按钮消息响应函数对交会定点计算函数的调用:

    iMethod==1,调用前方交会函数;

    iMethod==2,调用测边交会函数;

    iMethod==3,调用后方交会函数。

    二、主要代码

    // function.h : 头文件
    #pragma once
    struct Point2D    //定义二维坐标结构体
    {
            double X;
            double Y;
    };
    double Dms_TO_Rad(double dDms);                         //度分秒转化为弧度制
    double Distance(struct Point2D &A,struct Point2D &B);   //计算两点间的距离
    double Azimuth(struct Point2D &A,struct Point2D &B);    //计算方位角
    struct Point2D ForwardIntersection(struct Point2D &A,struct Point2D &B,double Alpha,double Betta);
    //前方交会定点计算
    struct Point2D LinearIntersection(struct Point2D &A,struct Point2D &B,double Dap,double Dbp);
    //测边交会定点计算
    struct Point2D ResctionIntersection(struct Point2D &A,struct Point2D &B,struct Point2D &C,double Alpha,double Betta,double Gamma);
    //后方交会定点计算
    // function.cpp : 实现文件
    #include "StdAfx.h"
    #include "Function.h"
    #include "math.h"
     
    const double PI=4.0*atan(1.0);
     
     
    double Dms_TO_Rad(double dDms)
    {
            double dRad;
            int iDegree,iMin;
            double dSec;
            iDegree=int(dDms);
            iMin=int((dDms-iDegree)*100);
            dSec=((dDms-iDegree)*100-iMin)*100;
            dDms=iDegree+double(iMin)/60+dSec/3600;
            dRad=dDms/180*PI;
            return dRad;
    }
    double Distance(struct Point2D &A,struct Point2D &B)
    {
            double dDistance=0;
            dDistance=sqrt((A.X-B.X)*(A.X-B.X)+(A.Y-B.Y)*(A.Y-B.Y));
            return dDistance;
    }
    double Azimuth(struct Point2D &A,struct Point2D &B)
    {
            // TODO: 在此添加控件通知处理程序代码
     
            double dAzimuth=0;                     //方位角
            double dx=B.X-A.X;                     //X坐标增量
            double dy=B.Y-A.Y;                     //Y坐标增量
     
            if (dx>0)
            {
                    if (dy>0)
                    {
                            dAzimuth=atan(dy/dx);
                    }
                    else if(dy<0)
                    {
                            dAzimuth=atan(dy/dx)+PI*2;
                    }
                    else
                            dAzimuth=0;
            }
            else if (dx<0)
            {
                    if(dy>0)
                    {
                            dAzimuth=atan(dy/dx)+PI;
                    }
                    else if(dy<0)
                    {
                            dAzimuth=atan(dy/dx)+PI;
                    }
                    else
                            dAzimuth=PI;
            }
            else
            {
                    if(dy>0)
                            dAzimuth=PI/2;
                    else if(dy<0)
                            dAzimuth=3*PI/2;
            }
            return dAzimuth;
    }
    struct Point2D ForwardIntersection(struct Point2D &A,struct Point2D &B,double Alpha,double Betta)
    {
            struct Point2D P;
            //将度分秒转化为弧度制
            double dA=Dms_TO_Rad(Alpha);
            double dB=Dms_TO_Rad(Betta);
            //计算A,B角余切值
            double dCotA=1/tan(dA);
            double dCotB=1/tan(dB);
            //计算前方交会定位值
            P.X=((A.X*dCotB+B.X*dCotA)+(B.Y-A.Y))/(dCotA+dCotB);
            P.Y=((A.Y*dCotB+B.Y*dCotA)+(A.X-B.X))/(dCotA+dCotB);
            return P;
    }
    struct Point2D LinearIntersection(struct Point2D &A,struct Point2D &B,double Dap,double Dbp)
    {
            struct Point2D P;
            double Dab;                            //三角形AB边长
            Dab=Distance(A,B);
            double alpha_ab;                       //三角形AB方位角
            alpha_ab=Azimuth(A,B);
            double Angle_A;                        //三角形内角A
            Angle_A=acos((Dab*Dab+Dap*Dap-Dbp*Dbp)/(2*Dab*Dap));
            double alpha_ap=alpha_ab-Angle_A;
     
            P.X=A.X+Dap*cos(alpha_ap);
            P.Y=A.Y+Dap*sin(alpha_ap);
            return P;
    }
    struct Point2D ResctionIntersection(struct Point2D &A,struct Point2D &B,struct Point2D &C,double Alpha,double Betta,double Gamma)
    {
            struct Point2D P;
            double dA;                             //三角形内角A
            double dB;                             //三角形内角B
            double dC;                             //三角形内角C
            double a;                              //三角形边长a
            double b;                              //三角形边长b
            double c;                              //三角形边长c
     
            a=Distance(C,B);
            b=Distance(C,A);
            c=Distance(B,A);
            dA=acos((b*b+c*c-a*a)/(2*b*c));
            dB=acos((a*a+c*c-b*b)/(2*a*c));
            dC=acos((b*b+a*a-c*c)/(2*b*a));
     
            double C_TEMP;                        //定义判断是否在危险圆附近的判断变量
            double B_TEMP;
            double A_TEMP;
     
            C_TEMP=Alpha+Betta+dC;
            B_TEMP=Alpha+Gamma+dB;
            A_TEMP=Betta+Gamma+dA;
     
            if ((C_TEMP>170)&&C_TEMP<190||
                    (B_TEMP>170)&&B_TEMP<190||
                    (B_TEMP>170)&&B_TEMP<190)
            {
                    P.X=-1;
                    P.Y=-1;
                    return P;                         //坐标是(-1,-1)则在危险圆附近
            }
     
            double Alpha_Rad=Dms_TO_Rad(Alpha);   //度分秒转化为弧度
            double Betta_Rad=Dms_TO_Rad(Betta);
            double Gamma_Rad=Dms_TO_Rad(Gamma);
     
            double pa=(tan(Alpha_Rad)*tan(dA))/(tan(Alpha_Rad)-tan(dA));
            double pb=(tan(Betta_Rad)*tan(dB))/(tan(Betta_Rad)-tan(dB));
            double pc=(tan(Gamma_Rad)*tan(dC))/(tan(Gamma_Rad)-tan(dC));
     
            P.X=(pa*A.X+pb*B.X+pc*C.X)/(pa+pb+pc);
            P.Y=(pa*A.Y+pb*B.Y+pc*C.Y)/(pa+pb+pc);
            return P;
    }
    // Exercise_4Dlg.h : 头文件
    // CExercise_4Dlg 对话框
    class CExercise_4Dlg : public CDialog
    {
    public:
            //前方交会控件变量
            double dXA_F;
            double dYA_F;
            double dXB_F;
            double dYB_F;
            double dAlpha_F;
            double dBetta_F;
            double dXP_F;
            double dYP_F;
            //测边交会控件变量
            double dXA_L;
            double dYA_L;
            double dXB_L;
            double dYB_L;
            double dDap_L;
            double dDbp_L;
            double dXP_L;
            double dYP_L;
            //后方交会空间变量
            double dXA_R;
            double dYA_R;
            double dXB_R;
            double dYB_R;
            double dXC_R;
            double dYC_R;
            double dXP_R;
            double dYP_R;
            double dAlpha_R;
            double dBetta_R;
            double dGamma_R;
            int iMethod;                          //交会定点类型代号
            afx_msg void OnBnClickedBtnCal();     //计算按钮
            afx_msg void OnBnClickedBtnClear();   //重置编辑框数据
            afx_msg void OnBnClickedCancel();     //退出对话框
            void ForwardIntersection(double XA,double YA,double XB,double YB,double Alpha,double Betta,double &XP,double &YP);
            //前方交会定点计算
            void LinearIntersection(double XA,double YA,double XB,double YB,double Dap,double Dbp,double &XP,double &YP);
            //测边交会定点计算
            void ResctionIntersection(double XA,double YA,double XB,double YB,double XC,double YC,double Alpha,double Betta,double Gamma,double &XP,double &YP);
            //后方交会定点计算
    };
    // Exercise_4Dlg.cpp : 实现文件
    #include "stdafx.h"
    #include "Exercise_4.h"
    #include "Exercise_4Dlg.h"
    #include "Function.h"
    CExercise_4Dlg::CExercise_4Dlg(CWnd* pParent /*=NULL*/)
            : CDialog(CExercise_4Dlg::IDD, pParent)
            , dXA_F(116.942)     //构造函数中给定变量初始值
            , dYA_F(683.295)
            , dXB_F(522.909)
            , dYB_F(794.647)
            , dAlpha_F(59.1042)
            , dBetta_F(56.3254)
            , dXP_F(398.151)
            , dYP_F(413.249)
            , dXA_L(1807.041)
            , dYA_L(719.853)
            , dXB_L(1646.382)
            , dYB_L(830.660)
            , dDap_L(105.983)
            , dDbp_L(159.648)
            , dXP_L(1805.957)
            , dYP_L(825.830)
            , dXA_R(840.134)
            , dYA_R(844.422)
            , dXB_R(1001.542)
            , dYB_R(1620.616)
            , dXC_R(659.191)
            , dYC_R(1282.629)
            , dXP_R(503.702)
            , dYP_R(1500.075)
            , dAlpha_R(-68.0237)
            , dBetta_R(-8.2414)
            , dGamma_R(76.2651)
            , iMethod(1)
    {
            m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
    }
    void CExercise_4Dlg::OnBnClickedBtnCal()
    {
            // TODO: 在此添加控件通知处理程序代码
            UpdateData(TRUE);
     
            struct Point2D A,B,C,P;
            //定义个二维结构体变量存放A,B,C,P坐标
     
            if (iMethod==1)
            {
                    A.X=dXA_F;
                    A.Y=dYA_F;
                    B.X=dXB_F;
                    B.Y=dYB_F;
                    P=ForwardIntersection(A,B,dAlpha_F,dBetta_F);
                    //调用前方交会定点计算函数
                    dXP_F=P.X;
                    dYP_F=P.Y;
            }
            else if (iMethod==2)
            {
                    A.X=dXA_L;
                    A.Y=dYA_L;
                    B.X=dXB_L;
                    B.Y=dYB_L;
                    P=LinearIntersection(A,B,dDap_L,dDbp_L);
                    //调用测边交会定点计算函数
                    dXP_L=P.X;
                    dYP_L=P.Y;
            }
            else if (iMethod==3)
            {
                    A.X=dXA_R;
                    A.Y=dYA_R;
                    B.X=dXB_R;
                    B.Y=dYB_R;
                    C.X=dXC_R;
                    C.Y=dYC_R;
                    P=ResctionIntersection(A,B,C,dAlpha_R,dBetta_R,dGamma_R);
                    //调用后方交会定点计算函数
     
                    dXP_R=P.X;
                    dYP_R=P.Y;
                    if ((dXP_R==-1)&&(dXP_R=-1))
                    {
                            MessageBox(_T("在危险圆附近"));
                    }
            }
            else
                    MessageBox(_T("代号输入有误"));
     
            UpdateData(FALSE);
    }

    三、运行结果

    展开全文
  • 基于VC6.0开发的摄影测量空间前方交会空间后方交会程序。
  • C#+VS2017编写,完整项目,打开即用。包括前方交会和后方交会,可视化界面,表格存放数据。
  • 利用C语言实现摄影测量中的前方交会与后方交会,而且非常好用,
  • 前方交会、测边交会、后方交会计算程序,供参考,如有不足之处还请指出。
  • 本程序给出了摄影测量学中常见后方交会前方交会的MFC实现方法,除此之外还给出了特征点提取和同名点匹配过程。希望对你们学习有所启发和帮助!
  • 这是一个用于交会测量计算的小软件,可以尽心前方交会、后方交会、测边交会的计算。
  • 近景摄影测量空间后方交会空间前方交会解法程序,matlab
  • 可以实现摄影测量中的前方交会,后方交会等功能
  • MATLAB语言编写空间后方交会-空间前方交会程序
  • 摄影测量最基础的程序代码,包括后方交会前方交会等等。
  • 本程序为近景摄影测量控制场定标的一部分,提供加密...函数重用率比较高,能用于航摄的后方交会前方交会! 另外带有批量处理算法! 相关文档: http://hi.baidu.com/yiyiyis/blog/item/69858ec4628919da38db4934.html
  • 空间后方交会前方交会 MFC实现 CSU 摄影测量学

    千次阅读 多人点赞 2019-08-02 16:29:50
    空间前方-后方交会 C++实现 CSU空间前方-后方交会 C++实现 CSU一、 实验目的二、实验内容与要求三、设计与实现:3.1设计思路3.1.1 基本框架:3.1.2 成员关系:3.2 界面设计及属性:3.3 主要代码:3.3.1 文件:空间...

    空间后方交会前方交会 MFC 实现 CSU 摄影测量学

    摄影测量学基础

    (工具:VS2010)

    一、 实验目的

    • 通过对数字影像空间后交前交的程序设计实验,要求我们进一步理解和掌握影像外方位元素的有关理论、原理和方法。掌握VC++.net中创建类
    • 利用计算机程序设计语言编写摄影测量空间交会软件进行快速确定影像的外方位元素及其精度,然后通过求得的外方位元素求解未知点的地面摄影测量坐标,达到通过摄影测量量测地面地理数据的目的。

    二、实验内容与要求

    要求:用C、VB或者Matlab编写空间后方交会一前方交会计算机程序。
    ➢提交实验报告:程序框图,程序源代码、计算结果及体会。
    ➢计算结果:地面点坐标、外方位元素及精度。
    实验数据:
    在这里插入图片描述

    三、设计与实现:

    3.1设计思路

    3.1.1 基本框架:

    在这里插入图片描述

    3.1.2 成员关系:

    在这里插入图片描述

    3.2 界面设计及属性:

    在这里插入图片描述

    3.3 主要代码:

    3.3.1 文件:Matrix.h

    /***************************************************************************
    *  文件名:<Matrix.h>                                                      *
    *                                                                          *
    *  描述:矩阵类                                                             *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月19日        引用              ***                       *
    *                                                                          *
    *  外部过程:                                                              *
    *                                                                          *
    /**************************************************************************/
    
    
    #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);//调整数组的大小,原有数据不变(未测试)
    
    	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();//生成单位矩阵
    };
    
    

    3.3.2文件:Math.cpp

    #include "StdAfx.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;
    }
    //重载预算符-
    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;
    		 }
    	 }
    }
    
    

    3.3.3文件:Support.h

    /***************************************************************************
    *  文件名:<Support.h>                                                      *
    *                                                                          *
    *  描述:封装所有函数的类、调用矩阵CMatrix类                                *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月19日        创建              ***                       *
    *                                                                          *
    *  外部过程:                                                              *
    *                                                                          *
    /**************************************************************************/
    
    #pragma once
    #include"Matrix.h"
    
    class CSupport
    {
    	double dLx;
    	double dLy;
    	double dRx;
    	double dRy;
    	double X;
    	double Y;
    	double Z;
    public:
    	int SplitStringArray(CString str, char split, CStringArray& aStr);
    	double length(double x1,double y1,double x2,double y2);
    	CSupport(void);
    	~CSupport(void);
    	void R1(CMatrix &R,CMatrix X );
    	void read(CString &strExample,CSupport *GCP);
    	int SpaceBack(CSupport *GCP,int xx,int yy,CMatrix &Left,bool _L,double &Left_a0,CMatrix &LeftPhotoes_Pre);
    	void SpaceFront(CMatrix &Left,CMatrix &Right,CMatrix &B,CMatrix &R_Left,CMatrix &R_Right,CSupport *GCP);
    	void out(double Left_a0,double Right_a0,CMatrix LeftPhotoes,CMatrix LeftPhotoes_Pre,CMatrix RightPhotoes,CMatrix RightPhotoes_Pre,CSupport *GCP ,CString &strResult,int t1,int t2);
    	void main(CString &strExample,CString &strResult);
    };
    
    

    3.3.3文件:Support.cpp

    #include "StdAfx.h"
    #include "math.h"
    #include "Matrix.h"
    #include "Support.h"
    
    
    /***************************************************************************
    *  名字:double CSupport::SplitStringArray(CString str, char split, CStringArray& aStr)    *
    *                                                                          *
    *  描述:字符串分割函数                                                    *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年3月20日       引用                ***                      *
    *  参数: 1.CString str                                                    *
    *         2.char split                                                     *
    *         3.CStringArray& aStr                                             *
    *  返回值:int类型数据 返回n                                               *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    int CSupport::SplitStringArray(CString str, char split, CStringArray& aStr)
    {
    	int startIdx = 0;
    	int idx = str.Find(split, startIdx);
    	aStr.RemoveAll();//先清空
    
    	while (-1 != idx)
    	{
    		CString sTmp = str.Mid(startIdx, idx - startIdx);
    		aStr.Add(sTmp);
    		startIdx = idx + 1;
    		idx = str.Find(split, startIdx);
    	}
    	CString sTmp = str.Right(str.GetLength() - startIdx);
    	if (! sTmp.IsEmpty())
    		aStr.Add(sTmp);
    	return aStr.GetSize();
    }
    /***************************************************************************
    *  名字:double length(double x1,double y1,double x2,double y2)            *
    *                                                                          *
    *  描述:由(x1,y1)和(x2,y2)计算两点之间距离 长度                         *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年3月20日       创建该函数          ***                      *
    *  参数: 1.double x1                                                      *
    *         2.double y1                                                      *
    *         3.double x2                                                      *
    *         4.double y2                                                      *
    *  返回值:double类型数据 返回距离                                         *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    double CSupport::length(double x1,double y1,double x2,double y2)
    {
    	double tmp=((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
    	return sqrt(tmp);
    }
    
    
    /***************************************************************************
    *  名字:void CSupport::read(CString &strExample,CSupport *GCP)            *
    *                                                                          *
    *  描述:读取数据,将数据存至GCP开辟的动态数组中                           *
    *                   调用函数SplitStringArray()                           *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月20日       创建该函数          ***                      *
    *  参数: 1.CString &strExample                                            *
    *         2.CSupport *GCP                                                  *
    *  返回值:无                                                              *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    void CSupport::read(CString &strExample,CSupport *GCP)
    {
    	//将dx dy dz
        int iLine;
    	CStringArray aStrLine;
    	iLine=SplitStringArray(strExample,13,aStrLine);
    	if(iLine==0)
    	{
    		AfxMessageBox(_T("请输入数据!"));
    	}
    	
    	CStringArray aStrTmp;
    	int n;
    
    	n=SplitStringArray(aStrLine[0],',',aStrTmp);
    	int temp1=_tstof(aStrTmp[0]);
    	int temp2=_tstof(aStrTmp[1]);
    	
    	for(int i=0;i<temp1;i++)
    	{
    		n=SplitStringArray(aStrLine[i+1],',',aStrTmp);
    		GCP[i].dLx=_tstof(aStrTmp[0])/1000;
    		GCP[i].dLy=_tstof(aStrTmp[1])/1000;
    		GCP[i].dRx=_tstof(aStrTmp[2])/1000;
    		GCP[i].dRy=_tstof(aStrTmp[3])/1000;
    		GCP[i].X=_tstof(aStrTmp[4]);
    		GCP[i].Y=_tstof(aStrTmp[5]);
    		GCP[i].Z=_tstof(aStrTmp[6]);
    	}
    	for(int i=temp1;i<temp2+temp1;i++)
    	{
    		n=SplitStringArray(aStrLine[i+1],',',aStrTmp);
    		GCP[i].dLx=_tstof(aStrTmp[0])/1000;
    		GCP[i].dLy=_tstof(aStrTmp[1])/1000;
    		GCP[i].dRx=_tstof(aStrTmp[2])/1000;
    		GCP[i].dRy=_tstof(aStrTmp[3])/1000;
    	}
    }
    
    
    /***************************************************************************
    *  名字:int CSupport::SpaceBack()                                        *
    *                                                                          *
    *  描述:后方交会函数,调用函数length()、R1()两个函数                       *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月20日       创建该函数          ***                      *
    *  参数: 1.CSupport *GCP           //主要元素在这里                       *
    *         2.int xx                                                         *
    *         3.int yy                                                         *
    *         4.CMatrix &Left          //外方位6个元素                         *
    *         5.bool _L                //用来标识                              *
    *         6.double &Left_a0        //中误差                                *
    *         7.CMatrix &LeftPhotoes_Pre     //精度                            *
    *  返回值:int t        //返回迭代次数                                     *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    int CSupport::SpaceBack(CSupport *GCP,int xx,int yy,CMatrix &Left,bool _L,double &Left_a0,CMatrix &LeftPhotoes_Pre)
    {
    	int temp=1;
    	int t=0;
    	int n=4;//定义n个已知点
    	double f=0.15;//需修改
    	double m;
    	if(_L)
    	{
    	 m=(length(GCP[xx].X,GCP[xx].Y,GCP[yy].X,GCP[yy].Y)/length(GCP[xx].dLx,GCP[xx].dLy,GCP[yy].dLx,GCP[yy].dLy));
    	}
    	else 
    	{
    		m=(length(GCP[xx].X,GCP[xx].Y,GCP[yy].X,GCP[yy].Y)/length(GCP[xx].dRx,GCP[xx].dRy,GCP[yy].dRx,GCP[yy].dRy));
    	}
    	//m=10000;//m赋值
    	double H=m*f;
    	CMatrix A;
    	A.SetSize(2*n,6);
    	CMatrix X;
    	X.SetSize(6,1);
    	CMatrix L;
    	L.SetSize(2*n,1);
    	CMatrix _x;
    	_x.SetSize(6,1);
    
    	
    	//矩阵X的值初始化
    	
    	X(0,0)=(GCP[xx].X+GCP[yy].X)/2;//XS
    	X(1,0)=(GCP[xx].Y+GCP[yy].Y)/2;//YS
    	X(2,0)=((GCP[xx].Z+GCP[yy].Z)/2+m*f);//ZS
    	/*
    	X(0,0)=(GCP[1].X+GCP[2].X+GCP[3].X+GCP[4].X)/4;//XS
    	X(1,0)=(GCP[1].Y+GCP[2].Y+GCP[3].Y+GCP[4].Y)/4;;//YS
    	X(2,0)=((GCP[1].Z+GCP[2].Z+GCP[3].Z+GCP[4].Z)/4+m*f);//ZS*/
    	for (int i=3;i<6;i++)
    	{
    			X(i,0)=0;
    	}
    
    
    	CMatrix R;
    	R.SetSize(3,3);
    	CMatrix Lz;
    	Lz.SetSize(4,1);
    
    	do
    	{
         R1(R,X);
    
        for(int k=0;k<n;k++)
    	{	
    	Lz(k,0)=R(0,2)*(GCP[k].X-X(0,0))+R(1,2)*(GCP[k].Y-X(1,0))+R(2,2)*(GCP[k].Z-X(2,0));
    	if(_L==1){
    	A(2*k,0)=(R(0,0)*f+R(0,2)*GCP[k].dLx)/Lz(k,0);
    	A(2*k,1)=(R(1,0)*f+R(1,2)*GCP[k].dLx)/Lz(k,0);
    	A(2*k,2)=(R(2,0)*f+R(2,2)*GCP[k].dLx)/Lz(k,0);
    	A(2*k,3)=GCP[k].dLy*sin(X(4,0))-(GCP[k].dLx*(GCP[k].dLx*cos(X(5,0))/f-GCP[k].dLy*sin(X(5,0)))+f*cos(X(5,0)))*cos(X(4,0));
    	A(2*k,4)=-f*sin(X(5,0))-GCP[k].dLx*(GCP[k].dLx*sin(X(5,0))+GCP[k].dLy*cos(X(5,0)))/f;
    	A(2*k,5)=GCP[k].dLy;
    	A(2*k+1,0)=(R(0,1)*f+R(0,2)*GCP[k].dLy)/Lz(k,0);
    	A(2*k+1,1)=(R(1,1)*f+R(1,2)*GCP[k].dLy)/Lz(k,0);
    	A(2*k+1,2)=(R(2,1)*f+R(2,2)*GCP[k].dLy)/Lz(k,0);
    	A(2*k+1,3)=-GCP[k].dLx*sin(X(4,0))-(GCP[k].dLy*(GCP[k].dLx*cos(X(5,0))/f+GCP[k].dLy*sin(X(5,0)))-f*sin(X(5,0)))*cos(X(4,0));
        A(2*k+1,4)=-f*cos(X(5,0))-GCP[k].dLy*(GCP[k].dLx*sin(X(5,0))+GCP[k].dLy*cos(X(5,0)))/f;
    	A(2*k+1,5)=-GCP[k].dLx; 
    	
    	
    	double _1=R(0,0)*(GCP[k].X-X(0,0))+R(1,0)*(GCP[k].Y-X(1,0))+R(2,0)*(GCP[k].Z-X(2,0));
    	double _2=R(0,1)*(GCP[k].X-X(0,0))+R(1,1)*(GCP[k].Y-X(1,0))+R(2,1)*(GCP[k].Z-X(2,0));
    	double _3=R(0,2)*(GCP[k].X-X(0,0))+R(1,2)*(GCP[k].Y-X(1,0))+R(2,2)*(GCP[k].Z-X(2,0));
    	L(2*k,0)=GCP[k].dLx+(f*_1)/_3;
    	L(2*k+1,0)=GCP[k].dLy+(f*_2)/_3;
    	}
    	else
    	{
    		
    	A(2*k,0)=(R(0,0)*f+R(0,2)*GCP[k].dRx)/Lz(k,0);
    	A(2*k,1)=(R(1,0)*f+R(1,2)*GCP[k].dRx)/Lz(k,0);
    	A(2*k,2)=(R(2,0)*f+R(2,2)*GCP[k].dRx)/Lz(k,0);
    	A(2*k,3)=GCP[k].dRy*sin(X(4,0))-(GCP[k].dRx*(GCP[k].dRx*cos(X(5,0))/f-GCP[k].dRy*sin(X(5,0)))+f*cos(X(5,0)))*cos(X(4,0));
    	A(2*k,4)=-f*sin(X(5,0))-GCP[k].dRx*(GCP[k].dRx*sin(X(5,0))+GCP[k].dRy*cos(X(5,0)))/f;
    	A(2*k,5)=GCP[k].dRy;
    	A(2*k+1,0)=(R(0,1)*f+R(0,2)*GCP[k].dRy)/Lz(k,0);
    	A(2*k+1,1)=(R(1,1)*f+R(1,2)*GCP[k].dRy)/Lz(k,0);
    	A(2*k+1,2)=(R(2,1)*f+R(2,2)*GCP[k].dRy)/Lz(k,0);
    	A(2*k+1,3)=-GCP[k].dRx*sin(X(4,0))-(GCP[k].dRy*(GCP[k].dRx*cos(X(5,0))/f+GCP[k].dRy*sin(X(5,0)))-f*sin(X(5,0)))*cos(X(4,0));
        A(2*k+1,4)=-f*cos(X(5,0))-GCP[k].dRy*(GCP[k].dRx*sin(X(5,0))+GCP[k].dRy*cos(X(5,0)))/f;
    	A(2*k+1,5)=-GCP[k].dRx; 
    	
    	
    	double _1=R(0,0)*(GCP[k].X-X(0,0))+R(1,0)*(GCP[k].Y-X(1,0))+R(2,0)*(GCP[k].Z-X(2,0));
    	double _2=R(0,1)*(GCP[k].X-X(0,0))+R(1,1)*(GCP[k].Y-X(1,0))+R(2,1)*(GCP[k].Z-X(2,0));
    	double _3=R(0,2)*(GCP[k].X-X(0,0))+R(1,2)*(GCP[k].Y-X(1,0))+R(2,2)*(GCP[k].Z-X(2,0));
    	L(2*k,0)=GCP[k].dRx+(f*_1)/_3;
    	L(2*k+1,0)=GCP[k].dRy+(f*_2)/_3;
    	}
    	
    	}
    
    	_x=(~A*A).Inv()*~A*L;
    	X=X+_x;
    	++t;
    	
    	}while(fabs(_x(0,0))>1e-4||fabs(_x(1,0))>1e-4||fabs(-_x(2,0))>1e-4||fabs(_x(3,0))>4.481e-6||fabs(_x(4,0))>4.481e-6||fabs(_x(5,0))>4.481e-6);
    	//平定精度
    	CMatrix Temp;
    	Temp=~(_x)*(_x);
    	double temp1;
    	temp1=Temp(0,0);
    	Left_a0=temp1/(2*n-6);
    	Left_a0=sqrt(Left_a0);
    
    	CMatrix Qxx;
    	Qxx.SetSize(6,6);
    	Qxx=(~A*A).Inv();
    	for(int i=0;i<6;i++)
    	{
    		LeftPhotoes_Pre(i,0)=Left_a0*sqrt(Qxx(i,i));
    	}
    	Left=X;
    	return t;
    }
    
    
    /***************************************************************************
    *  名字:void CSupport::R1(CMatrix &R,CMatrix X )                          *
    *                                                                          *
    *  描述:R矩阵赋值                                                         *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月20日       创建该函数          ***                      *
    *  参数: 1.CMatrix &R           //R矩阵 (3,3)                          *
    *         2.CMatrix X            // X矩阵(6,1)                          *
    *  返回值:无                                                              *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    void CSupport::R1(CMatrix &R,CMatrix X )
    {
    	R(0,0)=cos(X(3,0))*cos(X(5,0))-sin(X(3,0))*sin(X(4,0))*sin(X(5,0));
    	R(0,1)=-cos(X(3,0))*sin(X(5,0))-sin(X(3,0))*sin(X(4,0))*cos(X(5,0));
    	R(0,2)=-sin(X(3,0))*cos(X(4,0));
    	R(1,0)=cos(X(4,0))*sin(X(5,0));
    	R(1,1)=cos(X(4,0))*cos(X(5,0));
    	R(1,2)=-sin(X(4,0));
    	R(2,0)=sin(X(3,0))*cos(X(5,0))+cos(X(3,0))*sin(X(4,0))*sin(X(5,0));
    	R(2,1)=-sin(X(3,0))*sin(X(5,0))+cos(X(3,0))*sin(X(4,0))*cos(X(5,0));
    	R(2,2)=cos(X(3,0))*cos(X(4,0));
    }
    
    /***************************************************************************
    *  名字:void CSupport::SpaceFront()                                       *
    *                                                                          *
    *  描述:前方交会函数                                                      *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月20日       创建该函数          ***                      *
    *  参数: 1.CMatrix &Left           //左片6个外方位元素                    *
    *         2.CMatrix &Right          //右片6个外方位元素                    *
    *         3.CMatrix &B              //摄影基线3个元素                      *
    *         4.CMatrix &R_Left          //左片旋转矩阵(3,3)                *
    *         5.CMatrix &R_Right         //右片旋转矩阵(3,3)                *
    *         6.CSupport *GCP           //主要元素                             *
    *  返回值:无                                                              *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    void CSupport::SpaceFront(CMatrix &Left,CMatrix &Right,CMatrix &B,CMatrix &R_Left,CMatrix &R_Right,CSupport *GCP)
    {
    	double f=0.15;//需修改
    	int n=4;//从第五行起开始为计算值
    	int m=9;//5--9行为数据
    	B(0,0)=Right(0,0)-Left(0,0);
    	B(1,0)=Right(1,0)-Left(1,0);
    	B(2,0)=Right(2,0)-Left(2,0);
    
    	R1(R_Left,Left);
    	R1(R_Right,Right);
    	//开辟一个数组
    	
    	for(int i=n;i<m;i++)
    	{
    		CMatrix L_xyf;
    		L_xyf.SetSize(3,1);
    		L_xyf(0,0)=GCP[i].dLx;
    		L_xyf(1,0)=GCP[i].dLy;
    		L_xyf(2,0)=-f;
    
    		CMatrix R_xyf;
    		R_xyf.SetSize(3,1);
    		R_xyf(0,0)=GCP[i].dRx;
    		R_xyf(1,0)=GCP[i].dRy;
    		R_xyf(2,0)=-f;
    
    		
    		CMatrix L;
    		L.SetSize(3,1);
    		CMatrix R;
    		R.SetSize(3,1);
    		L=R_Left*L_xyf;
    		R=R_Right*R_xyf;
    
    		double N1,N2;
    		N1=((B(0,0)*R(2,0)-B(2,0)*R(0,0))/(L(0,0)*R(2,0)-R(0,0)*L(2,0)));
    		N2=((B(0,0)*L(2,0)-B(2,0)*L(0,0))/(L(0,0)*R(2,0)-R(0,0)*L(2,0)));
    
    		GCP[i].X=Right(0,0)+N2*R(0,0);
    		GCP[i].Y=(Left(1,0)+Right(1,0)+N1*L(1,0)+N2*R(1,0))/2;
    		GCP[i].Z=Right(2,0)+N2*R(2,0);
    	}
    }
    /***************************************************************************
    *  名字:void CSupport::out()                                              *
    *                                                                          *
    *  描述:输出函数                                                          *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月20日       创建该函数          ***                      *
    *  参数: 1.double Left_a0          //左片单位权中误差                     *
    *         2.double Right_a0          //左片单位权中误差                    *
    *         3.CMatrix LeftPhotoes             //左片外方位元素               *
    *         4.CMatrix LeftPhotoes_Pre          //左片外方位元素精度          *
    *         5.CMatrix RightPhotoes            //右片外方位元素               *
    *         6.CMatrix RightPhotoes_Pre         //右片外方位元素精度          *
    *         7.CSupport *GCP                  //主要元素                      *
    *         8.CString &strResult                //输出结果CString            *
    *         9.int t1                           //左片迭代次数                *
    *         10.int t2                          //右片迭代次数                *
    *  返回值:无                                                              *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    void CSupport::out(double Left_a0,double Right_a0,CMatrix LeftPhotoes,CMatrix LeftPhotoes_Pre,CMatrix RightPhotoes,CMatrix RightPhotoes_Pre,CSupport *GCP ,CString &strResult,int t1,int t2)
    {
    		strResult.Format(_T("%s\t%s\t%s\t\r\n"),
    		_T("-------------------------------------------------------- 摄影外方位元素  -------------------------------------------------------- "),
    		_T(""),
    		_T("")
    		);
    	
    	CString strOutput;
    	strResult=strResult+strOutput;
    	strOutput.Format(_T("%s%f%s%d\r\n%s\r\n"),
    		_T("左片单位权中误差:"),
    		Left_a0,
    		_T("                                  迭代次数:"),
    		t1,
    		_T("左片元素:依次是 Xs  Ys  Zs   Alfa   w   k:                                      精度:")
    		);
    	strResult=strResult+strOutput;
    		for(int z=0;z<6;z++)
    	{
    		if(z<3)
    		{
    			strOutput.Format(_T("%f\t\t\t\t\t\t%f\r\n"),LeftPhotoes(z,0),LeftPhotoes_Pre(z,0));
    		}
    		else
    		{
    			strOutput.Format(_T("%f\t\t\t\t\t\t\t%f\r\n"),LeftPhotoes(z,0),LeftPhotoes_Pre(z,0));
    		}
    		strResult=strResult+strOutput;
    	}
    	strOutput.Format(_T("%s%f%s%d\r\n%s\r\n"),
    		_T("右片单位权中误差:"),
    		Right_a0,
    		_T("                                  迭代次数:"),
    		t2,
    		_T("右片元素:依次是 Xs  Ys  Zs   Alfa   w   k:                                      精度:")
    		);
    	strResult=strResult+strOutput;
    	for(int z=0;z<6;z++)
    	{
    		if(z<3)
    		{
    		strOutput.Format(_T("%f\t\t\t\t\t\t%f\r\n"),RightPhotoes(z,0),RightPhotoes_Pre(z,0));
    		}
    		else
    		{
    		strOutput.Format(_T("%f\t\t\t\t\t\t\t%f\r\n"),RightPhotoes(z,0),RightPhotoes_Pre(z,0));
    		}
    		strResult=strResult+strOutput;
    
    	}
    	strOutput.Format(_T("%s\r\n%s\r\n"),
    		_T("全部元素:"),
    		_T("左片x     左片y      右片x      右片y              X                       Y                     Z    ")
    		);
    	strResult=strResult+strOutput;	
    	for(int i=0;i<9;i++)
    	{
    		strOutput.Format(_T("%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t\r\n"),GCP[i].dLx*1000,GCP[i].dLy*1000,GCP[i].dRx*1000,GCP[i].dRy*1000,GCP[i].X,GCP[i].Y,GCP[i].Z);
    		strResult=strResult+strOutput;
    	}
    
    }
    
    
    /***************************************************************************
    *  名字:void CSupport::main(CString &strExample,CString &strResult)       *
    *                                                                          *
    *  描述:主要函数 ,调用函数read()、SpqceBack()、SpaceFront()、out()等函数 *
    *                    等4个函数                                             *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月20日       创建该函数          ***                      *
    *  参数: 1.CString &strExample          //读入数据                        *
    *         2.CString &strResult          //输出数据                         *
    *  返回值:无                                                              *
    *                                                                          *
    *  注:                                                                    *
    /**************************************************************************/
    void CSupport::main(CString &strExample,CString &strResult)
    {
    	
    	CSupport GCP[9];//每一行数据各个元素提取
    	read(strExample,GCP);//读取函数
    	
    	
    	CMatrix LeftPhotoes;//左片外方位元素
    	LeftPhotoes.SetSize(6,1);
    	double Left_a0;//中误差
    	CMatrix LeftPhotoes_Pre;//左片外方位元素
    	LeftPhotoes_Pre.SetSize(6,1);
    
    	CMatrix RightPhotoes;//右片外方位元素
    	RightPhotoes.SetSize(6,1);
    	double Right_a0;//中误差
    	CMatrix RightPhotoes_Pre;//左片外方位元素
    	RightPhotoes_Pre.SetSize(6,1);
    	
    	double LeftPhotoes_t1=SpaceBack(GCP,0,2,LeftPhotoes,1,Left_a0,LeftPhotoes_Pre);//左片后方交会
    	double RightPhotoes_t2=SpaceBack(GCP,1,3,RightPhotoes,0,Right_a0,RightPhotoes_Pre);//右片后方交会
    
    	
    
    
    	CMatrix B;
    	B.SetSize(3,1);//Bx By Bz
    	CMatrix R_Left;
    	R_Left.SetSize(3,3);//左片旋转矩阵R
    	CMatrix R_Right;                                                     
    	R_Right.SetSize(3,3);//右片旋转矩阵R
    	SpaceFront(LeftPhotoes,RightPhotoes,B,R_Left,R_Right,GCP);
    
    	out(Left_a0,Right_a0,LeftPhotoes,LeftPhotoes_Pre,RightPhotoes,RightPhotoes_Pre,GCP,strResult,LeftPhotoes_t1,RightPhotoes_t2);
    }
    
    CSupport::CSupport(void)
    {
    }
    
    
    CSupport::~CSupport(void)
    {
    }
    
    

    3.3.5文件:RS_110_Z_SpaceIntersectionDlg.cpp (仅摘取部分)

    /***************************************************************************
    *  文件名:<RS_110_Z_SpaceIntersectionDlg.cpp>                            *
    *                                                                          *
    *  描述: Dialog对话框使用                                                  *
    *                                                                          *
    *  历史:**日期**         **理由**            **签名**                     *
    *      2019年4月19日        创建              ***                       *
    *                                                                          *
    *  外部过程:                                                              *
    *                                                                          *
    /**************************************************************************/
    //计算按钮
    void CRS_110_Z_SpaceIntersectionDlg::OnBnClickedOk()
    {
    
    	// TODO: 在此添加控件通知处理程序代码
    	//CDialogEx::OnOK();
    	UpdateData(true);
    	CSupport k;
    	k.main(strExample,strResult);
    	UpdateData(false);
    }
    
    
    //取消按钮
    void CRS_110_Z_SpaceIntersectionDlg::OnBnClickedCancel()
    {
    	// TODO: 在此添加控件通知处理程序代码
    	CDialogEx::OnCancel();
    }
    
    

    3.4运行结果

    在这里插入图片描述
    在这里插入图片描述

    3.5 设计技巧

     创建多个类,层次分明。调用了Cmatrix矩阵类,省了好多力气。
     使用指针动态开辟数组,较为方便。
     提前做好函数模块,结构化,有助于专注于编程。

    代码虽多不要贪杯~

    展开全文
  • 本代码提供了摄影测量中单像空间后方交会解算外方位元素和立体像对前方交会解算地面点坐标的详细方法,对于摄影测量学专业的学生来说无疑是一个不可或缺的好资料,希望对大家有所帮助!
  • 摄影测量后方交会-前方交会(C#)

    万次阅读 多人点赞 2015-09-05 20:39:45
    解算过程分为两个阶段:后方交会前方交会。 思路为利用后方交会计算外方位元素,利用前方交会解算待定地面点坐标。 后方交会的解算过程为: (1)获取已知数据:比例尺分母m,内方位元素x0,y0,f,以及控制点的...

    双像解析摄影测量中的后交-前交解法,常用于已知像片的外方位元素、需确定少量待定点坐标的情况。解算过程分为两个阶段:后方交会、前方交会。
    思路为利用后方交会计算外方位元素,利用前方交会解算待定地面点坐标。
    后方交会的解算过程为:
    (1)获取已知数据:比例尺分母m,内方位元素x0,y0,f,以及控制点的地面摄影测量坐标X、Y、Z。
    (2)量测控制点的像点坐标:得到像点坐标x、y。
    (3)确定未知数的初始值:在竖直摄影情况下,外方位角元素的初始值为0,即phi0=omiga0=kappa0=0;
    线元素中,Zs0=H=mf,Xs0、Ys0的取值用四个控制点坐标的平均值,即Xs0=1/4*(X1+X2+X3+X4);Ys0=1/4 *(Y1+Y2+Y3+Y4).
    (4)计算旋转矩阵R:利用角元素的近似值计算方向余弦值,组成R阵。
    (5)逐点计算像点坐标的近似值:利用未知数的近似值按共线方程计算控制点像点坐标的近似值(x),(y)。
    (6)组成误差方程式:逐点计算误差方程式的系数和常数项。
    (7)组成法方程式:计算法方程式的系数矩阵ATA与常数项ATL。
    (8)解求外方位元素:根据法方程,解求外方位元素改正数,并与相应的近似值求和,得到外方位元素新的近似值。
    (9)检查计算是否收敛:将求得的外方位元素的改正数与规定的限差比较,小于限差则计算终止,否则用新的近似值重复第4至第8步骤的计算,直到满足要求为止。

    前方交会的解算过程:
    (1)按照后方交会过程,计算出左右双片的外方位元素,用各片的外方位角元素计算左右片的方向余弦值,组成旋转矩阵R1和R2。
    (2)逐点计算像点的像空间辅助坐标(u1,v1,w1)及(u2,v2,w2)。
    (3)根据外方位线元素计算基线分量(Bu,Bv,Bw)。
    (4)计算投影系数(N1,N2)。
    (5)计算待定点在各自的像空间辅助坐标中的坐标(U1,V1,W1)及(U2,V2,W2)。
    (6)最后计算待定点的地面摄影测量坐标(X,Y,Z)。

    以下为C#实现代码:

    using System;
    using System.Collections.Generic;
    using System.ComponentModel;
    using System.Data;
    using System.Drawing;
    using System.Linq;
    using System.Text;
    using System.Windows.Forms;
    
    namespace Resection
    {
        public partial class Form1 : Form
        {
            public Form1()
            {
                InitializeComponent();
            }
            public TextBox[] mytextbox = new TextBox[41];
    
            //矩阵打包成类,矩阵为m * n,直接调用
            public class Matrix
            {
                double[,] A;
                int m, n;
                string name;
                public Matrix(int am, int an)
                {
                    m = am;
                    n = an;
                    A = new double[m, n];
                    name = "Result";
                }
                public Matrix(int am, int an, string aName)
                {
                    m = am;
                    n = an;
                    A = new double[m, n];
                    name = aName;
                }
    
                public int getM
                {
                    get { return m; }
                }
                public int getN
                {
                    get { return n; }
                }
                public double[,] Detail
                {
                    get { return A; }
                    set { A = value; }
                }
                public string Name
                {
                    get { return name; }
                    set { name = value; }
                }
            }
    
            /***********矩阵通用操作打包*************/
    
            class MatrixOperator
            {
                //矩阵加法
                public static Matrix MatrixAdd(Matrix Ma, Matrix Mb)
                {
                    int m = Ma.getM;
                    int n = Ma.getN;
                    int m2 = Mb.getM;
                    int n2 = Mb.getN;
    
                    if ((m != m2) || (n != n2))
                    {
                        Exception myException = new Exception("数组维数不匹配");
                        throw myException;
                    }
    
                    Matrix Mc = new Matrix(m, n);
                    double[,] c = Mc.Detail;
                    double[,] a = Ma.Detail;
                    double[,] b = Mb.Detail;
    
                    for (int i = 0; i < m; i++)
                        for (int j = 0; j < n; j++)
                            c[i, j] = a[i, j] + b[i, j];
                    return Mc;
                }
    
                //矩阵减法
                public static Matrix MatrixSub(Matrix Ma, Matrix Mb)
                {
                    int m = Ma.getM;
                    int n = Ma.getN;
                    int m2 = Mb.getM;
                    int n2 = Mb.getN;
                    if ((m != m2) || (n != n2))
                    {
                        Exception myException = new Exception("数组维数不匹配");
                        throw myException;
                    }
                    Matrix Mc = new Matrix(m, n);
                    double[,] c = Mc.Detail;
                    double[,] a = Ma.Detail;
                    double[,] b = Mb.Detail;
    
                    for (int i = 0; i < m; i++)
                        for (int j = 0; j < n; j++)
                            c[i, j] = a[i, j] - b[i, j];
                    return Mc;
                }
    
                //矩阵乘法
                public static Matrix MatrixMulti(Matrix Ma, Matrix Mb)
                {
                    int m = Ma.getM;
                    int n = Ma.getN;
                    int m2 = Mb.getM;
                    int n2 = Mb.getN;
    
                    if (n != m2)
                    {
                        Exception myException = new Exception("数组维数不匹配");
                        throw myException;
                    }
    
                    Matrix Mc = new Matrix(m, n2);
                    double[,] c = Mc.Detail;
                    double[,] a = Ma.Detail;
                    double[,] b = Mb.Detail;
    
                    for (int i = 0; i < m; i++)
                        for (int j = 0; j < n2; j++)
                        {
                            c[i, j] = 0;
                            for (int k = 0; k < n; k++)
                                c[i, j] += a[i, k] * b[k, j];
                        }
                    return Mc;
    
                }
    
                //矩阵数乘
                public static Matrix MatrixSimpleMulti(double k, Matrix Ma)
                {
                    int m = Ma.getM;
                    int n = Ma.getN;
                    Matrix Mc = new Matrix(m, n);
                    double[,] c = Mc.Detail;
                    double[,] a = Ma.Detail;
    
                    for (int i = 0; i < m; i++)
                        for (int j = 0; j < n; j++)
                            c[i, j] = a[i, j] * k;
                    return Mc;
                }
    
                //矩阵转置
                public static Matrix MatrixTrans(Matrix MatrixOrigin)
                {
                    int m = MatrixOrigin.getM;
                    int n = MatrixOrigin.getN;
                    Matrix MatrixNew = new Matrix(n, m);
                    double[,] c = MatrixNew.Detail;
                    double[,] a = MatrixOrigin.Detail;
                    for (int i = 0; i < n; i++)
                        for (int j = 0; j < m; j++)
                            c[i, j] = a[j, i];
                    return MatrixNew;
                }
    
                //矩阵求逆(伴随矩阵法)
                public static Matrix MatrixInvByCom(Matrix Ma)
                {
                    double d = MatrixOperator.MatrixDet(Ma);
                    if (d == 0)
                    {
                        Exception myException = new Exception("没有逆矩阵");
                        throw myException;
                    }
                    Matrix Ax = MatrixOperator.MatrixCom(Ma);
                    Matrix An = MatrixOperator.MatrixSimpleMulti((1.0 / d), Ax);
                    return An;
                }
                //对应行列式的代数余子式矩阵
                public static Matrix MatrixSpa(Matrix Ma, int ai, int aj)
                {
                    int m = Ma.getM;
                    int n = Ma.getN;
                    if (m != n)
                    {
                        Exception myException = new Exception("矩阵不是方阵");
                        throw myException;
                    }
                    int n2 = n - 1;
                    Matrix Mc = new Matrix(n2, n2);
                    double[,] a = Ma.Detail;
                    double[,] b = Mc.Detail;
    
                    //左上
                    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 Mc;
    
                }
    
                //矩阵的行列式,矩阵必须是方阵
                public static double MatrixDet(Matrix Ma)
                {
                    int m = Ma.getM;
                    int n = Ma.getN;
                    if (m != n)
                    {
                        Exception myException = new Exception("数组维数不匹配");
                        throw myException;
                    }
                    double[,] a = Ma.Detail;
                    if (n == 1) return a[0, 0];
    
                    double D = 0;
                    for (int i = 0; i < n; i++)
                    {
                        D += a[1, i] * MatrixDet(MatrixSpa(Ma, 1, i));
                    }
                    return D;
                }
    
                //矩阵的伴随矩阵
                public static Matrix MatrixCom(Matrix Ma)
                {
                    int m = Ma.getM;
                    int n = Ma.getN;
                    Matrix Mc = new Matrix(m, n);
                    double[,] c = Mc.Detail;
                    double[,] a = Ma.Detail;
    
                    for (int i = 0; i < m; i++)
                        for (int j = 0; j < n; j++)
                            c[i, j] = MatrixDet(MatrixSpa(Ma, j, i));
    
                    return Mc;
                }
            }
    
            private void button1_Click(object sender, EventArgs e)      //计算左片后方交会
            {
                double[] x = { Convert.ToDouble(t1.Text), Convert.ToDouble(t6.Text), Convert.ToDouble(t11.Text), Convert.ToDouble(t16.Text) };
                double[] y = { Convert.ToDouble(t2.Text), Convert.ToDouble(t7.Text), Convert.ToDouble(t12.Text), Convert.ToDouble(t17.Text) };
                double[] X = { Convert.ToDouble(t3.Text), Convert.ToDouble(t8.Text), Convert.ToDouble(t13.Text), Convert.ToDouble(t18.Text) };
                double[] Y = { Convert.ToDouble(t4.Text), Convert.ToDouble(t9.Text), Convert.ToDouble(t14.Text), Convert.ToDouble(t19.Text) };
                double[] Z = { Convert.ToDouble(t5.Text), Convert.ToDouble(t10.Text), Convert.ToDouble(t15.Text), Convert.ToDouble(t20.Text) };
                double f = Convert.ToDouble(t22.Text) / 1000, _m = Convert.ToDouble(t21.Text);
                for (int i = 0; i < 4; i++)    //单位换算
                {
                    x[i] = x[i] / 1000;
                    y[i] = y[i] / 1000;
                }
    
                /////定义外方位元素,并附初值
                double Xs, Ys, Zs, phi = 0, omiga = 0, kappa = 0;
                Xs = (X[0] + X[1] + X[2] + X[3]) / 4.0;
                Ys = (Y[0] + Y[1] + Y[2] + Y[3]) / 4.0;
                Zs = _m * f;
    
                /////定义x,y近似值,即计算值
                double[] _x = new double[4];
                double[] _y = new double[4];
    
                /////定义共线方程中的分子分母项,便于计算
                double[] _X = new double[4];
                double[] _Y = new double[4];
                double[] _Z = new double[4];
    
                /////定义旋转矩阵R
                double[,] R = new double[3, 3];
                double[,] L = new double[8, 1];//误差方程中的常数项
                double[,] XX = new double[6, 1];//X向量
    
                /////开始迭代
                do
                {
                    /////计算旋转矩阵
                    R[0, 0] = Math.Cos(phi) * Math.Cos(kappa) - Math.Sin(phi) * Math.Sin(omiga) * Math.Sin(kappa);//a1
                    R[0, 1] = -Math.Cos(phi) * Math.Sin(kappa) - Math.Sin(phi) * Math.Sin(omiga) * Math.Cos(kappa);//a2
                    R[0, 2] = -Math.Sin(phi) * Math.Cos(omiga);//a3
                    R[1, 0] = Math.Cos(omiga) * Math.Sin(kappa);//b1
                    R[1, 1] = Math.Cos(omiga) * Math.Cos(kappa);//b2
                    R[1, 2] = -Math.Sin(omiga);//b3
                    R[2, 0] = Math.Sin(phi) * Math.Cos(kappa) + Math.Cos(phi) * Math.Sin(omiga) * Math.Sin(kappa);//c1
                    R[2, 1] = -Math.Sin(phi) * Math.Sin(kappa) + Math.Cos(phi) * Math.Sin(omiga) * Math.Cos(kappa);//c2
                    R[2, 2] = Math.Cos(phi) * Math.Cos(omiga);//c3
    
                    for (int i = 0; i < 4; i++)
                    {
                        //用共线方程计算 x,y 的近似值 ,即计算值       
                        _X[i] = R[0, 0] * (X[i] - Xs) + R[1, 0] * (Y[i] - Ys) + R[2, 0] * (Z[i] - Zs);
                        _Y[i] = R[0, 1] * (X[i] - Xs) + R[1, 1] * (Y[i] - Ys) + R[2, 1] * (Z[i] - Zs);
                        _Z[i] = R[0, 2] * (X[i] - Xs) + R[1, 2] * (Y[i] - Ys) + R[2, 2] * (Z[i] - Zs);
    
                        _x[i] = -f * _X[i] / _Z[i];
                        _y[i] = -f * _Y[i] / _Z[i]; 
                    }
    
                    Matrix B = new Matrix(8, 6, "B");//4个控制点,每个是2行6列,4个是8行6列
                    int n = B.getN;
                    int m = B.getM;
                    double[,] b = B.Detail;
                    for (int i = 0; i < 4; i++)
                    {
                        //计算系数矩阵
                        b[2 * i, 0] = (R[0, 0] * f + R[0, 2] * x[i]) / _Z[i];
                        b[2 * i, 1] = (R[1, 0] * f + R[1, 2] * x[i]) / _Z[i];
                        b[2 * i, 2] = (R[2, 0] * f + R[2, 2] * x[i]) / _Z[i];
                        b[2 * i, 3] = y[i] * Math.Sin(omiga) - ((x[i] / f) * (x[i] * Math.Cos(kappa) - y[i] * Math.Sin(kappa)) + f * Math.Cos(kappa)) * Math.Cos(omiga);
                        b[2 * i, 4] = -f * Math.Sin(kappa) - (x[i] / f) * (x[i] * Math.Sin(kappa) + y[i] * Math.Cos(kappa));
                        b[2 * i, 5] = y[i];
    
                        b[2 * i + 1, 0] = (R[0, 1] * f + R[0, 2] * y[i]) / _Z[i];
                        b[2 * i + 1, 1] = (R[1, 1] * f + R[1, 2] * y[i]) / _Z[i];
                        b[2 * i + 1, 2] = (R[2, 1] * f + R[2, 2] * y[i]) / _Z[i];
                        b[2 * i + 1, 3] = -x[i] * Math.Sin(omiga) - ((x[i] / f) * (x[i] * Math.Cos(kappa) - y[i] * Math.Sin(kappa)) - f * Math.Sin(kappa)) * Math.Cos(omiga);
                        b[2 * i + 1, 4] = -f * Math.Cos(kappa) - (y[i] / f) * (x[i] * Math.Sin(kappa) + y[i] * Math.Cos(kappa));
                        b[2 * i + 1, 5] = -x[i];
    
                        //计算常数项
                        L[2 * i, 0] = x[i] - _x[i];
                        L[2 * i + 1, 0] = y[i] - _y[i];
                    }
                    Matrix C = MatrixOperator.MatrixTrans(B);          //系数矩阵的转置矩阵
                    C.Name = "C";
                    Matrix D = MatrixOperator.MatrixMulti(C, B);       //系数矩阵与其转置矩阵相乘
                    D.Name = "C*B";
                    Matrix dn = MatrixOperator.MatrixInvByCom(D);      //系数矩阵与其转置矩阵积的逆矩阵
                    dn.Name = "dn";
                    Matrix dn2 = MatrixOperator.MatrixMulti(dn, C);       //ATA的逆阵乘以A的转置
                    dn2.Name = "dn2";
                    double[,] ATARAT = dn2.Detail;
    
                    ////计算ATARAT * L,存在XX中
                    for (int i = 0; i < 6; i++)
                        for (int j = 0; j < 1; j++)
                        {
                            XX[i, j] = 0;
                            for (int l = 0; l < 8; l++)
                                XX[i, j] += ATARAT[i, l] * L[l, 0];
                        }
    
    
                    ////计算外方位元素值
                    Xs += XX[0, 0];
                    Ys += XX[1, 0];
                    Zs += XX[2, 0];
                    phi += XX[3, 0];
                    omiga += XX[4, 0];
                    kappa += XX[5, 0];
    
    
                }
                while (Math.Abs(XX[0, 0]) >= 0.000029 || Math.Abs(XX[1, 0]) >= 0.000029 || Math.Abs(XX[2, 0]) >= 0.000029 || Math.Abs(XX[3, 0]) >= 1000 * 0.000029 || Math.Abs(XX[4, 0]) >= 1000 * 0.000029 || Math.Abs(XX[5, 0]) >= 1000 * 0.000029);
    
                t23.Text = Convert.ToString(Xs);                         //输出到文本框
                t24.Text = Convert.ToString(Ys);
                t25.Text = Convert.ToString(Zs);
                t26.Text = Convert.ToString(phi);
                t27.Text = Convert.ToString(omiga);
                t28.Text = Convert.ToString(kappa);    
            }   
    
            private void button2_Click(object sender, EventArgs e)     //清空数据,并输入右片数据
            {
                for (int i = 0; i < 20; i++)
                {
                    mytextbox[i].Text = "";
                }
            }
    
            private void Form1_Load(object sender, EventArgs e)       //文本框编译成数组,方便调用
            {
                int i = 0, j = 0;
                TextBox t;
                foreach (Control c in this.Controls)                 
                {
                    if (c is TextBox)
                    {
                        mytextbox[i] = (TextBox)c;
                        i++;
                    }
                }
                for (i = 0; i < 41; i++)
                {
                    for (j = 0; j < 41; j++)
                    {
                        if (mytextbox[j].Name == "t" + Convert.ToString(i + 1))
                        {
                            t = mytextbox[i];
                            mytextbox[i] = mytextbox[j];
                            mytextbox[j] = t;
                           // mytextbox[i].Text = Convert.ToString(i);
                        }
                    }
                    //mytextbox[i].Text = " ";
                }
            }
    
            private void button3_Click(object sender, EventArgs e)     //计算右片后方交会
            {
                double[] x = { Convert.ToDouble(t1.Text), Convert.ToDouble(t6.Text), Convert.ToDouble(t11.Text), Convert.ToDouble(t16.Text) };
                double[] y = { Convert.ToDouble(t2.Text), Convert.ToDouble(t7.Text), Convert.ToDouble(t12.Text), Convert.ToDouble(t17.Text) };
                double[] X = { Convert.ToDouble(t3.Text), Convert.ToDouble(t8.Text), Convert.ToDouble(t13.Text), Convert.ToDouble(t18.Text) };
                double[] Y = { Convert.ToDouble(t4.Text), Convert.ToDouble(t9.Text), Convert.ToDouble(t14.Text), Convert.ToDouble(t19.Text) };
                double[] Z = { Convert.ToDouble(t5.Text), Convert.ToDouble(t10.Text), Convert.ToDouble(t15.Text), Convert.ToDouble(t20.Text) };
                double f = Convert.ToDouble(t22.Text) / 1000, _m = Convert.ToDouble(t21.Text);
                for (int i = 0; i < 4; i++)    //单位换算
                {
                    x[i] = x[i] / 1000;
                    y[i] = y[i] / 1000;
                }
    
                /////定义外方位元素,并附初值
                double Xs, Ys, Zs, phi = 0, omiga = 0, kappa = 0;
                Xs = (X[0] + X[1] + X[2] + X[3]) / 4.0;
                Ys = (Y[0] + Y[1] + Y[2] + Y[3]) / 4.0;
                Zs = _m * f;
    
                /////定义x,y近似值,即计算值
                double[] _x = new double[4];
                double[] _y = new double[4];
    
                /////定义共线方程中的分子分母项,便于计算
                double[] _X = new double[4];
                double[] _Y = new double[4];
                double[] _Z = new double[4];
    
                /////定义旋转矩阵R
                double[,] R = new double[3, 3];
                double[,] L = new double[8, 1];//误差方程中的常数项
                double[,] XX = new double[6, 1];//X向量
    
                /////开始迭代
                do
                {
                    /////计算旋转矩阵
                    R[0, 0] = Math.Cos(phi) * Math.Cos(kappa) - Math.Sin(phi) * Math.Sin(omiga) * Math.Sin(kappa);//a1
                    R[0, 1] = -Math.Cos(phi) * Math.Sin(kappa) - Math.Sin(phi) * Math.Sin(omiga) * Math.Cos(kappa);//a2
                    R[0, 2] = -Math.Sin(phi) * Math.Cos(omiga);//a3
                    R[1, 0] = Math.Cos(omiga) * Math.Sin(kappa);//b1
                    R[1, 1] = Math.Cos(omiga) * Math.Cos(kappa);//b2
                    R[1, 2] = -Math.Sin(omiga);//b3
                    R[2, 0] = Math.Sin(phi) * Math.Cos(kappa) + Math.Cos(phi) * Math.Sin(omiga) * Math.Sin(kappa);//c1
                    R[2, 1] = -Math.Sin(phi) * Math.Sin(kappa) + Math.Cos(phi) * Math.Sin(omiga) * Math.Cos(kappa);//c2
                    R[2, 2] = Math.Cos(phi) * Math.Cos(omiga);//c3
    
                    for (int i = 0; i < 4; i++)
                    {
                        //用共线方程计算 x,y 的近似值 ,即计算值       
                        _X[i] = R[0, 0] * (X[i] - Xs) + R[1, 0] * (Y[i] - Ys) + R[2, 0] * (Z[i] - Zs);
                        _Y[i] = R[0, 1] * (X[i] - Xs) + R[1, 1] * (Y[i] - Ys) + R[2, 1] * (Z[i] - Zs);
                        _Z[i] = R[0, 2] * (X[i] - Xs) + R[1, 2] * (Y[i] - Ys) + R[2, 2] * (Z[i] - Zs);
    
                        _x[i] = -f * _X[i] / _Z[i];
                        _y[i] = -f * _Y[i] / _Z[i];
                    }
    
                    Matrix B = new Matrix(8, 6, "B");//4个控制点,每个是2行6列,4个是8行6列
                    int n = B.getN;
                    int m = B.getM;
                    double[,] b = B.Detail;
                    for (int i = 0; i < 4; i++)
                    {
                        //计算系数矩阵
                        b[2 * i, 0] = (R[0, 0] * f + R[0, 2] * x[i]) / _Z[i];
                        b[2 * i, 1] = (R[1, 0] * f + R[1, 2] * x[i]) / _Z[i];
                        b[2 * i, 2] = (R[2, 0] * f + R[2, 2] * x[i]) / _Z[i];
                        b[2 * i, 3] = y[i] * Math.Sin(omiga) - ((x[i] / f) * (x[i] * Math.Cos(kappa) - y[i] * Math.Sin(kappa)) + f * Math.Cos(kappa)) * Math.Cos(omiga);
                        b[2 * i, 4] = -f * Math.Sin(kappa) - (x[i] / f) * (x[i] * Math.Sin(kappa) + y[i] * Math.Cos(kappa));
                        b[2 * i, 5] = y[i];
    
                        b[2 * i + 1, 0] = (R[0, 1] * f + R[0, 2] * y[i]) / _Z[i];
                        b[2 * i + 1, 1] = (R[1, 1] * f + R[1, 2] * y[i]) / _Z[i];
                        b[2 * i + 1, 2] = (R[2, 1] * f + R[2, 2] * y[i]) / _Z[i];
                        b[2 * i + 1, 3] = -x[i] * Math.Sin(omiga) - ((x[i] / f) * (x[i] * Math.Cos(kappa) - y[i] * Math.Sin(kappa)) - f * Math.Sin(kappa)) * Math.Cos(omiga);
                        b[2 * i + 1, 4] = -f * Math.Cos(kappa) - (y[i] / f) * (x[i] * Math.Sin(kappa) + y[i] * Math.Cos(kappa));
                        b[2 * i + 1, 5] = -x[i];
    
                        //计算常数项
                        L[2 * i, 0] = x[i] - _x[i];
                        L[2 * i + 1, 0] = y[i] - _y[i];
                    }
                    Matrix C = MatrixOperator.MatrixTrans(B);          //系数矩阵的转置矩阵
                    C.Name = "C";
                    Matrix D = MatrixOperator.MatrixMulti(C, B);       //系数矩阵与其转置矩阵相乘
                    D.Name = "C*B";
                    Matrix dn = MatrixOperator.MatrixInvByCom(D);      //系数矩阵与其转置矩阵积的逆矩阵
                    dn.Name = "dn";
                    Matrix dn2 = MatrixOperator.MatrixMulti(dn, C);       //ATA的逆阵乘以A的转置
                    dn2.Name = "dn2";
                    double[,] ATARAT = dn2.Detail;
    
                    ////计算ATARAT * L,存在XX中
                    for (int i = 0; i < 6; i++)
                        for (int j = 0; j < 1; j++)
                        {
                            XX[i, j] = 0;
                            for (int l = 0; l < 8; l++)
                                XX[i, j] += ATARAT[i, l] * L[l, 0];
                        }
    
    
                    ////计算外方位元素值
                    Xs += XX[0, 0];
                    Ys += XX[1, 0];
                    Zs += XX[2, 0];
                    phi += XX[3, 0];
                    omiga += XX[4, 0];
                    kappa += XX[5, 0];
    
    
                }
                while (Math.Abs(XX[0, 0]) >= 0.000029 || Math.Abs(XX[1, 0]) >= 0.000029 || Math.Abs(XX[2, 0]) >= 0.000029 || Math.Abs(XX[3, 0]) >= 1000 * 0.000029 || Math.Abs(XX[4, 0]) >= 1000 * 0.000029 || Math.Abs(XX[5, 0]) >= 1000 * 0.000029);
    
                t29.Text = Convert.ToString(Xs);                         //输出到文本框
                t30.Text = Convert.ToString(Ys);
                t31.Text = Convert.ToString(Zs);
                t32.Text = Convert.ToString(phi);
                t33.Text = Convert.ToString(omiga);
                t34.Text = Convert.ToString(kappa);
            }   
    
            private void button5_Click(object sender, EventArgs e)     //清空待定地面点像点坐标
            {
                t35.Text = "";
                t36.Text = "";
                t37.Text = "";
                t38.Text = "";
                t39.Text = "";
                t40.Text = "";
                t41.Text = "";
    
            }
    
            private void button4_Click(object sender, EventArgs e)    //计算前方交会
            {
                double Xs1, Ys1, Zs1, phi_1, omiga_1, kappa_1, Xs2, Ys2, Zs2, phi_2, omiga_2, kappa_2;
                double N1, N2, X, Y, Z;
                double Bu, Bv, Bw;
                Xs1 = Convert.ToDouble(t23.Text);               //十二个外方位元素
                Ys1 = Convert.ToDouble(t24.Text);
                Zs1 = Convert.ToDouble(t25.Text);
                phi_1 = Convert.ToDouble(t26.Text);
                omiga_1 = Convert.ToDouble(t27.Text);
                kappa_1 = Convert.ToDouble(t28.Text);
                Xs2 = Convert.ToDouble(t29.Text);
                Ys2 = Convert.ToDouble(t30.Text);
                Zs2 = Convert.ToDouble(t31.Text);
                phi_2 = Convert.ToDouble(t32.Text);
                omiga_2 = Convert.ToDouble(t33.Text);
                kappa_2 = Convert.ToDouble(t34.Text);
                Matrix tt1 = new Matrix(3, 1, "tt1");
                Matrix tt2 = new Matrix(3, 1, "tt2");
                double[,] a1 = tt1.Detail;
                double[,] a2 = tt2.Detail;
                int n = tt1.getN;
                int m = tt1.getM;
    
                a1[0, 0] = Convert.ToDouble(t35.Text);            //地面点A相应的像点a1、a2的像空间坐标系
                a1[1, 0] = Convert.ToDouble(t36.Text);
                a1[2, 0] = -Convert.ToDouble(t22.Text);
                a2[0, 0] = Convert.ToDouble(t37.Text);
                a2[1, 0] = Convert.ToDouble(t38.Text);
    
                Matrix rr1 = new Matrix(3, 3, "rr1");
                Matrix rr2 = new Matrix(3, 3, "rr2");
                double[,] rrr1 = rr1.Detail;
                double[,] rrr2 = rr2.Detail;
    
                /////计算左片旋转矩阵
                rrr1[0, 0] = Math.Cos(phi_1) * Math.Cos(kappa_1) - Math.Sin(phi_1) * Math.Sin(omiga_1) * Math.Sin(kappa_1);//a1
                rrr1[0, 1] = -Math.Cos(phi_1) * Math.Sin(kappa_1) - Math.Sin(phi_1) * Math.Sin(omiga_1) * Math.Cos(kappa_1);//a2
                rrr1[0, 2] = -Math.Sin(phi_1) * Math.Cos(omiga_1);//a3
                rrr1[1, 0] = Math.Cos(omiga_1) * Math.Sin(kappa_1);//b1
                rrr1[1, 1] = Math.Cos(omiga_1) * Math.Cos(kappa_1);//b2
                rrr1[1, 2] = -Math.Sin(omiga_1);//b3
                rrr1[2, 0] = Math.Sin(phi_1) * Math.Cos(kappa_1) + Math.Cos(phi_1) * Math.Sin(omiga_1) * Math.Sin(kappa_1);//c1
                rrr1[2, 1] = -Math.Sin(phi_1) * Math.Sin(kappa_1) + Math.Cos(phi_1) * Math.Sin(omiga_1) * Math.Cos(kappa_1);//c2
                rrr1[2, 2] = Math.Cos(phi_1) * Math.Cos(omiga_1);//c3
    
                /////计算右片旋转矩阵
                rrr2[0, 0] = Math.Cos(phi_2) * Math.Cos(kappa_2) - Math.Sin(phi_2) * Math.Sin(omiga_2) * Math.Sin(kappa_2);//a1
                rrr2[0, 1] = -Math.Cos(phi_2) * Math.Sin(kappa_2) - Math.Sin(phi_2) * Math.Sin(omiga_2) * Math.Cos(kappa_2);//a2
                rrr2[0, 2] = -Math.Sin(phi_2) * Math.Cos(omiga_2);//a3
                rrr2[1, 0] = Math.Cos(omiga_2) * Math.Sin(kappa_2);//b1
                rrr2[1, 1] = Math.Cos(omiga_2) * Math.Cos(kappa_2);//b2
                rrr2[1, 2] = -Math.Sin(omiga_2);//b3
                rrr2[2, 0] = Math.Sin(phi_2) * Math.Cos(kappa_2) + Math.Cos(phi_2) * Math.Sin(omiga_2) * Math.Sin(kappa_2);//c1
                rrr2[2, 1] = -Math.Sin(phi_2) * Math.Sin(kappa_2) + Math.Cos(phi_2) * Math.Sin(omiga_2) * Math.Cos(kappa_2);//c2
                rrr2[2, 2] = Math.Cos(phi_2) * Math.Cos(omiga_2);//c3
    
                Matrix Ss1 = MatrixOperator.MatrixMulti(rr1, tt1);  //地面点A在像空间辅助坐标系的矩阵
                Matrix Ss2 = MatrixOperator.MatrixMulti(rr2, tt2);
                double[,] S1 = Ss1.Detail;
                double[,] S2 = Ss2.Detail;
                Bu = Xs2 - Xs1;                                     //摄影基线B的分量
                Bv = Ys2 - Ys1;
                Bw = Zs2 - Zs1;
    
                N1 = (Bu * S2[2, 0] - Bw * S2[0, 0]) / (S1[0, 0] * S2[2, 0] - S2[0, 0] * S1[2, 0]);
                N2 = (Bu * S1[2, 0] - Bw * S1[0, 0]) / (S1[0, 0] * S2[2, 0] - S2[0, 0] * S1[2, 0]);
    
                X = Xs1 + N1 * S1[0, 0];
                Y = 0.5 * ((Ys1 + N1 * S1[1, 0]) + (Ys2 + N2 * S2[1, 0]));
                Z = Zs1 + N1 * S1[2, 0];
                t39.Text = Convert.ToString(X);       //输出地面坐标
                t40.Text = Convert.ToString(Y);
                t41.Text = Convert.ToString(Z);
            }
        }
    
    }

    程序界面如下:
    界面

    程序源代码链接:后方交会_前方交会 C#

    展开全文
  • 摄影测量常用程序,包括空间后方交会前方交会,以及区域网平差程序等,界面不错
  • 在学习摄影测量课程过程中,老师要求我们采用 C#语言编写程序实现 前方交会和后方交会数据梳理,绝对定向和相对定向处理以及光束法 。 经过一段时间努力 是顺利完成程序的编写。写完以后,一直存在电脑中没有理会。...

           在学习摄影测量课程过程中,老师要求我们采用 C#语言编写程序实现 前方交会和后方交会数据梳理,绝对定向和相对定向处理以及光束法 。 经过一段时间努力 是顺利完成程序的编写。写完以后,一直存在电脑中没有理会。前段时间 电脑系统崩溃,以为程序不见了。今天,在U盘中发现了摄影测量程序,分享给大家,希望可以给大家一点帮助。

         1 后方交会与前方交会

     

     

      2 相对定向与绝对定向

     

     

     

     3光束法

     

     

     

      总结

       以上为三个程序的界面,程序计算结果与书籍对比没有明显差距。 程序仅做学习参考。

    程序源码以及数据下载: 链接:http://pan.baidu.com/s/1i5QVgOH 密码:23hl

     

    转载于:https://www.cnblogs.com/hlx-blogs/p/7760943.html

    展开全文

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 8,392
精华内容 3,356
关键字:

前方交会