untitled4.m

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

close all
clear 
clc
warning off MATLAB:griddata:DuplicateDataPoints

CCDStartR=1;CCDStartC=1;%CCD相位展开起始点
ProjectorStartR=1;ProjectorStartC=1;%Projector相位展开起始点
CCDRow=512;CCDColumn=512;
ProjectorRow=512;ProjectorColumn=512;

%确定反向条纹写入路径
[FileNameInv,PathNameInv]=uiputfile({'*.bmp';'*.*'},'Writen the Inverse Image to disk');

%----------------------水平条纹--------------------------------
%读入拍摄到的投影到物体表面上的水平条纹图(PMP四步相移)
a = [32 31 30 28 24 16];
p=512./a;
M = 512;
N = 512;
% x=meshgrid(-M/2:(M/2-1),-N/2:(N/2-1));

thea = 3/5;
obj = zeros(M,N);
obj= peaks(M);
obj = obj*thea;

x= meshgrid(-M/2:M/2-1,-N/2:N/2-1);
y= meshgrid(-N/2:N/2-1,-M/2:M/2-1);

for i=1:6
h(:,:,1,i) = 127.5+127.5*cos(2*pi*(x+obj)/p(i)); 
h(:,:,2,i) = 127.5+127.5*cos(2*pi*(x+obj)/p(i)+pi/2);
h(:,:,3,i)= 127.5+127.5*cos(2*pi*(x+obj)/p(i)+pi);
h(:,:,4,i) = 127.5+127.5*cos(2*pi*(x+obj)/p(i)+3*pi/2);
Pw(:,:,a(i))=atan2(h(:,:,4,i)-h(:,:,2,i),h(:,:,1,i)-h(:,:,3,i));
end

v1= 127.5+127.5*cos(2*pi*x/p(1)); 
v2 = 127.5+127.5*cos(2*pi*x/p(2)+pi/2);
v3= 127.5+127.5*cos(2*pi*x/p(3)+pi);
v4 = 127.5+127.5*cos(2*pi*x/p(4)+3*pi/2);
wrappedv=atan2((v4-v2),(v1-v3));
 RUnwrappedPhaseVertical=(unwrap(unwrap(wrappedv)'))';

s = 5;% 2^s 为最大的投影条纹数目
for i = 1:s
    t1 = 2^(i-1);%%t1=1 2 4 8 16
    if 2^s-2*t1 == 0 %%t1 = 16;
        Pw(:,:,1) = 0;
        deta_tPu(:,:,1) = unwrap2(wrap0(Pw(:,:,2^s-t1)-Pw(:,:,1)),deta_Pu(:,:,2^s-t1));
        deta_Pu(:,:,1) =  deta_tPu(:,:,1)+deta_Pu(:,:,2^s-t1);
    else 
        deta_Pu(:,:,2^s-1) = wrap0(Pw(:,:,2^s)-Pw(:,:,2^s-1));
        deta_tPu(:,:,2^s-2*t1) = unwrap2(wrap0(Pw(:,:,2^s-t1)-Pw(:,:,2^s-2*t1)),deta_Pu(:,:,2^s-t1));
        deta_Pu(:,:,2^s-2*t1) =  deta_tPu(:,:,2^s-2*t1)+deta_Pu(:,:,2^s-t1);
    end 
end
% clear Pw;
Pu32 = deta_Pu(:,:,1);
clear deta_Pu;
Pu16 = deta_tPu(:,:,1);
Pu24 = deta_tPu(:,:,16)+Pu16;
Pu28 = deta_tPu(:,:,24)+Pu24;
Pu30 = deta_tPu(:,:,28)+Pu28;
Pu31 = deta_tPu(:,:,30)+Pu30;
clear deta_tPu;
%————————————最小二乘法拟合出斜率————————————————
ws = (32*Pu32+31*Pu31+30*Pu30+28*Pu28+24*Pu24+16*Pu16)/(32^2+31^2+30^2+28^2+24^2+16^2);
Pufit = ws*32;
UnwrappedPhaseVertical=Pufit;
p=ProjectorColumn/32;

RUnwrappedPhaseVertical=RUnwrappedPhaseVertical-32*pi;
UnwrappedPhaseVertical=UnwrappedPhaseVertical-RUnwrappedPhaseVertical(1,1)+2*pi/p;
RUnwrappedPhaseVertical=RUnwrappedPhaseVertical-RUnwrappedPhaseVertical(1,1)+2*pi/p;

for k1=1:CCDRow
    for j1=1:CCDColumn
%         if NN(k1,j1)==1
           x1(k1,j1)=UnwrappedPhaseVertical(k1,j1)*p/(2*pi);
%         end
    end
end
for k1=1:CCDRow
    for j1=1:CCDColumn
%         if NN(k1,j1)==1
           if round(x1(k1,j1))<1
              x1(k1,j1)=1;
           elseif round(x1(k1,j1))>ProjectorColumn
              x1(k1,j1)=ProjectorColumn;
           end
%        end
    end
end


%----------------------竖直条纹--------------------------------
%读入拍摄到的投影到物体表面上的竖直条纹图(PMP四步相移)
a = [32 31 30 28 24 16];
p=512./a;

p1= 127.5+127.5*cos(2*pi*y/p(1)); 
p2 = 127.5+127.5*cos(2*pi*y/p(2)+pi/2);
p3= 127.5+127.5*cos(2*pi*y/p(3)+pi);
p4 = 127.5+127.5*cos(2*pi*y/p(4)+3*pi/2);
wrappedp=atan2((p4-p2),(p1-p3));
 RUnwrappedPhaseParallel=(unwrap(unwrap(wrappedp)'))';

for i=1:6
h(:,:,1,i) = 127.5+127.5*cos(2*pi*(y+obj)/p(i)); 
h(:,:,2,i) = 127.5+127.5*cos(2*pi*(y+obj)/p(i)+pi/2);
h(:,:,3,i) = 127.5+127.5*cos(2*pi*(y+obj)/p(i)+pi);
h(:,:,4,i) = 127.5+127.5*cos(2*pi*(y+obj)/p(i)+3*pi/2);
Pw(:,:,a(i))=atan2(h(:,:,4,i)-h(:,:,2,i),h(:,:,1,i)-h(:,:,3,i));
end

s = 5;% 2^s 为最大的投影条纹数目
for i = 1:s
    t1 = 2^(i-1);%%t1=1 2 4 8 16
    if 2^s-2*t1 == 0 %%t1 = 16;
        Pw(:,:,1) = 0;
        deta_tPu(:,:,1) = unwrap2(wrap0(Pw(:,:,2^s-t1)-Pw(:,:,1)),deta_Pu(:,:,2^s-t1));
        deta_Pu(:,:,1) =  deta_tPu(:,:,1)+deta_Pu(:,:,2^s-t1);
    else 
        deta_Pu(:,:,2^s-1) = wrap0(Pw(:,:,2^s)-Pw(:,:,2^s-1));
        deta_tPu(:,:,2^s-2*t1) = unwrap2(wrap0(Pw(:,:,2^s-t1)-Pw(:,:,2^s-2*t1)),deta_Pu(:,:,2^s-t1));
        deta_Pu(:,:,2^s-2*t1) =  deta_tPu(:,:,2^s-2*t1)+deta_Pu(:,:,2^s-t1);
    end 
end
% clear Pw;
Pu32 = deta_Pu(:,:,1);
clear deta_Pu;
Pu16 = deta_tPu(:,:,1);
Pu24 = deta_tPu(:,:,16)+Pu16;
Pu28 = deta_tPu(:,:,24)+Pu24;
Pu30 = deta_tPu(:,:,28)+Pu28;
Pu31 = deta_tPu(:,:,30)+Pu30;
clear deta_tPu;
%————————————最小二乘法拟合出斜率————————————————
ws = (32*Pu32+31*Pu31+30*Pu30+28*Pu28+24*Pu24+16*Pu16)/(32^2+31^2+30^2+28^2+24^2+16^2);
Pufit = ws*32;
UnwrappedPhaseParallel=Pufit;



p=ProjectorRow/32;

RUnwrappedPhaseParallel=RUnwrappedPhaseParallel-32*pi;
UnwrappedPhaseParallel=UnwrappedPhaseParallel-RUnwrappedPhaseParallel(1,1)+2*pi/p;
RUnwrappedPhaseParallel=RUnwrappedPhaseParallel-RUnwrappedPhaseParallel(1,1)+2*pi/p;

for k1=1:CCDRow
    for j1=1:CCDColumn
%         if NN(k1,j1)==1
           x2(k1,j1)=UnwrappedPhaseParallel(k1,j1)*p/(2*pi);
%         end
    end
end
for k1=1:CCDRow
    for j1=1:CCDColumn
%         if NN(k1,j1)==1
           if round(x2(k1,j1))<1
              x2(k1,j1)=1;
           elseif round(x2(k1,j1))>ProjectorRow
              x2(k1,j1)=ProjectorRow;
           end
%        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[x,y]=meshgrid(1:CCDColumn,1:CCDRow);
% xf2=zeros(ProjectorRow,ProjectorColumn);%列坐标的反向传递关系
xf1=zeros(ProjectorRow,ProjectorColumn);%行坐标的反向传递关系
[x1Pro,x2Pro]=meshgrid(1:ProjectorColumn,1:ProjectorRow);
xf1=griddata(x1,x2,x,x1Pro,x2Pro,'cubic');
% xf1=griddata(x2,x1,y,x1Pro,x2Pro,'cubic');
% xf2=xf2.*xcClosest;
% xf1=xf1.*xcClosest;
%--------------------------------------------------------------------------
%--------------------------产生满调制度反向条纹----------------------------
%-------------------------------------------------------------------------
p=16;
for i=1:4
    x4l=zeros(ProjectorRow,ProjectorColumn);
    for k1=1:ProjectorRow
        for j1=1:ProjectorColumn
%             if xf2(k1,j1)~=0
                x4l(k1,j1)=(0.5+0.5*cos(2*pi*xf1(k1,j1)/p+(i-1)*pi/2));
%             end
        end
    end
    File4=fullfile(PathNameInv,strcat('C_CubicInverseL',num2str(i),'.bmp'));
    imwrite(x4l,File4);
end






% wpH = atan2((h4-h2),(h1-h3));
% [uwpH] = zlf_dunwrap_m(wpH,ones(M,N),[1 1]);
% idisp(uwpH),title('水平方向上的展开相位分布');
% 
% %----------------------竖直条纹--------------------------------
% %读入拍摄到的投影到物体表面上的竖直条纹图(PMP四步相移)
% v1 = 127.5+127.5*cos(2*pi*(y+obj)/p); 
% v2 = 127.5+127.5*cos(2*pi*(y+obj)/p+pi/2);
% v3 = 127.5+127.5*cos(2*pi*(y+obj)/p+pi);
% v4 = 127.5+127.5*cos(2*pi*(y+obj)/p+3*pi/2);
% 
% wpV = atan2((v4-v2),(v1-v3));
% [uwpV] = zlf_dunwrap_m(wpV,ones(M,N),[1 1]);
% idisp(uwpV),title('垂直方向上的展开相位分布');
% 
% %----------------------计算反向条纹--------------------------------
% %从物体表面的相位分布计算要投影的反向条纹;
% l = p*uwpV/(2*pi)+1;
% m = p*uwpH/(2*pi)+1;
% 
% ii = (meshgrid(1:M,1:N))';
% jj = meshgrid(1:N,1:M);
% 
% lM = round(max(l(:)));
% lm = round(min(l(:)));
% mM = round(max(m(:)));
% mm = round(min(m(:)));
% [x1,y1] = meshgrid(lm:lM,mm:mM);
% z1 = griddata(l,m,jj,x1,y1,'cubic');      %cubic,nearest,linear
% idisp(z1);
% InvFriV = 127.5 + 127.5*cos(2*pi*z1/p);
% idisp(InvFriV);title('垂直反向条纹计算结果');grid
% idisp(255*(v1-min(v1(:)))/(max(v1(:))-min(v1(:))));title('初始投影垂直条纹计算结果');grid
% 
% z2 = griddata(l,m,ii,x1,y1,'cubic');      %cubic,nearest,linear
% idisp(z2);
% InvFriH = 127.5 + 127.5*cos(2*pi*z2/p);
% idisp(InvFriH);title('水平反向条纹计算结果');grid
% idisp(255*(h1-min(h1(:)))/(max(h1(:))-min(h1(:))));title('初始投影水平条纹计算结果');grid
% 
% 
% d = 127.5+127.5*cos(2*pi*(z1+obj)/p);
% idisp(d); title('投影效果');grid
% idisp(d-v1); title('Difference');grid
% 
% aa = 2*pi*(z2-x).^2;
% an = aa(find(~isnan(aa)));
% [Num,temp] = size(an);
% max(max(sqrt(an)))
% sum(sqrt(an))/Num
% %平均误差
% Di = sqrt(sum(an)/(Num-1))
% return
% imwrite(InvFirV,gray(256),'InverseImageProV.bmp');
% imwrite(InvFirH,gray(256),'InverseImageProH.bmp');
% 
%