摄影测量实习作业—. 后方交会解算程序。

（1）原程序：A=input('请输入地面点坐标:')B=input('请输入影像坐标点:')/1000f=input('请输入主距:')/1000C=input('请输入像主点坐标初始值:')E=input('请输入外方位角初始值:')i=size(B,1);B=[B;B(1,:)];A=[A;A(1,:)];Xs=0;Ys=0;for j=1:im=sqrt(((A(j,1)-A(j+1,1)).^2+(A(j,2)-A(j+1,2)).^2)/((B(j,1)-B(j+1,1)).^2+(B(j,2)-B(，2)).^2))/(i-1);Xs=Xs+A(j,1)/i;Ys=Ys+A(j,2)/i;endZs=m*f;for j=1:1000R(1,1)=cos(E(1))*cos(E(3))-sin(E(1))*sin(E(2))*sin(E(3));R(1,2)=-cos(E(1))*sin(E(3))-sin(E(1))*sin(E(2))*cos(E(3));R(1,3)=-sin(E(1))*cos(E(2));R(2,1)=cos(E(2))*sin(E(3));R(2,2)=cos(E(2))*cos(E(3));R(2,3)=-sin(E(2));R(3,1)=sin(E(1))*cos(E(3))+cos(E(1))*sin(E(2))*sin(E(3));R(3,2)=-sin(E(1))*sin(E(3))+cos(E(1))*sin(E(2))*cos(E(3));R(3,3)=cos(E(1))*cos(E(2));for n=1:iXp(n)=R(1,1)*(A(n,1)-Xs)+R(2,1)*(A(n,2)-Ys)+R(3,1)*(A(n,3)-Zs);Yp(n)=R(1,2)*(A(n,1)-Xs)+R(2,2)*(A(n,2)-Ys)+R(3,2)*(A(n,3)-Zs);Zp(n)=R(1,3)*(A(n,1)-Xs)+R(2,3)*(A(n,2)-Ys)+R(3,3)*(A(n,3)-Zs);xc(n)=-f*Xp(n)/Zp(n);yc(n)=-f*Yp(n)/Zp(n);D(2*n-1,1)=(R(1,1)*f+R(1,3)*xc(n))/Zp(n);D(2*n-1,2)=(R(2,1)*f+R(2,3)*xc(n))/Zp(n);D(2*n-1,3)=(R(3,1)*f+R(3,3)*xc(n))/Zp(n);D(2*n,1)=(R(1,2)*f+R(1,3)*yc(n))/Zp(n);D(2*n,2)=(R(2,2)*f+R(2,3)*yc(n))/Zp(n);D(2*n,3)=(R(3,2)*f+R(3,3)*yc(n))/Zp(n);D(2*n-1,4)=yc(n)*sin(E(2))-(xc(n)/f*(xc(n)*cos(E(3))-yc(n)*sin(E(3)))+f*cos(E(3))) *cos(E(2));D(2*n-1,5)=-f*sin(E(3))-xc(n)*(xc(n)*sin(E(3))+yc(n)*cos(E(3)))/f;D(2*n-1,6)=yc(n);D(2*n,4)=-xc(n)*sin(E(2))-(yc(n)/f*(xc(n)*cos(E(3))-yc(n)*sin(E(3)))-f*sin(E(3)))*cos(E(2));D(2*n,5)= -f*cos(E(3))-yc(n)*(xc(n)*sin(E(3))+yc(n)*cos(E(3)))/f;D(2*n,6)=-xc(n);L(2*n-1,1)=Y(n,1)-xc(n);L(2*n,1)=Y(n,2)-yc(n);endX=inv(D'*D)*D'*L;Xs=Xs+X(1);Ys=Ys+X(2);Zs=Zs+X(3);E(1)=E(1)+X(4);E(2)=E(2)+X(5);E(3)=E(3)+X(6);M=max(abs(X));if M<=0.001XsYsZsRbreak;endend二．运行实例（课本p-44 习题24）27.已知四队点的影像坐标和地点坐标如下：试计算近似垂直摄影情况下空间后方交会的解。

假设内方位元素已知：f=153.24mm,x0=y0=0.（3）运行结果：请输入地面点坐标:[36589.41 25273.32 2195.17;37631.08 31324.51 728.69;39100.97 24934.98 2386.50;40426.54 30319.91 757.31]D =1.0e+004 *3.6589 2.5273 0.21953.7631 3.1325 0.07293.9101 2.4935 0.23874.0427 3.0320 0.0757请输入影像坐标点:[-86.15 -68.99;-53.40 82.21;-14.78 -76.63;10.46 64.43] Y =-0.0862 -0.0690-0.0534 0.0822-0.0148 -0.07660.0105 0.0644请输入主距:[153.24]f =0.1532请输入像主点坐标初始值:[0 0 0]C =0 0 0请输入外方位角初始值:[0 0 0]J =0 0 0Xs =3.9795e+004Ys =2.7477e+004Zs =7.5728e+003R =0.9977 0.0675 0.0039-0.0675 0.9977 -0.0021-0.0041 0.0018 1.0000二．Moravec点特征提取算法程序。

（1）原程序G=input('请输入搜索区灰度矩阵：')w=input('请输入窗口大小（象元表示;）：')s=size(G);if max(s)<wdisp('搜索区太小或窗口太大!');pause;return;endk=fix(w/2);for r=k+1:s(1)-kfor c=k+1:s(2)-kv1=sum(diff(G(r,c-k:c+k)).^2);v2=sum(diff(diag(G(r-k:r+k,c-k:c+k))).^2);v3=sum(diff(G(r-k:r+k,c)).^2);v4=sum(diff(diag(flipud(G(r-k:r+k,c-k:c+k)))).^2);iv(r-k,c-k)=min([v1 v2 v3 v4]);endendL=nan*ones(s);L(k+1:s(1)-k,k+1:s(2)-k)=ivy=input('请输入经验阈值:')[r,c]=find(L>y)（2）运行实例（自由编辑一组灰度值）[45 56 87 98 98 26 34 54 55;23 54 67 87 42 32 33 65 44;43 87 98 44 22 44 33 55 54 ;26 34 76 43 75 83 84 98 23;25 65 87 87 33 32 43 54 63]（3）运行结果：请输入搜索区灰度矩阵：[45 56 87 98 98 26 34 54 55;23 54 67 87 42 32 33 65 44;43 87 98 44 22 44 33 55 54 ;26 34 76 43 75 83 84 98 23;25 65 87 87 33 32 43 54 63]G =45 56 87 98 98 26 34 54 5523 54 67 87 42 32 33 65 4443 87 98 44 22 44 33 55 5426 34 76 43 75 83 84 98 2325 65 87 87 33 32 43 54 63请输入窗口大小（象元表示;）：4w =4L =Columns 1 through 5NaN NaN NaN NaN NaNNaN NaN NaN NaN NaNNaN NaN 1966 1405 2481NaN NaN NaN NaN NaNNaN NaN NaN NaN NaNColumns 6 through 9NaN NaN NaN NaNNaN NaN NaN NaN1573 1090 NaN NaNNaN NaN NaN NaNNaN NaN NaN NaN 请输入经验阈值:1500y =1500r =333c =356二．Marover算子的点特征提取原程序:clear;G=input('请输入搜索区灰度矩阵：')w=input('请输入窗口大小（象元表示;（如5*5，则输入5）：') s=size(G);if max(s)<wdisp('搜索区太小或窗口太大!');pause;return;endk=fix(w/2);for r=k+1:s(1)-kfor c=k+1:s(2)-kv1=sum(diff(G(r,c-k:c+k)).^2);v2=sum(diff(diag(G(r-k:r+k,c-k:c+k))).^2);v3=sum(diff(G(r-k:r+k,c)).^2);v4=sum(diff(diag(flipud(G(r-k:r+k,c-k:c+k)))).^2);iv(r-k,c-k)=min([v1 v2 v3 v4]);endendL=nan*ones(s);L(k+1:s(1)-k,k+1:s(2)-k)=ivy=input('请输入阈值:')[r,c]=find(L>y);for i=1:size(r)Z(i,:)=[r(i),c(i)];enddisp('特征点坐标表示（如点（2，3）则表示为2 3）如下：'); Z运行实例：已知一灰度距阵（计算中的边界点可视为2）如下：2 2 2 2 22 2 2 2 22 10 2 6 22 2 2 2 22 2 2 2 2试用morver算子提取其中的特征点。

运行结果：>请输入搜索区灰度矩阵：[2 2 2 2 2;2 2 2 2 2;2 10 2 6 2;2 2 2 2 2;2 2 2 2 2] G =2 2 2 2 22 2 2 2 22 10 2 6 22 2 2 2 22 2 2 2 2请输入窗口大小（象元表示;（如5*5，则输入5）：3w =3L =NaN NaN NaN NaN NaNNaN 0 0 0 NaNNaN 128 0 32 NaNNaN 0 0 0 NaNNaN NaN NaN NaN NaN请输入阈值:32y =32特征点坐标表示（如点（2，3）则表示为2 3）如下：Z =3 2三Wong-Trinder圆点定位算子原程序:W=input('请输入灰度矩阵:')s=size(W);T=(min(min(W))+mean(mean(W)))/2;G=W>T;M=sum(sum(G.*W));M=sum(sum(G.*W));I=repmat(1:s(2),s(1),1);J=repmat((1:s(1))',1,s(2));x=round(sum(sum(I.*G.*W))/M);y=round(sum(sum(J.*G.*W))/M);disp('目标重心坐标(x,y)');disp([x y]);运行实例：运用MATLAB软件自动生成一个9*8的距阵，作为目标矩阵。