插值分类查找mod.m

这是在MATLAB下生成的彩色条纹,对于傅里叶变换很有用,可以扩大实验的研究范围,更好的设计实验

p=32;%条纹周期
q=32;%等效周期
i1=imread('D:\matlab\work\cy2\4-0.bmp');
i1=double(i1);
i2=imread('D:\matlab\work\cy2\4-1.bmp');
i2=double(i2);
i3=imread('D:\matlab\work\cy2\4-2.bmp');
i3=double(i3);
i4=imread('D:\matlab\work\cy2\4-3.bmp');
i4=double(i4);%读入测量图象
t1=atan2((i4-i2),(i1-i3));
t2=myunwrap(t1);
t2(1,:)=t2(2,:);%消除边缘毛刺
[row,column]=size(i1);%提取图像大小，为处理方便
for i=1:row
r1(i,:)=i1(3,:);%初始条纹图，从变形条纹图中提取
r2(i,:)=i2(3,:);
r3(i,:)=i3(3,:);
r4(i,:)=i4(3,:);
end
r1=atan2((r4-r2),(r1-r3));
r0=myunwrap(r1);
clear i3;clear i4;clear t1;%清除不用数据，减少内存占用
clear r2;clear r3;clear r4;
%--------利用解方程的方法求出条纹调制度和背景光强------
% ri--I0, ra--I0*Y
for i=1:row
for j=1:column
a=[1 cos(t2(i,j));
1 -sin(t2(i,j))];
c=[i1(i,j); i2(i,j)];
x=a\c;
ri(i,j)=x(1,1);
ra(i,j)=x(2,1);
end
end  
ra1=medfilt2(ra,[5,5]);
ri1=medfilt2(ri,[5,5]);%用中值滤波消除方程求解过程中出现的奇异点
%---------------  双三次插值  -------------------
%插值点数不宜过大，不然会造成数据量太大，占用内存过多
%而计算速度缓慢
x1=zeros(row,column);
x2=zeros(row,column);
[x,y]=meshgrid(1:0.25:column,1:0.5:row);
[krow,kcolumn]=size(x);
t3=interp2(t2,x,y,'bicubic');%位相插值
rii=interp2(ri1,x,y,'bicubic');%背景
raa=interp2(ra1,x,y,'bicubic');%调制度
%--------------------------------------建立几何反转关系----------------------------
%为了减少计算时间，考虑到条纹变形不会超出一定的范围，所以只需比较象素点周围一定区域的位相值
%用zone_x,zone_y分别表示区域的范围，由于条纹近似垂直，zone_y可大大小与zone_x
%先对y方向分三种情况，每种情况又包含x方向的三种情况
zone_x=round(kcolumn/column)*50;
zone_y=round(krow/row)*10;
mul_x=4;%列的扩大倍数
mul_y=2;%行
offset_x=3;%列的位移量
offset_y=1;%行
%-------------
for k1=1:row
for j1=1:column
b=abs(r0(k1,j1)-t2(k1,1));
x1(k1,j1)=k1;
x2(k1,j1)=1;
if k1<10+1
%．．．．．y的第一种情况的三种分情况．．．．．．
if j1<50+1
for k2=1:k1*mul_y-offset_y+zone_y
  for j2=1:j1*mul_x-offset_x+zone_x
    if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));%找出差值最小的点
       x1(k1,j1)=k2;%存储横坐标
       x2(k1,j1)=j2;%存储纵坐标
    end
       if b<0.01
          break;%如果位相差值满足条件表明已找到，跳出行的比较
       end
  end
       if b<0.01
          break;%如果位相差值满足条件表明已找到，跳出列的比较
       end
end
elseif j1>column-50
for k2=1:k1*mul_y-offset_y+zone_y
for j2=j1*mul_x-offset_x-zone_x:kcolumn
if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
else 
    for k2=1:k1*mul_y-offset_y+zone_y
    for j2=j1*mul_x-offset_x-zone_x:j1*mul_x-offset_x+zone_x
if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
end
%．．．．．y的第2种情况的三种分情况．．．．．．
elseif k1>row-10
if j1<50+1
    for k2=k1*mul_y-offset_y-zone_y:krow
for j2=1:j1*mul_x-offset_x+zone_x
    if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
elseif j1>column-50
for k2=k1*mul_y-offset_y-zone_y:krow
for j2=j1*mul_x-offset_x-zone_x:kcolumn
if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
else 
    for k2=k1*mul_y-offset_y-zone_y:krow
    for j2=j1*mul_x-offset_x-zone_x:j1*mul_x-offset_x+zone_x
if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
end
%．．．．．y的第3种情况的三种分情况．．．．．．
else
    if j1<50+1
    for k2=k1*mul_y-offset_y-zone_y:k1*mul_y-offset_y+zone_y
    for j2=1:j1*mul_x-offset_x+zone_x
    if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
elseif  j1>column-50
    for k2=k1*mul_y-offset_y-zone_y:k1*mul_y-offset_y+zone_y
    for j2=j1*mul_x-offset_x-zone_x:kcolumn
if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
else 
    for k2=k1*mul_y-offset_y-zone_y:k1*mul_y-offset_y+zone_y
    for j2=j1*mul_x-offset_x-zone_x:j1*mul_x-offset_x+zone_x
if b>=abs(r0(k1,j1)-t3(k2,j2))
       b=abs(r0(k1,j1)-t3(k2,j2));
       x1(k1,j1)=k2;
       x2(k1,j1)=j2;
    end
       if b<0.01
          break;
       end
  end
       if b<0.01
          break;
       end
end
end
end
end
end
%----------------期望条纹图------------------------
x3=zeros(row,column);
for i=1:column
x3(:,i)=0.5*(1+cos(2*pi*i/p));
end
for i=1:row
for j=1:column
x22(i,j)=x(x1(i,j),x2(i,j));
end
end
%----------------产生反向条纹------------------------
x4=zeros(row,column);
x4=zeros(row,column);
for k1=1:row
for j1=1:column
x4(k1,j1)=(128*(1+cos(2*pi*x22(k1,j1)/p))-rii(x1(k1,j1),x2(k1,j1))+raa(x1(k1,j1),x2(k1,j1)))./(2*raa(x1(k1,j1),x2(k1,j1))+eps);
end
end 
figure,subplot(2,2,1),imshow(i1/256);
title('MASTER的测量条纹');
subplot(2,2,2),imshow(x3);
title('期望得到的条纹');
subplot(2,2,3),imshow(x4/2);
title('投影的反向条纹');