%xy方向的插值比較，標定了，查找范圍改進后應該會有結果
%------行方向坐標－－－－－－
p=64;%條紋周期
q=64;%等效周期
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));
ty2=myunwrappro(t1,1,146,539);
ty2(1,:)=ty2(2,:);%消除邊緣毛刺
[row,column]=size(i1);%提取圖像大小，為處理方便
r1=imread('D:\matlab\work\cy2\i1.bmp');
r1=double(r1);
r2=imread('D:\matlab\work\cy2\i2.bmp');
r2=double(r2);
r3=imread('D:\matlab\work\cy2\i3.bmp');
r3=double(r3);
r4=imread('D:\matlab\work\cy2\i4.bmp');
r4=double(r4);
r1=atan2((r4-r2),(r1-r3));
ry0=myunwrappro(r1,1,154,526);
ry0=ry0+ty2(146,539)-ry0(154,526);
k=1:0.2:column;
for i=1:row
ty3(i,:)=interp1(1:column,ty2(i,:),k,'cubic');
end
[krow1,kcolumn]=size(k);
clear i3;clear i4;clear t1;%清除不用數據，減少內存占用
clear r2;clear r3;clear r4;
%----------列方向坐標－－－－－－－－
p=32;
q=32;
i1=imread('D:\matlab\work\cy2\iy1.bmp');
i1=double(i1);
i2=imread('D:\matlab\work\cy2\iy2.bmp');
i2=double(i2);
i3=imread('D:\matlab\work\cy2\iy3.bmp');
i3=double(i3);
i4=imread('D:\matlab\work\cy2\iy4.bmp');
i4=double(i4);%讀入測量圖象
t1=atan2((i4-i2),(i1-i3));
tx2=myunwrappro(t1,1,146,539);
tx2(:,1)=tx2(:,2);%消除邊緣毛刺
r1=imread('D:\matlab\work\cy2\r1.bmp');
r1=double(r1);
r2=imread('D:\matlab\work\cy2\r2.bmp');
r2=double(r2);
r3=imread('D:\matlab\work\cy2\r3.bmp');
r3=double(r3);
r4=imread('D:\matlab\work\cy2\r4.bmp');
r4=double(r4);
r1=atan2((r4-r2),(r1-r3));
rx0=myunwrappro(r1,1,154,526);
rx0=rx0+tx2(146,539)-rx0(154,526);
kx=1:0.2:row;
kx=rot90(rot90(rot90(kx)));
[krow,kcolumn1]=size(kx);
for i=1:column
tx3(:,i)=interp1(1:row,tx2(:,i),kx,'cubic');
%雙三次插值
end
for i=1:kcolumn
    ty4(:,i)=interp1(1:row,ty3(:,i),kx,'cubic');
end
clear ty2;clear ty3;
for i=1:krow
tx4(i,:)=interp1(1:column,tx3(i,:),k,'cubic');
end
clear i3;clear i4;clear t1;%清除不用數據，減少內存占用
clear r2;clear r3;clear r4;clear tx2;
clear tx3;
%————————坐標反轉------------------
x1=zeros(row,column);
x2=zeros(row,column);
for k1=1:row
    for j1=1:column
        x1(k1,j1)=k1;
        x2(k1,j1)=j1;
        a=abs(rx0(k1,j1)-tx4(k1,j1));
        b=abs(ry0(k1,j1)-ty4(k1,j1));
        for k2=max(1,k1*5-4-200):min(k1*5+200,row*5-8)
            for j2=max(400,j1*5-4+200-200):min(j1+200+200,column*5-8)
                if a>=abs(rx0(k1,j1)-tx4(k2,j2))&b>=abs(ry0(k1,j1)-ty4(k2,j2))
                    a=abs(rx0(k1,j1)-tx4(k2,j2));
                    b=abs(ry0(k1,j1)-ty4(k2,j2));
                    x1(k1,j1)=k2;
                    x2(k1,j1)=j2;
                end
                if a<0.01&b<0.01
                    break;
                end
            end
            if a<0.01&b<0.01
                break;
            end
        end
    end
end