摄影测量实习空间后方交会

摄影测量实习空间后方交会—.后方交会解算程序。（1）原程序：D=input('请输入地面点坐标:')Y=input('请输入影像坐标点:')/1000f=input('请输入主距:')/1000C=input('请输入像主点坐标初始值:')J=input('请输入外方位角初始值:')i=size(Y,1);Y=[Y;Y(1,:)];D=[D;D(1,:)];Xs=0;Ys=0;forj=1:im=sqrt(((D(j,1)-D(j+1,1)).^2+(D(j,2)-D(j+1,2)).^2)/((Y(j,1)-Y(j+1,1)).^2+(Y(j,2)-Y(，2)).^2))/(i-1);Xs=Xs+D(j,1)/i;Ys=Ys+D(j,2)/i;endZs=m*f;forj=1:1000R(1,1)=cos(J(1))*cos(J(3))-sin(J(1))*sin(J(2))*sin(J(3));R(1,2)=-cos(J(1))*sin(J(3))-sin(J(1))*sin(J(2))*cos(J(3));R(1,3)=-sin(J(1))*cos(J(2));R(2,1)=cos(J(2))*sin(J(3));R(2,2)=cos(J(2))*cos(J(3));R(2,3)=-sin(J(2));R(3,1)=sin(J(1))*cos(J(3))+cos(J(1))*sin(J(2))*sin(J(3));R(3,2)=-sin(J(1))*sin(J(3))+cos(J(1))*sin(J(2))*cos(J(3));R(3,3)=cos(J(1))*cos(J(2));forn=1:iXp(n)=R(1,1)*(D(n,1)-Xs)+R(2,1)*(D(n,2)-Ys)+R(3,1)*(D(n,3)-Zs);Yp(n)=R(1,2)*(D(n,1)-Xs)+R(2,2)*(D(n,2)-Ys)+R(3,2)*(D(n,3)-Zs);Zp(n)=R(1,3)*(D(n,1)-Xs)+R(2,3)*(D(n,2)-Ys)+R(3,3)*(D(n,3)-Zs);xc(n)=-f*Xp(n)/Zp(n);yc(n)=-f*Yp(n)/Zp(n);A(2*n-1,1)=(R(1,1)*f+R(1,3)*xc(n))/Zp(n);A(2*n-1,2)=(R(2,1)*f+R(2,3)*xc(n))/Zp(n);A(2*n-1,3)=(R(3,1)*f+R(3,3)*xc(n))/Zp(n);A(2*n,1)=(R(1,2)*f+R(1,3)*yc(n))/Zp(n);A(2*n,2)=(R(2,2)*f+R(2,3)*yc(n))/Zp(n);A(2*n,3)=(R(3,2)*f+R(3,3)*yc(n))/Zp(n);A(2*n-1,4)=yc(n)*sin(J(2))-(xc(n)/f*(xc(n)*cos(J(3))-yc(n)*sin(J(3)))+f*cos(J(3)))*cos(J(2));A(2*n-1,5)=-f*sin(J(3))-xc(n)*(xc(n)*sin(J(3))+yc(n)*cos(J(3)))/f;A(2*n-1,6)=yc(n);A(2*n,4)=-xc(n)*sin(J(2))-(yc(n)/f*(xc(n)*cos(J(3))-yc(n)*sin(J(3)))-f*sin(J(3)))*cos(J(2));A(2*n,5)=-f*cos(J(3))-yc(n)*(xc(n)*sin(J(3))+yc(n)*cos(J(3)))/f;A(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(A'*A)*A'*L;Xs=Xs+X(1);Ys=Ys+X(2);Zs=Zs+X(3);J(1)=J(1)+X(4);J(2)=J(2)+X(5);J(3)=J(3)+X(6);M=max(abs(X));ifM=0.001XsYsZsRbreak;endend二．运行实例（课本p-44习题24）27.已知四队点的影像坐标和地点坐标如下：影像坐标地面坐标x(mm)y(mm0X(mm)Y(mm)Z(mm)1—86.15—68.9936589.4125273.322195.172—53.4082.2137631.0831324.51728.693—14.78—76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31试计算近似垂直摄影情况下空间后方交会的解。假设内方位元素已知：f=153.24mm,x0=y0=0.（3）运行结果：请输入地面点坐标:[36589.4125273.322195.17;37631.0831324.51728.69;39100.9724934.982386.50;40426.5430319.91757.31]D=1.0e+004*3.65892.52730.21953.76313.13250.07293.91012.49350.23874.04273.03200.0757请输入影像坐标点:[-86.15-68.99;-53.4082.21;-14.78-76.63;10.4664.43]Y=-0.0862-0.0690-0.05340.0822-0.0148-0.07660.01050.0644请输入主距:[153.24]f=0.1532请输入像主点坐标初始值:[00]C=000请输入外方位角初始值:[000]J=000Xs=3.9795e+004Ys=2.7477e+004Zs=7.5728e+003R=0.99770.06750.0039-0.06750.9977-0.0021-0.00410.00181.0000二．Moravec点特征提取算法程序。（1）原程序G=input('请输入搜索区灰度矩阵：')w=input('请输入窗口大小（象元表示;）：')s=size(G);ifmax(s)wdisp('搜索区太小或窗口太大!');pause;return;endk=fix(w/2);forr=k+1:s(1)-kforc=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([v1v2v3v4]);endendL=nan*ones(s);L(k+1:s(1)-k,k+1:s(2)-k)=ivy=input('请输入经验阈值:')[r,c]=find(Ly)（2）运行实例（自由编辑一组灰度值）[455687989826345455;235467874232336544;438798442244335554;263476437583849823;256587873332435463]（3）运行结果：请输入搜索区灰度矩阵：[455687989826345455;235467874232336544;438798442244335554;263476437583849823;256587873332435463]G=455687989826345455235467874232336544438798442244335554263476437583849823256587873332435463请输入窗口大小（象元表示;）：4w=4L=Columns1through5NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN196614052481NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNColumns6through9NaNNaNNaNNaNNaNNaNNaNNaN15731090NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN请输入经验阈值:1500y=1500r=333c=356二．Marover算子的点特征提取原程序:clear;G=input('请输入搜索区灰度矩阵：')w=input('请输入窗口大小（象元表示;（如5*5，则输入5）：')s=size(G);ifmax(s)wdisp('搜索区太小或窗口太大!');pause;return;endk=fix(w/2);forr=k+1:s(1)-kforc=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([v1v2v3v4]);endendL=nan*ones(s);L(k+1:s(1)-k,k+1:s(2)-k)=ivy=input('请输入阈值:')[r,c]=find(Ly);fori=1:size(r)Z(i,:)=[r(i),c(i)];enddisp('特征点坐标表示（如点（2，3）则表示为23）如下：');Z运行实例：已知一灰度距阵（计算中的边界点可视为2）如下：22222222222102622222222222试用morver算子提取其中的特征点。运行结果：请输入搜索区灰度矩阵：[22222;22222;210262;22222;22222]G=22222222222102622222222222请输入窗口大小（象元表示;（如5*5，则输入5）：3w=3L=NaNNaNNaNNaNNaNNaN000NaNNaN128032NaNNaN000NaNNaNNaNNaNNaNNaN请输入阈值:32y=32特征点坐标表示（如点（2，3）则表示为23）如下：Z=32三Wong-Trinder圆点定位算子原程序:W=input('请输入灰度矩阵:')s=size(W);T=(min(min(W))+mean(mean(W)))/2;G=WT;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([xy]);运行实例：运用MATLAB软件自动生成一个9*8的距阵，作为目标矩阵。运行结果：;请输入灰度矩阵:fix(rand(9,8)*20)W=313401513111214867832131013816163216751515813102143731491011716417163175138151915990141718106271587617101目标重心坐标(x,y)45四移动曲面曲面原程序：clear;D=xlsread('hu.xlsx')%DEM高程数据C=input('请输入待求高程点坐标[xy]:')%待求高程坐标w=input('请输入拟合点下标(例如输入[1234567]):')%按半径大小选取所需点D1=[D(:,1)-C(1)D(:,2)-C(2)D(:,3)];%坐标重心化W1=D1(:,1).^2+D1(:,2).^2;%求距离平方大小并排序[W2,index]=sort(W1);fori=1:size(index)W3(i,1:3)=D1(index(i),1:3);%经排序后的坐标endk=0;fori=1:size(D,1)ifsize(find(w==i),2)==0%未代入的点即检验