二维图像的离轴数字全息及显示

二维图像的离轴数字全息及重建

本文主要参考<衍射计算及数字全息>${}^{[1-3]}$一书中相关内容，若有不足的地方，恳请批评指正。
全息摄像主要经历光学全息，数字全息，和计算全息这样几个过程，其基本原理都是通过干涉纪录及衍射重建。当两束光波进行干涉时，纪录物体的振幅和相位信息，通过衍射计算，重构物体的振幅和相位信息。原则上，全息面任一一个点，或者一小块都包含原物体的所有信息，因而，当全息图损坏时，可以通过一小块物体重构原物体。若已知在物平面$y_0$，实际模拟运算时，可另物平面为一幅图像，物平面经过衍射距离$d$后，在衍射观测面$y_1$的复振幅分布可以以下公式获得。
$$\begin{gathered} U(p\Delta x,q\Delta y)=\frac{\exp (jkd)}{j\lambda d}\exp [\frac{jk}{2d}(p\Delta x^2)+q\Delta y^2] \\ \times ϝ[U_0(m\Delta x_0,n\Delta y_0)\times \exp (\frac{jk}{2d}(m\Delta x_0{)}^2,(n\Delta y_0{)}^2){]}_{\frac{p\Delta x}{\lambda d},\frac{p\Delta y}{\lambda d}} \end{gathered}$$
具体实现公式1的过程如下：
(1)物平面函数$U_0(m\Delta x_0,n\Delta y_0)与二次相位因子\exp (\frac{jk}{2d}(m\Delta x_0{)}^2,(n\Delta y_0{)}^2)]$相乘，对相乘的结果做傅里叶变换；
(2)对上述的结果与$U(p\Delta x,q\Delta y)=\frac{\exp (jkd)}{j\lambda d}\exp [\frac{jk}{2d}(p\Delta x^2)+q\Delta y^2]$相乘，即可得到衍射观测面的光波厂复振幅。因为整个计算过程只进行了一次快速傅里叶变换，因而被称为一次傅里叶变换(Singe-FFT)
物面观测面衍射尺寸分析：
之所以要分析物面尺寸关系，是因为不论是数字全息还是计算全息，纪录的相机如CCD，重建的显示设备如空间光调制器LCOS，其靶面尺寸都是有限的。因CCD或者LCOS的靶面尺寸是固定的，因而通过该尺寸反推物面尺寸。选用的CCD的像素尺寸为0.00465mm，像素数为1024则，靶面尺寸为$L$ =0.00465×1024；采样数为$N$ , 根据S-FFT物平面$L_0$和衍射观测面(观察面)的关系，可以
$$L_0=\frac{N\lambda d}{L}$$
在衍射重构时，全息面作为物平面； MATLAB模拟程序如下

%           h ——波长（mm）；          I ——数字全息图；
%           L ——全息图宽度（mm）；    z0——记录全息图的距离(mm)；
XRGB=imread('baboon.jpg');
X0=rgb2gray(XRGB);
figure,imshow(X0,[]);
[M0,N0]=size(X0);
N1=min(M0,N0);
N=1024;          %模拟形成的全息图取样数
X1=imresize(X0,N/4/N1);
[M1,N1]=size(X1);
X=zeros(N,N);
X(N/2-M1/2+1:N/2+M1/2,N/2-N1/2+1:N/2+N1/2)=X1(1:M1,1:N1);
h=0.632e-3;      %波长(mm), 可按需要修改
k=2*pi/h;
pix=0.00465     ;%CCD像素宽度(mm), 可按需要修改
L=N*pix;         %CCD宽度(mm)
z0=1000;         %衍射距离(mm), 可按需要修改
L0=h*N*z0/L;     %物平面宽度(mm)
Y=double(X);
a=ones(N,N);
b=rand(N,N)*2*pi;
f=Y.*exp(i.*b);  %叠加随机相位噪声,形成振幅正比于图像的初始场复振幅

figstr=strcat('初始物平面宽度=',num2str(L0),'mm');
figure,imshow(X,[]),colormap(gray); xlabel(figstr);title('物平面图像');
%---------------菲涅耳衍射的S-FFT计算开始
n=1:N;
x=-L0/2+L0/N*(n-1);	 					
y=x;
[yy,xx] = meshgrid(y,x); %网格化数据平面
Fresnel=exp(i*k/2/z0*(xx.^2+yy.^2));
f2=f.*Fresnel;
Uf=fft2(f2,N,N);
Uf=fftshift(Uf);
x=-L/2+L/N*(n-1);%CCD宽度取样(mm) 					
y=x;
[yy,xx] = meshgrid(y,x); 
phase=exp(i*k*z0)/(i*h*z0)*exp(i*k/2/z0*(xx.^2+yy.^2));%菲涅耳衍射积分前方的相位因子
Uf=Uf.*phase;
%---------------S-FFT计算结束
figstr=strcat('模拟CCD宽度=',num2str(L),'mm');
figure,imshow(abs(Uf),[]),colormap(gray); xlabel(figstr);title('到达CCD平面的物光振幅分布');
%---------------形成0-255灰度级的数字全息图
fex=N/L;
Qx=(4-2.5)*L0/8/z0;          %按照优化设计定义参考光方向余弦
Qy=Qx
x=[-L/2:L/N:L/2-L/N];
y=x;
[X,Y]=meshgrid(x,y);
Ar=max(max(abs(Uf)));        %按物光场振幅最大值定义参考光振幅
Ur=Ar*exp(i*k*(X.*Qx+Y.*Qy));%参考光复振幅
Uh=Ur+Uf;                    %物光与参考光干涉
Wh=Uh.*conj(Uh);             %干涉场强度
Imax=max(max(Wh));
I=uint8(Wh./Imax*255);      %形成0-255灰度级的数字全息图
imwrite(I,'D:\Ih.tif');     %形成数字全息图文件
figstr=strcat('全息图宽度=',num2str(L),'mm');
figure,imshow(Ih,[]),colormap(gray);xlabel(figstr);title('模拟形成的数字全息图');

%基于S-FFT的全息重建
f0 = double(Ih);
[N1,N2]=size(f0);
N=min(N1,N2);                      
h=0.000632;             %波长(mm)              
z0=1000;
L=N*pix;                %CCD宽度(mm)                        
In(1:N,1:N)=f0(1:N,1:N);
%-----------------------------1-FFT重建开始
n=1:N;
x=-L/2+L/N*(n-1);	 					
y=x;
[yy,xx] = meshgrid(y,x); 
k=2*pi/h;	
Fresnel=exp(i*k/2/z0*(xx.^2+yy.^2));
f2=In.*Fresnel;
Uf=fft2(f2,N,N);
Uf=fftshift(Uf);
L0=h*z0*N/L;
x=-L0/2+L0/N*(n-1);	 					
y=x;
[yy,xx] = meshgrid(y,x); 
phase=exp(i*k*z0)/(i*h*z0)*exp(i*k/2/z0*(xx.^2+yy.^2));
U0=Uf.*phase;%积分运算结果乘积分号前方相位因子
%-----------------------------1-FFT重建结束
If=U0*conj(U0);
Gmax=max(max(abs(U0)));
Gmin=min(min(abs(U0)));
figstr=strcat('重建物平面宽度=',num2str(L0),'mm');
figure(2),
imshow(abs(U0),[Gmin Gmax/1]),colormap(gray); xlabel(figstr);title('1-FFT物平面重建图像');
p=10;
while p
    figure(3);
    imshow(abs(U0),[Gmin Gmax/p]),colormap(gray);xlabel(figstr);title('1-FFT物平面重建图像');
    p=input('Gmax/p,p=10?');
end;
Up=fft2(U0);
figure,imshow(abs(Up),[])
figure,plot(abs(Up(round(512)+1,:)));
hold on;

[1] 古德曼 ，傅里叶光学导论[M]，北京：电子工业出版社，2006，第3版
[2]李俊昌，衍射计算及数字全息[M]，北京：科学出版社，2014
[3]陈家壁，苏显渝.光学信息技术原理及应用[M].北京：高等教育出版社，2001