zuijinlin.m

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

close all
clear 
clc

 Range=1;%插值邻域范围，1为4邻域
% 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);
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

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

%载入参考条纹位相减少计算时间
[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);
            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]%%%第一步：解出每套条纹测量后的截断相位图
% 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));
        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));
        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);
       
       %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=ProjectorColumn/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))>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    

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
%--------------------------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=x1(k1,j1);
            CountX=x2(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
%          end
    end
end
%-------------处理四领域范围没处理到的projector象素点-----------------
if Range>0
    %--------------
    for k1=1:CCDRow
        for j1=1:CCDColumn
%             if  NN(k1,j1)==1
                CountY=x1(k1,j1);
                CountX=x2(k1,j1);%与CCD象素点对应的DLP亚象素点
                NearestY=round(CountY);
                NearestX=round(CountX);
                for k2=max(1,NearestX-Range):min(ProjectorRow,NearestX+Range)
                    for j2=max(1,NearestY-Range):min(ProjectorColumn,NearestY+Range)
                        if xc(k2,j2)==0
                            DistanceY=Range+1.1-abs(j2-CountY);
                            DistanceX=Range+1.1-abs(k2-CountX);
                            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
                end
            end
%         end
    end
    %--------
end
%--------处理奇异点-------------
for k1=1:ProjectorRow
    for j1=1:ProjectorColumn
        if xf2(k1,j1)~=0
            if xf2(k1,j1)>CCDColumn
                xf2(k1,j1)=CCDColumn;
            elseif xf2(k1,j1)<1
                xf2(k1,j1)=1;
            end
        end
        if xf1(k1,j1)~=0
            if xf1(k1,j1)>CCDRow
                xf1(k1,j1)=CCDRow;
            elseif xf1(k1,j1)<1
                xf1(k1,j1)=1;
            end
        end
    end
end
%-------------------
xf1Closest=xf1;
xf2Closest=xf2;
xcClosest=xc;
%--------------------------------------------------------------------------
%--------------------------产生满调制度反向条纹----------------------------
%-------------------------------------------------------------------------
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*xf2(k1,j1)/p+(i-1)*pi/2))/255;
            end
        end
    end
    File=fullfile(PathNameInv,strcat('A_ClosestInverseL',num2str(i),'.bmp'));
    imwrite(x4l,File);
end