distortion.m

图像模板匹配计算

clc;
tic;
Ideal=imread('ideal.bmp');
Real=imread('real_d.bmp');
%figure;
%imshow(Real);           %显示实际图像
%figure;
%imshow(Ideal);            %显示理想图像

P_ideal=Ideal(1:1200,1:1600,1)*0.11+Ideal(1:1200,1:1600,2)*0.30+Ideal(1:1200,1:1600,3)*0.59;
P_real=Real(1:1200,1:1600,1)*0.11+Real(1:1200,1:1600,2)*0.30+Real(1:1200,1:1600,3)*0.59;
i=1;

while(i<=1200)
    j=1;
 while(j<=1600)
        if P_real(i,j)>200
            P_real(i,j)=255;
        else
            P_real(i,j)=0;
        end
        if P_ideal(i,j)>150
            P_ideal(i,j)=255;
        else
            P_ideal(i,j)=0;
        end
        j=j+1;
 end
    i=i+1;
end



%figure;
%imshow(P_ideal);           %显示实际图像
%figure;
%imshow(P_real);            %显示理想图像
p=0;
q=0;                       %控制点计数器
tic;
x=0;
y=0;

i=80;

while( (i>=80) & (i<1200-100))              %基于图像特征的，控制点寻找确定算法
       flag=0;
       j=30;
    while ( (j>=30) & (j<1600-60))
                        nu1=0;
                        nu2=0;
                        nu3=0;
                        nu4=0;
                       if length(find(P_ideal(i-(1:20),j)<1))>15
                           nu1=1;
                       end
                       if length(find(P_ideal(i+(1:20),j)<1))>15
                           nu2=1;
                       end
                       if length(find(P_ideal(i,j-(1:20))<1))>15
                           nu3=1;
                       end
                       if length(find(P_ideal(i,j+(1:20))<1))>15
                           nu4=1;
                       end
                       if (nu1==1 | nu2==1)&(nu3==1 | nu4==1)
                                         p=p+1;
                                         x(p)=j;
                                         y(p)=i;
                                         flag=flag+1;
                                         j=j+40;
                         else
                           j=j+1;
                       end
    end
    if flag==25
        i=i+40;
    else
        i=i+1;
    end
end
i=80;
j=30;
u(1:425)=0;
v(1:425)=0;
flag=0;
k=0;
while( (i>=80) & (i<1200-40))              %基于图像特征的，控制点寻找确定算法
     j=21;
    while ( (j>=21) & (j<1600-20))
        r_nu1=0;
        r_nu2=0;
        r_nu3=0;
        r_nu4=0;
      if length(find(P_real(i-(1:20),j)<1))>15
           r_nu1=1;
       end
       if length(find(P_real(i+(1:20),j)<1))>15
           r_nu2=1;
       end
       if length(find(P_real(i,j-(1:20))<1))>15
           r_nu3=1;
       end
       if length(find(P_real(i,j+(1:20))<1))>15
           r_nu4=1;
       end
       if (r_nu1==1 | r_nu2==1)&(r_nu3==1 | r_nu4==1)
              if length(find(u(q*25+1:(q+1)*25)))<25
                  w=uint16((j+30)/60-0.1);
                  if w>13
                      w=w-1;
                  end
                        if u(q*25+w)==0
                                u(q*25+w)=j;
                                v(q*25+w)=i;
                                flag=flag+1;
                                 if flag==25
                                  k=1;
                                 end
                              else
                                  if u(q*25+w)>j
                                      u(q*25+w)=j;
                                      v(q*25+w)=i;
                                  end
                              end
                             
                end
            if k==1
                j=1600;
            else
              j=j+40;
            end
       else
           j=j+1;
       end
    end
    if k==1
        i=i+30;
        q=q+1;
        flag=0;
        k=0;
       
    else
        i=i+1;
   end
end
%length(x)
%length(y)
%length(u)
%length(v)
%x(1:425)
%u(1:425)

%P_real=Real(1:1200,1:1600,1);
%for L=1:425
%P_real(v(L)-5:v(L)+5,u(L)-5:u(L)+5)=255;
%end
%figure;
%imshow(P_real);

pw=4;                                      %多项式次数
k=(pw+1)*pw/2;                             %系数向量的长度
T(1:k,1:k)=0;
T0(1:k,1:k)=0;
T1(1:k,1:k)=0;
X(1:k)=0;
Y(1:k)=0;
X0(1:k)=0;
Y0(1:k)=0;
delat=0;
L=1;
while L<=425
        i=1;
        j=1;
        k=1;
            while k<=pw
                p=1;
                while p<=pw-k+1
                   m=1;
                    while m<=pw
                       n=1;
                        while n<=pw-m+1
                         T0(i,j)=u(L)^(k+m-2)*v(L)^(p+n-2);
                         j=j+1;
                         n=n+1;
                         end
                         m=m+1;
                    end
                    j=1;
                    i=i+1;
                    p=p+1;
                end
               k=k+1;
            end
            j=1;
      m=1;     
     while m<=pw
         n=1;
         while n<=pw-m+1
              X0(j)=u(L)^(m-1)*v(L)^(n-1);
               j=j+1;
              n=n+1;
         end
         m=m+1;
     end  
   delat=delat+(x(L)-u(L))^2+(y(L)-v(L))^2;
    X=X+X0*x(L);
    Y=Y+X0*y(L);
    T=T+T0;
  L=L+1; 
end
%T=T/(1.0e+011)
%X=X/(1.0e+011)
%Y=Y/(1.0e+011)
delat=(1/425)*delat^(0.5)
%T1=pinv(T)
T1=inv(T);
T2=pinv(T);
if pw>=4
a=T1*X';
b=T1*Y';
else
a=T2*X';
b=T2*Y';
end
a=a';
b=b';
delat=0;
P_real=Real(1:1200,1:1600,1);
n=1;
p=1;
q=1;
 while n<=1200
    m=1;
   while m<=1600
        m1=0;
        n1=0;
        k=1;
        i=1;
        while i<=pw
        j=1;
        while j<=pw-i+1
              m1=m1+a(k)*m^(i-1)*n^(j-1);
              n1=n1+b(k)*m^(i-1)*n^(j-1);
              k=k+1;
             j=j+1;
         end
        i=i+1;
        end  
        m2=uint16(m1);
        n2=uint16(n1);
       if m2>m1
          p1=m2-m1;
          m1=uint16(m1-0.5);
       else
         p1=m1-m2;
         m1=uint16(m1);
       end
       if n2>n1
         q1=n2-n1;
         n1=uint16(n1-0.5);
       else
          q1=n1-n2;
          n1=uint16(n1);
       end
      
    %  uint8(q1);
    %  uint8(p1);
%    P_ideal_real(n1,m1)=(1-double(q1))*(double(q1)*P_real(n+1,m+1)+(1-double(q1))*P_real(n,m+1))+double(q1)*(double(q1)*P_real(n+1,m)+(1-double(q1))*P_real(n,m));
      %P_ideal_real(n1,m1)=(1-uint8(q1))*(uint8(q1)*P_real(n+1,m+1)+(1-uint8(q1))*P_real(n,m+1))+uint8(q1)*(uint8(q1)*P_real(n+1,m)+(1-uint8(q1))*P_real(n,m));
      %P_ideal_real(n1,m1)=(1-q1)*(q1*P_real(n+1,m+1)+(1-q1)*P_real(n,m+1))+q1*(q1*P_real(n+1,m)+(1-q1)*P_real(n,m));
    
  %       if m1<1
  %          m1=1;
  %      end
 %       if n1<1
 %           n1=1;
 %       end
 %       if m1>1600
 %           m1=1600;
 %       end
 %       if n1>1200
 %           n1=1200;
 %       end
        P_ideal_real(n1,m1)=P_real(n,m);
     if (u(p)==m)&(v(p)==n)
            delat=delat+(m1-x(p))^2+(n1-v(p))^2;
            p=p+1;
        end
      m=m+1;
   end
 n=n+1;
 end
delat=double(delat);
delat=(1/425)*delat^(0.5)
%figure;
%imshow(P_ideal_real);


border1_left=min(u);
border2_right=max(u);
border1_up=min(v);
border2_down=max(v);
sector=5;
i=border1_up;
i1=1;
while (i<=border2_down)
    j=border1_left;
    k=1;
    while(j<=border2_right)
      Pin(i1:i1+sector,k:k+sector)=P_ideal_real(i,j);
      j=j+1;
      k=k+sector+1;
    end
    i=i+1;
    i1=i1+sector+1;
end
size(Pin)
imwrite(Pin,'Pin.jpg','jpg');
t=toc