www.pudn.com > fc.rar > ed3333.m, change:2015-05-18,size:19941b

```close all;
clear all;
clc
figure(1);plot(mtr(:,:,1),mtr(:,:,2),'b*');
gd=22;        % define thr grid of the dots in the corrected image
for ix=1:19
for iy=1:13
mtr1(ix,iy,1)=(ix-1)*gd+15; % x-coordinate
mtr1(ix,iy,2)=(iy-1)*gd+15; % y-coordinate
end
end
figure(2);plot(mtr1(:,:,1),mtr1(:,:,2),'r+');
figure(3);imshow(I);
[m,n]=size(I);
J=zeros(m,n);
M=[];
for i2=1:1:m
for j2=1:1:n
J(i2,j2)=255;
end
end
for i1=1:1:339
for j1=1:1:444 %(i1 j1)为要进行矫正的点
for i=1:1:18
for j=1:1:12
k11=(mtr(i,j,2)-mtr(i+1,j,2))/(mtr(i,j,1)-mtr(i+1,j,1));%指定方框上下两边的直线方程，其中k11,b11为下面的直线
b11=(mtr(i,j,1)*mtr(i+1,j,2)-mtr(i+1,j,1)*mtr(i,j,2))/(mtr(i,j,1)-mtr(i+1,j,1));
k22=(mtr(i,j+1,2)-mtr(i+1,j+1,2))/(mtr(i,j+1,1)-mtr(i+1,j+1,1));
b22=(mtr(i,j+1,1)*mtr(i+1,j+1,2)-mtr(i+1,j+1,1)*mtr(i,j+1,2))/(mtr(i,j+1,1)-mtr(i+1,j+1,1));
if ((abs(mtr(i,j,1)-j1)<0.01)&&(abs(mtr(i,j,2)-i1)<0.01))
K=0;
KK=0;
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);
break;
end
if ((abs(mtr(i+1,j,1)-j1)<0.01)&&(abs(mtr(i+1,j,2)-i1)<0.01))
K=0;
KK=1;
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);
break;
end
if ((abs(mtr(i,j+1,1)-j1)<0.01)&&(abs(mtr(i,j+1,2)-i1)<0.01))
K=1;
KK=0;
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);
break;
end
if ((abs(mtr(i+1,j+1,1)-j1)<0.01)&&(abs(mtr(i+1,j+1,2)-i1)<0.01))
K=1;
KK=1;
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);
break;
end
if (abs(mtr(i,j,1)-mtr(i,j+1,1))<0.01)&&(abs(mtr(i+1,j,1)-mtr(i+1,j+1,1))<0.01) %k1 k2为左右两条直线的斜率，此时k1 k2为无穷大，即左右两直线都垂直于x轴
x1=mtr(i,j,1);
x2=mtr(i+1,j,1);%由于左右两直线都垂直于x轴，所以过(j1,i1)与两直线的交点横坐标便知
if ((k11*j1+b11)<=i1)&&((k22*j1+b22)>=i1) %将已知点的横坐标分别带入上、下边方程，如果下边方程得到的纵坐标小于已知点纵坐标i1，说明该点在下方直线的上侧，同理若
%上边方程的值大于已知点纵坐标i1，说明该点在上方直线的下侧，若上述两者同时满足且结合上面的条件语句可知找到了该点所在的方框，并将该方框的左下角点坐标记录(i,j)
%figure(1);hold on;plot(mtr(i,j,1),mtr(i,j,2),'r+');%将找到的已知点所在的方框的左下角坐标画出来标记
d=1;
h=[];
for k0=-10:0.01:10
b0=i1-j1*k0;
y1=k0*x1+b0;
if  double(y1)<mtr(i,j,2)||double(y1)>mtr(i,j+1,2)
continue;
end
y2=k0*x2+b0;%(x1,y1)(x2,y2)为要求的直线与左右两条直线的交点
if  double(y2)<mtr(i+1,j,2)||double(y2)>mtr(i+1,j+1,2)
continue;
end
%%%求距离
a=sqrt((y1-mtr(i,j,2))^2);
b=sqrt((y2-mtr(i+1,j,2))^2);
A=sqrt((mtr(i,j,1)-mtr(i,j+1,1))^2+(mtr(i,j,2)-mtr(i,j+1,2))^2);%求已知两点之间的距离
B=sqrt((mtr(i+1,j,1)-mtr(i+1,j+1,1))^2+(mtr(i+1,j,2)-mtr(i+1,j+1,2))^2);
hh=a/A-b/B;
h=[h;k0,hh];
if sign(hh)~=sign(d)
break;
end
d=hh;
end
K=a/A;
%%%用来求已知点与上下两边的交点(x11,y11)(x22,y22)
dd=-1;
kkk=[];
for k00=-50:0.1:50
b00=i1-j1*k00;
x11=((b00-b11)/(k11-k00));
if  double(x11)<mtr(i,j,1)||double(x11)>mtr(i+1,j,1)
continue;
end
y11=(k11*b00-k00*b11)/(k11-k00);
x22=(b00-b22)/(k22-k00);
y22=(k22*b00-k00*b22)/(k22-k00);%(x11,y11)(x22,y22)为要求的直线与已知直线的下、上两个交点
%%%求距离
aa=sqrt((x11-mtr(i,j,1))^2+(y11-mtr(i,j,2))^2);
bb=sqrt((x22-mtr(i,j+1,1))^2+(y22-mtr(i,j+1,2))^2);
AA=sqrt((mtr(i,j,1)-mtr(i+1,j,1))^2+(mtr(i,j,2)-mtr(i+1,j,2))^2);%求已知两点之间的距离
BB=sqrt((mtr(i,j+1,1)-mtr(i+1,j+1,1))^2+(mtr(i,j+1,2)-mtr(i+1,j+1,2))^2);
kk=aa/AA-bb/BB;
kkk=[kkk;k00,kk];
%用来判断k发生符号变化
if sign(aa/AA-bb/BB)~=sign(dd)
break;
end
dd=kk;
end
KK=aa/AA;
%%求左右两条边的成等比的比例值
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);
break;
end
else if (abs(mtr(i,j,1)-mtr(i,j+1,1))<=0.01)&&(abs(mtr(i+1,j,1)-mtr(i+1,j+1,1))>0.01)%第一个点和第三个点横坐标相等且第二个和第四个横坐标不相等
x1=mtr(i,j,1);
if ((k11*j1+b11)<=i1)&&((k22*j1+b22)>=i1) %将已知点的横坐标分别带入上、下边方程，如果下边方程得到的纵坐标小于已知点纵坐标j1，说明该点在下方直线的上侧，
d=1;
h=[];
for k0=-10:0.01:10
b0=i1-j1*k0;
%%%求交点坐标(x1,y1)(x2,y2)
y1=k0*x1+b0;
if  double(y1)<mtr(i,j,2)||double(y1)>mtr(i,j+1,2)
continue;
end
x2=(b0-b2)/(k2-k0);
y2=(k2*b0-k0*b2)/(k2-k0);%(x1,y1)(x2,y2)为要求的直线与左右两条直线的交点
if  double(y2)<mtr(i+1,j,2)||double(y2)>mtr(i+1,j+1,2)
continue;
end
%%%求距离
a=sqrt((y1-mtr(i,j,2))^2);
b=sqrt((x2-mtr(i+1,j,1))^2+(y2-mtr(i+1,j,2))^2);
A=sqrt((mtr(i,j,1)-mtr(i,j+1,1))^2+(mtr(i,j,2)-mtr(i,j+1,2))^2);%求已知两点之间的距离
B=sqrt((mtr(i+1,j,1)-mtr(i+1,j+1,1))^2+(mtr(i+1,j,2)-mtr(i+1,j+1,2))^2);
hh=a/A-b/B;
h=[h;k0,hh];
if sign(hh)~=sign(d)
break;
end
d=hh;
end
K=a/A;
%%%用来求已知点与上下两边的交点(x11,y11)(x22,y22)
dd=-1;
kkk=[];
for k00=-50:0.1:50
b00=i1-j1*k00;
x11=((b00-b11)/(k11-k00));
if  double(x11)<mtr(i,j,1)||double(x11)>mtr(i+1,j,1)
continue;
end
y11=(k11*b00-k00*b11)/(k11-k00);
x22=(b00-b22)/(k22-k00);
y22=(k22*b00-k00*b22)/(k22-k00);%(x11,y11)(x22,y22)为要求的直线与已知直线的下、上两个交点
%%%求距离
aa=sqrt((x11-mtr(i,j,1))^2+(y11-mtr(i,j,2))^2);
bb=sqrt((x22-mtr(i,j+1,1))^2+(y22-mtr(i,j+1,2))^2);
AA=sqrt((mtr(i,j,1)-mtr(i+1,j,1))^2+(mtr(i,j,2)-mtr(i+1,j,2))^2);%求已知两点之间的距离
BB=sqrt((mtr(i,j+1,1)-mtr(i+1,j+1,1))^2+(mtr(i,j+1,2)-mtr(i+1,j+1,2))^2);
kk=aa/AA-bb/BB;
kkk=[kkk;k00,kk];
if sign(aa/AA-bb/BB)~=sign(dd)
break;
end
dd=kk;
end
KK=aa/AA;
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);
break;
end
else if (abs((mtr(i,j,1)-mtr(i,j+1,1))>0.01)&&(abs(mtr(i+1,j,1)-mtr(i+1,j+1,1))<=0.01))%第一个点横坐标不等于第三个点横坐标且第二个点横坐标和第四个点横坐标相等
x2=mtr(i+1,j,1);
if ((k11*j1+b11)<=i1)&& ((k22*j1+b22)>=i1)
d=1;
h=[];
for k0=-10:0.01:10
b0=i1-j1*k0;
%%%求交点坐标(x1,y1)(x2,y2)
x1=((b0-b1)/(k1-k0));
y1=(k1*b0-k0*b1)/(k1-k0);
if  double(y1)<mtr(i,j,2)||double(y1)>mtr(i,j+1,2)
continue;
end
y2=k0*x2+b0;%(x1,y1)(x2,y2)为要求的直线与左右两条直线的交点
if  double(y2)<mtr(i+1,j,2)||double(y2)>mtr(i+1,j+1,2)
continue;
end
%%%求距离
a=sqrt((x1-mtr(i,j,1))^2+(y1-mtr(i,j,2))^2);
b=sqrt((y2-mtr(i+1,j,2))^2);
A=sqrt((mtr(i,j,1)-mtr(i,j+1,1))^2+(mtr(i,j,2)-mtr(i,j+1,2))^2);%求已知两点之间的距离
B=sqrt((mtr(i+1,j,1)-mtr(i+1,j+1,1))^2+(mtr(i+1,j,2)-mtr(i+1,j+1,2))^2);
hh=a/A-b/B;
h=[h;k0,hh];
if sign(hh)~=sign(d)
break;
end
d=hh;
end
K=a/A;
%%%用来求已知点与上下两边的交点(x11,y11)(x22,y22)
dd=-1;
kkk=[];
for k00=-50:0.1:50
b00=i1-j1*k00;
x11=((b00-b11)/(k11-k00));
if  double(x11)<mtr(i,j,1)||double(x11)>mtr(i+1,j,1)
continue;
end
y11=(k11*b00-k00*b11)/(k11-k00);
x22=(b00-b22)/(k22-k00);
y22=(k22*b00-k00*b22)/(k22-k00);%(x11,y11)(x22,y22)为要求的直线与已知直线的下、上两个交点
%%%求距离
aa=sqrt((x11-mtr(i,j,1))^2+(y11-mtr(i,j,2))^2);
bb=sqrt((x22-mtr(i,j+1,1))^2+(y22-mtr(i,j+1,2))^2);
AA=sqrt((mtr(i,j,1)-mtr(i+1,j,1))^2+(mtr(i,j,2)-mtr(i+1,j,2))^2);%求已知两点之间的距离
BB=sqrt((mtr(i,j+1,1)-mtr(i+1,j+1,1))^2+(mtr(i,j+1,2)-mtr(i+1,j+1,2))^2);
kk=aa/AA-bb/BB;
kkk=[kkk;k00,kk];
if sign(aa/AA-bb/BB)~=sign(dd)
break;
end
dd=kk;
end
KK=aa/AA;
%%求左右两条边的成等比的比例值
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);
break;
end
else k1=(mtr(i,j,2)-mtr(i,j+1,2))/(mtr(i,j,1)-mtr(i,j+1,1));%%第一个点横坐标不等于第三个点横坐标且第二个点横坐标和第四个点横坐标也不相等
if k1==0
break;
end
b1=(mtr(i,j,1)*mtr(i,j+1,2)-mtr(i,j+1,1)*mtr(i,j,2))/(mtr(i,j,1)-mtr(i,j+1,1));
k2=(mtr(i+1,j,2)-mtr(i+1,j+1,2))/(mtr(i+1,j,1)-mtr(i+1,j+1,1));
if k2==0
break;
end
b2=(mtr(i+1,j,1)*mtr(i+1,j+1,2)-mtr(i+1,j+1,1)*mtr(i+1,j,2))/(mtr(i+1,j,1)-mtr(i+1,j+1,1));%左右两条直线的参数，其中k1,b1为左面的直线，另一组为右面的直线
if (i1-b1)/k1>j1||(i1-b2)/k2<j1  %将已知点的纵坐标分别带入左、右边方程，如果左边方程得到的横坐标大于已知点横坐标i1，说明该点不在左方直线的右侧，同理若右边方程的值
continue;                     %小于已知点横坐标i1，说明该点不在右方直线的左侧，若上述两者有一者不满足则跳出此层循环
end
if ((k11*j1+b11)<=i1)&& ((k22*j1+b22)>=i1)
d=1;
h=[];
for k0=-10:0.01:10
b0=i1-j1*k0;
%%%求交点坐标(x1,y1)(x2,y2)
x1=((b0-b1)/(k1-k0));
y1=(k1*b0-k0*b1)/(k1-k0);
if  double(y1)<mtr(i,j,2)||double(y1)>mtr(i,j+1,2)
continue;
end
x2=(b0-b2)/(k2-k0);
y2=(k2*b0-k0*b2)/(k2-k0);%(x1,y1)(x2,y2)为要求的直线与左右两条直线的交点
if  double(y2)<mtr(i+1,j,2)||double(y2)>mtr(i+1,j+1,2)
continue;
end
%%%求距离
a=sqrt((x1-mtr(i,j,1))^2+(y1-mtr(i,j,2))^2);
b=sqrt((x2-mtr(i+1,j,1))^2+(y2-mtr(i+1,j,2))^2);
A=sqrt((mtr(i,j,1)-mtr(i,j+1,1))^2+(mtr(i,j,2)-mtr(i,j+1,2))^2);%求已知两点之间的距离
B=sqrt((mtr(i+1,j,1)-mtr(i+1,j+1,1))^2+(mtr(i+1,j,2)-mtr(i+1,j+1,2))^2);
hh=a/A-b/B;
h=[h;k0,hh];
if sign(hh)~=sign(d)
break;
end
d=hh;
end
K=a/A;
%%%用来求已知点与上下两边的交点(x11,y11)(x22,y22)
dd=-1;
kkk=[];
for k00=-50:0.1:50
b00=i1-j1*k00;
x11=((b00-b11)/(k11-k00));
if  double(x11)<mtr(i,j,1)||double(x11)>mtr(i+1,j,1)
continue;
end
y11=(k11*b00-k00*b11)/(k11-k00);
x22=(b00-b22)/(k22-k00);
y22=(k22*b00-k00*b22)/(k22-k00);%(x11,y11)(x22,y22)为要求的直线与已知直线的下、上两个交点
%%%求距离
aa=sqrt((x11-mtr(i,j,1))^2+(y11-mtr(i,j,2))^2);
bb=sqrt((x22-mtr(i,j+1,1))^2+(y22-mtr(i,j+1,2))^2);
AA=sqrt((mtr(i,j,1)-mtr(i+1,j,1))^2+(mtr(i,j,2)-mtr(i+1,j,2))^2);%求已知两点之间的距离
BB=sqrt((mtr(i,j+1,1)-mtr(i+1,j+1,1))^2+(mtr(i,j+1,2)-mtr(i+1,j+1,2))^2);
kk=aa/AA-bb/BB;
kkk=[kkk;k00,kk];
%用来判断k发生符号变化
if sign(aa/AA-bb/BB)~=sign(dd)
break;
end
dd=kk;
end
KK=aa/AA;
y=round((mtr1(i+1,j,1)-mtr1(i,j,1))*KK+mtr1(i,j,1));
x=round((mtr1(i,j+1,2)-mtr1(i,j,2))*K+mtr1(i,j,2));
M(i,j,1)=x;
M(i,j,2)=y;
J(x,y)=I(i1,j1);%将畸变图上的灰度值移植到在矫正后图像的对应点上
break;
end
end
end
end
end
end
end
end
J=(J>30);
figure(4);imshow(J);

```