test1.m

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

%          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);
                            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
%----------------
%---------------------------------------------------------------------------------------------B_BilinearInverseL----------------------
%-------------Projector象素对应的CCD亚像素点计算---------------------------
%---------------------即反向坐标关系的建立---------------------------------
%-----------------双线性projector座标平面插值------------------------------
%--------------------------------------------------------------------------
%----------------------四领域坐标线性插值-------------------
xa=zeros(ProjectorRow,ProjectorColumn);
xb=zeros(ProjectorRow,ProjectorColumn);%各点权重和
xc=zeros(ProjectorRow,ProjectorColumn);%标志位
%处理遗漏点的标志位
for k1=1:CCDRow-1
    for j1=1:CCDColumn-1
        if  (NN(k1,j1)==1)&&(NN(k1+1,j1)==1)&&(NN(k1,j1+1)==1)&&(abs(x2(k1,j1)-x2(k1+1,j1))<10)&&(abs(x2(k1,j1)-x2(k1,j1+1))<10)...
                &&(abs(x2(k1,j1+1)-x2(k1+1,j1))<10) && (abs(x1(k1,j1)-x1(k1+1,j1))<10)&&(abs(x1(k1,j1)-x1(k1,j1+1))<10)...
                &&(abs(x1(k1,j1+1)-x1(k1+1,j1))<10)
            CountY=(x2(k1,j1)+x2(k1+1,j1)+x2(k1,j1+1))/3;
            CountX=(x1(k1,j1)+x1(k1+1,j1)+x1(k1,j1+1))/3;
            LowestY=floor(CountY);
            LowestX=floor(CountX);
            for k2=max(1,LowestX):min(ProjectorRow,LowestX+1)
                for j2=max(1,LowestY):min(ProjectorColumn,LowestY+1)
                    xa(k2,j2)=xa(k2,j2)+1.1-abs(k2-CountX);
                    xb(k2,j2)=xb(k2,j2)+1.1-abs(j2-CountY);%计算每个点权重之和
                    xc(k2,j2)=1;
                end
            end
        end
    end
end

xf1=zeros(ProjectorRow,ProjectorColumn);
xf2=zeros(ProjectorRow,ProjectorColumn);

for k1=1:CCDRow-1
    for j1=1:CCDColumn-1
        if  (NN(k1,j1)==1)&&(NN(k1+1,j1)==1)&&(NN(k1,j1+1)==1)&&(abs(x2(k1,j1)-x2(k1+1,j1))<10)&&(abs(x2(k1,j1)-x2(k1,j1+1))<10)...
                &&(abs(x2(k1,j1+1)-x2(k1+1,j1))<10) && (abs(x1(k1,j1)-x1(k1+1,j1))<10)&&(abs(x1(k1,j1)-x1(k1,j1+1))<10)...
                &&(abs(x1(k1,j1+1)-x1(k1+1,j1))<10)
            CountY=(x2(k1,j1)+x2(k1+1,j1)+x2(k1,j1+1))/3;
            CountX=(x1(k1,j1)+x1(k1+1,j1)+x1(k1,j1+1))/3;
            LowestY=floor(CountY);
            LowestX=floor(CountX);
            A1=[j1+1-j1 0;0 k1+1-k1];
            A2=[-x1(k1,j1)+x1(k1+1,j1) x2(k1,j1)-x2(k1+1,j1);-x1(k1,j1+1)+x1(k1,j1) x2(k1,j1+1)-x2(k1,j1)];
            A3=(x2(k1,j1)*(-x1(k1+1,j1)+x1(k1,j1+1))+x2(k1,j1+1)*(-x1(k1,j1)+x1(k1+1,j1))+x2(k1+1,j1)*(-x1(k1,j1+1)+x1(k1,j1)));
            if A3~=0
                A=1/A3*A1*A2;
                C=[j1;k1]-A*[x2(k1,j1);x1(k1,j1)];
                for k2=max(1,LowestX):min(ProjectorRow,LowestX+1)
                    for j2=max(1,LowestY):min(ProjectorColumn,LowestY+1)
                        B=A*[j2;k2]+C;
                        xf2(k2,j2)=B(1,1)*(1.1-abs(j2-CountY))/xb(k2,j2)+xf2(k2,j2);
                        xf1(k2,j2)=B(2,1)*(1.1-abs(k2-CountX))/xa(k2,j2)+xf1(k2,j2);
                    end
                end
            else
                for k2=max(1,LowestX):min(ProjectorRow,LowestX+1)
                    for j2=max(1,LowestY):min(ProjectorColumn,LowestY+1)
                        xf2(k2,j2)=j1*(1.1-abs(j2-CountY))/xb(k2,j2)+xf2(k2,j2);
                        xf1(k2,j2)=k1*(1.1-abs(k2-CountX))/xa(k2,j2)+xf1(k2,j2);
                    end
                end
            end
        end
    end
end


%-------------处理四领域范围没处理到的projector象素点-----------------
if Range>0
    %--------------
    for k1=1:CCDRow-1
        for j1=1:CCDColumn-1
            if  (NN(k1,j1)==1)&&(NN(k1+1,j1)==1)&&(NN(k1,j1+1)==1)&&(abs(x2(k1,j1)-x2(k1+1,j1))<10)&&(abs(x2(k1,j1)-x2(k1,j1+1))<10)...
                    &&(abs(x2(k1,j1+1)-x2(k1+1,j1))<10) && (abs(x1(k1,j1)-x1(k1+1,j1))<10)&&(abs(x1(k1,j1)-x1(k1,j1+1))<10)...
                    &&(abs(x1(k1,j1+1)-x1(k1+1,j1))<10)
                CountY=(x2(k1,j1)+x2(k1+1,j1)+x2(k1,j1+1))/3;
                CountX=(x1(k1,j1)+x1(k1+1,j1)+x1(k1,j1+1))/3;
                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
                            xa(k2,j2)=xa(k2,j2)+Range+1.1-abs(k2-CountX);
                            xb(k2,j2)=xb(k2,j2)+Range+1.1-abs(j2-CountY);
                        end
                    end
                end
            end
        end
    end

    for k1=1:CCDRow-1
        for j1=1:CCDColumn-1
            if  (NN(k1,j1)==1)&&(NN(k1+1,j1)==1)&&(NN(k1,j1+1)==1)&&(abs(x2(k1,j1)-x2(k1+1,j1))<10)&&(abs(x2(k1,j1)-x2(k1,j1+1))<10)...
                    &&(abs(x2(k1,j1+1)-x2(k1+1,j1))<10) && (abs(x1(k1,j1)-x1(k1+1,j1))<10)&&(abs(x1(k1,j1)-x1(k1,j1+1))<10)...
                    &&(abs(x1(k1,j1+1)-x1(k1+1,j1))<10)
                CountY=(x2(k1,j1)+x2(k1+1,j1)+x2(k1,j1+1))/3;
                CountX=(x1(k1,j1)+x1(k1+1,j1)+x1(k1,j1+1))/3;
                LowestY=floor(CountY);
                LowestX=floor(CountX);
                A1=[j1+1-j1 0;0 k1+1-k1];
                A2=[-x1(k1,j1)+x1(k1+1,j1) x2(k1,j1)-x2(k1+1,j1);-x1(k1,j1+1)+x1(k1,j1) x2(k1,j1+1)-x2(k1,j1)];
                A3=(x2(k1,j1)*(-x1(k1+1,j1)+x1(k1,j1+1))+x2(k1,j1+1)*(-x1(k1,j1)+x1(k1+1,j1))+x2(k1+1,j1)*(-x1(k1,j1+1)+x1(k1,j1)));
                if A3~=0
                    A=1/A3*A1*A2;
                    C=[j1;k1]-A*[x2(k1,j1);x1(k1,j1)];
                    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  %&&(xb(k2,j2)~=0)
                                B=A*[j2;k2]+C;
                                xf2(k2,j2)=B(1,1)*(Range+1.1-abs(j2-CountY))/xb(k2,j2)+xf2(k2,j2);
                                xf1(k2,j2)=B(2,1)*(Range+1.1-abs(k2-CountX))/xa(k2,j2)+xf1(k2,j2);
                            end
                        end
                    end
                else
                    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  %&&(xb(k2,j2)~=0)

                                xf2(k2,j2)=j1*(Range+1.1-abs(j2-CountY))/xb(k2,j2)+xf2(k2,j2);
                                xf1(k2,j2)=k1*(Range+1.1-abs(k2-CountX))/xa(k2,j2)+xf1(k2,j2);
                            end
                        end
                    end
                end
            end
        end
    end
    %--------
end
%--------处理奇异点-------------
for k1=1:ProjectorRow
    for j1=1:ProjectorColumn
        if xf1(k1,j1)>CCDRow
            xf1(k1,j1)=CCDRow;
        elseif (xf1(k1,j1)<1)&&(xf1(k1,j1)>0)
            xf1(k1,j1)=1;
        else
            if (xf1(k1,j1)==0)&&(xf1Closest(k1,j1)~=0)
                xf1(k1,j1)=xf1Closest(k1,j1);
            end
        end

        if xf2(k1,j1)>CCDColumn
            xf2(k1,j1)=CCDColumn;
        elseif (xf2(k1,j1)<1)&&(xf2(k1,j1)>0)
            xf2(k1,j1)=1;
        else
            if (xf2(k1,j1)==0)&&(xf2Closest(k1,j1)~=0)
                xf2(k1,j1)=xf2Closest(k1,j1);
            end
        end
    end
end
%------------------
%--------------------------------------------------------------------------
%--------------------------产生满调制度反向条纹----------------------------
%-------------------------------------------------------------------------
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('B_BilinearInverseL',num2str(i),'.bmp'));
    imwrite(x4l,File);
end
%---------------------------------------------------------------------------------------------C_CubicInverseL----------------------
%-------------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);
InverseAfterProjectedPhase=zeros(CCDRow,CCDColumn);
xf2=griddata(x2,x1,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*xf2(k1,j1)/p+(i-1)*pi/2))/255;
            end
        end
    end
    File=fullfile(PathNameInv,strcat('C_CubicInverseL',num2str(i),'.bmp'));
    imwrite(x4l,File);
end