一、程序流程图
二、部分代码 两个结构体,分别用来保存外方位元素和原始数据。 - struct EOEO // elements of exterior orientation
- {
- double Xs;
- double Ys;
- double Zs;
- double omega;
- double phi;
- double kappa;
- };
复制代码- 几个主要函数
- //////////////////////////////////////////////////////////////////////////
- // 函数功能:初始化坐标数据
- // 参数说明:sd:保存原始数据的结构体数组,PhotographicScale:摄影比例尺,focus:摄影机主距
- // filename:保存坐标数据的文件名
- //////////////////////////////////////////////////////////////////////////
- void InitData(SourceData* sd, char* filename)
- {
- //打开文件以备读取数据
- cout <<"开始读取数据..."<<endl;
- ifstream datafile(filename,ios::in | ios::nocreate);
- if( !datafile)
- {
- cout<<"打开文件失败, 文件可能不存在"<<endl;
- system("pause");
- _exit(1);
- }
- memset(sd, 0, N*sizeof(SourceData));
-
- int i;
- for (i = 0; i != N; i ++)
- {
- datafile >> sd[i].x >> sd[i].y >> sd[i].X >> sd[i].Y >> sd[i].Z;
- }
- datafile.close();
- cout <<"读取数据完毕..."<<endl;
- }
- //////////////////////////////////////////////////////////////////////////
- // 函数功能:检查改正数是否已达精度要求
- // 参数说明:deta:保存改正数的数组
- //////////////////////////////////////////////////////////////////////////
- bool CheckPrecision(double* deta)
- {
- //2.9088820866572159615394846141477e-5 即 0.1′的弧度表示
- bool ret;
- ret = (fabs(deta[0])<0.000001 && fabs(deta[1])<0.000001 && fabs(deta[2])<0.000001 && \
- fabs(deta[3])<2.90888208656e-5 && fabs(deta[4])<2.90888208656e-5 && \
- fabs(deta[5])<2.90888208656e-5);
- return ret;
- }
- //////////////////////////////////////////////////////////////////////////
- // 函数功能:迭代器,计算的主体部分
- // 参数说明:sd:保存原始数据的结构体数组,PhotographicScale:摄影比例尺,focus:摄影机主距
- //////////////////////////////////////////////////////////////////////////
- void Iterator(SourceData sd[N], double PhotographicScale, double Focus)
- {
- double phi, omega, kappa, Xs, Ys, Zs;
- phi = omega = kappa =0.0;
- Xs = Ys = Zs = 0.0;
- for (int k=0;k<N;k++)
- {
- sd[k].x /= 1000.0;
- sd[k].y /= 1000.0;
- Xs += sd[k].X;
- Ys += sd[k].Y;
- }
- Xs /= N;
- Ys /= N;
-
- double f = Focus / 1000.0;
- double m = PhotographicScale;
- Zs =m*f;
-
- cout<<endl<<"六元素初始化值:"<<endl;
- cout <<"Xs: "<<Xs<<'\t'<<"Ys: "<<Ys<<'\t'<<"Zs: "<<Zs<<endl;
- cout <<"phi: "<<phi<<"\t\t"<<"omega: "<<omega<<"\t"<<"kappa: "<<kappa<<endl<<endl;
- //声明并初始化六元素改正数矩阵
- double deta[6] = {1,1,1,1,1,1};
- double x0(0);
- double y0(0);//内方位元素
- double X0[N] = {0.0};
- double Y0[N] = {0.0};
- double Z0[N] = {0.0};
- //声明旋转矩阵
- double R[9];
- double A[2*6] = {0.0};
- double AT[6*2] = {0.0};
- double Buf1[36] = {0.0};
- double Buf2[36] = {0.0};//ATA累加
- double Buf3[6] = {0.0};
- double Buf4[6] = {0.0};//ATL累加
-
- double Buf5[8*6] = {0.0};//存储8×6的A矩阵,没办法
- double Buf6[8*1] = {0.0};//存储8×1的L矩阵,同上
- double V[8*1] = {0.0};
- double ATA[36] = {0.0};
- double ATL[6] = {0.0};
- double L[2] = {0.0};
- int iCount = 0;//迭代次数
- cout << "开始迭代计算..."<<endl;
- while(!CheckPrecision(deta))
- {
- cout <<endl<<"第 "<< ++iCount<<"次迭代:"<<endl;
- if (iCount == MAXITERATION)
- {
- cout << ">>>迭代次数超限,可能不收敛"<<endl;
- break;
- }
- //每次迭代之前必须清空两个保存累加值的矩阵ATA与ATL
- for (int i = 0; i != 36; i ++)
- {
- ATA[i] = 0.0;
- if (i < 6)
- {
- ATL[i] = 0.0;
- }
- }
- //计算旋转矩阵
- R[0] = cos(phi)*cos(kappa) - sin(phi)*sin(omega)*sin(kappa);
- R[1] = -cos(phi)*sin(kappa) - sin(phi)*sin(omega)*cos(kappa);
- R[2] = -sin(phi)*cos(omega);
- R[3] = cos(omega)*sin(kappa);
- R[4] = cos(omega)*cos(kappa);
- R[5] = -sin(omega);
- R[6] = sin(phi)*cos(kappa) + cos(phi)*sin(omega)*sin(kappa);
- R[7] = -sin(phi)*sin(kappa) + cos(phi)*sin(omega)*cos(kappa);
- R[8] = cos(phi)*cos(omega);
- for (i = 0; i != N; i ++)
- {
- Z0[i] = R[2]*(sd[i].X-Xs) + R[5]*(sd[i].Y-Ys) + R[8]*(sd[i].Z-Zs);
- X0[i] = x0 - f*(R[0]*(sd[i].X-Xs) + R[3]*(sd[i].Y-Ys) + R[6]*(sd[i].Z-Zs)) / Z0[i];
- Y0[i] = y0 - f*(R[1]*(sd[i].X-Xs) + R[4]*(sd[i].Y-Ys) + R[7]*(sd[i].Z-Zs)) / Z0[i];
- A[0] = ((R[0]*f + R[2]*(sd[i].x-x0)) / Z0[i]);
- A[1] = ((R[3]*f + R[5]*(sd[i].x-x0)) / Z0[i]);
- A[2] = ((R[6]*f + R[8]*(sd[i].x-x0)) / Z0[i]);
- A[3] = ((sd[i].y-y0)*sin(omega) - ((sd[i].x-x0)*((sd[i].x-x0)*cos(kappa)\
- - (sd[i].y-y0)*sin(kappa))/f + f*cos(kappa))*cos(omega));
- A[4] = (-f*sin(kappa) - (sd[i].x-x0)*((sd[i].x-x0)*sin(kappa) + (sd[i].y-y0)*cos(kappa))/f);
- A[5] = (sd[i].y - y0);
- A[6] = ((R[1]*f + R[2]*(sd[i].y-y0)) /Z0[i]);
- A[7] = ((R[4]*f + R[5]*(sd[i].y-y0)) /Z0[i]);
- A[8] = ((R[7]*f + R[8]*(sd[i].y-y0)) /Z0[i]);
- A[9] = (-(sd[i].x-x0)*sin(omega) - ((sd[i].y-y0)*((sd[i].x-x0)*cos(kappa)\
- - (sd[i].y-y0)*sin(kappa))/f - f*sin(kappa))*cos(omega));
- A[10] = (-f*cos(kappa) - (sd[i].y-y0)*((sd[i].x-x0)*sin(kappa) + (sd[i].y-y0)*cos(kappa))/f);
- A[11] = (-(sd[i].x-x0));
- //该循环保存A矩阵,最后评定精度用
- for (int l=0;l<12;l++)
- {
- Buf5[12*i+l] = A[l];
- }
-
- //所谓的逐步法化,即要在循环内部就将ATA计算出来累加,下面的L矩阵类似
- MatrixTranspose(A,AT,2,6);
- MatrixMultiply(AT,A,Buf1,6,2,6);
- MatrixCopy(ATA, Buf2, 36);
- MatrixAdd(Buf1,Buf2,ATA,36); // 为逐步法化后的ATA矩阵累加
- L[0] = (sd[i].x - X0[i]);
- L[1] = (sd[i].y - Y0[i]);
- //保存L矩阵,最后评定精度用
- for (l=0;l<2;l++)
- {
- Buf6[2*i+l] = L[l];
- }
- MatrixMultiply(AT,L,Buf3,6,2,1);
- MatrixCopy(ATL, Buf4, 6);
- MatrixAdd(Buf3,Buf4,ATL,6);
- }//for
- //“逐步法化”的另一处不同,出循环即可直接计算ATA逆乘ATL
- MatrixInversion(ATA,6);
- MatrixMultiply(ATA,ATL,deta,6,6,1);
- //deta即为改正数
- Xs += deta[0];
- Ys += deta[1];
- Zs += deta[2];
- phi += deta[3];
- omega += deta[4];
- kappa += deta[5];
- cout<<"改正数值:"<<endl;
- for (i=0; i != 6; i ++)
- {
- cout <<"deta["<<i<<"]: "<<deta[i]<<endl;
- }
- cout << endl<<"六元素值:"<<endl;
- cout <<"Xs: "<<Xs<<endl<<"Ys: "<<Ys<<endl<<"Zs: "<<Zs<<endl;
- cout <<"kappa: "<<kappa<<endl<<"omega: "<<omega<<endl<<"phi: "<<phi<<endl<<endl;
- }//while
- EOEO eoeo;
- eoeo.kappa = kappa;
- eoeo.omega = omega;
- eoeo.phi = phi;
- eoeo.Xs = Xs;
- eoeo.Ys = Ys;
- eoeo.Zs = Zs;
- cout << ">>>正常退出迭代"<<endl<<endl<<endl;
- //精度评定
- double Q[6] = {0.0};
- for (int h=0;h<6;h++)
- {
- Q[h] = ATA[h*6+h];
- }
- MatrixMultiply(Buf5,deta,V,8,6,1);//V=Ax-L
- MatrixMinus(V,Buf6,V,8);
- double m0(0);//单位权中误差
- double VSum(0);//[vv],即平方和
- int i;
- for (i=0;i<8;i++)
- {
- VSum += V[i]*V[i];
- }
- m0=sqrt(VSum / (2*N - 6));//中误差m0
- double M[6] = {0.0};//保存六个值的中误差
- for (i = 0; i != 6; i ++)
- {
- M[i] = m0 * sqrt(Q[i]);
- if (i >= 3)
- {
- M[i] = M[i]*180*3600/PI;
- }
- }
- OutputResult(&eoeo, R, M, m0);
- cout<<endl<<"解算全部完成"<<endl<<endl;
- }
复制代码三、计算结果
四、本程序的两种拓展形式 1、 将空间后方交会计算封装成类 为了实现空间后交计算的移植性,这里我将该过程封装成一个类CResection,并定义了相应的接口。 该类接受原始数据的接口为带参构造函数,其原型为: CResection (); 输出接口为两个类成员函数: OutputResult(); //该成员适用于windows console application 程序调用时调用iostream库输出结果。 SaveResult(); //该成员适用于任何情况下将计算结果保存为文件,用户指定文件路径与文件名。 以上三个函数类型声明为public以满足从外部环境声明类对象调用的需要。另包含3个私有成员函数Iterator (), CheckPrecision(),InitData() 进行内部计算,编码与运行时对用户不可见。 详细请参看工程CResection代码。 用户调用该类时,只需将CResection.h与CResection.cpp添加进目标工程并在文件适当位置包含该类头文件即可。 使用实例: #include “CResection.h” //………………… CResection res (50000,154.23, “SourceData.txt”); res.SaveResult (“Results.txt”); res.OutputResult(); //仅在控制台程序可用 2、 将该功能封装成动态链接库 封装成动态链接库,只要将类CResection稍加修改即可,对于代码能实现可靠的隐蔽性。 调用该动态链接库,只需在工程设置中link标签页下添加CResection.lib, 并在需要调用该库的文件中包含CResection.h,而后即可声明该类对象进行空间后交计算。 详细请参看工程CResection_dll代码 上位机VC MFC程序开发精典实例大全源码与视频讲解配套下载408例 经历1年的编程与录制点击进入查看 如果您认可,可联系功能定制! 如果您着急,充值会员可直接联系发您资料!
|