www.pudn.com > fc.rar > ed3333.asv, change:2015-05-15,size:19982b


close all; 
clear all; 
clc 
load mtr mtr; 
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+'); 
%I=imread('E:\小英article\查表法\un6\useful\fig1.jpg'); 
%C:\Users\Administrator\Desktop\fig1.jpg 
I=imread('fig1.jpg'); 
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);