直接线性变换 Matlab 实现的程序源代码 function re=DLT(A,B) %imco为像方坐标，输入单位是像素 imco=A; %此处为控制点像方坐标，格式为2×n，单位：像素 %obco为物方坐标，输入单位是毫米 obco=B; imco_be=[];B=[];M=[]; for i=1:size(imco,2) %此处为控制点物方坐标，格式为n×3单位：毫米 imco_be=[imco_be;imco(:,i)]; end for i=1:size(imco,2) A1=[obco(i,:),1,0,0,0,0]; A2=[0,0,0,0,obco(i,:),1]; M=[M;A1;A2]; B1=obco(i,:).*imco_be(2*i-1); B2=obco(i,:).*imco_be(2*i); B=[B;B1;B2]; end M=[M,B]; N=M(1:11,:); L=N\(-imco_be(1:11,:)); X0=-((L(1)*L(9)+L(2)*L(10)+L(3)*L(11))/(L(9)*L(9)+L(10)*L(10)+L(11)*L(11))); Y0=-((L(5)*L(9)+L(6)*L(10)+L(7)*L(11))/(L(9)*L(9)+L(10)*L(10)+L(11)*L(11))); L=[L;0];M3=[];W=[]; for i=1:size(imco,2) xyz=obco(i,:); A=xyz(1)*L(9)+xyz(2)*L(10)+xyz(3)*L(11)+1; r2=(imco_be(2*i-1)-X0)*(imco_be(2*i-1)-X0)+(imco_be(2*i)-Y0)*(imco_be(2*i)-Y 0); M1=[A*(imco_be(2*i-1)-X0)*r2;A*(imco_be(2*i)-Y0)*r2]; M2=-[M(2*i-1:2*i,:),M1]/A; M3=[M3;M2]; W=[W;-[imco_be(2*i-1);imco_be(2*i)]/A]; end WP=M3'*W; NBBN=inv(M3'*M3); LP=-NBBN*WP; v=M3*LP+W; imco_be=imco_be+v; X0=-(LP(1)*LP(9)+LP(2)*LP(10)+LP(3)*LP(11))/(LP(9)*LP(9)+LP(10)*LP(10)+LP (11)*LP(11)); Y0=-(LP(5)*LP(9)+LP(6)*LP(10)+LP(7)*LP(11))/(LP(9)*LP(9)+LP(10)*LP(10)+LP (11)*LP(11)); 1 

FX=(-X0*X0+(LP(1)*LP(1)+LP(2)*LP(2)+LP(3)*LP(3))/(LP(9)*LP(9)+LP(10)*LP( 10)+LP(11)*LP(11)))^0.5; V1(1)=FX; V0(1)=v'*v; for J=1:1:8%由此控制迭代平差次数 M3=[]; W=[]; for i=1:size(imco,2) xyz=obco(i,:); A=xyz(1)*LP(9)+xyz(2)*LP(10)+xyz(3)*LP(11)+1; r2=(imco_be(2*i-1)-X0)*(imco_be(2*i-1)-X0)+(imco_be(2*i)-Y0)*(imco_be(2*i)-Y 0); M1=[A*(imco_be(2*i-1)-X0)*r2;A*(imco_be(2*i)-Y0)*r2]; M2=-[M(2*i-1:2*i,:),M1]/A; M3=[M3;M2]; W=[W;-[imco_be(2*i-1);imco_be(2*i)]/A]; end WP=M3'*W; NBBN=inv(M3'*M3); LP=-NBBN*WP; v=M3*LP+W; imco_be=imco_be+v; X0=-(LP(1)*LP(9)+LP(2)*LP(10)+LP(3)*LP(11))/(LP(9)*LP(9)+LP(10)*LP(10)+LP (11)*LP(11)); Y0=-(LP(5)*LP(9)+LP(6)*LP(10)+LP(7)*LP(11))/(LP(9)*LP(9)+LP(10)*LP(10)+LP (11)*LP(11)); FX=(-X0*X0+(LP(1)*LP(1)+LP(2)*LP(2)+LP(3)*LP(3))/(LP(9)*LP(9)+LP(10)*LP( 10)+LP(11)*LP(11)))^0.5; V1(J+1,:)=FX; V2(J,:)=V1(J+1)-V1(J); V0(J+1,:)=v'*v; %由此记录每次计算后改正数的大小 end FX=(-X0*X0+(LP(1)*LP(1)+LP(2)*LP(2)+LP(3)*LP(3))/(LP(9)*LP(9)+LP(10)*LP( 10)+LP(11)*LP(11)))^0.5; FY=(-Y0*Y0+(LP(5)*LP(5)+LP(6)*LP(6)+LP(7)*LP(7))/(LP(9)*LP(9)+LP(10)*LP( 10)+LP(11)*LP(11)))^0.5; %F为求得的焦距 format F=(FX+FY)/2; F X0 Y0 

LP