姚建云，應曉霖，許富洋，吳瓊，張?；?，李勇

（1 浙江師范大學信息光學研究所，浙江金華，321004）

（2 浙江省光信息檢測與顯示技術(shù)研究重點實驗室，浙江金華，321004）

0 引言

全息術(shù)由GABOR D 在1948年提出［1］，三維顯示是全息術(shù)的重要應用領域之一［2-7］。全息圖能夠同時記錄物光場的振幅和相位信息，其再現(xiàn)像可以提供人眼觀察所需的所有深度線索，因此全息技術(shù)被認為是理想的三維顯示技術(shù)。隨著計算機技術(shù)與電子技術(shù)的不斷發(fā)展，計算全息（Computer Generated Hologram，CGH）應運而生［8-9］。它通過計算機數(shù)值計算代替光學記錄過程來獲得全息圖，再將全息圖加載到空間光調(diào)制器或輸出到記錄介質(zhì)上實現(xiàn)三維物體的再現(xiàn)。這一技術(shù)的出現(xiàn)使得虛實場景結(jié)合的全息圖制作成為可能，并且大大推動了動態(tài)全息三維顯示的發(fā)展。

龐大的計算量是實現(xiàn)高質(zhì)量計算機制全息三維顯示的瓶頸之一。研究者提出了一系列的全息圖快速算法。常見的CGH 算法有點元法、面元法及層析法等等。點元法將物體分解為一系列離散發(fā)光點，每個點發(fā)出球面波。全息圖由這些球面波在全息面上疊加后的物光與參考光干涉形成。該方法適合計算復雜場景的全息圖。為提高點元法的計算速度，國內(nèi)外學者提出了諸多算法，如查表法［10-13］（Look-up Table，LUT）、差分法［14］、分離變量法［15］及波前記錄面［16-17］（Wave Recording Plane，WRP）法等。面元法主要分為基于采樣的面元法［18-19］和解析面元法［20-21］。同一物體的面元數(shù)量遠小于點元數(shù)量，因此面元法計算量比點元法少，但是該方法計算紋理復雜物體的全息圖依然耗時較多。層析法［22-24］將三維物體分為多個相互平行的平面，利用快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）求解每一層物體衍射到全息面上的復振幅，然后疊加并與參考光干涉，最終得到全息圖。由于層析法采用了FFT，有效地提高了計算速度。但采用FFT 計算物面到全息面上的衍射時，對數(shù)據(jù)的插值操作降低了計算速度，尤其是物面或全息面的長寬比大時。

當計算機制全息圖的分辨率在數(shù)億像素時，通常采用并行計算或光場與計算全息結(jié)合的方法進行全息圖的計算。并行計算方法利用全息圖計算中的獨立性，采用多臺計算設備同時計算，設備較為昂貴，計算時電能消耗大。文獻［25］提出采用服務器集群的并行計算方法實現(xiàn)高加速比的全息圖計算。文獻［26］提出采用GPU 集群的并行算法，在由16 個節(jié)點構(gòu)成、每個節(jié)點有16 個GPU 的超大規(guī)模并行機上實現(xiàn)了44 000×44 000 像素菲涅爾全息圖的快速計算。光場與計算全息結(jié)合的方法利用一系列二維光場圖像來計算全息圖，其運動視差是離散的。增加運動視差的平滑性需要密集的光場圖像，全息圖的計算時間也隨之明顯增加。文獻［27］提出采用頻域拼接的方法實現(xiàn)了94 340×94 340 像素的彩色彩虹全息快速計算。文獻［28］提出采用光場圖像與全息透鏡結(jié)合的方法實現(xiàn)了200 000×200 000 像素的全視差全息圖快速計算。文獻［29］提出采用角頻率切片的由光場圖像計算全息圖的方法，實現(xiàn)了25 000×25 000 像素的全息圖快速計算。本文通過分析層析法CGH 計算過程，提出了一種由三維體數(shù)據(jù)計算CGH 的快速算法。采用積分項周期延拓消除了冗余計算，同時利用人眼的低通濾波特性省去了圖像插值中的數(shù)字低通濾波。

1 原理

1.1 離散菲涅爾衍射公式分析

圖1 為層析法計算菲涅爾全息圖的流程示意圖，R表示參考光，H表示全息面。首先對物體分層，然后對每一層的數(shù)據(jù)計算其在全息面上的光場復振幅分布，將得到的所有復振幅分布求和并與參考光干涉得到要記錄的全息圖。

圖1 層析法菲涅爾全息圖計算Fig.1 Layer-based calculation model of Fresnel hologram

在計算每一層數(shù)據(jù)在全息面上的衍射場分布時需要通過菲涅爾衍射積分公式來完成，即

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

積分項式（2）為二維傅里葉變換，采用計算機數(shù)值計算時對空域（物面）和空頻域（全息面）進行采樣離散化后，利用FFT 算法實現(xiàn)。在水平、垂直方向采樣間隔相同時，采用Δo和Δh分別表示物面和全息面的采樣間隔。對物面離散化后式（2）表示為

由（3）式可知，物面離散化后全息面上的積分項是一個周期函數(shù)，是連續(xù)情況下全息面上積分項的周期延拓。離散傅里葉變換（Discrete Fourier Transformation，DFT）是取式（3）的一個周期并采樣離散化，F(xiàn)FT 是DFT 的快速算法。因此FFT 實際隱含了周期性。當全息面或物面水平和垂直方向的采樣點數(shù)分別為M和N時，根據(jù)FFT 原理，物面和全息面的尺寸必須滿足

則式（3）水平和垂直方向的周期分別為MΔh和NΔh。

全息圖尺寸與觀察視角有關。如圖2所示，在記錄物體大小一定時，人眼透過全息圖觀察再現(xiàn)像，只有部分與瞳孔對應區(qū)域的衍射光進入人眼。人眼移動時，對應的全息圖區(qū)域也隨之移動，則觀察到不同視角處的再現(xiàn)像?？梢姡谖矬w大小、全息面與物面距離確定且全息面采樣間隔滿足要求時，全息圖越大，其可觀察視角越大。

圖2 在不同位置透過全息圖觀察再現(xiàn)像Fig.2 Observe the reconstructed image in different positions through the hologram

由式（4）可知，全息圖越大，要求物面的采樣間隔越小。人眼是一個低通系統(tǒng)，分辨率有限，物面采樣間隔小到一定程度，超過了人眼分辨率，會造成信息冗余，增加不必要的存儲和計算量。因此，在大視角全息圖計算時，以人眼分辨率為基準確定物面的采樣間隔。物面大小確定后，這一采樣間隔對應的只是某一視角下與瞳孔對應的小部分全息圖。增加視角，全息圖尺寸必須增加。由式（4）知，必須減小物面采樣間隔。在數(shù)字信號處理中，整數(shù)倍插值用以下模型描述：在原有樣點之間插零，然后進行數(shù)字低通濾波。設采樣后的物面復振幅為o（pΔo，qΔo），在物面相鄰樣點間插kx、ky個零，數(shù)據(jù)采樣間隔變?yōu)棣/kx和Δo/ky，則插零后的物面復振幅表示為o（m′Δo/kx，n′Δo/ky）。根據(jù)式（3），積分項變?yōu)?/p>

這是一個周期為（λzkx/Δo，λzkx/Δo）的周期函數(shù)。插零前后的物面采樣點關系可表示為

則式（5）也可以表示為

式（7）是周期為（λz/Δo，λz/Δo）的周期函數(shù)，周期為式（5）的1/kx和1/ky。因此，插零后積分項的一個周期內(nèi)有kx×ky個周期的原積分項。綜上所述，物面插零積分項周期變大，變大部分是原積分項的周期延拓。插零后積分項一個周期的分布可表示為

如圖2所示，人眼透過全息圖觀察再現(xiàn)像時，只有部分全息圖上衍射的光進入人眼。相當于進行了如下運算

式中，（xc，yc）是眼睛對應的全息圖窗口中心，Rect 表示二維矩形窗函數(shù)，X、Y分別是窗口x、y方向的寬度。當窗口寬度等于積分項周期時，一個完整的周期信息進入人眼。此時式（9）的逆傅里葉變換就是內(nèi)插公式形式的采樣定理，而全息圖衍射再現(xiàn)就實現(xiàn)了逆傅里葉變換。這樣，在全息面積分項周期延拓再利用人眼的低通濾波特性實現(xiàn)物面的插值運算，減少了大量冗余計算并且省去了數(shù)字插值中的數(shù)字低通濾波運算。

由上述分析可知，對于大視角計算全息圖的積分項，在采樣間隔滿足人眼分辨率時，只需計算其中一個周期，而其他部分直接周期延拓即可。設要計算的全息圖大小為M×N，將已知參量代入約束條件式（4）中計算出全息面積分項一個周期的像素數(shù)L×K。此時x，y方向上積分項的周期數(shù)為和。在傳統(tǒng)算法中，需要根據(jù)計算的全息圖大小和衍射距離等參量對物體的分層數(shù)據(jù)進行重采樣或補零到與全息圖相同大小，再對每一層數(shù)據(jù)做菲涅爾衍射計算求出全息面上的衍射場。此過程中利用FFT 計算積分項的計算量為。新算法中用FFT 計算積分項的計算量變?yōu)榭梢?，當需要計算的全息圖越大（即包含的周期數(shù)越多）時，新算法的計算速度優(yōu)勢越明顯。通常物體與全息面水平和垂直方向的采樣間隔相同，要求積分區(qū)域必須是正方形。對于非正方形全息圖，通常通過物面補零到正方形計算后再裁剪來達到這一要求。在長寬比大的全息圖計算時，本算法能夠進一步減少冗余計算量。

1.2 積分項系數(shù)特點

用uc（x，y）表示式（1）的系數(shù)（后續(xù)表述中統(tǒng)稱為系數(shù)），即

在計算某一層物體的衍射場時，z為定值則為常數(shù)。顯然可以對x、y進行變量分離，表示為

式中，A表示常數(shù)。計算M×N大小的全息圖時對整個全息面計算uc(x，y)需要進行2M×N次三角函數(shù)與4M×N次乘法計算。式（11）分離變量后，系數(shù)uc(x，y)在全息面上的分布分解成行列相互獨立的分量ucx和ucy。因此只需要計算長度分別為1×M和1×N的兩個向量存儲在內(nèi)存中，將積分項周期性延拓后乘上對應位置處的行、列系數(shù)即可得到全息面上的衍射場復振幅分布。此時對系數(shù)uc(x，y)的計算變?yōu)?M+N次三角函數(shù)與2(M×N+M+N) +M次乘法計算。由于三角函數(shù)計算耗時較長，采用行列分解后明顯提高了計算速度。

1.3 算法流程

根據(jù)以上分析，將菲涅爾衍射計算分為系數(shù)和積分項兩部分。系數(shù)部分采用行、列正交分解再合成計算；積分項部分采用計算一個周期再周期延拓；最后兩部分合成得到全息面上的物光復振幅分布。快速算法具體流程為：

1）讀入點云數(shù)據(jù)，進行分層；

2）計算一層物體在全息面上積分項的一個周期，將求得的積分項存入內(nèi)存；

3）根據(jù)最終要計算的全息圖大小計算該層數(shù)據(jù)對應系數(shù)行、列分量存入內(nèi)存；

4）將該層的積分項周期延拓，根據(jù)對應位置的行號列號獲取系數(shù)分量相乘得到完整系數(shù)，再將系數(shù)與積分項相乘得到該位置的物光復振幅，并與前一次計算好的復振幅疊加；

5）重復步驟2）～4）直到所有層計算完成；

6）將疊加后的復振幅與參考光干涉，經(jīng)編碼得到全息圖。

2 實驗結(jié)果與分析

為驗證提出算法的正確性并分析計算速度提高情況，設計了程序進行實驗。對比了新算法與傳統(tǒng)層析法在再現(xiàn)像清晰度、視差效果、深度再現(xiàn)效果的差異；研究了不同全息圖尺寸下，計算速度的提高情況。如圖3所示，實驗所選用的三維模型點云數(shù)據(jù)是由斯坦福大學三維點云數(shù)據(jù)庫提供，模型包含35 947 個物點。實驗中物體在x、y、z三個方向上的實際尺寸分別設置為46.8 mm、46.4 mm、36.3 mm。將點云數(shù)據(jù)分為40層，為了更好地觀察視差效果，在第41 層處放置了棋盤格。中心層到全息面的距離為310 mm。實驗中全息圖計算采用Matlab 編程實現(xiàn)，在PC 機（CUP：Intel（R）Core（TM）i7-10700K@3.80GHz，內(nèi)存：16GB）上運行。

圖3 三維物體模型Fig.3 Model of three-dimensional object

2.1 實驗與結(jié)果

為觀察采用本算法所制作全息圖的再現(xiàn)效果，計算了像素數(shù)為54 179×399 087 的全息圖。之后用自行研制的全息打印機將全息圖輸出到光刻膠版上。打印機的分辨率為318 nm。物空間內(nèi)有效樣點數(shù)為520點×515 點×41 層。物體的采樣間隔為0.09 mm，積分項一個周期的采樣點數(shù)為4927×4927。將打印好的全息圖放在如圖4所示的再現(xiàn)光路中觀察再現(xiàn)像，光源選擇波長為473 nm 的激光，通過擴束和濾波后得到球面波照射到全息圖上。拍攝實驗結(jié)果所用相機為華為榮耀20 手機所搭載的相機。

圖4 再現(xiàn)光路圖Fig.4 The diagram of optical setup for reconstruction

圖5（a）為全息圖的局部，圖5（b）和（c）分別為傳統(tǒng)算法與快速算法的再現(xiàn)效果，未觀察到二者的差異。圖6 為相機在三個不同視角下拍攝到的再現(xiàn)像。紅色框標示區(qū)域中，物體與棋盤格之間的遮擋關系正確。當多角度連續(xù)拍攝或人眼直接平移觀察時可以看到連續(xù)的視差變化。

圖5 部分全息圖及兩種算法的再現(xiàn)結(jié)果Fig.5 Partial hologram and reconstruction results of two algorithms

圖6 不同視角下的再現(xiàn)像Fig.6 Reconstructed image from different viewing-angle

觀察不同深度的再現(xiàn)像時，在圖4 光路基礎上在全息圖右側(cè)加入焦距為30 cm 的透鏡，將再現(xiàn)像成像為等大的實像并用毛玻璃承接。經(jīng)過測量得知加入透鏡后整個成像系統(tǒng)在實像位置附近的景深約為14 mm，因此選擇間隔大于景深，且特征較為明顯的兩個深度進行拍攝。圖7 為兩個不同深度處的再現(xiàn)像。圖中深度為換算成虛像的深度，再現(xiàn)的深度信息準確。

圖7 不同深度再現(xiàn)像對比Fig.7 Comparison of the reconstructed image in different depth

用層析法計算三維物體的全息圖時，再現(xiàn)像的質(zhì)量與數(shù)據(jù)分層間隔大小密切相關。層間隔滿足人眼分辨率的情況下所觀察到的再現(xiàn)像在各個視角下都是連續(xù)的，當不滿足人眼分辨率時，大角度傾斜觀察容易出現(xiàn)再現(xiàn)像的分層。在其他參數(shù)不變的前提下，將三維物體拉近到距離全息面100 mm 附近計算、制作全息圖，在右邊大角度傾斜觀察得到如圖8所示結(jié)果。由于觀察距離變小，此位置上人眼分辨率更高，再現(xiàn)像在大視角下觀察出現(xiàn)了明顯的分層。要避免這一現(xiàn)象的出現(xiàn)需要根據(jù)人眼分辨率、物體的三維形貌、觀察視角等選擇合理的層間隔。

圖8 再現(xiàn)像分層現(xiàn)象Fig.8 The slice phenomenon of reconstructed image

為研究快速算法的計算速度提高情況，選用上述三維數(shù)據(jù)模型在同一PC 機上用兩種算法計算其菲涅爾全息圖，取50 次運行時間的平均值。在全息面采樣間隔318 nm 時，積分項一個周期的采樣點數(shù)為4 927×4 927，當周期數(shù)在3 以上后傳統(tǒng)算法計算過程中會由于計算機物理內(nèi)存不足而需要讀取硬盤，計算速度會明顯降低，此時的計算速度不具可比性。因此采用較小的積分項周期進行計算，全息面的采樣間隔為4.03 μm（實驗室中SLM 的像素尺寸），物體距離全息面900 mm。根據(jù)1.1 節(jié)所述求出全息面上積分項一個周期內(nèi)包含的像素數(shù)為1 174×1 174。計算了九組周期數(shù)從kx=2，ky=2 到kx=10，ky=10 的全息圖。每一組分別采用快速算法周期延拓加系數(shù)行列分解（方法1）、快速算法僅積分項周期延拓（方法2）、及傳統(tǒng)算法（方法3）三種方法計算，得到如表1所示的結(jié)果。圖9 為不同算法計算時間與積分項周期數(shù)關系曲線。由結(jié)果可見系數(shù)行列分解后的快速算法速度最快，僅周期延拓的快速計算次之，且均快于傳統(tǒng)算法。整體來看所計算的全息圖越大時快速算法的速度優(yōu)勢越明顯，當全息圖分辨率達到11 740×11 740 像素時計算速度達到傳統(tǒng)算法的13 倍。

圖9 計算時間與周期數(shù)關系Fig.9 The relationship between the computing time and number of period

表1 不同算法計算時間對比Table 1 Comparison of computational time of different algorithm

為探究數(shù)據(jù)層數(shù)對速度提高倍數(shù)的影響，進行了層數(shù)從10 到50，以5 層為增量的全息圖計算速度提高倍數(shù)對比實驗。每組實驗用兩種算法計算周期數(shù)從kx=2，ky=2 到kx=8，ky=8 的全息圖以及提速倍數(shù)，得到了分層數(shù)對快速算法的速度提高倍數(shù)無影響的結(jié)論。表2 給出了其中5 組數(shù)據(jù)。

表2 層數(shù)與速度提高倍數(shù)關系Table 2 Relationship between number of layers and speed increase multiple

2.2 分析與討論

由實驗結(jié)果可知，提出的算法在保證良好再現(xiàn)效果的同時顯著提高了計算機制菲涅爾全息圖的計算速度。在大尺寸全息圖的計算上有著良好的表現(xiàn)。另外需要說明的是在表1 傳統(tǒng)算法用時的數(shù)據(jù)中，全息圖像素數(shù)增加到9 392×9 392 后計算用時會發(fā)生激增。這是由于傳統(tǒng)層析法計算全息圖時每一層數(shù)據(jù)都要插值或補零到與全息圖相同的像素。計算量隨著全息圖尺寸的增加而呈平方式的增長，從而導致計算用時激增。相比于傳統(tǒng)算法，快速算法在計算完成每一層在全息面上一個周期的積分項后，即可通過周期延拓來得到整個全息面的積分項。計算量隨著全息圖尺寸增大而平緩增長。

所提算法在計算完成積分項的一個周期后，后續(xù)計算互相獨立，使得采用并行計算進一步提高全息圖的計算速度得以實現(xiàn)。優(yōu)化算法，采用多核CPU 或GPU 并行計算加速全息圖的制作將是下一步的研究目標。

3 結(jié)論

計算全息數(shù)據(jù)量巨大。本文提出的快速算法，從離散菲涅爾衍射計算模型出發(fā)分析，考慮人眼分辨率，從系數(shù)和積分項兩方面改進。將系數(shù)項行列正交分解再合成；只計算部分區(qū)域的積分項，再周期延拓。該算法所得全息圖再現(xiàn)像質(zhì)量與傳統(tǒng)算法相當，同時減少了大量的冗余計算，計算速度隨全息圖尺寸的增大而快速提高。