test1.m

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

close all
clear 
clc

Range=1;%插值邻域范围，1为4邻域
% ThreshLCD=30;%有效投影区的阈值

%------------------CCD象素对应的Projector亚像素点计算----------------------
%---------||读入投影仪180背景图，计算投影区域||----------
%-------------LCD的180背景--------------
[FileNameLCD,PathNameLCD]=uigetfile({'*.bmp','Bitmap files(*.bmp)';...
        '*.*','All Files(*.*)'},'get the LCD Background images');   %读入投影仪180背景图的文件名和路径
if ~isequal([FileNameLCD,PathNameLCD],[0,0])%isequal若两组数相同则为真(if:如果后面这个语句为真)
        File=fullfile(PathNameLCD,strcat(num2str(180),'.bmp'));%fullfile使用指定部分建立文件全名
                                                               %strcat结成长串
        LCDBack=imread(File);
end
LCDBack=rgb2gray(LCDBack);
% LCDBack=double(LCDBack);
%----------确定LCD的投影区域--------------
[CCDRow,CCDColumn]=size(LCDBack);%CCDRow为多少行，如（480乘640的一个矩阵，CCDRow=480,CCDColumn=640）
% 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
%水平变形条纹
[FileNameParal,PathNameParal]=uigetfile({'*.bmp','Bitmap files(*.bmp)';...
        '*.*','All Files(*.*)'},'get the LCD Parallel Deformed fringes');%%24幅水平变形条纹图
%载入参考条纹位相减少计算时间

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[FileNameOrig,PathNameOrig]=uigetfile({'*.mat','Mat files(*.mat)';...
        '*.*','All Files(*.*)'},'Load the Original Phase');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

if ~isequal([FileNameParal,PathNameParal],[0,0])
    for i=1:6
        for j=1:4
            File=fullfile(PathNameParal,strcat(num2str(FringeNumber(i)),'TPUV',num2str(j),'.bmp'));
            T=imread(File);
            File=rgb2gray(T);
            I(:,:,j,i)=File;%%%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]%%%第一步：解出每套条纹测量后的截断相位图
% 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 = 16;32-2^4
        Pw(:,:,1) = 0;
        deta_tPu(:,:,1) = unwrap2(Pw(:,:,2^s-t1)-Pw(:,:,1),deta_Pu(:,:,2^s-t1),1);
        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));
        deta_tPu(:,:,2^s-2*t1) = unwrap2(Pw(:,:,2^s-t1)-Pw(:,:,2^s-2*t1),deta_Pu(:,:,2^s-t1),1);
        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);
%竖向对应点计算:摄像机像素映射到投影器坐标系--RUnwrappedPhaseParallel=RUnwrappedPhaseParallel-32*pi;
p=ProjectorRow/32;

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))>ProjectorRow
              x1(k1,j1)=ProjectorRow;
           end
%        end
    end
end
%----------------------------||垂直变形条纹||---------------------------
    for i=1:6
        for j=1:4
            File=fullfile(strrep(PathNameParal,'TPUV','TPU'),strcat(num2str(FringeNumber(i)),'TPUV',num2str(j),'.bmp'));
              T1=imread(File);
              T1=rgb2gray(T1);
            I(:,:,j,i)=T1;
        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),1);
        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),1);
        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=ProjectorColumn/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))>ProjectorColumn
              x2(k1,j1)=ProjectorColumn;
           end
%        end
    end
end

%--------------------------A_ClosestInverseL----------------------
%-------------Projector象素对应的CCD亚像素点计算---------------------------
%---------------------即反向坐标关系的建立---------------------------------
%-----------------最近邻projector座标平面插值------------------------------
%--------------------------------------------------------------------------

%----------------------四领域坐标线性插值-------------------
xa=zeros(ProjectorRow,ProjectorColumn);%行坐标的各点权重和
xb=zeros(ProjectorRow,ProjectorColumn);%列坐标的各点权重和
xc=zeros(ProjectorRow,ProjectorColumn);%标志位
xf1=zeros(ProjectorRow,ProjectorColumn);%行坐标的反向传递关系
xf2=zeros(ProjectorRow,ProjectorColumn);%列坐标的反向传递关系
for k1=1:CCDRow
    for j1=1:CCDColumn
%         if  NN(k1,j1)==1%只计算projector有效照明区域
            CountY=x2(k1,j1);
            CountX=x1(k1,j1);%与CCD象素点对应的DLP亚象素点
            LowestY=floor(CountY);
            LowestX=floor(CountX);
            for k2=max(1,LowestX):min(ProjectorRow,LowestX+1)
                for j2=max(1,LowestY):min(ProjectorColumn,LowestY+1)
                    DistanceX=1.1-abs(k2-CountX);
                    DistanceY=1.1-abs(j2-CountY);
                    if DistanceY>xb(k2,j2)
                        xb(k2,j2)=DistanceY;
                        xf2(k2,j2)=j1;
                        xc(k2,j2)=1;
                    end
                    if DistanceX>xa(k2,j2)
                        xa(k2,j2)=DistanceX;
                        xf1(k2,j2)=k1;
                        xc(k2,j2)=1;
                    end
                end
            end