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
2、#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);void transpose(double *m1,double *m2,int m,int n);void inverse(double *a,int n);void multi(double *mat1,doub
3、le * 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;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
4、!=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) 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
5、+=mat1i*b+k*mat2k*c+j; return;double *readdata() 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,c
6、ordata+i*5+j); printf(控制点坐标数据读取成功！n); return cordata;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) FILE *fp; char *file1=结算数据.txt; fp=fopen(file1,w); fprintf(fp,-原始坐标数据为-：n); for (int i=0;iha
7、ng;i+) for (int j=0;j5;j+) fprintf(fp,%7.4lf ,datai*5+j); fprintf(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+) fprin
8、tf(fp,%7.5lf ,faxishui*5+j); fprintf(fp,n); fprintf(fp,-n); fprintf(fp,-误差方程常数项为:-：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); fprin
9、tf(fp, fai= %lf, oumiga=%lf, kapa=%lfn,fai,oumiga,kapa); fprintf(fp,-n); fclose(fp); return;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; do
10、uble correct61=0.0,AT68=0.0,ATA66=0.0,ATL61=0.0; int row; /row用于存放坐标行数 double *controlpoint; controlpoint=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,&x
11、0); /double类型数据要用%lf printf(y0=); y0/=1000.0; scanf(%lf,&y0); printf(f=); scanf(%lf,&f); f=f/1000.0; /*-*/ /-确定未知数初始值- for(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
12、)-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)*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);
13、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,diedainumber); 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); pri
14、ntf( Ys= %lfn,Ys); printf( Zs= %lfn,Zs); printf( fai= %lfn,fai); printf( oumiga= %lfn,oumiga); printf( kapa= %lfn,kapa); savedata(row,controlpoint,A0,ATA0,L0,diedainumber,Xs,Ys,Zs,fai,oumiga,kapa); printf(-解算结束！-n); system(pause);解算结果：2. 后方交会-前方交会求解地面点坐标已知左右像片外方位元素，给出像点坐标：左像点坐标：右像点坐标：x（/m）y（/m）x（/m）
15、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 double *readdata();void savedata(int hang,double *data);double *readdata() FILE *fp; int i,j,k; int number; char datacatolog100; char leftdata300; /s
16、canf(%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+) fread(&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
17、,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+) fprintf(fp,%7.4lf ,datai*3+j); fprintf(fp,nn); fprintf(fp,-); fclose(fp); return;void main() double *imagepoint; int row;
18、 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,oumiga2=0.0178,kapa2=0.1489,Ys2=397367.234,Xs2=3445959.266,Zs2=1490.096; / printf(请输入左像片的外方位元素：n); /printf
19、(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); /printf(请输入右像片的外方位元素：n); /printf(Xs2= ); /scanf(%lf,&Xs2); /printf(Ys2= ); /scanf(%lf,&Ys
20、2); /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,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+)
