内容发布更新时间 : 2024/5/17 9:57:55星期一 下面是文章的全部内容请认真阅读。

单像空间后方交会程序

西南交通大学 土木工程学院 测绘工程 张慧鑫 学号：20030593 输入文件形式如下：

C++源程序如下：

#include #include #include #include #include using namespace std; const int n=6;

void inverse (double c[n][n]);

templatevoid transpose (T1*mat1,T2*mat2,int a,int b); templatevoid multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c);

templatevoid input (T*mat,int a,int b);

templatevoid output(T*mat,char*s,int a,int b); int main() {

ofstream outFile; cout.precision(5);

double x0=0.0, y0=0.0; double fk=0.15324; //内方位元素 double m=39689; //估算比例尺

double B[4][5]={0.0},R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];

input (B,4,5); //从文件中读取控制点的影像坐标和地面坐标，存入数组B double Xs=0.0, Ys=0.0, Zs=0.0,Q=0.0,W=0.0,K=0.0; double X,Y,Z,L[8][1],A[8][6]; //确定未知数的出始值

for(int i=0;i<4;i++) {Xs=Xs+B[i][2]; Ys=Ys+B[i][3]; Zs=Zs+B[i][4]; }

Xs=Xs/4; Ys=Ys/4; Zs=Zs/4+m*fk; int f=0;

do//迭代计算 {f++;

//组成旋转矩阵

R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K); R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K); R[0][2]=-sin(Q)*cos(W); R[1][0]=cos(W)*sin(K); R[1][1]=cos(W)*cos(K); R[1][2]=-sin(W);

R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K); R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K); R[2][2]=cos(Q)*cos(W); //计算系数阵和常数项

for(int i=0,k=0,j=0;i<=3;i++,k++,j++) { X=R[0][0]*(B[i][2]-Xs)+R[1][0]*(B[i][3]-Ys)+R[2][0]*(B[i][4]-Zs); Y=R[0][1]*(B[i][2]-Xs)+R[1][1]*(B[i][3]-Ys)+R[2][1]*(B[i][4]-Zs); Z=R[0][2]*(B[i][2]-Xs)+R[1][2]*(B[i][3]-Ys)+R[2][2]*(B[i][4]-Zs); L[j][0]=B[i][0]-(x0-fk*X/Z); L[j+1][0]=B[i][1]-(y0-fk*Y/Z); j++; A[k][0]=(R[0][0]*fk+R[0][2]*(B[i][0]-x0))/Z; A[k][1]=(R[1][0]*fk+R[1][2]*(B[i][0]-x0))/Z; A[k][2]=(R[2][0]*fk+R[2][2]*(B[i][0]-x0))/Z;

A[k][3]=(B[i][1]-y0)*sin(W)-((B[i][0]-x0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk+fk*cos(K))*cos(W);

A[k][4]=-fk*sin(K)-(B[i][0]-x0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk; A[k][5]=B[i][1]-y0; A[k+1][0]=(R[0][1]*fk+R[0][2]*(B[i][1]-y0))/Z; A[k+1][1]=(R[1][1]*fk+R[1][2]*(B[i][1]-y0))/Z; A[k+1][2]=(R[2][1]*fk+R[2][2]*(B[i][1]-y0))/Z;

A[k+1][3]=-(B[i][0]-x0)*sin(W)-((B[i][1]-y0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk-fk*sin(K))*cos(W);

A[k+1][4]=-fk*cos(K)-(B[i][1]-y0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk; A[k+1][5]=-(B[i][0]-x0);

k++; }

transpose(A,AT,6,8); multi(AT,A,ATA,6,8,6); inverse(ATA);

multi(AT,L,ATL,6,8,1); multi(ATA,ATL,XG,6,6,1);

Xs=Xs+XG[0][0]; Ys=Ys+XG[1][0]; Zs=Zs+XG[2][0]; Q=Q+XG[3][0]; W=W+XG[4][0]; K=K+XG[5][0];

}while(XG[3][0]>=6.0/206265.0||XG[4][0]>=6.0/206265.0||XG[5][0]>=6.0/206265.0);

cout<<\迭代次数为：\ //精度评定

double AXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6]; multi(A,XG,AXG,8,6,1);

for( i=0;i<8;i++) //计算改正数 V[i][0]=AXG[i][0]-L[i][0]; transpose (V,VT,1,8); multi(VT,V,VTV,1,8,1); m0=VTV[0][0]/2; for(i=0;i<6;i++) for(int j=0;j<6;j++) D[i][j]=m0*ATA[i][j];

//屏幕输出误差方程系数阵、常数项、改正数 output(A,\误差方程系数阵A为：\ output(L,\常数项L为：\ output(XG,\改正数为：\

outFile.open(\ //打开并添加aim.txt文件 outFile.precision(10);

//以文件的形式输出像片外方位元素、旋转矩阵、方差阵 outFile<<\一、像片的外方位元素为：\

outFile<

<

outFile<

\旁向倾角为：\像片旋角为：\

outFile<

outFile<

outFile<

for(i=0;i<6;i++) {for(int j=0;j<6;j++) outFile<

outFile.close(); return 0; }

templatevoid transpose(T1*mat1,T2*mat2,int a,int b) { int i,j;

for(i=0;i

templatevoid multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) { int i,j,k;

for(i=0;i

return; }

template void input (T*mat,int a,int b) { ifstream inFile;

inFile.open(\控制点坐标.txt\ while(!inFile.eof()) {for (int i=0;i>mat[i][j]; }

inFile.close(); return; }

templatevoid output(T*mat,char*s,int a,int b) { cout<

cout<

cout<

return; }

void inverse(double c[n][n]) { int i,j,h,k; double p;

double q[n][12];

for(i=0;i

q[i][j]=0;}

for(h=k=0;k

p=q[k][h]/q[i][h]; for(j=0;j<12;j++) { q[i][j]*=p;

q[i][j]-=q[k][j]; } }

for(h=k=n-1;k>0;k--,h--) // 消去对角线以上的数据 for(i=k-1;i>=0;i--) {if(q[i][h]==0) continue;

p=q[k][h]/q[i][h]; for(j=0;j<12;j++) {q[i][j]*=p;

q[i][j]-=q[k][j];}}

for(i=0;i

for(i=0;i

程序的结果输出如下：（包括文本输出结果和荧屏输出中间数据）