面向真彩色三維顯示的分層角譜算法和Gerchberg-Saxton算法研究?

范爽 張亞萍 王帆 高云龍 錢曉凡 張永安 許蔚 曹良才

1)(昆明理工大學理學院,昆明 650500)

2)(昆明理工大學建筑工程學院,昆明 650500)

3)(清華大學精密儀器系,精密測試技術及儀器國家重點實驗室,北京 100084)

1 引 言

計算全息較之傳統(tǒng)光學全息省去了光源、記錄介質(zhì)和光學平臺等硬件,噪聲低,可復制性高,可計算虛擬物體,具有極大的靈活性[1?6],可為理想的真三維顯示技術提供技術支持.計算三維物體的全息圖的基本算法主要有:點源法[7]、面元法[8,9]、層析法[10]以及在頻域中進行計算的角譜法[11,12]等.目前的研究主要致力于面向復雜形貌三維物體提高再現(xiàn)像質(zhì)量和計算速度.

在改善全息圖重建像質(zhì)量方面,2004年,北京理工大學謝敬輝等[10]提出了菲涅耳波帶板掃描全息術,成功地計算了三維物體的全息圖,并通過數(shù)字高通濾波抑制了背景噪聲,提高了成像質(zhì)量;2005年,華沙理工大學Michal等[13]提出用迭代菲涅耳算法計算三平面的相位型全息圖,計算簡單物體,迭代10000次,再現(xiàn)像存在2.5倍噪聲;2014年,北京郵電大學曹雪梅等[14]提出利用二維彩色圖像和深度圖生成菲涅耳計算全息圖,強度疊加80次,成功地進行了復雜彩色三維場景的全息記錄和再現(xiàn);2015年,東南大學Chang等[15]提出像面上的振幅和相位被理想的振幅和統(tǒng)一的相位所取代的雙重約束的Gerchberg-Saxton(GS)算法,與傳統(tǒng)的GS算法相比,相同迭代次數(shù)下,同一物體重建像的散斑噪聲對比度降低41.7%,提高了成像質(zhì)量;2017年,中國科學院光電技術研究所Pang等[16]在此基礎上提出在物面和全息面之間加入兩倍物面大小的虛擬媒介平面的分層迭代算法,計算簡單物體的再現(xiàn)像,其相位的均方根誤差(RMSE)在迭代800次后可近似收斂為0,再現(xiàn)像質(zhì)量得到了大幅提升.在提高計算速度方面,2015年,清華大學Zhao等[17]提出分層角譜算法,同一硬件條件下計算同一復雜三維物體,與傳統(tǒng)點源法相比,計算負載大大降低,計算速度提高約860倍.

本文基于分層角譜算法分析了在真彩色物體計算方面的多波長采樣,實現(xiàn)了彩色物體的分層角譜算法,并得到清晰的再現(xiàn)像;同時,為了消除算法實施過程中嚴重的散斑噪聲,引入GS算法[18],提出基于分層角譜的GS算法,有效地抑制散斑噪聲,提高再現(xiàn)像質(zhì)量,并通過計算引入GS算法前后再現(xiàn)像的RMSE和峰值信噪比(PSNR)加以證明.同時,相比于雙重約束的GS算法[15]和分層迭代算法[16],該算法編碼復雜度更低,對復雜形貌三維物體的計算可有效提高算法速度.

2 真彩色物體的分層角譜算法

2.1 分層角譜算法的數(shù)學模型推導

分層角譜算法的基本物理模型是將三維物體沿深度方向離散成s個相互平行的平面,各層物面需附加隨機相位以防止物光波頻譜信息丟失,再將每一層物面通過平面波角譜理論衍射到達全息面后進行疊加,即可得到整個三維物體的物光場,進而求得全息圖,物理過程如圖1所示,其中hj是第j層物面的全息圖.

圖1 分層角譜算法的原理[17]Fig.1.Principle of angular-spectrum layer-oriented method[17].

假設三維物光場的復振幅為

其中,(x0,y0,z)是物光場的坐標,j是層的序號,s是層的總數(shù),aj為第j層的振幅,zj為第j層物面的深度值,k是波數(shù),φj是第j層物面的隨機相位.

平面波角譜算法的傳遞函數(shù)為

在平面波角譜理論中,平面波在自由空間中傳播時,其波面形狀不發(fā)生改變,僅產(chǎn)生一個與傳播距離有關的相移[19].每一層衍射物面的角譜函數(shù)乘以一個對應的與傳播距離相關的傳遞函數(shù),得到該衍射物面在全息面上的角譜函數(shù)為

將進行傅里葉逆變換得到該衍射物面的全息圖,所有物面的全息圖疊加得到整個三維物光場的全息圖,其復振幅為

其中,(x1,y1)是全息面的坐標,?表示二維傅里葉變換,??1是二維傅里葉逆變換.

2.2 真彩色分層角譜算法的仿真計算

分層角譜算法主要根據(jù)三維物體深度圖的灰度數(shù)據(jù)進行分層,首先將實際采集或軟件生成的三維物體,通過計算機圖形學方法得到其強度圖和深度圖;或三維物體的二維平面圖像通過立體匹配技術得到圖像的視差,再利用視差圖得到相應的深度圖[14].深度圖的灰度變化表示三維物體的深度變化,8 bit深度圖像即存在28=256個灰度值,物體可分為256層,每一層均采用平面波衍射角譜理論進行計算得到子全息圖,將子全息圖疊加成整個物體的全息圖[17],最后加入?yún)⒖脊庵亟ㄈS物光場.

在實現(xiàn)真彩色分層角譜算法時,由于紅、綠、藍三色光的波長不同,為了避免混頻現(xiàn)象,采樣首先要滿足奈奎斯特采樣定理.

一維傳遞函數(shù)為

其相位分布的局域空間頻率表示為[17]

采樣間隔應滿足奈奎斯特采樣定理[12],得到

設物面尺寸L=?x·N,?x為采樣間隔,N為采樣點數(shù),當最大頻率fxmax=N/(2L)時,有效衍射距離為

根據(jù)(8)式的計算可知,當各層物面的采樣間隔?x和采樣點數(shù)N都相同且確定時,各層物面的尺寸L也相同且確定,這時為避免混頻現(xiàn)象,有效衍射距離的取值應為λ=λR=625 nm時的值.

利用計算機圖形學繪制一套彩色茶具,并通過渲染技術得到茶具強度圖和深度圖,分別對應于茶具的振幅信息和深度信息,如圖2(a)和圖2(b)所示.

首先,根據(jù)深度圖的灰度值范圍將茶壺和茶杯分為兩個大的層面,每層物面都附加隨機相位,當深度圖和強度圖的?x=8μm,N=1920時,L=15.36 mm,所以設遠處茶杯的近處茶壺的然后,設紅、綠、藍三個波長分別為λR=625 nm,λG=532 nm,λB=473 nm.將強度圖分解為R,G,B三基色圖像,如圖3(a)—(c)所示.

將三基色圖分別采用分層角譜算法計算對應的相位型全息圖,如圖4所示.

圖2 彩色茶具 (a)強度圖;(b)深度圖Fig.2.Colorful tea set:(a)Intensity map;(b)depth map.

圖3 茶具的三基色圖 (a)R基色圖;(b)G基色圖;(c)B基色圖Fig.3.Three color images of tea set:(a)R-based color map;(b)G-based color map;(c)B-based color map.

最后,選擇垂直平面參考光R=1,分別對每個相位全息圖進行模擬重構,生成對應顏色分量的重構結果,使用圖像處理的方法將其合成彩色圖像[20],其數(shù)值模擬結果如圖5所示.

由圖5(a)和圖5(b)所示的重構結果得出該方法可成功地對復雜彩色三維場景進行全息記錄與再現(xiàn).相較于文獻[14]中的方法而言,該方法不用多次進行強度疊加即可得到較為清晰的再現(xiàn)像,提高了計算速度;同時,平面波角譜理論的應用避免了單次菲涅耳衍射過程中波長不同導致的衍射角大小不同,最終影響單色全息重構像的尺寸大小,也避免了菲涅耳衍射存在的傍軸近似和大數(shù)值孔徑系統(tǒng)中較為嚴重的計算誤差.但是,三維物體深度圖的準確性會直接影響再現(xiàn)像的準確性[21],且在實驗中仍要考慮波長不同引起的色差,一般采用消色差傅里葉透鏡解決這類問題.從圖5(a)和圖5(b)可以看出,由于隨機不穩(wěn)定的相位分布導致臨近采樣點的相位起伏波動大,再現(xiàn)像存在一些散斑噪聲.

為抑制散斑噪聲可以將物光和會聚光結合,但生成的全息圖中會有邊緣振蕩效應,分辨率低[16].所以,本文進一步引入GS算法,提出基于分層角譜的GS算法的模型,用于計算復雜形貌的三維物體,有效地抑制散斑噪聲,降低算法復雜度.

圖4 茶具的三基色圖對應的全息圖 (a)R通道全息圖;(b)G通道全息圖;(c)B通道全息圖Fig.4.Corresponding holograms of three color images of tea set:(a)R-channel hologram;(b)G-channel hologram;(c)B-channel hologram.

圖5 彩色分層角譜算法的數(shù)值模擬結果 (a)在180 mm處再現(xiàn)像;(b)在240 mm處再現(xiàn)像Fig.5.Numerical simulation results of colorful angular-spectrum layer-oriented algorithm:(a)The reconstruction of image at 180 mm;(b)the reconstruction of image at 240 mm.

3 基于分層角譜算法的GS算法

傳統(tǒng)GS算法[22?24]的目標物體是二維平面圖像.首先,通過傅里葉變換生成目標圖像的全息圖.然后,全息圖的相位不變,振幅用常數(shù)調(diào)制得到調(diào)制后的全息圖,該全息圖通過逆傅里葉變換得到再現(xiàn)像;再現(xiàn)像的相位不變,振幅用目標圖像振幅調(diào)制得到調(diào)制后的再現(xiàn)像,該再現(xiàn)像再通過傅里葉變換得到全息圖,如此構成一個循環(huán).最后,當?shù)螖?shù)達到設定的值時,輸出全息圖的相位,再現(xiàn)原物體.

三維GS算法是在傳統(tǒng)GS算法的基礎上,將一個二維物平面變?yōu)橛刹煌疃鹊亩S物平面組成的三維物體,因為三維物體涉及深度信息,所以將迭代中的傅里葉變換和逆傅里葉變換改成與深度有關的菲涅耳衍射和逆菲涅耳衍射[25,26],如圖6所示.

根據(jù)三維GS算法的原理,將其與分層角譜算法結合,光在全息面與物面之間傳播可以用平面波角譜衍射代替菲涅耳衍射,分層角譜GS算法的流程圖如圖7所示.

圖6 三維GS算法Fig.6.Three-dimensional GS algorithm.

圖7 基于分層角譜法的GS算法Fig.7.GS algorithm based on angular-spectrum layer-oriented method.

首先設第n次循環(huán)中全息圖的復振幅為Un(x1,y1),當n=1時,U1(x1,y1)=U(x1,y1).取第n次循環(huán)中全息圖的相位θn:

其中,angle{·}表示求相位的操作.用常數(shù)1替換第n次循環(huán)中全息圖的振幅,得到調(diào)制后的第n次循環(huán)的全息圖:

再將調(diào)制后的全息圖通過平面波角譜衍射傳到各層物面,得到第n次循環(huán)中各層物面的復振幅,其中第j層物面的復振幅為

其中,是第n次循環(huán)中第j層物面的振幅,為第n次循環(huán)中第j層物面的相位.用原始物體的第j層物面的振幅aj替換得到調(diào)制后的第n次循環(huán)中第j層物面的復振幅:

最后,每層物面再通過平面波角譜衍射傳到全息面并疊加,得到第n+1次循環(huán)中的全息圖,其復振幅表示為

當n=N時終止循環(huán),N為迭代次數(shù).

將三維火車頭通過計算機渲染技術得到強度圖和深度圖,如圖8(a)和圖8(b)所示.

基于圖8(a)和圖8(b),首先將三維火車頭的強度圖分解為R,G,B三基色圖,如圖9(a)—(c)所示.

然后,分別用傳統(tǒng)的分層角譜算法和基于分層角譜的GS算法計算三維火車頭的三基色全息圖,得到圖10所示結果,其中基于分層角譜的GS算法迭代5000次.

根據(jù)圖10所示的全息圖分別重建三維火車頭,結果如圖11所示.

圖8 三維火車頭 (a)強度圖;(b)深度圖Fig.8.Three-dimensional locomotive:(a)Intensity map;(b)depth map.

圖9 火車頭的三基色圖 (a)R基色圖;(b)G基色圖;(c)B基色圖Fig.9.Three color images of three-dimensional locomotive:(a)R-based color map;(b)G-based color map;(c)B-based color map.

圖10 采用兩種算法計算火車頭的三基色圖對應的全息圖 (a)R通道全息圖(傳統(tǒng)的分層角譜算法);(b)G通道全息圖(傳統(tǒng)的分層角譜算法);(c)B通道全息圖(傳統(tǒng)的分層角譜算法);(d)R通道全息圖(基于分層角譜算法的GS算法);(e)G通道全息圖(基于分層角譜算法的GS算法);(f)B通道全息圖(基于分層角譜算法的GS算法)Fig.10.Corresponding holograms of three color images of tea set:(a)R-channel hologram(traditional angularspectrum layer-oriented);(b)G-channel hologram(traditional angular-spectrum layer-oriented);(c)B-channel hologram(traditional angular-spectrum layer-oriented);(d)R-channel hologram(the GS algorithm based on angularspectrum layer-oriented method,iteration 5000 times);(e)G-channel hologram(the GS algorithm based on angularspectrum layer-oriented method,iteration 5000 times);(f)B-channel hologram(the GS algorithm based on angularspectrum layer-oriented method,iteration 5000 times).

圖11 傳統(tǒng)分層角譜算法與基于分層角譜的GS算法的再現(xiàn)像的比較 (a)在180 mm處再現(xiàn)像(分層角譜);(b)在180 mm處再現(xiàn)像(基于分層角譜的GS算法);(c)在210 mm處再現(xiàn)像(分層角譜);(d)在210 mm處再現(xiàn)像(基于分層角譜的GS算法);(e)在240 mm處再現(xiàn)像(分層角譜);(f)在240 mm處再現(xiàn)像(基于分層角譜的GS算法)Fig.11.Comparison of reconstructive images between angular-spectrum layer-oriented method and the GS algorithm based on angular-spectrum layer-oriented method:(a)The reconstruction of image at 180 mm(angularspectrum layer-oriented);(b)the reconstructionof image at 180 mm(the GS algorithm based on angular-spectrum layer-oriented method);(c)the reconstruction of image at 210 mm(angular-spectrum layer-oriented);(d)the reconstruction of image at 210 mm(the GS algorithm based on angular-spectrum layer-oriented method);(e)the reconstruction of image at 240 mm(angular-spectrum layer-oriented);(f)the reconstruction of image at 240 mm(the GS algorithm based on angular-spectrum layer-oriented method).

從圖11(a)與圖11(b)、圖11(c)與圖11(d)、圖11(e)與圖11(f)之間的比較可以看出,傳統(tǒng)的分層角譜算法和基于分層角譜的GS算法計算得到的再現(xiàn)像較之原物體都存在一定的顏色失真,但傳統(tǒng)的分層角譜算法計算得到的再現(xiàn)像,其背景及物體本身的散斑噪聲較大,較為模糊;基于分層角譜的GS算法的再現(xiàn)像散斑噪聲較少,更為清晰.為直觀地證明兩種算法再現(xiàn)像的差異,計算兩種算法的RMSE和PSNR來進行質(zhì)量評估,RMSE越小,再現(xiàn)像質(zhì)量越好;PSNR越大,再現(xiàn)像質(zhì)量越好.根據(jù)文獻[16],RMSE的公式為

其中是第j層重建像的振幅分布,是第j層原始物體的振幅分布.根據(jù)文獻[17],PSNR的公式為

其中,M,N分別是水平和豎直方向的像素數(shù);I0,Ir分別為原始圖像和再現(xiàn)像的光強分布.繪制出RMSE和PSNR隨著迭代次數(shù)的變化曲線,如圖12所示.

圖12 分層角譜算法與基于分層角譜的GS算法的誤差比較(a)RMSE;(b)PSNRFig.12.Error comparison between angular-spectrum layeroriented method and GS algorithm based on angularspectrum layer-oriented method:(a)RMSE;(b)PSNR.

在圖12中,當?shù)螖?shù)為0即未進行GS迭代時,RMSE和PSNR是利用傳統(tǒng)的分層角譜算法計算得到的,當開始進行GS迭代時,利用基于分層角譜的GS算法計算再現(xiàn)像的RMSE和PSNR.在5000次迭代內(nèi)RMSE和PSNR最終都趨近于收斂,且圖12(a)中RMSE隨著迭代次數(shù)的增加而減小,在前100次迭代內(nèi)迅速減小,圖12(b)中PSNR隨著迭代次數(shù)的增加而增大,在前100次迭代內(nèi)增加速度快,該結果證明了基于分層角譜的GS算法收斂速度快,再現(xiàn)像質(zhì)量更好.但圖12(a)中的RMSE和圖12(b)中的PSNR在某個幾百次迭代內(nèi)均有小的起伏,RMSE的起伏范圍是±0.02,PSNR的起伏范圍±0.08,這是GS算法的局部收斂特性引起.基于分層角譜的GS算法計算5000次迭代內(nèi)火車頭的再現(xiàn)像,其頭部的RMSE近似收斂到0.77,PSNR近似收斂到65.7,火車頭中部和尾部的RMSE近似收斂到0.68,PSNR近似收斂到70.0,即火車頭頭部再現(xiàn)像質(zhì)量明顯低于中部和尾部;同傳統(tǒng)分層角譜算法的再現(xiàn)像相比,其再現(xiàn)像的頭部、中部和尾部的RMSE分別減小了約0.11,0.40和0.41,PSNR分別增大了約1.15,5.70和4.13,證明了基于分層角譜的GS算法可以對細節(jié)信息多的復雜物體進行高質(zhì)量地重建,有效地抑制散斑噪聲.

本文提出的基于分層角譜的GS算法的模型結合了分層角譜算法和GS算法的優(yōu)點.算法的復雜度只與分層層數(shù)相關,與物體的復雜程度無關,可以進行復雜三維形貌物體的計算;避免了傍軸近似帶來的計算誤差;且收斂速度快,有效地抑制散斑噪聲,提高成像質(zhì)量.

4 結 論

本文分析了分層角譜算法在計算真彩色物體全息圖時的多波長采樣問題,并分別采用分層角譜算法和基于分層角譜的GS算法計算真彩色物體的全息圖和再現(xiàn)像.通過計算RMSE和PSNR評估再現(xiàn)像的質(zhì)量,驗證本文提出的基于分層角譜的GS算法可有效地抑制散斑噪聲,提高再現(xiàn)像質(zhì)量.在本文算法中,迭代次數(shù)越多,成像質(zhì)量越好,散斑越少,但計算時間增加.所以,在實際高精度的全息計算過程中,可采用GPU等硬件[27?29],并行處理模式、大量的線程運算等結合[30,31]的方法,進行高速計算全息圖并高質(zhì)量成像.

本文研究工作承蒙美國弗吉尼亞理工暨州立大學(Virginia Tech)Ting-Chung Poon教授指導修改完成,謹致謝意！同時,研究工作得到安徽大學韋穗教授的指點和幫助,在此一并表示感謝!

參考文獻

[1]Goodman J W,Lawrence R W 1967Appl.Phys.Lett.11 77

[2]Adam K,David L K 1965Appl.Opt.4 387

[3]Qian X F 2015Information Optical Digital Laboratory(Beijing:Science Press)p200(in Chinese)[錢曉凡 2015信息光學數(shù)字實驗室(北京:科學出版社)第200頁]

[4]Gabor D 1948Nature161 777

[5]Leith E N,Upatnieks J 1963J.Opt.Soc.Am.53 1377

[6]Tan F 2012Opt.Instrum.34 16(in Chinese)[覃芳 2012光學儀器34 16]

[7]Waters J P 1966Appl.Phys.Lett.9 405

[8]Matsushima K,Nakahara S 2009Appl.Opt.48 H54

[9]Pan Y J,Wang Y T,Liu J,Li X,Jia J 2014Appl.Opt.53 1354

[10]Sun P,Xie J H,Zhou Y L 2004Acta Opt.Sin.24 110(in Chinese)[孫萍,謝敬輝,周元林2004光學學報24 110]

[11]Li J C 2014Di ff raction Calculations and Digital Holography(Vol.1)(Beijing:Science Press)p261(in Chinese)[李俊昌2014衍射計算及數(shù)字全息(上冊)(北京:科學出版社)第261頁]

[12]Poon T C,Liu J P 2014Introduction to Modern Digital Holography with MATLAB(London:Cambridge University Press)p5

[13]Michal M,Maciei S,Andrzei K,Grzegorz M 2005Opt.Eng.44 125805

[14]Cao X M,Sang X Z,Chen Z D,Leng J M,Zhang M,Guo N,Yu C X,Xu D X 2014Chin.J.Lasers41 0609002(in Chinese)[曹雪梅,桑新柱,陳志東,冷俊敏,張明,郭南,余重秀,徐大雄2014中國激光41 0609002]

[15]Chang C L,Xia J,Yang L,Lei W,Yang Z M,Chen J H 2015Appl.Opt.54 6994

[16]Pang H,Wang J Z,Cao A X,Zhang M,Shi L F,Deng Q L 2017IEEE Photon.J.9 1

[17]Zhao Y,Cao L C,Zhang H,Kong D Z,Jin G F 2015Opt.Express23 25440

[18]Gerchberg R W,Saxton W O 1972Optik35 237

[19]Xie J H,Liao N F,Cao L C 2007Fundamentals of Fourier Optics and Contemporary Optics(Beijing:Beijing Polytechnic University Press)p79(in Chinese)[謝景輝,廖寧放,曹良才2007傅里葉光學與現(xiàn)代光學基礎(北京:北京理工大學出版社)第79頁]

[20]Shen C,Wei S,Liu K F,Zhang F,Li H,Wang Y 2014Laser Optoelectr.Prog.51 030005(in Chinese)[沈川,韋穗,劉凱峰,張芬,李浩,王岳2014激光與光電子學進展51 030005]

[21]Piao Y L,Kwon K C,Jeong J S,Kim N 20163D Image Acquisition and Display:Technology,Perception and ApplicationsHeidelberg,Germany,July 25–28,2016 JW4A.29

[22]Peng J M,Du S J,Jiang P Z 2013High Power Laser and Particle Beams25 315(in Chinese)[彭金錳,杜少軍,蔣鵬志2013強激光與粒子束25 315]

[23]Gu X,Xu K S 2000J.Fudan Univ.(Natural Science)39 205(in Chinese)[顧翔,徐克璹2000復旦學報(自然科學版)39 205]

[24]Wang J Y,Cao J H 2015J.Tianjin Univ.Technol.Education25 36(in Chinese)[王金洋,曹繼華2015天津職業(yè)技術師范大學學報25 36]

[25]Li F,Bi Y,Wang H,Sun M Y,Kong X X 2012Chin.J.Lasers39 1009001

[26]Zhou P H,Bi Y,Sun M Y,Wang H,Li F,Qi Y 2014Appl.Opt.53 G209

[27]Pan Y C,Xu X W,Liang X N 2013Appl.Opt.52 6562

[28]Liu J,Ma X,Wang Y T,Jia J,Zhang Y X 2015 CN104281490A(in Chinese)[劉娟,馬曉,王涌天,賈甲,張迎曦2015 CN104281490A]

[29]Kwon M W,Kim S C,Yoon S E,Kim E S 2015Opt.Express23 2101

[30]Chen H R,Fu S H,Wang Y Q 2014Opto-Electron.Eng.41 48(in Chinese)[陳慧榮,付勝豪,王元慶2014光電工程41 48]

[31]Fu S H,Wang Y Q,Bao X L,Fan K F 2013Electron.Opt.Control20 61(in Chinese)[付勝豪,王元慶,鮑緒良,范科峰2013電光與控制20 61]