面向瞬态成像仿真的遮挡目标逆向重构方法

李海瑞， 吴　鑫， 张建奇

(西安电子科技大学 物理与光电工程学院，陕西 西安 710071)

瞬态成像(Transient Imaging，TI)是一种新颖的成像方式．瞬态成像源于上世纪70年代末由文献[1]完成的“飞光记录”(Light-in-Flight)研究，其提出了利用皮秒级激光照明场景并进行全息记录，能够对某一特定时刻的波前进行光学重建．尽管重建场景的复杂度受全息设备性能的限制，但其他研究者已经利用这一方法进行了可见场景三维几何结构测量的工作[2]，形成了瞬态成像的概念．在瞬态成像过程中，短脉冲激光在遍历完场景并达到光能分布平衡之前，就会以很高的时间分辨率被探测到．因此，得到的瞬态图像是一系列表示场景对脉冲响应的时序图像．

随着超快成像器件的高速发展，一个更简单高效的成像方式被提了出来，被称为瞬态成像框架(Transient Imaging Framework)[3]．该框架的成像设备由一个飞秒级激光器和一个条纹相机组成，其中飞秒级激光器作为主动脉冲照明源，条纹相机用于捕获场景对脉冲响应的时序图像．这一项新的成像方式催生出了很多新的应用，如重建被遮挡场景的几何结构、获得物体的表面反射率、复杂场景中光传输的可视化[4]以及透明物体的深度探测等．遮挡重构是一个重要的研究领域，图像的形成过程十分复杂，既包含可见物体的直接反射，还包括遮蔽物体的间接散射．如何把这些反射分离开来，是一个巨大的挑战．目前的重构方案仍然是以文献[5-7]所提出的方法为基础，他们利用一个严格控制的实验室设备演示了遮挡重构的过程．尽管重构结果没有使用到场景先验假设，但是却依赖于一个昂贵的定制设备，从而具有很大的局限性．近年来，其他研究者相继提出了两种更廉价的瞬态成像设备[8-9]，尽管可能牺牲一定的时间和空间分辨率，却对推动瞬态成像技术的商业化应用具有重要的意义．由于瞬态成像的基础理论尚未完全建立，目前的工作中仍然存在一些问题难以回答．文献[10-11]利用信号处理方面的理论对重构过程进行了详细分析，指出可以利用信号处理方面的理论和方法提高重构质量．

图1 激光脉冲的多次反射传输模型

笔者利用文献[12]提出的一个瞬态渲染框架获得场景响应的仿真数据．针对瞬态成像中实际硬件水平的不足和实验的复杂性，改进了原有框架的逆向重构算法．利用软件仿真的手段将瞬态成像过程进行可视化分析，提高了目标重构的效果．

1　逆向重构算法的简化

图1是描述了一束激光在场景中传输的物理模型．首先，由激光器B发出的一束脉冲光到达散射墙面R上的点L；然后，光从点L处向四周散射，其中有一小部分沿路径rl到达了未知表面S上的点s；接着，光从点s处再次向四周散射，其中有一小部分沿路径rc到达了散射墙面R上的点w；最后，光从点w又一次向四周散射，其中一小部分到达条纹相机C中的像素p上，并被探测到．逆向重构的任务是通过对反射光线的计算，获得被遮挡表面S的几何信息，又称为绕射探测或遮挡成像．

1.1　逆向重构算法

逆向重构算法包括两个计算过程: 工作空间体素化和逆向投影．假设在空间中存在一个完全包围未知表面S的长方体空间，记为工作空间V，把工作空间V均匀地划分为一系列边长相等的立方体，并称之为体素，记为v．记激光器发射激光脉冲为时间原点，t=0．由图1可知，在t时刻，只有满足以下条件的体素v才对点w处有光强贡献:

ct=|L-B|+|v-L|+|w-v|+|p-w|．

(1)

由于散射墙面R是可见场景，条纹相机C和激光器B的设置已知，因此，变量ct，|L-B|，|p-w|都是可知的．所以有

|v-L|+|w-v|=d,d=ct-|L-B|-|p-w|．

(2)

这是一个以点L和点w为焦点、长轴为d的椭球面，称满足该式关系的p,w,v为互相贡献关系．

对于工作空间V中的每一个体素v，按如下方式计算一个关于v的像素累积函数H:

(3)

其中，p是所有与体素v为互相贡献关系的像素；Ip为像素p处探测到的光强；a为一常数因子，通常取a=1．

1.2　逆向重构算法的简化

将像素累积函数式(3)简化为

(4)

即对于每个体素v，其所有的贡献像素p对函数H的贡献值相同，皆取值为1．这意味着重构过程中Ip的大小并不重要，只需记录体素在光的传播过程有无响应．因此，在重构过程中就不必反复计算探测到光强的大小，而只需要判断能否探测到光强即可，从而可以极大地简化重构过程的复杂度，提高了算法的处理效率．

从式(2)可以看出，在t时刻，对点w有光强贡献的体素不是惟一的．式(2)确定了一个可能体素的候选范围，即由式(2)确定的一个椭球面．假设光的衰减只与其传播距离有关，即满足Ie=Is/r的关系．图2(a)描述了当点s在由点L和点w确定的椭圆轨迹上移动时，光经过光路rl和rc的衰减程度 |rl| |rc| 是不同的．但是点L处的光强IL却是相同的，因而到达点w的光强Ip=IL/ (|rl| |rc|) 就会不同．利用式(3)计算，恰好消除了经过不同光路衰减程度不同的问题，得到相同贡献的光强IL．因此，在这一点上，式(4)和式(3)是完全等价的．

图2　距离积分布图

但是实际上由式(2)确定的椭球面上只有一个点s对点w有光强贡献，即光路是惟一的，光强的衰减程度是惟一的．但当按照式(3)对式(2)确定的椭球面上所有的可能点s计算时，同样大小的Ip却因为体素v的空间位置不同，而对H的贡献大小不同．在xy平面中取两个定点F1,F2，其坐标分别为 (-c, 0)，(c,0)，这里取c=1．设点f为平面内的任意一点．按下式可以得到一个距离等积曲线:

|f-F1| |f-F2|=a1，

(5)

其中，a1为一个常数．当a1从0逐渐增大时，距离等积曲线的分布如图2(b)所示．由图2(b)可以看出，距离等积曲线的分布不是各个方向均匀的．当a1≈c时，靠近线段F1F2的中垂线的体素v对函数H的贡献更大; 当a1≫c时，逐渐接近均匀．这说明利用式(3)进行重构时，可能会使结果具有一定的空间差异性．

当把式(3)简化为式(4)后，重构问题就会完全归结为一个概率统计的问题．对重构问题的分析和研究会变得更加简单．

2　逆向重构算法的补充

在对逆向重构算法进行简化的基础上，笔者提出了一种与其相适应的取样模式，可以加快重构结果的收敛速度．在描述取样模式之前，首先利用检测概率[13]的概念给出逆向重构算法的概率意义．设v为工作空间V内的任意一个体素，w为散射墙面R上的任意一点，定义事件A(v,w)为在点w处检测到体素v．在这里，由式(2)和式(4)可得

(6)

设W为相机在散射墙面R上的可探测范围，w∈W，对应的函数H记为HW．对HW进行归一化处理:

(7)

图3　取样过程示意图

采用如图3所示的取样模式．

直接利用上面描述的逆向重构算法，可以得到一个目标物体的粗略重构结果，进而对目标物体的中心位置O0做出一个估计．调整图3中的圆心O，使直线OO0垂直于散射墙面R．要求 |w-L|≪r，接下来描述时，把点w和点L简化为一点，即点W．规定点W只能取图3中圆心O和圆周上的位置．以点(r,0)为第1个n等分点，按逆时针方向，顺序得到圆周上第i(i=1，2，…，n)个n等分点(要求n≥3)．记在圆周的第i个n等分点位置处按式(4)得到的函数H为Hi．特别地，圆心O处得到的函数H记为H0．

改进后的逆向重构算法可以按以下步骤实现:

(8)

(3) 阈值分割．经过汇总之后，Hcircle已经收敛了，为了获得比较好的显示效果，可以选取一个较小的阈值对结果进行分割，阈值的选取如下:

Hthresh=β·max{Hcircle(v),v∈V},

(9)

其中，β是一个常数因子．取β=0.001～0.010，可以得到较好的结果．

图4　仿真场景布局

3　重构结果与分析

图5　重构结果的三维显示

图6　重构结果的二维显示

3.1　重构结果的可视化显示

分别以汉字“西”和“电”作为待重构目标的重构结果的三维显示如图5所示．灰色点代表重构点，白色点是利用其他方法得到的原模型的点云，作为对比．重构结果的二维显示如图6所示．图6(a)中的灰度大小表示式(8)中计算的Hcircle的相对大小．图6(b)为图6(a)的二值化显示．

3.2　对重构结果的分析

记实际包含待重构物体表面点的体素点为“真点”，否则为“假点”．由式(2)可知，虽然该式把“真点”的范围从整个工作空间限制到了一个椭球面上，但椭球面上包含“真点”的同时，必然包含很多“假点”．因此，仅仅依靠式(6)虽然可以重构出“真点”，但必然引入大量的噪声．值得注意的是，由式(2)引入的“假点”受椭球面的位置和形状的影响远远大于“真点”．因此可以通过改变椭球面的形状和位置得到多组重构结果，然后结合这些结果达到抑制噪声的目的．上文的取样模式正是基于这一点提出的．由重构结果可知，通过该方法确实可以有效地抑制噪声，从而加快重构处理过程．需要指出的是，笔者提出的取样模式在有效地抑制噪声的同时，也会对一部分“真点”进行抑制．理论上，这种抑制几乎不会造成目标物体的信息丢失; 但实际上，由于取样量少的原因，每组均会有部分“真点”遗漏，最终造成汇总结果损失部分信息．因此，每组取样的量必须足够大，特别是如图3中圆周上的位置，需要进行多次重复取样，以减少信息丢失．

4　结 束 语

笔者改进了瞬态成像中的逆向重构算法，基于文献[12]提出的瞬态渲染框架，通过软件仿真的方法实现了被遮挡物体几何结构的重构．仿真结果显示，改进后的逆向重构算法不仅大大简化了处理过程，而且取得了更好的重构效果．文中的仿真方法与结果有助于深入探索光的传输过程，并且对超高速成像系统的设计与测试也有指导意义．

参考文献：

[1] ABRAMSON N. Light-in-flight Recording by Holography[J]. Optics Letters, 1978, 3(4): 121-123.

[2]NILSSON B, CARLSSON T E. Direct Three-dimensional Shape Measurement by Digital Light-in-flight Holography[J]. Applied Optics, 1998, 37(34): 7954-7959.

[3]KIRMANI A, HUTCHISON T, DAVIS J, et al. Looking around the Corner Using Ultrafast Transient Imaging[J]. International Journal of Computer Vision, 2011, 95(1): 13-28.

[4]VELTEN A, WU D, MASIA B, et al. Imaging the Propagation of Light through Scenes at Picosecond Resolution[J]. Communications of the ACM, 2016, 59(9): 79-86.

[5]VELTEN A, WILLWACHER T, GUPTA O, et al. Recovering Three-dimensional Shape around a Corner Using Ultrafast Time-of-flight Imaging[J]. Nature Communications, 2012, 3(2): 745.

[6]GUPTA O, WILLWACHER T, VELTEN A, et al. Reconstruction of Hidden 3D Shapes Using Diffuse Reflections[J]. Optics Express, 2012, 20(17): 19096-19108.

[7]RASKAR R, VELTEN A. Methods and Apparatus for Imaging of Occluded Objects from Scattered Light: US9148649[P]. 2015.

[8]HEIDE F, HULLIN M B, GREGSON J, et al. Low-budget Transient Imaging Using Photonic Mixer Devices[J]. ACM Transactions on Graphics, 2013, 32(4): 1-10.

[9]BUTTAFAVA M, ZEMAN J, Tosi A, et al. Non-line-of-sight Imaging Using a Time-gated Single Photon Avalanche Diode[J]. Optics Express, 2015, 23(16): 20997-21011.

[10]KADAMBI A, ZHAO H, SHI B, et al. Occluded Imaging with Time-of-flight Sensors[J]. ACM Transactions on Graphics, 2016, 35(2): 2836164.

[11]BHANDARI A, RASKAR R. Signal Processing for Time-of-flight Imaging Sensor: an Introduction to Inverse Problems in Computational 3-D Imaging[J]. IEEE Signal Processing Magazine, 2016, 33(5): 45-58.

[12]JARABO A, MARCO J, MUNOZ A, et al. A Framework for Transient Rendering[J]. ACM Transactions on Graphics, 2014, 33(6): 1-10.

[13]徐军, 周翔, 梁昌洪. 一种计算红外点目标检测概率的方法[J]. 西安电子科技大学学报, 2001, 28(1): 18-21.

XU Jun, ZHOU Xiang, LIANG Changhong. Calculation of the Detecting Probability of an Infrared Point Target[J]. Journal of Xidian University, 2001, 28(1): 18-21.