《摄影测量程序汇总(后方交会前方交会单模型光束法平差).doc》由会员分享，可在线阅读，更多相关《摄影测量程序汇总(后方交会前方交会单模型光束法平差).doc（21页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、-/程序运行环境为Visual Studio2010.运行前请先将坐标数据放在debug下。1. 单像空间后方交会原始数据：内方位元素x0(/mm)y0(/mm)主距f(/mm)比例尺分母00153.2450000像点坐标(/mm)物点坐标(/m)-86.15-68.9936589.4125273.322195.17-53.482.2137631.0831324.51728.69-14.78-76.6339100.9724934.982386.510.4664.4340426.5430319.81757.31C语言程序：#include #include #include double *re
2、addata();void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa);void transpose(double *m1,double *m2,int m,int n);void inverse(double *a,int n);void multi(double *mat1,double * mat2,double * result,in
3、t a,int b,int c);void inverse(double *a,int n)/*正定矩阵求逆*/ int i,j,k;for(k=0;kn;k+)for(i=0;in;i+)if(i!=k)*(a+i*n+k)=-*(a+i*n+k)/(*(a+k*n+k);*(a+k*n+k)=1/(*(a+k*n+k);for(i=0;in;i+)if(i!=k)for(j=0;jn;j+)if(j!=k)*(a+i*n+j)+=*(a+k*n+j)* *(a+i*n+k);for(j=0;jn;j+)if(j!=k)*(a+k*n+j)*=*(a+k*n+k);void transpos
4、e(double *m1,double *m2,int m,int n) /矩阵转置 int i,j;for(i=0;im;i+)for(j=0;jn;j+)m2j*m+i=m1i*n+j;return;void multi(double *mat1,double *mat2,double * result,int a,int b,int c) int i,j,k; for(i=0;ia;i+) for(j=0;jc;j+) resulti*c+j=0; for(k=0;kb;k+) resulti*c+j+=mat1i*b+k*mat2k*c+j; return;double *readda
5、ta()FILE *fp;int i,j;int number;char datacatolog100;/scanf(%s,datacatolog);if (fp=fopen(控制点坐标.txt,r)=NULL)printf(读取数据出错!n);return false; fscanf(fp,%d,&number);double *cordata=new doublenumber*5;for (i=0;inumber;i+)for (j=0;j5;j+) fscanf(fp,%lf,cordata+i*5+j);printf(控制点坐标数据读取成功！n);return cordata;void
6、 savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa)FILE *fp;char *file1=结算数据.txt;fp=fopen(file1,w);fprintf(fp,-原始坐标数据为-：n);for (int i=0;ihang;i+)for (int j=0;j5;j+)fprintf(fp,%7.4lf ,datai*5+j);fprintf
7、(fp,n);fprintf(fp,-n);fprintf(fp,-误差方程系数阵为：-：n);for (int i=0;ihang*2;i+)for (int j=0;j6;j+)fprintf(fp,%7.4lf ,xishuarrayi*5+j);fprintf(fp,n);fprintf(fp,-n);fprintf(fp,-法方程系数阵为：-：n);for (int i=0;i6;i+)for (int j=0;j6;j+)fprintf(fp,%7.5lf ,faxishui*5+j);fprintf(fp,n);fprintf(fp,-n);fprintf(fp,-误差方程常数项
8、为:-：n);for (int i=0;ihang*2;i+)fprintf(fp,%lf ,li);fprintf(fp,n);fprintf(fp,-n);fprintf(fp,-迭代次数为:-：n);fprintf(fp,%dn,i);fprintf(fp,-n);fprintf(fp,-外方位元素为：-n);fprintf(fp, Xs= %lf, Ys=%lf, Zs=%lfn,xs,ys,zs);fprintf(fp, fai= %lf, oumiga=%lf, kapa=%lfn,fai,oumiga,kapa);fprintf(fp,-n);fclose(fp);return;
9、void main()int i,j;int ii,jj;int diedainumber=0;double x0=0.0,y0=0.0,f=0.0;double m=50000; /估算比例尺double fai=0,oumiga=0,kapa=0,Xs=0,Ys=0,Zs=0;double R33=0.0;double X=0.0,Y=0.0,Z=0.0,L81=0.0,A86=0.0; double correct61=0.0,AT68=0.0,ATA66=0.0,ATL61=0.0;int row; /row用于存放坐标行数double *controlpoint;controlpoi
10、nt=readdata();row=sizeof(controlpoint);for (i=0;irow;i+)for (j=0;j5;j+)printf(%3.3lf ,controlpointi*5+j);printf(n);/*-内方位元素-*/printf(请输入像片的内方位元素(mm)：n);printf(x0=);x0/=1000.0;scanf(%lf,&x0); /double类型数据要用%lfprintf(y0=);y0/=1000.0;scanf(%lf,&y0);printf(f=);scanf(%lf,&f);f=f/1000.0;/*-*/-确定未知数初始值-for(
11、int i=0;irow;i+)Xs=Xs+controlpointi*5+2; Ys=Ys+controlpointi*5+3; Zs=Zs+controlpointi*5+4;Xs/=row;Ys/=row;Zs=Zs/row+m*f;/- do diedainumber+;/-组成旋转矩阵-R00=cos(fai)*cos(kapa)-sin(fai)*sin(oumiga)*sin(kapa);R01=-cos(fai)*sin(kapa)-sin(fai)*sin(oumiga)*cos(kapa);R02=-sin(fai)*cos(oumiga);R10=cos(oumiga)*
12、sin(kapa);R11=cos(oumiga)*cos(kapa);R12=-sin(oumiga);R20=sin(fai)*cos(kapa)+cos(fai)*sin(oumiga)*sin(kapa);R21=-sin(fai)*sin(kapa)+cos(fai)*sin(oumiga)*cos(kapa);R22=cos(fai)*cos(oumiga);/-/计算系数阵和常数项for(int i=0,k=0,j=0;i=6.0/206265.0|correct40=6.0/206265.0|correct50=6.0/206265.0); printf(迭代次数为：%dn,d
13、iedainumber); printf(-误差方程系数为：-n);for (i=0;i8;i+)for (j=0;j6;j+)printf(%4.4lf ,Aij);printf(n);printf(-n); printf(求解得到的外方位元素为：n);printf( Xs= %lfn,Xs);printf( Ys= %lfn,Ys);printf( Zs= %lfn,Zs);printf( fai= %lfn,fai);printf( oumiga= %lfn,oumiga);printf( kapa= %lfn,kapa);savedata(row,controlpoint,A0,ATA
14、0,L0,diedainumber,Xs,Ys,Zs,fai,oumiga,kapa);printf(-解算结束！-n);system(pause);解算结果：2. 后方交会-前方交会求解地面点坐标已知左右像片外方位元素，给出像点坐标：左像点坐标：右像点坐标：x（/m）y（/m）x（/m）y（/m）0.00530.00690.004820.00270.005560.002290.00521-0.001910.008150.004990.007730.00089-0.002470.00156-0.00279-0.00292C语言代码：#include #include #include doub
15、le *readdata();void savedata(int hang,double *data);double *readdata()FILE *fp;int i,j,k;int number;char datacatolog100;char leftdata300;/scanf(%s,datacatolog);if (fp=fopen(像点坐标数据.txt,r)=NULL)printf(读取数据出错!n);system(pause);exit(0);fscanf(fp,%d,&number);double *c=new doublenumber*4;for (k=0;k2;k+)fre
16、ad(&leftdata,14,1,fp);for (i=0;inumber;i+) for (j=0;j2;j+)fscanf(fp,%lf,c+k*2+i*4+j);fclose(fp);return c;void savedata(int hang,double *data)FILE *fp;char *file1=地面点坐标数据.txt;fp=fopen(file1,w);fprintf(fp,-像点对应地面点坐标为-：n);fprintf(fp,n);for (int i=0;ihang;i+)fprintf(fp,第%d点： ,i+1);for (int j=0;j3;j+)fpr
17、intf(fp,%7.4lf ,datai*3+j);fprintf(fp,nn);fprintf(fp,-);fclose(fp);return;void main()double *imagepoint;int row;int i,j;imagepoint=readdata();row=sizeof(imagepoint);/-double f=24;f/=1000;double fai1=-0.0061,oumiga1=0.0327,kapa1=0.1711,Ys1=397367.171,Xs1=3445820.098,Zs1=1486.212;double fai2=0.0063,ou
18、miga2=0.0178,kapa2=0.1489,Ys2=397367.234,Xs2=3445959.266,Zs2=1490.096;/ printf(请输入左像片的外方位元素：n);/printf(Xs1= );/scanf(%lf,&Xs1);/printf(Ys1= );/scanf(%lf,&Ys1);/printf(Zs1= );/scanf(%lf,&Zs1);/printf(fai1= );/scanf(%lf,&fai1);/printf(oumiga1= );/scanf(%lf,&oumiga1);/printf(kapa1= );/scanf(%lf,&kapa1)
19、;/printf(请输入右像片的外方位元素：n);/printf(Xs2= );/scanf(%lf,&Xs2);/printf(Ys2= );/scanf(%lf,&Ys2);/printf(Zs2= );/scanf(%lf,&Zs2);/printf(fai2= );/scanf(%lf,&fai2);/printf(oumiga2= );/scanf(%lf,&oumiga2);/printf(kapa2= );/scanf(%lf,&kapa2);double Bx=Xs2-Xs1,By=Ys2-Ys1,Bz=Zs2-Zs1;double N1=0,N2=0;double X1=0,
20、Y1=0,Z1=0,X2=0,Y2=0,Z2=0;double R133=0.0;double R233=0.0;double GEOdata43=0.0;for (i=0;irow;i+)/-组成左影像旋转矩阵-R100=cos(fai1)*cos(kapa1)-sin(fai1)*sin(oumiga1)*sin(kapa1);R101=-cos(fai1)*sin(kapa1)-sin(fai1)*sin(oumiga1)*cos(kapa1);R102=-sin(fai1)*cos(oumiga1);R110=cos(oumiga1)*sin(kapa1);R111=cos(oumig
21、a1)*cos(kapa1);R112=-sin(oumiga1);R120=sin(fai1)*cos(kapa1)+cos(fai1)*sin(oumiga1)*sin(kapa1);R121=-sin(fai1)*sin(kapa1)+cos(fai1)*sin(oumiga1)*cos(kapa1);R122=cos(fai1)*cos(oumiga1);/-/-组成右影像旋转矩阵-R200=cos(fai2)*cos(kapa2)-sin(fai2)*sin(oumiga2)*sin(kapa2);R201=-cos(fai2)*sin(kapa2)-sin(fai2)*sin(ou
22、miga2)*cos(kapa2);R202=-sin(fai2)*cos(oumiga2);R210=cos(oumiga2)*sin(kapa2);R211=cos(oumiga2)*cos(kapa2);R212=-sin(oumiga2);R220=sin(fai2)*cos(kapa2)+cos(fai2)*sin(oumiga2)*sin(kapa2);R221=-sin(fai2)*sin(kapa2)+cos(fai2)*sin(oumiga2)*cos(kapa2);R222=cos(fai2)*cos(oumiga2);/-像空辅系坐标-X1=R100*imagepoint
23、i*4+0+R101*imagepointi*4+1-R102*f;Y1=R110*imagepointi*4+0+R111*imagepointi*4+1-R112*f;Z1=R120*imagepointi*4+0+R121*imagepointi*4+1-R122*f;X2=R200*imagepointi*4+2+R201*imagepointi*4+3-R202*f;Y2=R210*imagepointi*4+2+R211*imagepointi*4+3-R212*f;Z2=R220*imagepointi*4+2+R221*imagepointi*4+3-R222*f;/-/-点投
24、影系数-N1=(Bx*Z2-Bz*X2)/(X1*Z2-Z1*X2);N2=(Bx*Z1-Bz*X1)/(X1*Z2-Z1*X2);/-/-计算地面点坐标-GEOdatai0=Xs1+N1*X1;GEOdatai1=Ys1+By+N2*Y2;GEOdatai2=Zs1+N1*Z1;/-/-for (i=0;i4;i+)printf(第%d个地面点坐标： ,i+1);for (j=0;j3;j+)printf(%lf ,GEOdataij);printf(nn);savedata(row,GEOdata0);system(pause);测试结果：3. 单模型光束法严密平差缺少已知数据进行验证，因
25、此如果有已知数据请代入已知数据进行验证。C语言代码：#include #include double *readdata();void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa,double xs2,double ys2,double zs2,double fai2,double oumiga2,double kapa2);void tran
26、spose(double *m1,double *m2,int m,int n);void inverse(double *a,int n);void multi(double *mat1,double * mat2,double * result,int a,int b,int c);void inverse(double *a,int n)/*正定矩阵求逆*/ int i,j,k;for(k=0;kn;k+)for(i=0;in;i+)if(i!=k)*(a+i*n+k)=-*(a+i*n+k)/(*(a+k*n+k);*(a+k*n+k)=1/(*(a+k*n+k);for(i=0;in
27、;i+)if(i!=k)for(j=0;jn;j+)if(j!=k)*(a+i*n+j)+=*(a+k*n+j)* *(a+i*n+k);for(j=0;jn;j+)if(j!=k)*(a+k*n+j)*=*(a+k*n+k);void transpose(double *m1,double *m2,int m,int n) /矩阵转置 int i,j;for(i=0;im;i+)for(j=0;jn;j+)m2j*m+i=m1i*n+j;return;void multi(double *mat1,double *mat2,double * result,int a,int b,int c)
28、 int i,j,k;for(i=0;ia;i+)for(j=0;jc;j+)resulti*c+j=0;for(k=0;kb;k+)resulti*c+j+=mat1i*b+k*mat2k*c+j;return;double *readdata()FILE *fp;int i,j,k;int number;char datacatolog100;char leftdata300;/scanf(%s,datacatolog);if (fp=fopen(控制点坐标.txt,r)=NULL)printf(读取数据出错!n);return false;fscanf(fp,%d,&number);do
29、uble *c=new doublenumber*7;for (k=0;k3;k+)fread(&leftdata,14,1,fp);for (i=0;i4;i+) if (k=2)for (j=0;j3;j+)fscanf(fp,%lf,c+k*2+i*7+j);elsefor (j=0;j2;j+)fscanf(fp,%lf,c+k*2+i*7+j);fclose(fp);return c;void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,doub
30、le ys,double zs,double fai,double oumiga,double kapa,double xs2,double ys2,double zs2,double fai2,double oumiga2,double kapa2)FILE *fp;char *file1=结算数据.txt;fp=fopen(file1,w);fprintf(fp,-原始坐标数据为-：n);for (int i=0;ihang;i+)for (int j=0;j7;j+)fprintf(fp,%7.4lf ,datai*7+j);fprintf(fp,n);fprintf(fp,-n);fprintf(fp,-误差方程系数阵为：-：n);for (int i=0;ihang*4;i+)for (int j=0;j12;j+)fprintf(fp,%4.3e ,xishuarrayi*12+j);fprintf(fp,n);fprintf(fp,-n);fprintf(fp,-法方程系数阵为：-：n);for (int i=0;i12;i+)for (int j=0;j12;j+)fprintf(fp,%4.3e ,faxishui*12+j);fprintf(fp,n);fprint