Matlab+实现直接线性变换
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
直接线性变换Matlab实现的程序源代码
function re=DLT(A,B)
%imco为像方坐标，输入单位是像素
imco=A; %此处为控制点像方坐标，格式为2×n，单位：像素
%obco为物方坐标，输入单位是毫米
obco=B; %此处为控制点物方坐标，格式为n×3单位：毫米
imco_be=[];B=[];M=[];
for i=1:size(imco,2)
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