%逐行逐列查找位相相同點(diǎn)，進(jìn)行了插值和查找范圍分類，程序較復(fù)雜，但只進(jìn)行了垂直條紋的投影，存在對(duì)準(zhǔn)問題
%仍沒進(jìn)行標(biāo)定
p=32;%期望條紋周期
q=32;%等效周期
%--------------讀入測(cè)量圖象----------------
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);
%---------------  雙三次插值  -------------------
%插值點(diǎn)數(shù)不宜過大，不然會(huì)造成數(shù)據(jù)量太大，占用內(nèi)存過多
%而計(jì)算速度緩慢
[x,y]=meshgrid(1:0.2:column,1:1:row);
[krow,kcolumn]=size(x);
t3=interp2(t2,x,y,'bicubic');
%--------------------------------------建立幾何反轉(zhuǎn)關(guān)系----------------------------
%為了減少計(jì)算時(shí)間，考慮到條紋變形不會(huì)超出一定的范圍，所以只需比較象素點(diǎn)周圍一定區(qū)域的位相值
%用zone_x,zone_y分別表示區(qū)域的范圍，由于條紋近似垂直，zone_y可大大小與zone_x
%先對(duì)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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
end
if b==0
break;%如果位相差值為零表明已找到，跳出行的比較
end
end
if b==0
break;%如果位相差值為零表明已找到，跳出列的比較
end
end
end
%．．．．．．y的點(diǎn)三種情況的分情況．．．．．．．
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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));%找出差值最小的點(diǎn)
x1(k1,j1)=k2;%存儲(chǔ)橫坐標(biāo)
x2(k1,j1)=j2;%存儲(chǔ)縱坐標(biāo)
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
%----------------產(chǎn)生反向條紋------------------------
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));%插值求反向條紋，沒考慮背景和調(diào)制度
end
end
subplot(2,2,1),imshow(i1/256);
title('MASTER的測(cè)量條紋');
subplot(2,2,2),imshow(x3);
title('期望得到的條紋');
subplot(2,2,3),imshow(x4/2);
title('投影的反向條紋');