OpenMP并行程序在数字全息三维重构中的应用

马 静，邸江磊，肖 锋

(1.西安工业大学 计算机科学与工程学院，陕西 西安 710032;2.西北工业大学 理学院，陕西 西安 710072)

0 引 言

数字全息术是顾德门[1]在1967年提出的一种全新的全息成像方法，以CCD等光电耦合器件取代传统的干版记录全息图，并由计算机以数字的形式对全息图进行再现。随着全息技术的发展，数字全息技术在显微三维物体表面形变测量、粒子场测试、流场测定[2-3]、生物样品显微测量等方面有了新的研究[4-5]。数字全息不需要成像透镜，灵敏度高，是一种比较理想的微小物体数字三维测量方法[6-7]。

数字全息将全息图记录的物光波和参考光波相干涉叠加时产生的干涉条纹数字化，然后在计算机中通过再现算法进行处理，从而得到物体的三维强度像和相位像。早期的单次傅里叶变换算法应用较为广泛，该算法使再现像素大小与再现距离呈正比例变化。在光学的标量衍射理论和傅里叶变换方法的基础上，又发展出了三种最常用的再现算法：菲涅耳变换法、卷积法和角谱法[8]。

文中基于数字全息技术，通过卷积算法，实现了数字全息三维重构，针对重构计算量大、重构时间过长等问题，利用OpenMP并行程序对数字全息三维重构过程进行了提速，并对其进行验证。

1 基于卷积的数字全息重构原理

数字全息重构是利用电荷耦合器件CCD获得干涉全息图后，模拟生成参考光，然后利用计算机技术和相关算法对全息图进行处理，恢复物体的振幅和相位信息，重构物体的三维模型。通过数字全息重构，可以直接得到被记录物体再现像的复振幅分布，从而获得被记录物体的表面亮度和形貌分布。而全息中采用的CCD曝光时间短，分辨率高，记录和再现过程都比经典光全息方法更快捷。

卷积法是基于快速傅里叶变换的数值再现方法[9]。利用卷积法重构图像的尺寸不随重建距离变化。衍射积分可以看作是物波函数与自由空间脉冲响应函数的卷积。

g(x',y',ξ,η)=

(1)

利用解卷积的方法就可以重构原物像，即

u'=F-1{F[h·r]·F[g]}

(2)

对于离轴全息而言，需要通过离轴物理光路，利用干涉原理将物体信息以干涉条纹的形式记录下来，通过CCD相机以图像的形式保存下来，形成全息图[10]。

设照射到全息平面上的物光波和参考光波分别为O(x,y)和C(x,y)，则全息图平面的干涉条纹强度为H(x,y)：

O(x,y)C*(x,y)+O*(x,y)C(x,y)

(3)

其中，C*(x,y)和O*(x,y)分别表示物光波和参考光波的共轭光波。

上述全息图需要通过CCD记录离散化的光强分布。设CCD的像素数为Nx×Ny，光敏面尺寸为Lx×Ly，则像元尺寸分别为Δx=Lx/Nx，Δy=Ly/Ny，CCD记录的离散图像为：

(4)

CCD记录的离散化的干涉强度即为数字全息图。

根据卷积理论，物光的复振幅U(x,y)可根据全息图平面的干涉条纹强度、参考光波以及脉冲响应函数表示为：

U(xi,yi)=F-1{F[H(x,y)C(x,y)]F[g(x,y,xi,yi)]}

(5)

其中，H(x,y)为全息图平面的干涉条纹强度；C(x,y)为参考光波；g(x,y,xi,yi)为脉冲响应函数。

脉冲响应函数与具体的物理光路有关，计算方法如下：

g(x,y,xi,yi)=

g(xi-x,yi-y)

(6)

其中，λ为波长。

通过上述方法可以得到物光的复振幅[11-13]。物光的复振幅中携带了物体的三维信息，如果物体起伏超出了波长范围，则需要进行解包裹，从而提取物体的三维分布信息。

2 数字全息三维重构过程

数字全息三维重构过程如图1所示。数字全息处理的对象为重构硬件系统采集的全息图。根据重构硬件系统组成的不同，对全息图的处理也不同。

图1 数字全息三维重构过程

对全息图的处理需要使用参考光，如果参考光是平面光，可直接用采集到的全息图进行后续处理，如果是非平面光，则需用参考光乘以全息图。处理后的图像要进行傅里叶变换，从而得到其频谱分布图。

频谱分布图包含零级区域和正负一级区域。其中物体信息的频谱分布在正负一级中心附近，而位于图像中心最强的零级则没有任何物体信息，因此需要进行选频和移频。选频是选择正一级或负一级频谱区域，移频是将正一级或负一级频谱中心移到图像中心。这一步需要用户参与，为了便于用户进行选频操作，文中根据频谱图的特征实现了正一级或负一级频谱中心及频谱区域的自动生成。

正一级或负一级频谱中心自动生成的基本思想就是令图像的极大值为零级中心，次极大值和第三极大值为正一级或负一级频谱中心。由于数字全息三维重构的物理系统中都使用了激光，因此全息图存在大量激光散斑，从而导致次极大值有可能出现在零级中心区域附近，而不是正一级或负一级中心。为此，首先采用均值滤波器对需要显示的频谱图像进行滤波处理，然后在寻找次极大值和第三极大值，这样可以有效消除激光散斑的影响，准确识别出正负一级中心。有了零级中心、正负一级中心，就可以生成正负一级区域。

文中生成的正负一级中心及其区域，只是为用户提供初始依据，用户还可以手动对其进行修正。另外需要说明的是，文中只是对显示给用户的频谱图进行了滤波处理，实际的频谱图并没有改变。

根据脉冲响应函数得到传递函数，将传递函数与移频后的图像相乘，对其乘积进行傅里叶逆变换。逆变换的图像是一个复数图像，其相位包含了物体信息，此时相位在0～2 pi之间，需要再进行解包裹处理。

解包裹采用最小Lp范数法，这种方法不依赖路径的全局算法，计算量大，但对误差点的控制较好。最小范数法是要得到解包裹相位Φ的相邻像素相位差和包裹相位Ψ的相邻像素相位差最小范数Lp。最小Lp范数解要使式(7)中的J取到极小值。

(7)

范数p的选取是关键。当p≫2时，解包裹曲面与真实梯度面相差很大；当p<2时，结果与局部梯度较匹配，但是权重的作用变大；p=2是目前使用最多的最小范数解包裹方法，也就是最小二乘法，可将解包裹过程转化为泊松方程，再用迭代方法求解[14-15]。

解包裹后的相位代表实际物体高程，需要根据标定值得到物体相位信息。

上述三维数字全息重构算法是一个复杂耗时的过程。在重构过程中，一般需要选择合适的频谱区域，并对重构距离、波长等重构参数进行设定。

3 三维重构的并行处理

数字全息三维重构过程要进行大量的图像处理工作，运算量非常大，采用并行处理算法可加快运算速度，从而提高数字全息三维重构的实时性。

OpenMP是一个在共享存储的多处理机上编写并行程序的应用编程接口。对于同步共享变量、合理分配负载等任务都提供了有效的支持。

OpenMP的共享存储模型最底层是处理器集合，这些处理器可以访问内存中的同一个共享位置，因而可以通过共享变量进行交互和同步。OpenMP可根据需要设置包含子句项,在没有其他约束条件下，子句可以无序，也可以任意地选择。#pragma omp parallel for[子句#]是最频繁使用的编译指导语句，可搭配firstprivate，if，lastprivate，private，reduction，schedule等子句使用。

在进行数字全息三维重构时，一般有两种应用情况，一种是针对单幅图像的重构，另一种是针对视频图像的重构。一种并行设计思想是针对不同任务进行并行处理，另一种思想是针对某一图像重构的不同任务进行并行处理。为了兼顾对图像和视频的重构任务，采用第二种并行设计思想。就是对每一幅图像数字全息三维重构的各个环节采用并行算法进行处理。

根据数字全息三维重构过程，采用并行处理算法的环节包括：全息图与参考光的相乘；全息图的傅里叶变换；频谱图移频；传递函数的生成；传递函数与移频后图像的相乘；相位解包裹；解包裹后相位转换为物高以及显示图像生成等。

上述任务的共同特点是都需要对图像进行处理，需要采用二维循环来完成。因此文中采用针对循环运算的并行处理方法，在需要进行并行运算的for语句前增加以下OpenMP编译指令：

#pragma omp parallel for schedule(kind,[size])

其中kind表示模式，共有三种模式，分别为static、dynamic和guided。static模式下为每个线程分配size次循环任务，按照线程数量一次性完成任务分配。如果size数值过大，则可能导致出现有的内核任务很重，有的内核没有任务的情况。

dynamic模式也是为每个内核分配size次循环任务，但是分配过程是动态的，哪个内核处于空闲状态就给它分配size次循环任务。

guided模式按照从大到小的方式给每个内核分配任务，先分配较多任务，然后逐渐减小任务量，最小任务量为size次循环。

为了检测并行算法以及各种并行运算模式的运行速度，在配有AMD Phenom II X6 1055T六盒处理器的台式机上，对软件运行时间进行了测试。所采用全息图的分辨率为1 024×768，测试时对全息图进行400次相同的重构，计算重构时间的均值及其方差，上述过程多次测量选择中值，结果如表1所示。

表1 测量结果

从表1可以看出，未采用并行处理算法时，每次重构的平均时间为0.369 s，在三种不同的并行方法处理下，重构时间缩减到0.23 s左右，单幅重构图像时间提高约1.56倍。不同参数下，三种并行方法重构的时间相差0.002～0.004 s，采用guided方式，方差相比而言略小。

4 处理结果

根据上述重构算法，用C++语言开发了数字全息三维重构软件。选取一幅通过CCD采集的分辨率板全息图，如图2所示，该幅分辨率板全息图像的主要参数如下：波长632.8 nm；像素尺寸4.65 μm；重构距离118.025 mm。

图2 分辨率板全息图

根据上述图像参数对分辨率板全息图进行重构，重构结果如图3所示。由于分辨率板为相位型物体，因此相位图就已经反映了分辨率板的形状。图中的三维重构图是没有进行解包裹的结果。

图3 分辨率板重构图像

从频谱图中可以看出，在零级周围有一个正一级中心，还有一个负一级中心，无论是正一级还是负一级都可以重构出相位信息，其对应的重构距离一个为正，另一个为负。

5 结束语

利用OpenMP技术，根据卷积的离轴数字全息重构的基本原理，设计并开发了数字全息三维重构软件。采用并行运算提高了三维重构的运行速度，并对OpenMP实现并行运算的方式进行了对比分析，结果表明各种并行模式的平均重构时间基本相同，但是guided模式的方差相对较小。

[1] GOODMAN J W,LWRENCE R W.Digital image formulation from electronically detected holograms[J].Applied Physics Letters,1967,11(3):77-79.

[2] FLEMING C P，ECKERT J，HALPERN E F，et al.Depth resolved detection of lipid using spectroscopic optical coherence tomography[J].Biomedical Optics Express,2013,4(8):1269-1284.

[3] LUCENTE M.Interactive three-dimensional holographic displays:seeing the future in depth[J].ACM SIGGRAPH Computer Graphics,1997,31(2):63-67.

[4] DI CAPRIO G,DARDANO P,COPPOLA G,et al.Digital holographic microscopy characterization of superdirectivebeam by metamaterial[J].Optics Letters,2012,37(7):1142-1144.

[5] ALIEVA T,BASTIAANS M J.On fractional Fourier transform moments[J].IEEE Signal Processing Letters,2000,7(11):320-323.

[6] LIEBLING M,BLU T,UNSER M.Fresnelets:new multriesolution wavelet bases for digital holography[J].IEEE Transactions on Image Processing,2003,12(1):29-43.

[7] 赵 洁,王大勇，李 艳，等.数字全息显微术应用于生物样品相衬成像的实验研究[J].中国激光,2010,37(11):2906-2911.

[8] 张志会，王华英，刘佐强，等.基于快速傅里叶变换的相位解包裹算法[J].激光与光电子学进展,2012,49(12):59-65.

[9] 王华英，刘佐强，廖 薇，等.基于最小范数的四种相位解包裹算法比较[J].中国激光,2014,41(2):122-127.

[10] 郭恒光，瞿 军．基于小波域双谱分析的磨粒图像多尺度形状特征提取[J]．计算机应用与软件，2016，33(9)：224-226．

[11] 赖建新,胡长军,赵宇迪,等.OpenMP任务调度开销及负载均衡分析[J].计算机工程，2006，32(18)：58-60.

[12] 马旭龙，林 峰.基于OpenMP的快速并行分层算法[J].计算机辅助设计与图形学学报，2015，27(4):747-753.

[13] 邹贤才,李建成,汪海洪,等.OpenMP并行计算在卫星重力数据处理中的应用[J].测绘学报,2010,39(6):636-641.

[14] 奎 因.并行程序设计C、MPI与OpenMP[M].北京:清华大学出版社,2005.

[15] 潘哲郎，李仕萍，钟金钢.用数字全息层析成像技术测量毛细管的内径及壁厚[J].光学精密工程，2013，21(7):1643-1650.