www.pudn.com > distor.rar > 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