插值查找（没标定）.m

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

%逐行逐列查找位相相同点，进行了插值和查找范围分类，程序较复杂，但只进行了垂直条纹的投影，存在对准问题
%仍没进行标定
p=32;%期望条纹周期
q=32;%等效周期
%--------------读入测量图象----------------
i1=imread('D:\matlab\work\cy\4-0.bmp');
i1=double(i1);
i2=imread('D:\matlab\work\cy\4-1.bmp');
i2=double(i2);
i3=imread('D:\matlab\work\cy\4-2.bmp');
i3=double(i3);
i4=imread('D:\matlab\work\cy\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);
x1=zeros(row,column);
x2=zeros(row,column);
%---------------  双三次插值  -------------------
%插值点数不宜过大，不然会造成数据量太大，占用内存过多
%而计算速度缓慢
[x,y]=meshgrid(1:0.2:column,1:1:row);
[krow,kcolumn]=size(x);
t3=interp2(t2,x,y,'bicubic');
%--------------------------------------建立几何反转关系----------------------------
%为了减少计算时间，考虑到条纹变形不会超出一定的范围，所以只需比较象素点周围一定区域的位相值
%用zone_x,zone_y分别表示区域的范围，由于条纹近似垂直，zone_y可大大小与zone_x
%先对y方向分三种情况，每种情况又包含x方向的三种情况
zone_x=round(kcolumn/column)*50;
zone_y=round(krow/row)*20;
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<zone_y+1
%．．．．．y的第一种情况的三种分情况．．．．．．
    if j1<zone_x+1
for k2=1:k1+zone_y
for j2=1:j1+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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
elseif j1>column-zone_x
for k2=1:k1+zone_y
for j2=j1-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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
else
for k2=1:k1+zone_y
for j2=j1-zone_x:j1+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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
end
%．．．．．．y的第二种情况三种分情况．．．．．．
elseif k1>row-zone_y
if j1<zone_x+1
for k2=k1-zone_y:krow
for j2=1:j1+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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
elseif j1>cloumn-zone_x
for k2=k1-zone_y:krow
for j2=j1-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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
else
for k2=k1-zone_y:krow
for j2=j1-zone_x:j1+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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
end
%．．．．．．y的点三种情况的分情况．．．．．．．
else   
if j1<zone_x+1
for k2=k1-zone_y:k1+zone_y
for j2=1:j1+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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
elseif j1>column-zone_x
for k2=k1-zone_y:k1+zone_y
for j2=j1-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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
break;%如果位相差值为零表明已找到，跳出列的比较
end
end
else
for k2=k1-zone_y:k1+zone_y
for j2=j1-zone_x:j1+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
break;%如果位相差值为零表明已找到，跳出行的比较
end
end
if b==0
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);
for k1=1:row
for j1=1:column
x4(k1,j1)=0.5*(1+cos(2*pi*x22(k1,j1)/p));%插值求反向条纹，没考虑背景和调制度
end
end
subplot(2,2,1),imshow(i1/256);
title('MASTER的测量条纹');
subplot(2,2,2),imshow(x3);
title('期望得到的条纹');
subplot(2,2,3),imshow(x4/2);
title('投影的反向条纹');