相对定向—绝对定向解法实验报告
1、
实验代码
1.1根据所给同名像点的像平面坐标进行相对定向，求解相对的相对定向元素 function xP= xiangduidingxiang( Lxy,Rxy,f )
%UNTITLED Summary of this function goes here % Detailed explanation goes here %设置相对定向元素是初始值 u=0; v=0; w=0; q=0; k=0;
bu=Rxy(1,1)-Lxy(1,1); while (1)
%求解余弦元素
a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k); a2=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k); a3=-sin(q)*cos(w); b1=cos(w)*sin(k); b2=cos(w)*cos(k); b3=-sin(w);
c1=sin(q)*cos(k)+cos(q)*sin(w)*sin(k); c2=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k); c3=cos(q)*cos(w);
R=[a1,a2,a3;b1,b2,b3;c1,c2,c3];
[n,m]=size(Lxy); u2=[];v2=[];w2=[]; for i=1:n
u2(i)=a1*Rxy(i,1)+a2*Rxy(i,2)-a3*f; v2(i)=b1*Rxy(i,1)+b2*Rxy(i,2)-b3*f; w2(i)=c1*Rxy(i,1)+c2*Rxy(i,2)-c3*f; end
for i=1:n
u1(i)=Lxy(i,1); v1(i)=Lxy(i,2);
w1(i)=-f; end
bv=bu*u; bw=bu*v;
for i=1:n
N1(i)=(bu*w2(i)-bw*u2(i))/(u1(i)*w2(i)-u2(i)*w1(i)); N2(i)=(bu*w1(i)-bw*u1(i))/(u1(i)*w2(i)-u2(i)*w1(i)); end
for i=1:n
a(i)=-u2(i)*v2(i)*N2(i)/w2(i);
b(i)=-(w2(i)+v2(i)*v2(i)/w2(i))*N2(i); c(i)=u2(i)*N2(i); d(i)=bu;
e(i)=-v2(i)*bu/w2(i);
l(i)=N1(i)*v1(i)-N2(i)*v2(i)-bv; end
%组成法方程系数阵A
A=zeros(n,5); %c个控制点，A：c行，5列 for i=1:n
A(i,1)=a(i); A(i,2)=b(i); A(i,3)=c(i); A(i,4)=d(i); A(i,5)=e(i); L(i,1)=l(i); end
%求解改正数X
X=inv((A')*A)*(A')*L;
q=q+X(1,1);w=w+X(2,1);k=k+X(3,1);u=u+X(4,1);v=v+X(5,1);
%求解改正数绝对值的最大项，判断最大项是否小于限差 Xabs=abs(X); aaa=max(Xabs);
if aaa<0.00003 %当改正数中绝对值最大的改正数小于限差0.00003 break; %后跳出循环，计算结果已经收敛 end
xP=[u,v,q,w,k ]; end
1.2根据所给控制点的像平面坐标，求解控制点的模型坐标 function M= calmodelcord(xP,Lxy,Rxy,f,m) h(1,:)=[0,0,0,0,0,0];
bu=(Lxy(1,1)-Rxy(1,1))*m; bv=bu*xP(1); bw=bu*xP(2);
h(2,:)=[bu,bv,bw,xP(3),xP(4),xP(5)];
M= qianfang(h,Lxy,Rxy,f); end
1.3利用控制点的地面摄影测量坐标和模型坐标求解相对立体模型的绝对定向元素 function [jP,Accuracy]= jueduidingxiang(M,G)
%UNTITLED Summary of this function goes here % Detailed explanation goes here %设置绝对定向元素是初始值 Xs=0; Ys=0; Zs=0; q=0; w=0; k=0; r=1;
[n,m]=size(G);
gt=sum(G)/n; gm=sum(M)/n;
for i=1:n
%Mg(i,:)=M(i,:)-gm; %Gg(i,:)=G(i,:)-gt; Mg(i,:)=M(i,:) ; Gg(i,:)=G(i,:) ; end
%组成法方程系数阵A
%A=zeros(3*n,4); %c个控制点，A：2c行，6列 A=zeros(3*n,7); for i=1:n
A(3*i-2,:)=[1,0,0,Mg(i,1),-Mg(i,3),0,-Mg(i,2)]; A(3*i-1,:)=[0,1,0,Mg(i,2),0,-Mg(i,3),Mg(i,1)]; A(3*i-0,:)=[0,0,1,Mg(i,3),Mg(i,1),Mg(i,2),0]; end
while(1) %求解余弦元素
a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k); a2=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k); a3=-sin(q)*cos(w); b1=cos(w)*sin(k); b2=cos(w)*cos(k); b3=-sin(w);
c1=sin(q)*cos(k)+cos(q)*sin(w)*sin(k); c2=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k); c3=cos(q)*cos(w);
R=[a1,a2,a3;b1,b2,b3;c1,c2,c3];
L=zeros(3*n,1); for i=1:n
l(i,:)=Gg(i,:)-r*Mg(i,:)*R'-[Xs,Ys,Zs]; L(3*i-2)=l(i,1); L(3*i-1)=l(i,2); L(3*i-0)=l(i,3); end
%求解改正数X
X=inv((A')*A)*(A')*L;
q=q+X(5,1);w=w+X(6,1);k=k+X(7,1);r=r+X(4,1);
Xs=Xs+X(1,1);Ys=Ys+X(2,1);Zs=Zs+X(3,1) ;
%q=q+X(1,1);w=w+X(2,1);k=k+X(3,1);r=r+X(4,1); %Xs=Xs+X(1,1);Ys=Ys+X(2,1);Zs=Zs+X(3,1) ;
%求解改正数绝对值的最大项，判断最大项是否小于限差 Xabs=abs(X);
%X2=X(1:3); X2=X(1:7); aaa=max(X2);
if aaa<0.00003 %当改正数中绝对值最大的改正数小于限差0.00003 break; %后跳出循环，计算结果已经收敛 end
V=A*X-L;
Qx=inv((A')*A);
m=sqrt(V'*V/(3*n-7)); mx=m*sqrt(Qx(1,1)); my=m*sqrt(Qx(2,2)); mz=m*sqrt(Qx(3,3)); mr=m*sqrt(Qx(4,4)); mq=m*sqrt(Qx(5,5)); mw=m*sqrt(Qx(5,5)); mk=m*sqrt(Qx(7,7));
Accuracy=[m,mx,my,mz,mr,mq,mw,mk]; jP=[Xs,Ys,Zs,q,w,k,r]; end
1.4根据同名像点在左右像片上的坐标，运用相对定向-绝对定向求解其对应的地面点在摄影测量坐标系中的坐标
function G= modeltoground(M,jP)
%UNTITLED Summary of this function goes here % Detailed exjPlanation goes here %设置绝对定向元素是初始值 Xs=jP(1); Ys=jP(2); Zs=jP(3); q=jP(4); w=jP(5); k=jP(6); r=jP(7);
%求解余弦元素
a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);