计算全息图的制作及数字再现

计算全息图的制作及其数字再现
物理科学与工程技术学院作者姓名：杨煦、杨康明
指导老师：蔡志岗教授
摘要：计算机制全息图是制作全息图的一种新技术，它是利用数字计算机来综合的全息图，它不需要物体的实际存在，而是把物波的数学描述输入计算机处理后，控制绘图仪输出或显示器显示二制成的全息图。

计算全息图的数字再现是利用计算机模拟光学全息的光路，仿真菲涅尔衍射、透镜傅里叶变换等光学过程从而在虚拟的观察屏上得到全息再现像。

关键词：计算全息数字再现
一、引言：
早在1965年，Kozman和Kelly就提出了计算机生成全息图（Computer Generated Holography，简称CGH）的概念，那时受计算机速度、容量和显示器分辨率等因素的约束，直到80年代中期以前计算机全息图的研究一直未取得大的进展。

国内对全息技术的研究主要集中在物理光学领域。

而目前由于计算机技术的发展以及计算机硬件的进步，已经可以制作空间带宽积很大的计算全息图，但是由于输出设备的精度问题，难以制作质量很高的全息图。

因此我们将以此为研究重点，希望从编码方法上有所突破，解决这个问题。

二、实验原理
计算全息图的制作和再现过程主要分为以下几个步骤：
1、抽样，得到物体或波面在离散样点上的值；
2、计算，计算物光波在全息平面上的光场分布；
3、编码，把全息平面上光波的复振幅分布编码成为全息图的透过率变化；
4、成图，在计算机控制下，将全息图的透过率变化绘制成图，如果绘图设备分辨率不够，则绘制一个较大的图，再缩版到得到使用的全息图；
5、再现，这一步骤与光学全息图的再现没有什么区别。

制作一个傅立叶变换全息图的典型流程如下：
（一）、抽样
抽样包括对输入图像的抽样和对全息图的抽样。

实际上，输入图像和全息图像的信号都是连续的。

而计算机只能对离散的数据进行处理，所以必须对物光和全息图像进行离散化，即抽样处理。

由空间带宽积的传递不变性可以知道，在全息图平面上的空间带宽积SW 应该和物体的空间带宽积SW 相等。

换句话说，全息图上的抽样数(可分辨单元数)要等于(或大于)输入图像的抽样数。

于是我们可方便地确定在计算机全息图制作时，输入图像和全息图平面上应该具
有的抽样点数总数，即2()N x v =∆∆。

其中2()x ∆是输入图像的面积，2()v ∆是其频带面
积。

另外，给出一幅图像，如何知道其带宽呢？由于成像系统是低通滤波器，其传递图像的最高频率是有限的，存在着截止频率，故得到的图像可近似地看作限带的，输入图像的带宽可近似地用光学通道(包括接收器)的截止频率来表示。

还有一种方法是用空间局部频率的概念来定义带宽：
由空间局部频率出发，可以来定义物波函数的带宽。

定义：物波函数在x 、y 方向的带宽分别为Bx ，By 。

max
2(,)x W x y B x λ⎛∂⎫= ⎪∂⎝⎭ max
2(,)y W x y B y λ⎛⎫∂= ⎪∂⎝⎭ 即用空间局部频率的最大值的两倍来表示物波函数的带宽。

这是一种有效的近似处理。

一般地，以图片格式存储在计算机里面的图像已经是图像抽样之后的离散数据了。

我们只需要知道图像的尺寸和其空间带宽积，这样便可以换算出原图像的抽样间隔和全息图的抽样间隔。

所以以图片格式存储的图像我们一般不必再作抽样处理。

如果图像是以其他连续形
式格式存储的话，就需要对图像进行抽样。

(二)、计算
计算图像的全息图的过程，是运用物理定律求出全息平面上的复振幅分布的过程，运用计算机模拟光的传播、衍射、干涉，得到全息平面上的复振幅分布。

具体计算时，常会用到傅立叶变换、菲涅尔衍射公式等等。

计算机全息中运用的是经典的标量波衍射理论。

这一理论把光波看作标量波来处理，能很好地解释全息过程。

理论里面涉及到的公式一般是连续形式的。

计算机要运用这些理论必须把输入量、输出量和公式三者都离散化。

输入和输出量的离散化由前面的抽样方法实现。

公式的离散化则是数值计算关心的问题。

然而得到离散化后的公式与原来的公式并不完全等价。

用计算机制作的全息图，经常需要用到傅立叶变换算法。

计算机里面常用的傅立叶算法是快速傅立叶变换算法(FFT)。

离散形式的傅立叶变换和连续傅立叶变换并不完全相同。

（三）、编码
由于前面计算得到的只是全息平面上的复振幅分布。

而实际上的全息图一般是透过率的分布(振幅型)、折射率的分布(相位型)和厚度分布(相位型)等等。

这些分布涉及到的参数都是正的实数，所以必须把复数的复振幅变换到对应的正实数参数，这就是编码过程。

在光全息中，用底片记录全息平面的光强就是把复振幅以光强的形式进行编码。

而在计算全息中，编码的方法则很灵活。

下面介绍一种传统的编码方法：罗曼III编码。

1、迂回相位效应
迂回相位效应是不规则光栅的衍射效应。

考虑一个规则的二元光栅，光栅周期为d，如图1(a)。

(a)规则二元光栅(b)某条纹发生移动后的二元光栅
图1 迂回相位效应示意图
对规则光栅某一级衍射相邻两条光线的光程差为。

若光栅的某个条纹发生p的位移，如图1(b)，则该处的光程差变为，比原来增加了光程差。

结合光栅方程，条纹变动处有相位变化。

这个附加的相位便是迂回相位。

可以看到，迂回相位跟波长无光，跟光的入射角出射角也无关，只跟衍射级、光栅常数、附加位移有关。

利用迂回相位效应，控制条纹的附加位移，就可以实现对波面的控制。

2、罗曼III编码的编码方法
罗曼III编码就是应用迂回相位效应来进行编码的。

具体编码过程如下：
通过前面“计算”这一个步骤，可以得到全息平面上的离散的复振幅分布。

现在把每个离散点作为一个全息像元，如图2。

像元的宽和高分别为du和dv(抽样间隔)。

像元内有一宽W 高H的矩形孔。

宽度方向矩形孔的中心偏离像元中心的距离为P，而高度方向两者中心同高。

除了矩形孔内部之外，像元的其它位置都不可以透光。

根据迂回相位效应，把每个矩形孔当作是一个条纹，用矩形孔的面积调制复振幅的绝对值(强度)，矩形孔的位置调制复振幅的复角(相位)，这便是罗曼III编码的方法。

首先求出全息平面上每个离散点的复振幅的绝对值，然后对其归一化得到归一化后的振幅值A，再求出复振幅的复角，接着用这些振幅值A和复角分别对每个像元内的矩形孔的高度和宽度进行调制。

假设要在第M级衍射光中得到最亮的像，就要使，，。

这样操作后，就得到了一张罗曼III编码的全息图。

图2 罗曼III编码示意图
在具体操作的时候，可能会遇到矩形孔超出像元的情况，此时要进行“模式溢出校正”。

“模式溢出校正法”就是把溢出部分的矩形孔循环移动到像元的另一侧(如图3)。

图3 模式溢出校正
这样操作即可以保证每个像元的透光面积保持不变，又可以不影响相位的调制，因为像元左右边界的相位差是。

由于矩形孔衍射的原因，这样编码得到的再现图的亮度并不均匀。

这种不均匀可以通过对输入图像的强度乘以一个sinc因子来改善。

还有一个问题需要注意。

前面迂回相位的推导并不严格。

如图4所示，光栅条纹移动后，实际上C和C’已经不是同一点了。

真正的附加光程并不是C与C’点光线的光程差，而应该是C’点光线与C’点在原来规则光栅同一位置处的光线的光程差。

图4 迂回效应的进一步说明图
也就是说，罗曼III编码中迂回相位不是相对于像元中心处的相位而附加上的量，而应该是相对于矩形孔移动后的中心处位置的相位而附加的量。

(四)、绘制和照相
得到了编码后的全息图像，可以用绘图仪或打印机复制到胶片上显示出来，然后把图像通过微缩照相缩小到合适的尺寸，以达到足够高的分辨率。

或者用高分辨率的空间光调制器直接调制图像。

（五）、再现
通过上面方法得到的计算机全息图，可以用光全息图的再现方法来再现。

除了光学再现外，还有一种再现——数字再现。

这方法是利用计算得到的全息图数据，通过计算机模拟光学再现系统，最后在计算机中计算出再现图像，实现全息图的再现。

（六）、全息图的数字再现
数字再现的实际上就是模拟光在光学系统里传播。

全息图的数字再现过程包括全息图的衍射和透镜的变换等过程。

在许多光学教科书上都可以找到菲涅尔衍射公式 ()221001()2200111111(,)(,)exp ikz i x y z e E x y e F E x y i x y i z z πλπλλ+⎧⎫⎡⎤⎪⎪=+⎨⎬⎢⎥⎪⎪⎣⎦⎩
⎭ 其中，F 代表傅立叶变换，是输入面复振幅，是输出面复振幅，是两平面的距离，k 为光波波矢，λ为光波波长。

还有透镜的透射系数公式
221111exp 2(,)0x y ik f t x y ⎧⎛⎫+-⎪ ⎪=⎨⎝⎭⎪⎩
孔径内孔径外 其中f 为透镜的焦距。

后文涉及到的全息图都是傅立叶变换全息图，所以其再现光学系统就只是一个傅立叶透镜。

把全息图放在透镜的前焦面，在后焦面上便得到其再现图像。

再现过程如图5。

图5 数字再现过程图
把前面编码后的全息图数据作为输入复振幅数据，经过一次菲涅尔衍射运算后，把得到的复振幅分布乘以透镜的透射系数，然后再经过一次菲涅尔衍射运算，在再现屏幕上便得到再现像的复振幅分布，取其模的平方，就可以得到再现图像。

数字再现的程序见附录1。

要强调的是，这里模拟的是一个理想的共轴光学系统，有许多实际因素(如：透镜的像差、离轴情况)并没有考虑。

三、实验制作计算机全息图及其数字再现
罗曼III 编码傅立叶变换全息图的计算的程序用Matlab 来编写的，具体程序代码见附录2。

程序的过程大致是：先读入图像，把图像变为灰度图，并且给图像加上随机相位。

这样做是原因是若直接对图像作傅里叶变换，会使得归一化的振幅变化很快。

读入的图像加上随机相
位处理(作用类似于光学实验中给光源加上毛玻璃)，能使图像经过傅里叶变换后，振幅的变化变缓慢，而再现的时候又不会影响光强的分布。

接着对图像进行二维的快速傅立叶变换。

然后用罗曼III编码算法进行编码。

编码过程只引入了“模式溢出校正”，没有引入其他的校正。

最后把全息图显示在屏幕上，得到的全息图保存为图片格式。

利用计算机模拟光学全息的光路，仿真菲涅尔衍射、透镜傅里叶变换等光学过程从而在虚拟的观察屏上得到全息再现像。

图6 计算全息的数字再现
实验中选用的对象是高宽都为32像素汉字：占。

选取这个字主要是因为他方向性强，仅含有横竖两种笔画，便于观察不同像元大小对全息再现的影响。

以下为调节像元大小为不同值时的全息编码图以及再现图。

图7 du=1，dv=1
图8 du=16，dv=16
图9 du=32，dv=32
图10 du=8，dv=36
附录：
1．全息图的数字再现Matlab程序
文件名：Fresnel.m
% 本程序是输入一幅图像,求z1距离后的面上菲涅尔衍射场的场分布
% 入射光波长lambada,单位m
% E1为输入图像代表的矩阵
% E0为输出图像代表的矩阵
% deltaX1为输入图像横向的宽度,单位m
% deltaY1为输入图像纵向的宽度,单位m
% deltaX0为输出图像横向的宽度,单位m
% deltaY0为输出图像纵向的宽度,单位m
% z1输入图像与衍射平面的距离,单位m
function [E0,deltaX0,deltaY0]=Fresnel(lambada,E1,deltaX1,deltaY1,z1) [ny,nx]=size(E1);
[x0,y0]=meshgrid(-(nx-1)/2:(nx-1)/2,-(ny-1)/2:(ny-1)/2);
dx1=deltaX1/nx;
dy1=deltaY1/ny;
dx0=lambada*z1/deltaX1;
dy0=lambada*z1/deltaY1;
x0=x0*dx0;
y0=y0*dy0;
c=exp(i*2*pi/lambada*z1)/(i*lambada*z1)*exp(i*pi*(x0.^2+y0.^2)/lambada/z1); % 相关的系数，参考菲涅尔衍射公式
for p=1:nx
for q=1:ny
E(q,p)=E1(q,p)*exp(i*pi/lambada/z1*((dx1*(p-(nx+1)/2))^2+(dy1*(q-(ny+1)/2))^2)) ; % 相关的计算，参考菲涅尔衍射公式
end
end
E=fftshift(E); % 把数据作“平移”到第一象限作处理(其实是把
图像作周期延拓，再取第一象限的那个周期)
Etemp=fftshift(fft2(E)); % 把傅立叶变换谱的零频率移动到中心
E0=c.*Etemp; % 参考菲涅尔衍射公式
deltaX0=dx0*nx;
deltaY0=dy0*ny;
文件名：lens.m
% 透镜输入输出面分别紧靠在透镜的两侧，认为输入和输出面之间的距离为零
% 透镜的中心在光轴上，透镜所在平面与光轴垂直
% 入射光波长lambada,单位m
% E1为输入图像代表的矩阵
% E0为输出图像代表的矩阵
% deltaX1为输入图像横向的宽度,单位m
% deltaY1为输入图像纵向的宽度,单位m
% F为透镜的焦距,单位m
% aperture为透镜的孔径，单位m
function E0=lens(lambada,E1,deltaX1,deltaY1,F,aperture)
[ny,nx]=size(E1); % 对一幅图片，纵向对应图片矩阵的行(ny)，横向对应矩阵的列(nx) dx=deltaX1/nx;
dy=deltaY1/ny;
[x,y]=meshgrid(-(nx-1)/2:(nx-1)/2,-(ny-1)/2:(ny-1)/2);
x=x*dx;
y=y*dy;
I=find(x.^2+y.^2>(aperture/2).^2);
x(I)=0;
y(I)=0;
E0=E1.*exp(-i*pi*(x.^2+y.^2)/lambada/F); % 透镜的复透射系数公式
文件名：main.m
clear;
lambada=632.8e-9;
p=imread('j:\hologram datas\pp.bmp'); %全息图所在的路径
p=rgb2gray(p);
%figure;
%imshow(p);
p0=double(p);
deltaX=1e-2; % 全息图宽(m)
deltaY=1e-2; % 全息图高(m)
F=27e-2; % 焦距(m)
aperture=10e-2; % 孔径(m)
obj=F;
img=F;
[p1,U,V]=Fresnel(lambada,p0,deltaX,deltaY,obj); % 第一次菲涅尔衍射
p2=lens(lambada,p1,U,V,F,aperture); % 经过透镜
[p3,U,V]=Fresnel(lambada,p2,U,V,img); % 第二次菲涅尔衍射
pout=p3;
pout=abs(pout);
pout=pout/max(max(pout));
figure
imshow(pout*100);
axis on
xlabel(U);
ylabel(V);
2．罗曼III编码制作二维傅立叶全息图的matlab程序
clear;
close;
%%%%%%%%%%%%%%%%%%% 参数输入区 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% p0=imread('c:\f.bmp'); % 读入图像
du=1; %全息图像素x方向的间隔
dv=1; %全息图像素y方向的间隔
M=1; %希望在第几级衍射光中得到图像(这里选择了第一级) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%p1=rgb2gray(p0); %图像变为灰度图像
p1=p0; %%%%%%%%%%%%%%%%%%%%%%%%%% 暂时使用 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure;
imshow(p1); %显示灰度图象
p2=double(p1);
[nx,ny]=size(p2); %得到图像的大小mx,ny，并且把每个像素作为一个
抽样点
p3=p2.*exp(i*2*pi*rand(nx,ny)); %对图片附加随机的相位
fftp3=fftshift(fft2(fftshift(p3))); %对图像进行二维傅里叶变换
amplitude=abs(fftp3); %得到复振幅的振幅值
phase=angle(fftp3); %得到复振幅的相位值
amplitude=amplitude/max(max(amplitude)); %归一化
H=amplitude*dv; %矩孔的高度
P=phase/(2*pi*M)*du; %矩孔中心与抽样区中心的偏移距离
W=0.5/M*du; %矩孔的宽度
normal=find(abs(P+sign(P).*W/2)<=du/2); %得到矩孔不用模式溢出校正的序号leftOF=find((P-W/2)<(-du/2)); %得到矩孔左边溢出的序号
rightOF=find((P+W/2)>(du/2)); %得到矩孔右边溢出的序号
uleft=zeros(nx,ny); %矩孔左侧坐标
uright=zeros(nx,ny); %矩孔右侧坐标
vup=zeros(nx,ny); %矩孔上侧坐标
vdown=zeros(nx,ny); %矩孔下侧坐标
mark=zeros(nx,ny);
mark(normal)=1; %标志那些位置没有溢出
%%% 先求出所有的矩孔四个边相对抽样区中心的偏移量
vup=H/2; %矩孔的上边界
vdown=-H/2; %矩孔的下边界
uleft(normal)=P(normal)-W/2; %没有溢出矩孔的左边界
uright(normal)=P(normal)+W/2; %没有溢出矩孔的右边界
%矩孔溢出,要分成两块,但是中间的白色块也是一个矩形,所以表示白色的矩形块
uleft(leftOF)=P(leftOF)+W/2; %白色矩形块左边界
uright(leftOF)=P(leftOF)-W/2+du; %白色矩形块右边界
uleft(rightOF)=P(rightOF)+W/2-du; %白色矩形块左边界
uright(rightOF)=P(rightOF)-W/2; %白色矩形块右边界
%%% 实现模式溢出的校正及画图
figure;
axis([0 nx 0 ny]);
hold on;。