摄影测量后方交会程序(共17页).doc

《摄影测量后方交会程序(共17页).doc》由会员分享，可在线阅读，更多相关《摄影测量后方交会程序(共17页).doc（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上摄影测量后方交会程序（c/c+）输入数据截图：结果截图：程序源代码（其中的矩阵求逆在前面已经有了，）：#include #include #include const double PRECISION=1e-5;typedef double DOUBLE5;int InputData(int &Num, DOUBLE *&Data,double &m,double &f);int Resection(const int &Num,const DOUBLE *&Data,const double &m,const double &f);int InverseMatrix

2、(double *matrix,const int &row);int main(int argc, char* argv)DOUBLE *Data=NULL;int Num;double f(0),m(0);if(InputData(Num,Data,m,f)if (Data!=NULL)delete Data;return 1;if(Resection(Num,Data,m,f)if (Data!=NULL)delete Data;return 1;if (Data!=NULL)delete Data;printf(解算完毕.n);do printf(计算结果保存于结果.txt文件中n请选

3、择操作（输入P打开结果数据，R打开原始数据，其它退出程序）：);fflush(stdin); /刷新输入流char order=getchar(); if (P=order | p=order) system(结果.txt);else if (R=order | r=order)system(data.txt);elsebreak;system(cls);while(1);system(PAUSE);return 0;/*函数名：InputData *函数介绍：从文件（data.txt）中读取数据,*文件格式如下： *点数 m(未知写作0)* 内方位元素(f x0 y0)*编号 x y X Y

4、 Z*下面是一个实例：4 0153.24 0 01 -86.15 -68.99 36589.41 25273.32 2195.172 -53.40 82.21 37631.08 31324.51 728.693 -14.78 -76.63 39100.97 24934.98 2386.504 10.46 64.43 40426.54 30319.81 757.31*参数：(in/out)Num(点数)，*(in/out)Data(存放数据),m,f,x0,y0*返回值：int ，0成功，1文件打开失败，2控制点个*数不足，3文件格式错误*作者：vcrs*完成时间：09-10-4*/int In

5、putData(int &Num, DOUBLE *&Data,double &m,double &f)double x0,y0;FILE *fp_input;if (!(fp_input=fopen(data.txt,r)return 1;fscanf(fp_input,%d%lf,&Num,&m);if (Num4)return 2;fscanf(fp_input,%lf%lf%lf,&f,&x0,&y0);f/=1000;if (m0 | f0)return 3;Data=new DOUBLENum;double *temp= new doubleNum-1;double scale=0

6、;int i;for (i=0;iNum;i+)/读取数据，忽略编号if(fscanf(fp_input,%*d%lf%lf%lf%lf%lf,&Datai0,&Datai1,&Datai2,&Datai3,&Datai4)!=5)return 3;/单位换算成mDatai0/=1000.0;Datai1/=1000.0;/如果m未知则归算其值if (0=m)for (i=0;iNum-1;i+)tempi=(Datai2-Datai+12)/(Datai0-Datai+10)+(Datai3-Datai+13)/(Datai1-Datai+11);scale+=tempi/2.0;m=sca

7、le/(Num-1);fclose(fp_input);delete temp;return 0;/*函数名：MatrixMul *函数介绍：求两个矩阵的积,*参数：Jz1(第一个矩阵),row(第一个矩阵行数),*Jz2(第二个矩阵),row(第二个矩阵列数),com(第一个*矩阵列数),(out)JgJz(存放结果矩阵)*返回值：void*作者：vcrs*完成时间：09-10-4*/void MatrixMul(double *Jz1,const int &row,double *Jz2,const int &line,const int &com,double *JgJz)for (in

8、t i=0;irow;i+)for (int j=0;jline;j+)double temp=0;for (int k=0;kcom;k+)temp+=*(Jz1+i*com+k)*(*(Jz2+k*line+j);*(JgJz+i*line+j)=temp; /*函数名：OutPut *函数介绍：向结果.txt文件输出数据*参数：Q协因数阵，m精度，m0单位权中误差，6个外*方位元素，旋转矩阵*返回值：int，0成功，1失败*作者：vcrs*完成时间：09-10-4*/int OutPut(const double *&Q,const double *&m,const double &m0

9、,const double &Xs,const double &Ys,const double &Zs,const double &Phi,const double &Omega,const double &Kappa,const double *R)FILE *fp_out;if (!(fp_out=fopen(结果.txt,w)return 1;FILE *fp_input;if (!(fp_input=fopen(data.txt,r)return 1;fprintf(fp_out,*n);fprintf(fp_out,n空间后方交会程序(CC+)n遥感信息工程学院n班级：00000n学

10、号：n姓名：vcrsnn);fprintf(fp_out,*n);fprintf(fp_out,已知数据：nn已知点数：);int num;double temp,x,y;fscanf(fp_input,%d%lf,&num,&temp);fprintf(fp_out,%dn,num);fprintf(fp_out,摄影比例尺（0表示其值位置）：);fprintf(fp_out,%10.0lfn,temp);fprintf(fp_out,内方位元素(f x0 y0)：);fscanf(fp_input,%lf%lf%lf,&temp,&x,&y);fprintf(fp_out,%10lft%1

11、0lft%10lfn,temp,x,y);for (int i=0;inum;i+)double temp5;fscanf(fp_input,%*d%lf%lf%lf%lf%lf,&temp0,&temp1,&temp2,&temp3,&temp4);fprintf(fp_out,%3dt%10lft%10lft%10lft%10lft%10lfn,i+1,temp0,temp1,temp2,temp3,temp4);fclose(fp_input);fprintf(fp_out,*n);fprintf(fp_out,计算结果如下：nn外方位元素：n);fprintf(fp_out,tXs=%

12、10lfn,Xs);fprintf(fp_out,tYs=%10lfn,Ys);fprintf(fp_out,tZs=%10lfn,Zs);fprintf(fp_out,tPhi=%10lfn,Phi);fprintf(fp_out,tOmega=%10lfn,Omega);fprintf(fp_out,tKappa=%10lfnn,Kappa);fprintf(fp_out,旋转矩阵：n);for (i=0;i3;i+)fprintf(fp_out,t);for (int j=0;j3;j+)fprintf(fp_out,%10lft,*(R+i*3+j);fprintf(fp_out,n)

13、;fprintf(fp_out,n单位权中误差：%10lfnn,m0);fprintf(fp_out,协因数阵：n);for (i=0;i6;i+)fprintf(fp_out,t);for (int j=0;j6;j+)fprintf(fp_out,%20lft,*(Q+i*6+j);fprintf(fp_out,n);fprintf(fp_out,n外方位元素精度：);for (i=0;i6;i+)fprintf(fp_out,%10lft,mi);fprintf(fp_out,n);fprintf(fp_out,*n);fclose(fp_out);return 0;/*函数名：Rese

14、ction *函数介绍：计算*参数：Num(点数)，Data(数据)，m,f(焦距),x0,y0*返回值：int，0成功，其它失败*作者：vcrs*完成时间：09-10-4*/int Resection(const int &Num,const DOUBLE *&Data,const double &m,const double &f)double Xs=0,Ys=0,Zs=0;int i,j;/设置初始值for (i=0;iNum;i+)Xs+=Datai2;Ys+=Datai3;Xs/=Num;Ys/=Num;Zs=m*f;double Phi(0),Omega(0),Kappa(0);d

15、ouble R33=0.0;double *L=new double2*Num;typedef double Double66;Double6 *A=new Double62*Num;double *AT=new double2*Num*6;double *ATA=new double6*6;double *ATL=new double6;double *Xg=new double6;/迭代计算do /旋转矩阵R00=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa);R01=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(

16、Omega)*cos(Kappa);R02=-sin(Phi)*cos(Omega);R10=cos(Omega)*sin(Kappa);R11=cos(Omega)*cos(Kappa);R12=-sin(Omega);R20=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa);R21=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa);R22=cos(Phi)*cos(Omega);for (i=0;iNum;i+)double X=R00*(Datai2-Xs)+R10*(Datai3-

17、Ys)+R20*(Datai4-Zs);double Y=R01*(Datai2-Xs)+R11*(Datai3-Ys)+R21*(Datai4-Zs);double Z=R02*(Datai2-Xs)+R12*(Datai3-Ys)+R22*(Datai4-Zs);double xxx,yyy;xxx=-f*X/Z;yyy=-f*Y/Z;/常数项L2*i=Datai0-(-f*X/Z);L2*i+1=Datai1-(-f*Y/Z);A2*i0=(R00*f+R02*(xxx)/Z;A2*i1=(R10*f+R12*(xxx)/Z;A2*i2=(R20*f+R22*(xxx)/Z;A2*i3=

18、(yyy)*sin(Omega)-(xxx)/f)*(xxx)*cos(Kappa)-(yyy)*sin(Kappa)+f*cos(Kappa)*cos(Omega);A2*i4=-f*sin(Kappa)-(xxx)/f)*(xxx)*sin(Kappa)+(yyy)*cos(Kappa);A2*i5=(yyy);A2*i+10=(R01*f+R02*(yyy)/Z;A2*i+11=(R11*f+R12*(yyy)/Z;A2*i+12=(R21*f+R22*(yyy)/Z;A2*i+13=-(xxx)*sin(Omega)-(yyy)/f)*(xxx)*cos(Kappa)-(yyy)*si

19、n(Kappa)-f*sin(Kappa)*cos(Omega);A2*i+14=-f*cos(Kappa)-(yyy)/f)*(xxx)*sin(Kappa)+(yyy)*cos(Kappa);A2*i+15=-(xxx);/求矩阵A的转置矩阵ATfor (i=0;i2*Num;i+)for (j=0;j=PRECISION |fabs(Xg1)=PRECISION |fabs(Xg2)=PRECISION |fabs(Xg3)=PRECISION | fabs(Xg4)=PRECISION | (Xg5)=PRECISION);/注：协因数阵，旋转矩阵等计算本应该使用最后外方位元素值，/由

20、于变换很小忽略double *Q=ATA;double *V=new double2*Num;MatrixMul(&A00,2*Num,Xg,1,6,V);double VTV=0;for(i=0;i2*Num;i+)Vi-=Li;VTV+=Vi*Vi;double m0=sqrt(VTV/(2*Num-6);double *mm=new double6;for (i=0;i6;i+)mmi=sqrt(*(Q+i*6+i)*m0;OutPut(Q,mm,m0,Xs,Ys,Zs,Phi,Omega,Kappa,&R00);delete L;delete A;delete AT;delete AT

21、A;delete ATL;delete Xg;delete mm;delete V;return 0;void swap(double &a,double &b)double temp=a;a=b;b=temp;/*函数名：InverseMatrix *函数介绍：求矩阵的逆（高斯-约当法） *输入参数：(in/out)matrix（矩阵首地址），*（in）row（矩阵阶数）*输出参数：matrix（原矩阵的逆矩阵）*返回值：int ，0成功，1失败*调用函数：swap(double&,double&)*作者：vcrs*完成时间：09-10-4*/int InverseMatrix(double

22、 *matrix,const int &row)double *m=new doublerow*row;double *ptemp,*pt=m;int i,j;ptemp=matrix;for (i=0;irow;i+)for (j=0;jrow;j+)*pt=*ptemp;ptemp+;pt+;int k;int *is=new introw,*js=new introw;for (k=0;krow;k+)double max=0;/全选主元/寻找最大元素for (i=k;irow;i+)for (j=k;jmax)max=*(m+i*row+j);isk=i;jsk=j;if (0 = m

23、ax)return 1;/行交换if (isk!=k)for (i=0;irow;i+)swap(*(m+k*row+i),*(m+isk*row+i);/列交换if (jsk!=k)for (i=0;irow;i+)swap(*(m+i*row+k),*(m+i*row+jsk);*(m+k*row+k)=1/(*(m+k*row+k);for (j=0;jrow;j+)if (j!=k)*(m+k*row+j)*=*(m+k*row+k);for (i=0;irow;i+)if (i!=k)for (j=0;jrow;j+)if(j!=k)*(m+i*row+j)-=*(m+i*row+k)*(m+k*row+j);for (i=0;i=0;r-)if (jsr!=r)for (j=0;jrow;j+)swap(*(m+r*row+j),*(m+jsr*row+j);if (isr!=r)for (i=0;irow;i+)swap(*(m+i*row+r),*(m+i*row+isr);ptemp=matrix;pt=m;for (i=0;irow;i+)for (j=0;jrow;j+)*ptemp=*pt;ptemp+;pt+;delete is;delete js;delete m;return 0;专心-专注-专业