interpcomparesimul.m

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

%--------------------------------------------------------------------------
%--------本程序为反向条纹产生中确立坐标几何反向传递关系的插值方法比较----------
%---------------包括最近邻插值,线性插值以及二元三次插值-----------------------
%-----------------------------计算机模拟------------------------------------
%--------------------------------------------------------------------------
close all
clear 
clc
p=8;q=8;
CCDStartR=1;CCDStartC=1;%CCD相位展开起始点
ProjectorStartR=1;ProjectorStartC=1;%Projector相位展开起始点
CCDRow=256;CCDColumn=256;
ProjectorRow=256;ProjectorColumn=256;
NN=ones(CCDRow,CCDColumn);%有效投影区域
Range=1;%插值邻域范围，1为4邻域
m=peaks(CCDRow);
m=flipud(m);
x=ones(CCDRow,CCDColumn);
y=ones(CCDRow,CCDColumn);%p是光栅周期，q是等效波长，m是代测面形
for n=1:CCDRow
    x(:,n)=n;
end

for n=1:CCDColumn
    y(n,:)=n;
end
%------
backm=ones(CCDRow,CCDColumn)+127;     
for n=1:CCDColumn
    backm(:,n)=backm(:,n)+n/2;
end
%------
modm=(m-min(min(m)))./max(max((m-min(min(m)))));
Mod=0.5*((384-backm)/256).*modm;
Back=0.5*((384-backm)/256);
clear modm backm;
%------
r1=0.5+0.5*cos(2*pi*x/p);
r2=0.5+0.5*cos(2*pi*x/p+pi/2);
r3=0.5+0.5*cos(2*pi*x/p+pi);
r4=0.5+0.5*cos(2*pi*x/p+3*pi/2);
r5=0.5+0.5*cos(2*pi*x/p+2*pi);%参考条纹

i1=Back+Mod.*cos(2*pi*x/p-2*pi*m/q);
i2=Back+Mod.*cos(2*pi*x/p-2*pi*m/q+pi/2);
i3=Back+Mod.*cos(2*pi*x/p-2*pi*m/q+pi);
i4=Back+Mod.*cos(2*pi*x/p-2*pi*m/q+3*pi/2);
i5=Back+Mod.*cos(2*pi*x/p-2*pi*m/q+2*pi);%变形条纹
%------------------
WrappedPhaseVertical=atan2(2*(i4-i2),(i1-2*i3+i5));
RWrappedPhaseVertical=atan2(2*(r4-r2),(r1-2*r3+r5));
UnwrappedPhaseVertical=myunwrappro(WrappedPhaseVertical,1,CCDStartR,CCDStartC);
RUnwrappedPhaseVertical=myunwrappro(RWrappedPhaseVertical,1,ProjectorStartR,ProjectorStartC);
clear WrappedPhaseVertical RWrappedPhaseVertical;
x1=zeros(CCDRow,CCDColumn);%竖向,即对应水平条纹
x2=zeros(CCDRow,CCDColumn);%横向,即对应垂直条纹
 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
%--------------水平条纹-----------------
r1=0.5+0.5*cos(2*pi*y/p);
r2=0.5+0.5*cos(2*pi*y/p+pi/2);
r3=0.5+0.5*cos(2*pi*y/p+pi);
r4=0.5+0.5*cos(2*pi*y/p+3*pi/2);
r5=0.5+0.5*cos(2*pi*y/p+2*pi);

i1=r1;i2=r2;i3=r3;i4=r4;
%--------
WrappedPhaseParallel=atan2(2*(i4-i2),(i1-2*i3+i5));
RWrappedPhaseParallel=atan2(2*(r4-r2),(r1-2*r3+r5));
UnwrappedPhaseParallel=myunwrappro(WrappedPhaseParallel,1,CCDStartR,CCDStartC);
RUnwrappedPhaseParallel=myunwrappro(RWrappedPhaseParallel,1,ProjectorStartR,ProjectorStartC);
  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
           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
clear r1 r2 r3 r4 r5  i2 i3 i4 i5 WrappedPhaseVertical WrappedPhaseParallel;
%------------------------A_ClosestInverseL----------------------
%-------------Projector象素对应的CCD亚像素点计算---------------------------
%---------------------即反向坐标关系的建立---------------------------------
%-----------------最近邻projector座标平面插值------------------------------
%--------------------------------------------------------------------------

%----------------------四领域坐标线性插值-------------------
xb=zeros(ProjectorRow,ProjectorColumn);%各点权重和
xc=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)
                    Distance=1.1-abs(j2-CountY);
                    if Distance>xb(k2,j2)
                        xb(k2,j2)=Distance;
                        xf2(k2,j2)=j1;
                        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
                            Distance=Range+1.1-abs(j2-CountY);
                            if Distance>xb(k2,j2)
                                xb(k2,j2)=Distance;
                                xf2(k2,j2)=j1;
                                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
    end
end
%-------------------产生反向条纹,计算反向条纹投影结果误差-------------------
InverseImageCloset=0.5+0.5*cos(2*pi*xf2/p);%反向条纹
figure,idisp(InverseImageCloset);
InverseAfterProjectedPhase=zeros(CCDRow,CCDColumn);%投影结果相位
for k1=1:CCDRow-1
    for j1=1:CCDColumn-1
        k2=x1(k1,j1);
        j2=x2(k1,j1);
        CountX=floor(k2);
        CountY=floor(j2);
        if (CountX*CountY)~=0
            InverseAfterProjectedPhase(k1,j1)=(xf2(CountX,CountY)*(1-abs(CountY-j2))+...
                xf2(CountX,CountY+1)*(1-abs(CountY+1-j2)))/...
                ((1-abs(CountY+1-j2))+(1-abs(CountY-j2)));
        end
    end
end

ImageAfterProjectedCloset=Back+Mod.*cos(2*pi*InverseAfterProjectedPhase/p);
%反向条纹投影结果
figure,idisp((ImageAfterProjectedCloset));
DifferenceShow=abs(InverseAfterProjectedPhase(3:253,3:253)-x(3:253,3:253));
  for i=1:CCDRow-6
     for j=1:CCDColumn-6
         if DifferenceShow(i,j)>1
             DifferenceShow(i,j)= DifferenceShow(i,j)/10;
         end
     end
 end

 for i=4:CCDRow-4
     for j=4:CCDColumn-4
         DifferenceJ(i,j)=((2*pi*InverseAfterProjectedPhase(i,j)/p-2*pi*x(i,j)/p)^2);  
     end
 end
figure,idisp(DifferenceShow);
DifferenceC= DifferenceJ(find(~isnan(DifferenceJ)));
[NumberDifference,temp]=size(DifferenceC);
sum(sqrt(DifferenceC))/NumberDifference
%平均误差
Difference=sqrt(sum(DifferenceC)/(NumberDifference-1))
%相位标准差

%---------------------------------------------------------------------------------------------B_BilinearInverseL----------------------
%-------------Projector象素对应的CCD亚像素点计算---------------------------
%---------------------即反向坐标关系的建立---------------------------------
%-----------------双线性projector座标平面插值------------------------------
%--------------------------------------------------------------------------
xb=zeros(ProjectorRow,ProjectorColumn);%各点权重和
xa=zeros(ProjectorRow,ProjectorColumn);
xc=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)
                    xc(k2,j2)=1;
                    xb(k2,j2)=xb(k2,j2)+1.1-abs(j2-CountY);%计算每个点权重之和
                    xa(k2,j2)=xa(k2,j2)+1.1-abs(k2-CountX);
                end
            end
        end
    end
end

xf2=zeros(ProjectorRow,ProjectorColumn);%列坐标的反向传递关系
xf1=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)
                    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


%-------------处理四领域范围没处理到的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
                            xb(k2,j2)=xb(k2,j2)+Range+1.1-abs(j2-CountY);
                            xa(k2,j2)=xa(k2,j2)+Range+1.1-abs(k2-CountX);%求权重
                        end
                    end
                end
            end
        end
    end

    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
                            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
%--------处理奇异点-------------
for k1=1:ProjectorRow
    for j1=1:ProjectorColumn
        if xf2(k1,j1)>CCDColumn
            xf2(k1,j1)=CCDColumn;
        elseif (xf2(k1,j1)<1)&&(xf2(k1,j1)>0)
            xf2(k1,j1)=1;
        end

        if xf1(k1,j1)>CCDRow
            xf2(k1,j1)=CCDRow;
        elseif (xf2(k1,j1)<1)&&(xf2(k1,j1)>0)
            xf2(k1,j1)=1;
        end
    end
end
%----------------------反向条纹的产生和投影结果相位差的计算-------------------
InverseImageLinear=0.5+0.5*cos(2*pi*xf2/p);%反向条纹
figure,idisp(InverseImageLinear);
InverseAfterProjectedPhase=zeros(CCDRow,CCDColumn);%反向条纹投影结果相位
for k1=1:CCDRow
     InverseAfterProjectedPhase(k1,:)=interp1(x(k1,:),xf2(k1,:),x2(k1,:),'cubic');
end


ImageAfterProjectedLinear=Back+Mod.*cos(2*pi*InverseAfterProjectedPhase/p);
%反向条纹投影结果
 figure,idisp((ImageAfterProjectedLinear));
 DifferenceShow=abs(InverseAfterProjectedPhase(3:253,3:253)-x(3:253,3:253));
  for i=1:CCDRow-6
     for j=1:CCDColumn-6
         if DifferenceShow(i,j)>1
             DifferenceShow(i,j)= DifferenceShow(i,j)/10;
         end
     end
 end

 for i=4:CCDRow-4
     for j=4:CCDColumn-4
         DifferenceJ(i,j)=((2*pi*InverseAfterProjectedPhase(i,j)/p-2*pi*x(i,j)/p)^2);  
     end
 end
figure,idisp(DifferenceShow);
DifferenceC= DifferenceJ(find(~isnan(DifferenceJ)));
[NumberDifference,temp]=size(DifferenceC);
max(max(sqrt(DifferenceC)))
sum(sqrt(DifferenceC))/NumberDifference
%平均误差
Difference=sqrt(sum(DifferenceC)/(NumberDifference-1))
%相位标准差
%--------------------------


%---------------------------------------------------------------------------------------------C_CubicInverseL----------------------
%-------------Projector象素对应的CCD亚像素点计算---------------------------
%---------------------即反向坐标关系的建立---------------------------------
%-----------------二元三次projector座标平面插值------------------------------
%--------------------------------------------------------------------------
xf2=zeros(CCDRow,CCDColumn);%列坐标的反向传递关系
xf1=zeros(CCDRow,CCDColumn);%行坐标的反向传递关系
[x1CCD,x2CCD]=meshgrid(1:CCDColumn,1:CCDRow);
InverseAfterProjectedPhase=zeros(CCDRow,CCDColumn);%反向条纹投影结果
xf2=griddata(x2,x1,x,x1CCD,x2CCD,'cubic');
xf1=griddata(x2,x1,y,x1CCD,x2CCD,'cubic');%二元三次插值
InverseImageCubic=0.5+0.5*cos(2*pi*xf2/p);%反向条纹
figure,idisp(InverseImageCubic);

InverseAfterProjectedPhase=griddata(x,y,xf2,x2,x1,'cubic');

ImageAfterProjectedCubic=Back+Mod.*cos(2*pi*InverseAfterProjectedPhase/p);
%反向条纹投影结果
 figure,idisp((ImageAfterProjectedCubic));
 DifferenceShow=abs(InverseAfterProjectedPhase(3:253,3:253)-x(3:253,3:253));
  for i=1:CCDRow-6
     for j=1:CCDColumn-6
         if DifferenceShow(i,j)>0.01
             DifferenceShow(i,j)= DifferenceShow(i,j)/100;
         end
     end
 end

 for i=4:CCDRow-4
     for j=4:CCDColumn-4
         DifferenceJ(i,j)=((2*pi*InverseAfterProjectedPhase(i,j)/p-2*pi*x(i,j)/p)^2);  
     end
 end
figure,idisp(DifferenceShow);
DifferenceC= DifferenceJ(find(~isnan(DifferenceJ)));
[NumberDifference,temp]=size(DifferenceC);
max(max(sqrt(DifferenceC)))
sum(sqrt(DifferenceC))/NumberDifference
%平均误差
Difference=sqrt(sum(DifferenceC)/(NumberDifference-1))
%相位标准差