计算全息通过计算机模拟三维物体在全息面上的复振幅分布(计算全息图)并加载到空间光调制器(spatial light modulator，SLM)上，由再现光在空间上呈现出三维物体真正的三维(3 dimention, 3D)像^[1]，是一种非常理想的真三维裸眼3D技术，近年来在理论与应用研究方面得到了广泛的重视^[2-3]。尽管计算全息在算法与实验设备等方面取得了很大进展，但离实用阶段仍相距甚远，其主要原因在于：三维物体信息量丰富，空间结构复杂，且难以用具体函数描述其物波分布, 导致其全息计算复杂度高，运算量巨大。全息图生成达不到计算全息三维实时显示要求，是亟待解决的关键问题之一^[4-6]。

计算全息图(computer generated hologram, CGH)的形成涉及物体表面光源衍射场的计算，目前主要方法可概略分为点源法和面元法两类^{[1, 7]}。点源法需将被记录物体离散成大量的点(光源)，并模拟计算所有点发出的光线在全息面上与参考光干涉叠加而获得全息图^[4]，运算量大，十分耗时。从实时处理角度出发，一系列基于点源法的计算全息优化算法相继被提出，其代表是采用“空间换时间”策略的查找表(look-up tables，LUT)及各种改进算法。通过离线计算单个物点在全息面的贡献, 并按特定结构(三维表格)预先存储，实际生成全息图时仅需根据物体表面位置与查找表间的关联进行简单计算即可^[8]，但需大量存储空间用于三维表格数据加载。沿着这一思路并通过对点源法计算公式的合理简化，将三维LUT拆分成两个分别存储水平、竖直方向调制因子的二维LUT^[9]，该分裂查找表(split look-up tables，S-LUT)方式不仅大幅降低了查找表的内存需求，还可提高计算速度，再现三维物体时只需调用两个LUT中的调制因子就能获得全息图；基于夫琅禾费近似原理提出的压缩查找表(compressed-LUT，C-LUT)可进一步节省存储空间^[10]，但易造成重建图像失真，其原因在于忽略深度信息降低了全息图水平和竖直方向的相位精度；文献[11]提出的三角函数查找表(trigonometric-LUT, T-LUT), 通过在查找表离线计算中引入深度信息可有效减少在线计算复杂度，若采用单精度数据存储，则实际内存占用还可降低一半。

面元法是基于特殊形状均匀面源的频谱有解析解及衍射的角谱理论而形成，其代表是平面波角谱及其改进算法^[12-14]。针对由二维表面组成的三维物体，文献[12]基于自定义的二维表面特性函数，利用快速傅里叶变换(fast Fourier transformation，FFT)单独计算每个表面发射光场，并在全息面上进行综合，从而获得全息图，因每个表面仅计算一次FFT，该方法可有效节省时间，但三维物体表面轮廓过于复杂时处理困难；文献[13]直接解析计算普通三角形在全息图平面的频率分布，可避免对每个小面片的FFT计算，算法复杂度仅取决于全息图分辨率和3D模型多边形数目，速度更快，但未考虑表面相互遮挡及隐藏面消除；文献[14]给出了有遮挡的三维表面物体辐射光场数学描述，全息图光场计算时，通过光线跟踪选取特定视角的可视面片，从而达到移除被遮挡面、解决遮挡目的。总体上，因计算时需对复杂物体表面进行面元优化分解、考虑面元遮挡及边界对重建物体形状影响，面元法具有较高的计算复杂度高。文献[15]对三维物体进行分层并计算每层图到全息面的菲涅尔衍射分布，叠加后再通过相干光的干涉得到全息图；文献[16]通过获取三维物体在不同视角的投影图，计算每个投影图对全息面像素的贡献值来获得计算全息图光场，其对光学系统要求较低，白光条件下即可获得真实物体的全息图，但再现像质量受投影图片数量影响大。

点源法三维全息计算原理简单，操作灵活，重建像质量高且具有真三维特性，本文通过改进S-LUT算法以进一步提升点源法的计算性能，为构建基于计算全息的三维显示系统提供参考。

1 S-LUT算法原理

如图 1所示，设被记录三维物体与全息图的坐标系分别为O-XYZ和o-xy，坐标轴OX//ox，OY//oy，若三维物体被离散为N个空间点$\left( {{X}_{j}}, {{Y}_{j}}, {{Z}_{j}} \right)$，$j\in \left[ 1, N \right]$, 则根据点源法可计算得到各物点发出的光在全息面(x, y)处的复振幅^[8]：

$ I\left( {{x}_{h}}, {{y}_{h}} \right)=\overset{N}{\mathop{\underset{j}{\mathop{\mathop{\sum }^{}}}\, }}\, {{A}_{j}}\text{exp}\left( \text{i}\frac{2\text{ }\!\!\pi\!\!\text{ }}{\lambda }R \right) $

(1)

式中，$I\left( {{x}_{h}}, {{y}_{h}} \right)$是物光波到达全息面的复振幅；R表示物点j到全息面(x, y点的距离；λ是记录光的波长；A_j为j点的振幅；N为物点总数。

由式(1)可看出，点源法存在大量距离(相位)重复计算，且主要为平方、开方、三角函数等耗时的符号运算，传统查找表算法在线运行时虽通过直接地址查找能减少运行时间，但需构建大规模LUT存储离线计算输出的相位信息，所需存储空间达GB量级^[1]。为有效降低存储要求，S-LUT算法采用了查找表分裂、降维策略，原理如下^[9]。

以规整的空间网格表征物体点，沿网格(空间坐标轴)方向的采样点数分别为N_X、N_Y、N_Z，物体空间点数为：$N={{N}_{X}}\times {{N}_{Y}}\times {{N}_{Z}}$，全息面尺寸为 W×H, 若$\left| {{x}_{h}}-{{X}_{j}} \right|\ll {{Z}_{j}}$，$\left| {{y}_{h}}-{{Y}_{j}} \right|\ll {{Z}_{j}}$，则式(1)经适当简化可改写为：

$ I\left( {{x}_{h}}, {{y}_{h}} \right)=\overset{N}{\mathop{\underset{j}{\mathop{\mathop{\sum\limits }^{}}}\, }}\, {{A}_{j}}H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right) $

(2)

其中，

式中，$H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)$为水平调制因子；$V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right)$为竖直调制因子。对于Y坐标相同的同列物点$\left( {{X}_{j}}, {{Y}_{j}}, {{Z}_{j}} \right)$，因$H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)$相同，它们在全息图上点$\left( {{x}_{h}}, {{y}_{h}} \right)$的作用可表示为：

$ I\left( {{x}_{h}}, {{y}_{h}} \right){{|}_{\left( {{X}_{j}}, {{Y}_{j}}, {{Z}_{j}} \right)}}=H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)\overset{{{N}_{Y}}}{\mathop{\underset{j=0}{\mathop{\mathop{\sum }^{}}}\, }}\, {{A}_{j}}V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right) $

(3)

因$\text{ }\!\!\Delta\!\!\text{ }x$和$\overset{N}{\mathop{\underset{j}{\mathop{\mathop{\sum }^{}}}\, }}\, {{A}_{j}}V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right)$相互独立，故式(3)可分两步进行：

$ S\left( {{y}_{h}} \right)=\overset{{{N}_{Y}}}{\mathop{\underset{j=0}{\mathop{\mathop{\sum }^{}}}\, }}\, {{A}_{j}}V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right) $

(4)

$ I\left( {{x}_{h}}, {{y}_{h}} \right){{|}_{\left( {{X}_{j}}, {{Y}_{j}}, {{Z}_{j}} \right)}}=H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)S\left( {{y}_{h}} \right) $

(5)

叠加各列物点复振幅获得最终全息图：

$ I\left( {{x}_{h}}, {{y}_{h}} \right)=\mathop{\sum }^{}I\left( {{x}_{h}}, {{y}_{h}} \right){{|}_{\left( {{X}_{j}}, {{Y}_{j}}, {{Z}_{j}} \right)}} $

(6)

由式(4)~(6)可知，S-LUT在线计算的总次数为${{N}_{Z}}\left[ {{N}_{X}}\left( {{N}_{Y}}H+WH \right)+WH \right]$，每次计算包含一次加法和乘法运算，其中，式(4)计算次数为N_YH，式(5)计算次数为WH。

2 算法加速

本文对S-LUT算法的改进主要有两个方面：(1)改变查找表存储内容，进一步减少在线计算次数，并对在线计算过程进行统一计算设备架构(compute unified device architecture，CUDA)并行加速；(2)实施相位迭代计算，提高离线查找表生成效率。

2.1 在线加速

为进一步减少在线计算次数，可考虑存储物点每一列对全息面每一行像素的贡献S(y_h)，即将式(4)也转化为离线计算，有：

$ S\left( {{y}_{h}}, k \right)=\overset{{{N}_{Y}}}{\mathop{\underset{j}{\mathop{\mathop{\sum }^{}}}\, }}\, {{A}_{j}}V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{z}_{j}} \right)|k~(k=0, 1\ldots {{N}_{x}}-1) $

(7)

相应地，以$V'\left( S\left( {{y}_{h}}, k \right), {{Z}_{j}} \right)$(称为新竖直调制因子)取代$V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right)$进行存储，$H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)$则保持不变，图 2~3分别给出了两种竖直调制因子的查找表数据存储结构。

分析可知，若全息图上每个像素所占空间为M，$V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right)$和$V'\left( S\left( {{y}_{h}}, k \right), {{Z}_{j}} \right)$所需存储空间分别为${{N}_{Y}}\times H\times {{N}_{Z}}\times M$和${{N}_{X}}\times H\times {{N}_{Z}}\times M$，若被记录物体沿坐标轴X、Y方向的采样点数N_X、N_Y相同，则两种竖直调制因子存储空间大小不变；若N_X>N_Y，则新竖直调制因子存储空间将增加，否则减少，但改进S-LUT的图形处理器(graphic processing unit，GPU)在线计算总次数减少为${{N}_{Z}}\times \left( {{N}_{X}}+1 \right)\times \left( W\times H \right)$，理论上效率提升为：

根据新查找表数据的存储内容、结构，本文对全息图在线计算部分进行CUDA并行加速，核函数设计为$W/b\times H/b$个线程块(线程块大小由GPU性能决定)，每个线程块负责大小为$b\times b$的线程(子全息图，共$W\times H$个)，单个线程负责子全息图内像素点复振幅计算。为减少全局存储器中表格数据重复读取次数，将两个表格均划分成维度为b的多个方阵，子全息图u_sub由各方阵乘积给出，过程见图 4。首先从全局存储器中载入表格方阵到共享存储器进行乘积运算，然后由各线程将像素复振幅计算结果汇总到寄存器中，最后写入全局存储器。这种方式可有效节省全局存储器带宽(表格分块后，仅需读取$W/b+H/b$次)，并充分利用速度较快的共享存储器降低延时。

2.2 离线加速

现有点源法在利用各种离线查找表进行计算全息加速时均未考虑查找表生成效率，从应用上讲应予以统筹兼顾，本文根据文献[17]的相位迭代计算思想，离线生成S-LUT查找表时以加、减替换耗时的平方、开根号运算，可显著提高其计算效率。因S-LUT水平、竖直因子计算形式一样，以$H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)$为例推导其递归计算公式，令${{\theta }_{H}}\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)$表示空间点$\left( {{X}_{j}}, {{Y}_{j}}, {{Z}_{j}} \right)$在全息面$\left( {{x}_{h}}, {{y}_{h}} \right)$处的衍射相位，即有：

$ \begin{align} & H\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)=\text{exp}\left( \text{i}k\left( {{Z}_{j}}+\frac{\text{ }\!\!\Delta\!\!\text{ }{{x}^{2}}}{2{{Z}_{j}}} \right) \right)= \\ & \ \ \ \ \ \ \ \ \text{exp}\left( \text{i}{{\theta }_{H}}\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right) \right) \\ \end{align} $

(8)

$ \begin{matrix} {{\theta }_{H}}\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)=k\left( {{Z}_{j}}+\frac{\text{ }\!\!\Delta\!\!\text{ }{{x}^{2}}}{2{{Z}_{j}}} \right)=2\text{ }\!\!\pi\!\!\text{ }\left( \frac{{{Z}_{j}}}{\lambda }+\frac{\text{ }\!\!\Delta\!\!\text{ }{{x}^{2}}}{2\lambda {{Z}_{j}}} \right)= \\ 2\text{ }\!\!\pi\!\!\text{ }\left( {{\theta }_{Z}}+{{\theta }_{X}} \right) \\ \end{matrix} $

(9)

式中，${{\theta }_{Z}}=\frac{{{Z}_{j}}}{\lambda };{{\theta }_{X}}\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)=\frac{\text{ }\!\!\Delta\!\!\text{ }{{x}^{2}}}{2\lambda {{Z}_{j}}}$。

令d表示全息面上像素水平移动距离，则如图 5所示，同一物点在全息面$\left( {{x}_{h}}+d, {{y}_{h}} \right)$处的衍射相位为：

$\begin{matrix} {{\theta }_{H}}\left( \text{ }\!\!\Delta\!\!\text{ }x+d, {{Z}_{j}} \right)={{\theta }_{Z}}+\frac{\text{ }\!\!\Delta\!\!\text{ }{{x}^{2}}}{2\lambda {{Z}_{j}}}+\frac{\left( 2d\text{ }\!\!\Delta\!\!\text{ }x+{{d}^{2}} \right)}{2\lambda {{Z}_{j}}}= \\ {{\theta }_{H}}\left( \text{ }\!\!\Delta\!\!\text{ }x, {{Z}_{j}} \right)+{{T}_{d}} \\ \end{matrix}$

(10)

则当d=1时，有：${{T}_{1}}=\frac{\left( 2\text{ }\!\!\Delta\!\!\text{ }x+{{1}^{2}} \right)}{2\lambda {{Z}_{j}}}$

当d=2时，有：

当d=n+1时：

以此类推，当d=n时：

${{T}_{n}}={{T}_{n-1}}+{{T}_{1}}+\left( n-1 \right)\Delta $

(11)

综合式(8)~(11)可知，同一物点对全息面上某一行像素的衍射相位可迭代计算得到，除首像素外，其余像素的相位计算均为简单的加、减及乘法，故能有效提高H(Δx,Z_j)的计算效率，即有：

$\begin{matrix} {{\theta }_{H}}\left( \text{ }\!\!\Delta\!\!\text{ }x+n, {{Z}_{j}} \right)=2\text{ }\!\!\pi\!\!\text{ }\left( {{\theta }_{Z}}+{{\theta }_{X}}\left( \text{ }\!\!\Delta\!\!\text{ }x+n, {{Z}_{j}} \right) \right)= \\ 2\text{ }\!\!\pi\!\!\text{ }\left( {{\theta }_{Z}}+{T_{n}} \right) \\ \end{matrix}$

(12)

$H\left( \text{ }\!\!\Delta\!\!\text{ }x+n, {{Z}_{j}} \right)=\text{exp}\left( \text{i}{{\theta }_{H}}\left( \text{ }\!\!\Delta\!\!\text{ }x+n, {{Z}_{j}} \right) \right)$

(13)

同理，可迭代计算得到竖直调节因子$V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right)$。需要说明的是，本文$V\text{ }\!\!'\!\!\text{ }\left( S\left( {{y}_{h}}, k \right), {{Z}_{j}} \right)$的计算是通过$V\left( \text{ }\!\!\Delta\!\!\text{ }y, {{Z}_{j}} \right)$的累加、乘法运算实现的，其计算效率也相应得到提升。

3 计算全息图实验与分析

本文在PC机(Windows 10操作系统，VS2012编译环境)上实现上述算法，算法验证所需两物体虚拟三维模型由离散点云构成(见图 6)，并按要求规整到空间网格：物体1(头像)的采样点数$\left( {{N}_{X}}, {{N}_{Y}}, {{N}_{Z}} \right)=\left( 300, \text{ }\!\!~\!\!\text{ }300, 300 \right)$，物体2(圣诞老人)采样点数$\left( {{N}_{X}}, {{N}_{Y}}, {{N}_{Z}} \right)=\left( 400, 400, 300 \right)$，全息面分辨率为$W\times H=1\text{ }\!\!~\!\!\text{ }024\times 768$像素，波长λ=655nm。算法实施平台性能主要参数见表 1。算法性能测试包含查找表离线生成效率以及全息图GPU在线计算效率。

根据物体采样点数、全息计算参数及§3给出的算法加速策略，本文先通过相位迭代离线计算出查找表数据并存储到硬盘，然后将查找表数据加载到GPU在线计算出振幅型全息图，表 2和表 3分别列出了两测试物体在不同层数下的正面计算全息图生成时间开销。定义查找表生成效率提升度R₁=原查找表生成时间/新查找表迭代生成时间－1, 则由表 2和表 3可知，本文算法R₁约为1~1.5, 同一物体不同层R₁保持稳定, 微小的数值波动源于不同层的物体空间点数存在差异。考虑到本文算法查找表是在原算法竖直调节因子基础上产生, 额外增加了计算开销，则相同表格下相位迭代计算带来的查找表生成效率提升度更高。

此外，需要指出的是，采用加、减运算生成查找表时更有利于其GPU加速(本文尚未实施)，若本文查找表计算GPU加速后能达到常规30倍提升，则其生成时间仅约占全息图在线计算时间的3%~5%，已无需进行离线计算，从而有利于节省存储空间。全息图GPU在线计算方面，定义其效率提升度R₂=(原算法GPU计算时间/本文算法GPU计算时间－1)×100%, ，则R₂与§2.1给出的理论值相吻合:物体1的R₂理论值为300/1 024≈30%，该物体在20层时R₂为26%, 随着层数增加，R₂提升幅度放缓, 当层数达到150层后趋于稳定在15%左右; 类似地, 物体2的R₂理论值为400/1 024≈40%，在20层时R₂为35.5%, 随着层数增加，R₂提升幅度放缓, 当层数达到150层后稳定在21%左右。

对比表 2和3还可发现, 相同全息面分辨率下, 物体空间点数越多, 改进算法效果越显著, 见图 7和图 8。

本文利用信噪比(signal-to-noise ratio，SNR)对全息图质量进行评估，令${{I}_{c}}\left( x, y \right)$、${{I}_{p}}\left( x, y \right)$分别表示原算法和本文算法生成的计算全息图，SNR计算公式为：

表 4列出了测试物体不同层数下的计算全息图SNR。文献[7]实施相位迭代计算获得的全息图SNR约为50dB，并认为当全息图SNR>30时, 其再现像视觉质量上的差异已难以识别，据此而言，本文算法两测试物体全息图SNR与文献[7]相接近, 且不同层SNR保持稳定(限于篇幅，仅列出50~200层结果)，表明本文查找表快速生成方式可保证计算全息图及其再现像的质量。图 9~10分别给出了50层时两测试物体的振幅型计算全息图及其再现像，可见两种方法输出结果在视觉上的差异已难以觉察。表 4同时给出了测试物体再现像SNR，可以看出其值在不同层保持稳定，但大小略低于相应的全息图SNR，原因在于物点空间再现时存在二次衍射噪声，故再现像SNR有所降低。综上所述，本文研究结果表明：

1) 采用相位迭代计算过程，以加、减替换耗时的平方、开根号运算可大幅提高查找表离线生成效率, 并随物体点数、层数增加而增大。

2) 通过改变竖直调制因子查找表内容为存储物点每一列对全息面每一行像素的贡献，可有效减少S-LUT的全息图在线计算次数，从而提升在线计算效率(不少于15%), 且CUDA并行设计更为简单。

4 结语

与利用人眼双目视差或视觉暂留效应的传统3D显示技术相比，计算全息三维显示能提供人眼视觉系统所需的全部深度信息，同时还具有传统光学全息所不具备的低噪声、高重复性及可获得虚拟物体全息图等方面的优点，是一种非常理想的真三维裸视技术，困难在于三维物体全息计算复杂度高、运算量巨大，达不到实时显示要求。本文从在线、离线计算两方面同时对S-LUT算法进行性能提升，在不影响计算全息图及其再现像质量的前提下显著提升计算效率。下一步将对改进的离线查找表生成算法进行GPU加速以使其能达到在线计算时间要求，进一步完善现有算法并将其应用于裸视3D显示系统应用开发之中。