积分项周期延拓快速计算大尺寸菲涅尔全息图

姚建云，应晓霖，许富洋，吴琼，张海花，李勇

（1 浙江师范大学信息光学研究所，浙江金华，321004）

（2 浙江省光信息检测与显示技术研究重点实验室，浙江金华，321004）

0 引言

全息术由GABOR D 在1948年提出［1］，三维显示是全息术的重要应用领域之一［2-7］。全息图能够同时记录物光场的振幅和相位信息，其再现像可以提供人眼观察所需的所有深度线索，因此全息技术被认为是理想的三维显示技术。随着计算机技术与电子技术的不断发展，计算全息（Computer Generated Hologram，CGH）应运而生［8-9］。它通过计算机数值计算代替光学记录过程来获得全息图，再将全息图加载到空间光调制器或输出到记录介质上实现三维物体的再现。这一技术的出现使得虚实场景结合的全息图制作成为可能，并且大大推动了动态全息三维显示的发展。

庞大的计算量是实现高质量计算机制全息三维显示的瓶颈之一。研究者提出了一系列的全息图快速算法。常见的CGH 算法有点元法、面元法及层析法等等。点元法将物体分解为一系列离散发光点，每个点发出球面波。全息图由这些球面波在全息面上叠加后的物光与参考光干涉形成。该方法适合计算复杂场景的全息图。为提高点元法的计算速度，国内外学者提出了诸多算法，如查表法［10-13］（Look-up Table，LUT）、差分法［14］、分离变量法［15］及波前记录面［16-17］（Wave Recording Plane，WRP）法等。面元法主要分为基于采样的面元法［18-19］和解析面元法［20-21］。同一物体的面元数量远小于点元数量，因此面元法计算量比点元法少，但是该方法计算纹理复杂物体的全息图依然耗时较多。层析法［22-24］将三维物体分为多个相互平行的平面，利用快速傅里叶变换（Fast Fourier Transform，FFT）求解每一层物体衍射到全息面上的复振幅，然后叠加并与参考光干涉，最终得到全息图。由于层析法采用了FFT，有效地提高了计算速度。但采用FFT 计算物面到全息面上的衍射时，对数据的插值操作降低了计算速度，尤其是物面或全息面的长宽比大时。

当计算机制全息图的分辨率在数亿像素时，通常采用并行计算或光场与计算全息结合的方法进行全息图的计算。并行计算方法利用全息图计算中的独立性，采用多台计算设备同时计算，设备较为昂贵，计算时电能消耗大。文献［25］提出采用服务器集群的并行计算方法实现高加速比的全息图计算。文献［26］提出采用GPU 集群的并行算法，在由16 个节点构成、每个节点有16 个GPU 的超大规模并行机上实现了44 000×44 000 像素菲涅尔全息图的快速计算。光场与计算全息结合的方法利用一系列二维光场图像来计算全息图，其运动视差是离散的。增加运动视差的平滑性需要密集的光场图像，全息图的计算时间也随之明显增加。文献［27］提出采用频域拼接的方法实现了94 340×94 340 像素的彩色彩虹全息快速计算。文献［28］提出采用光场图像与全息透镜结合的方法实现了200 000×200 000 像素的全视差全息图快速计算。文献［29］提出采用角频率切片的由光场图像计算全息图的方法，实现了25 000×25 000 像素的全息图快速计算。本文通过分析层析法CGH 计算过程，提出了一种由三维体数据计算CGH 的快速算法。采用积分项周期延拓消除了冗余计算，同时利用人眼的低通滤波特性省去了图像插值中的数字低通滤波。

1 原理

1.1 离散菲涅尔衍射公式分析

图1 为层析法计算菲涅尔全息图的流程示意图，R表示参考光，H表示全息面。首先对物体分层，然后对每一层的数据计算其在全息面上的光场复振幅分布，将得到的所有复振幅分布求和并与参考光干涉得到要记录的全息图。

图1 层析法菲涅尔全息图计算Fig.1 Layer-based calculation model of Fresnel hologram

在计算每一层数据在全息面上的衍射场分布时需要通过菲涅尔衍射积分公式来完成，即

式中，o(x0，y0)表示物面上的复振幅分布，u(x，y)表示全息面上物光波的复振幅分布，k=2π/λ，λ为波长。

积分项式（2）为二维傅里叶变换，采用计算机数值计算时对空域（物面）和空频域（全息面）进行采样离散化后，利用FFT 算法实现。在水平、垂直方向采样间隔相同时，采用Δo和Δh分别表示物面和全息面的采样间隔。对物面离散化后式（2）表示为

由（3）式可知，物面离散化后全息面上的积分项是一个周期函数，是连续情况下全息面上积分项的周期延拓。离散傅里叶变换（Discrete Fourier Transformation，DFT）是取式（3）的一个周期并采样离散化，FFT 是DFT 的快速算法。因此FFT 实际隐含了周期性。当全息面或物面水平和垂直方向的采样点数分别为M和N时，根据FFT 原理，物面和全息面的尺寸必须满足

则式（3）水平和垂直方向的周期分别为MΔh和NΔh。

全息图尺寸与观察视角有关。如图2所示，在记录物体大小一定时，人眼透过全息图观察再现像，只有部分与瞳孔对应区域的衍射光进入人眼。人眼移动时，对应的全息图区域也随之移动，则观察到不同视角处的再现像。可见，在物体大小、全息面与物面距离确定且全息面采样间隔满足要求时，全息图越大，其可观察视角越大。

图2 在不同位置透过全息图观察再现像Fig.2 Observe the reconstructed image in different positions through the hologram

由式（4）可知，全息图越大，要求物面的采样间隔越小。人眼是一个低通系统，分辨率有限，物面采样间隔小到一定程度，超过了人眼分辨率，会造成信息冗余，增加不必要的存储和计算量。因此，在大视角全息图计算时，以人眼分辨率为基准确定物面的采样间隔。物面大小确定后，这一采样间隔对应的只是某一视角下与瞳孔对应的小部分全息图。增加视角，全息图尺寸必须增加。由式（4）知，必须减小物面采样间隔。在数字信号处理中，整数倍插值用以下模型描述：在原有样点之间插零，然后进行数字低通滤波。设采样后的物面复振幅为o（pΔo，qΔo），在物面相邻样点间插kx、ky个零，数据采样间隔变为Δo/kx和Δo/ky，则插零后的物面复振幅表示为o（m′Δo/kx，n′Δo/ky）。根据式（3），积分项变为

这是一个周期为（λzkx/Δo，λzkx/Δo）的周期函数。插零前后的物面采样点关系可表示为

则式（5）也可以表示为

式（7）是周期为（λz/Δo，λz/Δo）的周期函数，周期为式（5）的1/kx和1/ky。因此，插零后积分项的一个周期内有kx×ky个周期的原积分项。综上所述，物面插零积分项周期变大，变大部分是原积分项的周期延拓。插零后积分项一个周期的分布可表示为

如图2所示，人眼透过全息图观察再现像时，只有部分全息图上衍射的光进入人眼。相当于进行了如下运算

式中，（xc，yc）是眼睛对应的全息图窗口中心，Rect 表示二维矩形窗函数，X、Y分别是窗口x、y方向的宽度。当窗口宽度等于积分项周期时，一个完整的周期信息进入人眼。此时式（9）的逆傅里叶变换就是内插公式形式的采样定理，而全息图衍射再现就实现了逆傅里叶变换。这样，在全息面积分项周期延拓再利用人眼的低通滤波特性实现物面的插值运算，减少了大量冗余计算并且省去了数字插值中的数字低通滤波运算。

由上述分析可知，对于大视角计算全息图的积分项，在采样间隔满足人眼分辨率时，只需计算其中一个周期，而其他部分直接周期延拓即可。设要计算的全息图大小为M×N，将已知参量代入约束条件式（4）中计算出全息面积分项一个周期的像素数L×K。此时x，y方向上积分项的周期数为和。在传统算法中，需要根据计算的全息图大小和衍射距离等参量对物体的分层数据进行重采样或补零到与全息图相同大小，再对每一层数据做菲涅尔衍射计算求出全息面上的衍射场。此过程中利用FFT 计算积分项的计算量为。新算法中用FFT 计算积分项的计算量变为可见，当需要计算的全息图越大（即包含的周期数越多）时，新算法的计算速度优势越明显。通常物体与全息面水平和垂直方向的采样间隔相同，要求积分区域必须是正方形。对于非正方形全息图，通常通过物面补零到正方形计算后再裁剪来达到这一要求。在长宽比大的全息图计算时，本算法能够进一步减少冗余计算量。

1.2 积分项系数特点

用uc（x，y）表示式（1）的系数（后续表述中统称为系数），即

在计算某一层物体的衍射场时，z为定值则为常数。显然可以对x、y进行变量分离，表示为

式中，A表示常数。计算M×N大小的全息图时对整个全息面计算uc(x，y)需要进行2M×N次三角函数与4M×N次乘法计算。式（11）分离变量后，系数uc(x，y)在全息面上的分布分解成行列相互独立的分量ucx和ucy。因此只需要计算长度分别为1×M和1×N的两个向量存储在内存中，将积分项周期性延拓后乘上对应位置处的行、列系数即可得到全息面上的衍射场复振幅分布。此时对系数uc(x，y)的计算变为2M+N次三角函数与2(M×N+M+N) +M次乘法计算。由于三角函数计算耗时较长，采用行列分解后明显提高了计算速度。

1.3 算法流程

根据以上分析，将菲涅尔衍射计算分为系数和积分项两部分。系数部分采用行、列正交分解再合成计算；积分项部分采用计算一个周期再周期延拓；最后两部分合成得到全息面上的物光复振幅分布。快速算法具体流程为：

1）读入点云数据，进行分层；

2）计算一层物体在全息面上积分项的一个周期，将求得的积分项存入内存；

3）根据最终要计算的全息图大小计算该层数据对应系数行、列分量存入内存；

4）将该层的积分项周期延拓，根据对应位置的行号列号获取系数分量相乘得到完整系数，再将系数与积分项相乘得到该位置的物光复振幅，并与前一次计算好的复振幅叠加；

5）重复步骤2）～4）直到所有层计算完成；

6）将叠加后的复振幅与参考光干涉，经编码得到全息图。

2 实验结果与分析

为验证提出算法的正确性并分析计算速度提高情况，设计了程序进行实验。对比了新算法与传统层析法在再现像清晰度、视差效果、深度再现效果的差异；研究了不同全息图尺寸下，计算速度的提高情况。如图3所示，实验所选用的三维模型点云数据是由斯坦福大学三维点云数据库提供，模型包含35 947 个物点。实验中物体在x、y、z三个方向上的实际尺寸分别设置为46.8 mm、46.4 mm、36.3 mm。将点云数据分为40层，为了更好地观察视差效果，在第41 层处放置了棋盘格。中心层到全息面的距离为310 mm。实验中全息图计算采用Matlab 编程实现，在PC 机（CUP：Intel（R）Core（TM）i7-10700K@3.80GHz，内存：16GB）上运行。

图3 三维物体模型Fig.3 Model of three-dimensional object

2.1 实验与结果

为观察采用本算法所制作全息图的再现效果，计算了像素数为54 179×399 087 的全息图。之后用自行研制的全息打印机将全息图输出到光刻胶版上。打印机的分辨率为318 nm。物空间内有效样点数为520点×515 点×41 层。物体的采样间隔为0.09 mm，积分项一个周期的采样点数为4927×4927。将打印好的全息图放在如图4所示的再现光路中观察再现像，光源选择波长为473 nm 的激光，通过扩束和滤波后得到球面波照射到全息图上。拍摄实验结果所用相机为华为荣耀20 手机所搭载的相机。

图4 再现光路图Fig.4 The diagram of optical setup for reconstruction

图5（a）为全息图的局部，图5（b）和（c）分别为传统算法与快速算法的再现效果，未观察到二者的差异。图6 为相机在三个不同视角下拍摄到的再现像。红色框标示区域中，物体与棋盘格之间的遮挡关系正确。当多角度连续拍摄或人眼直接平移观察时可以看到连续的视差变化。

图5 部分全息图及两种算法的再现结果Fig.5 Partial hologram and reconstruction results of two algorithms

图6 不同视角下的再现像Fig.6 Reconstructed image from different viewing-angle

观察不同深度的再现像时，在图4 光路基础上在全息图右侧加入焦距为30 cm 的透镜，将再现像成像为等大的实像并用毛玻璃承接。经过测量得知加入透镜后整个成像系统在实像位置附近的景深约为14 mm，因此选择间隔大于景深，且特征较为明显的两个深度进行拍摄。图7 为两个不同深度处的再现像。图中深度为换算成虚像的深度，再现的深度信息准确。

图7 不同深度再现像对比Fig.7 Comparison of the reconstructed image in different depth

用层析法计算三维物体的全息图时，再现像的质量与数据分层间隔大小密切相关。层间隔满足人眼分辨率的情况下所观察到的再现像在各个视角下都是连续的，当不满足人眼分辨率时，大角度倾斜观察容易出现再现像的分层。在其他参数不变的前提下，将三维物体拉近到距离全息面100 mm 附近计算、制作全息图，在右边大角度倾斜观察得到如图8所示结果。由于观察距离变小，此位置上人眼分辨率更高，再现像在大视角下观察出现了明显的分层。要避免这一现象的出现需要根据人眼分辨率、物体的三维形貌、观察视角等选择合理的层间隔。

图8 再现像分层现象Fig.8 The slice phenomenon of reconstructed image

为研究快速算法的计算速度提高情况，选用上述三维数据模型在同一PC 机上用两种算法计算其菲涅尔全息图，取50 次运行时间的平均值。在全息面采样间隔318 nm 时，积分项一个周期的采样点数为4 927×4 927，当周期数在3 以上后传统算法计算过程中会由于计算机物理内存不足而需要读取硬盘，计算速度会明显降低，此时的计算速度不具可比性。因此采用较小的积分项周期进行计算，全息面的采样间隔为4.03 μm（实验室中SLM 的像素尺寸），物体距离全息面900 mm。根据1.1 节所述求出全息面上积分项一个周期内包含的像素数为1 174×1 174。计算了九组周期数从kx=2，ky=2 到kx=10，ky=10 的全息图。每一组分别采用快速算法周期延拓加系数行列分解（方法1）、快速算法仅积分项周期延拓（方法2）、及传统算法（方法3）三种方法计算，得到如表1所示的结果。图9 为不同算法计算时间与积分项周期数关系曲线。由结果可见系数行列分解后的快速算法速度最快，仅周期延拓的快速计算次之，且均快于传统算法。整体来看所计算的全息图越大时快速算法的速度优势越明显，当全息图分辨率达到11 740×11 740 像素时计算速度达到传统算法的13 倍。

图9 计算时间与周期数关系Fig.9 The relationship between the computing time and number of period

表1 不同算法计算时间对比Table 1 Comparison of computational time of different algorithm

为探究数据层数对速度提高倍数的影响，进行了层数从10 到50，以5 层为增量的全息图计算速度提高倍数对比实验。每组实验用两种算法计算周期数从kx=2，ky=2 到kx=8，ky=8 的全息图以及提速倍数，得到了分层数对快速算法的速度提高倍数无影响的结论。表2 给出了其中5 组数据。

表2 层数与速度提高倍数关系Table 2 Relationship between number of layers and speed increase multiple

2.2 分析与讨论

由实验结果可知，提出的算法在保证良好再现效果的同时显著提高了计算机制菲涅尔全息图的计算速度。在大尺寸全息图的计算上有着良好的表现。另外需要说明的是在表1 传统算法用时的数据中，全息图像素数增加到9 392×9 392 后计算用时会发生激增。这是由于传统层析法计算全息图时每一层数据都要插值或补零到与全息图相同的像素。计算量随着全息图尺寸的增加而呈平方式的增长，从而导致计算用时激增。相比于传统算法，快速算法在计算完成每一层在全息面上一个周期的积分项后，即可通过周期延拓来得到整个全息面的积分项。计算量随着全息图尺寸增大而平缓增长。

所提算法在计算完成积分项的一个周期后，后续计算互相独立，使得采用并行计算进一步提高全息图的计算速度得以实现。优化算法，采用多核CPU 或GPU 并行计算加速全息图的制作将是下一步的研究目标。

3 结论

计算全息数据量巨大。本文提出的快速算法，从离散菲涅尔衍射计算模型出发分析，考虑人眼分辨率，从系数和积分项两方面改进。将系数项行列正交分解再合成；只计算部分区域的积分项，再周期延拓。该算法所得全息图再现像质量与传统算法相当，同时减少了大量的冗余计算，计算速度随全息图尺寸的增大而快速提高。