shiyan.m

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

close all
clear 
clc

% ThreshLCD=30;%有效投影区的阈值

%---------||读入投影仪180背景图，计算投影区域||----------
%-------------LCD的180背景--------------
[FileNameLCD,PathNameLCD]=uigetfile({'*.bmp','Bitmap files(*.bmp)';...
        '*.*','All Files(*.*)'},'get the LCD Background images');   
if ~isequal([FileNameLCD,PathNameLCD],[0,0])
        File1=fullfile(PathNameLCD,strcat(num2str(180),'.bmp'));%fullfile使用指定部分建立文件全名
        LCDBack=imread(File1);                                  %strcat结成长串
end 

LCDBack=rgb2gray(LCDBack);
LCDBack=double(LCDBack);
%----------确定LCD的投影区域--------------
[CCDRow,CCDColumn]=size(LCDBack);
NN=zeros(CCDRow,CCDColumn);
%未处理公共区域的有效投影区
% for i=1:CCDRow
%     for j=1:CCDColumn
%         if LCDBack(i,j)<ThreshLCD
%             LCDBack(i,j)=0;
%         else
%             NN(i,j)=1;
%         end
%     end
% end

%-----------------用时间相位展开求解绝对相位,条纹条数呈负指数增长---------
%----------------------------||垂直变形条纹||---------------------------
FringeNumber=[32 31 30 28 24 16]; %32 ,32-2^0 ,32-2^1 ,32-2^2 ,32-2^3 ,32-2^4

%%24幅垂直变形条纹图
[FileNameVertical,PathNameVertical]=uigetfile({'*.bmp','Bitmap files(*.bmp)';...
        '*.*','All Files(*.*)'},'get the Vertical Deformed fringes');

%载入参考条纹位相减少计算时间
[FileNameOrig,PathNameOrig]=uigetfile({'*.mat','Mat files(*.mat)';...
        '*.*','All Files(*.*)'},'Load the Vertical Original Phase');

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

if ~isequal([FileNameVertical,PathNameVertical],[0,0])
    for i=1:6
        for j=1:4
            File2=fullfile(PathNameVertical,strcat(num2str(FringeNumber(i)),'TPUV',num2str(j),'.bmp'));
            File2=imread(File2);
            File2=rgb2gray(File2);
            %keyboard
            I(:,:,j,i)=File2;%%%num2str(X):把非整数数组转化成串
        end                         %%%%strcat：结成长串
    end
end
%I(:,:,1,：)--I(:,:,2,：)--I(:,:,3,：)--I(:,:,4,：)六个不同周期的条纹图的四步相移
%I(:,:,:,1）到 I(:,:,:,6)是条纹数从32到16的四步相移

dIo(:,:,:) = double(I(:,:,1,:))-double(I(:,:,3,:)); % 分母
dIe(:,:,:) = double(I(:,:,4,:))-double(I(:,:,2,:));  % 分子
Pw(:,:,32) = atan2(dIe(:,:,1),dIo(:,:,1));%deat_Pu[6,0]%%%第一步：解出每套条纹测量后的截断相位图
%keyboard
% Pw(:,:,32) = Pw(:,:,32).*NN;
Pw(:,:,31) = atan2(dIe(:,:,2),dIo(:,:,2));%deat_Pu[7,0]
% Pw(:,:,31) = Pw(:,:,31).*NN;
Pw(:,:,30) = atan2(dIe(:,:,3),dIo(:,:,3));%deat_Pu[8,0]
% Pw(:,:,30) = Pw(:,:,30).*NN;
Pw(:,:,28) = atan2(dIe(:,:,4),dIo(:,:,4));%deat_Pu[7,0]
% Pw(:,:,28) = Pw(:,:,28).*NN;
Pw(:,:,24) = atan2(dIe(:,:,5),dIo(:,:,5));%deat_Pu[8,0]
% Pw(:,:,24) = Pw(:,:,24).*NN;
Pw(:,:,16) = atan2(dIe(:,:,6),dIo(:,:,6));%deat_Pu[8,0]
% Pw(:,:,16) = Pw(:,:,16).*NN;

% Pw(:,:,32) = atan2(dIo(:,:,1),dIe(:,:,1));%deat_Pu[7,0]
% Pw(:,:,31) = atan2(dIo(:,:,2),dIe(:,:,2));
% Pw(:,:,30) = atan2(dIo(:,:,3),dIe(:,:,3));
% Pw(:,:,28) = atan2(dIo(:,:,4),dIe(:,:,4));
% Pw(:,:,24) = atan2(dIo(:,:,5),dIe(:,:,5));
% Pw(:,:,16) = atan2(dIo(:,:,6),dIe(:,:,6));
clear  dIe dIo 
%keyboard

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(Pw(:,:,2^s-t1)-Pw(:,:,1),deta_Pu(:,:,2^s-t1));
        deta_Pu(:,:,1) =  deta_tPu(:,:,1)+deta_Pu(:,:,2^s-t1);
        % deta_tPu(:,:,1) = unwrap2(Pw(:,:,16)-Pw(:,:,1),deta_Pu(:,:,16),1);
        % deta_Pu(:,:,1) =  deta_tPu(:,:,1)+deta_Pu(:,:,16);
    else 
%         
        deta_Pu(:,:,2^s-1) = wrap0(Pw(:,:,2^s)-Pw(:,:,2^s-1));
        
        %keyboard
          % deta_Pu(:,:,2^s-1) = Pw(:,:,2^s)-Pw(:,:,2^s-1);
        deta_tPu(:,:,2^s-2*t1) = unwrap2(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);
%        deta_tPu(:,:,30) = unwrap2(Pw(:,:,31)-Pw(:,:,30),deta_Pu(:,:,31),1);
%        deta_tPu(:,:,28) = unwrap2(Pw(:,:,30)-Pw(:,:,28),deta_Pu(:,:,30),1);
%        deta_tPu(:,:,24) = unwrap2(Pw(:,:,28)-Pw(:,:,24),deta_Pu(:,:,28),1);
%        deta_tPu(:,:,16) = unwrap2(Pw(:,:,24)-Pw(:,:,16),deta_Pu(:,:,24),1);
%        
%        deta_Pu(:,:,31) =wrap0(Pw(:,:,32)-Pw(:,:,31));%%%%32和31两套条纹投影后相邻两幅图像同一点的展开相位差
%        deta_Pu(:,:,30) =  deta_tPu(:,:,30)+deta_Pu(:,:,31);
%       deta_Pu(:,:,28) =  deta_tPu(:,:,28)+deta_Pu(:,:,30);
%        deta_Pu(:,:,24) =  deta_tPu(:,:,24)+deta_Pu(:,:,28);
%        deta_Pu(:,:,16) =  deta_tPu(:,:,16)+deta_Pu(:,:,24);

    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;%%%%%%总的展开相位差

% clear Pu32 Pu16 Pu24 Pu30 Pu28 Pu31   Pufit ws;

%------------载入参考条纹位相减少计算时间-----------
if ~isequal([FileNameOrig,PathNameOrig],[0,0])
    RUnwrappedPhaseVertical=fullfile(PathNameOrig,FileNameOrig);
    load(RUnwrappedPhaseVertical);
end

clear PathNameOrig,FileNameOrig;
[ProjectorRow,ProjectorColumn]=size(RUnwrappedPhaseVertical);

p=ProjectorColumn/32;

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

x1=zeros(CCDRow,CCDColumn);
x2=zeros(CCDRow,CCDColumn);

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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------------------------||水平变形条纹||---------------------------

%%24幅水平变形条纹图
[FileNameParal,PathNameParal]=uigetfile({'*.bmp','Bitmap files(*.bmp)';...
        '*.*','All Files(*.*)'},'get the  Parallel Deformed fringes');
    
 if ~isequal([FileNameParal,PathNameParal],[0,0])
    for i=1:6
        for j=1:4
            File3=fullfile(PathNameParal,strcat(num2str(FringeNumber(i)),'TPU',num2str(j),'.bmp'));
            File3=imread(File3);
            File3=rgb2gray(File3);
            I(:,:,j,i)=File3;%%%num2str(X):把非整数数组转化成串
        end                         %%%%strcat：结成长串
    end
end   
%     for i=1:6
%         for j=1:4
%             File3=fullfile(strrep(PathNameParal,'TPUV','TPU'),strcat(num2str(FringeNumber(i)),'TPUV',num2str(j),'.bmp'));
%               File3=imread(File3);
%               File3=rgb2gray(File3);
%             I(:,:,j,i)=File3;
%         end
%     end
%     clear FileNameParal,PathNameParal;
dIo(:,:,:) = double(I(:,:,1,:))-double(I(:,:,3,:)); % 分母
dIe(:,:,:) = double(I(:,:,4,:))-double(I(:,:,2,:));  % 分子
clear I;
Pw(:,:,32) = atan2(dIe(:,:,1),dIo(:,:,1));%deat_Pu[6,0]
% Pw(:,:,32) = Pw(:,:,32).*NN;
Pw(:,:,31) = atan2(dIe(:,:,2),dIo(:,:,2));%deat_Pu[7,0]
% Pw(:,:,31) = Pw(:,:,31).*NN;
Pw(:,:,30) = atan2(dIe(:,:,3),dIo(:,:,3));%deat_Pu[8,0]
% Pw(:,:,30) = Pw(:,:,30).*NN;
Pw(:,:,28) = atan2(dIe(:,:,4),dIo(:,:,4));%deat_Pu[7,0]
% Pw(:,:,28) = Pw(:,:,28).*NN;
Pw(:,:,24) = atan2(dIe(:,:,5),dIo(:,:,5));%deat_Pu[8,0]
% Pw(:,:,24) = Pw(:,:,24).*NN;
Pw(:,:,16) = atan2(dIe(:,:,6),dIo(:,:,6));%deat_Pu[8,0]
% Pw(:,:,16) = Pw(:,:,16).*NN;
clear dIe dIo 

s = 5;% 2^s 为最大的投影条纹数目
for i = 1:s
    t1 = 2^(i-1);
    if 2^s-2*t1 == 0 %%t1 = 32;
        Pw(:,:,1) = 0;
        deta_tPu(:,:,1) = unwrap2(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(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;

clear Pu32 Pu16 Pu24 Pu30 Pu28 Pu31 Pufit Pufit ws;

%-----------------横向对应点计算:摄像机像素映射到投影器坐标系--------------------
p=ProjectorRow/32;

[FileNameOrig,PathNameOrig]=uigetfile({'*.mat','Mat files(*.mat)';...
        '*.*','All Files(*.*)'},'Load the Original P Phase');
    
 if ~isequal([FileNameOrig,PathNameOrig],[0,0])
    RUnwrappedPhaseParallel=fullfile(PathNameOrig,FileNameOrig);
    load(RUnwrappedPhaseParallel);
end
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


%-------------Projector象素对应的CCD亚像素点计算---------------------------
%---------------------即反向坐标关系的建立---------------------------------
%-----------------二元三次projector座标平面插值------------------------------
%--------------------------------------------------------------------------
[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)=(124+64*cos(2*pi*xf1(k1,j1)/p+(i-1)*pi/2))/255;
%             end
        end
    end
    File4=fullfile(PathNameInv,strcat('C_CubicInverseL',num2str(i),'.bmp'));
    imwrite(x4l,File4);
end