yuanshi.m

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

%--------------------------------------------------------------------------
%--------本程序为反向条纹产生中确立坐标几何反向传递关系的插值方法比较----------
%---------------包括最近邻插值,线性插值以及二元三次插值-----------------------
%--------------------------------------------------------------------------
%插值邻域范围，1为4邻域
Range=1;
%有效投影区的阈值
ThreshLCD=30;
%--------------------------------------------------------------------------
%------------------CCD象素对应的Projector亚像素点计算----------------------
%--------------------------------------------------------------------------

%---------||读入投影仪180背景图，计算投影区域||----------

%-------------LCD的180背景--------------
[FileNameLCD,PathNameLCD]=uigetfile({'*.bmp','Bitmap files(*.bmp)';...
        '*.*','All Files(*.*)'},'get the LCD Background images');
if ~isequal([FileNameLCD,PathNameLCD],[0,0])
        File=fullfile(PathNameLCD,strcat(num2str(180),'.bmp'));
        LCDBack=imread(File);
 end
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];
%水平变形条纹
[FileNameParal,PathNameParal]=uigetfile({'*.bmp','Bitmap files(*.bmp)';...
        '*.*','All Files(*.*)'},'get the LCD Parallel Deformed fringes');
%载入参考条纹位相减少计算时间
[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)),'TPU',num2str(j),'.bmp'));
            I(:,:,j,i)=imread(File);
        end
    end
end
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 = 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 M N  Pufit ws;

%------------载入参考条纹位相减少计算时间-----------
if ~isequal([FileNameOrig,PathNameOrig],[0,0])
    File=fullfile(PathNameOrig,FileNameOrig);
    load(File);
end
clear PathNameOrig,FileNameOrig;
[ProjectorRow,ProjectorColumn]=size(RUnwrappedPhaseVertical);
%------------------竖向对应点计算:摄像机像素映射到投影器坐标系-------------------
p=ProjectorRow/32;
RUnwrappedPhaseParallel=RUnwrappedPhaseParallel-32*pi;
UnwrappedPhaseParallel=UnwrappedPhaseParallel-RUnwrappedPhaseParallel(1,1)+2*pi/p;
RUnwrappedPhaseParallel=RUnwrappedPhaseParallel-RUnwrappedPhaseParallel(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)=UnwrappedPhaseParallel(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,'TPU','TPUV'),strcat(num2str(FringeNumber(i)),'TPU',num2str(j),'.bmp'));
            I(:,:,j,i)=imread(File);
        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;
UnwrappedPhaseVertical=Pufit;

clear Pu32 Pu16 Pu24 Pu30 Pu28 Pu31 Pufit Pufit ws;

%-----------------横向对应点计算:摄像机像素映射到投影器坐标系--------------------
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
           x2(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(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
        end
    end
end
%-------------处理四领域范围没处理到的projector象素点-----------------
if Range>0
    %--------------
    for k1=1:CCDRow
        for j1=1:CCDColumn
            if  NN(k1,j1)==1
                CountY=x2(k1,j1);
                CountX=x1(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);