基于光传播Monte Carlo模拟的三维模型特征提取方法

刘洪豪, 刘贤喜, 张开兴，2*, 卢 山， Lee Heow Pueh, 宋正河

1. 山东农业大学机械与电子工程学院， 山东 泰安 271018 2. 山东省农业装备智能化工程实验室， 山东 泰安 271018 3. Department of Mechanical Engineering, National University of Singapore, Singapore 117575, Singapore 4. 中国农业大学工学院， 北京 100083

引 言

三维模型特征提取是三维模型分类筛选的关键环节[1]， 高精度特征提取方法可为三维重构、 3D打印以及虚拟现实技术提供技术参考。 特征提取的核心是利用不同维度方法分析并挖掘模型的形状或拓扑结构， 建立与模型唯一对应且具有高区分度的代表特征。 传统的三维特征提取方法主要有数学统计法[2]、 模型投影法[3]和拓扑分析法[4]， 三种方法均不同程度提高了特征提取准确度， 但该三种方法的特征分析只停留在模型表面， 无法综合反映模型内部结构属性。

如何深入模型内部， 并提取能够反映三维模型表面属性和内部结构的综合特征信息已成为特征提取技术的研究焦点， 相关研究表明光具有较强的穿透特性， 光在不同介质中传播时， 不同光谱波段的光会与介质发生不同作用， 产生包括散射、 反射、 透射、 吸收等在内的多种物理现象[5]， 而这些现象反映了物体材料属性及物理形状。 因此， 在一定空间内模拟光与三维模型相互作用， 统计并分析能量吸收、 透射反射次数以及光传播路径等多种信息因子， 进而可获得一种代表性强、 区分度高的三维模型特征。

鉴于可见光独特的传播特性以及三维模型复杂多样性， 本文使用蒙特卡罗法[6]模拟光在三维模型中传播全过程， 将光谱分析理论与可见光传播特性分析拓展到三维模型特征提取技术中， 提出一种新的三维模型特征提取方法。

1 原理与方法

光在介质中传播时会发生多种物理现象， 图1显示了光在三维模型中产生各种光学现象的原理示意图， 不同波段可见光在不同形态的介质中产生光散射、 反射以及透射的概率具有随机性， 蒙特卡罗法能够较好解决随机性问题[7-8]。 使用蒙特卡罗法模拟光在三维模型中的传播并实现特征提取主要包括： 设置光子束的初始化状态， 根据光的波长特征和传播介质的光学属性确定光子移动步长， 分析并确定光散射、 界面反射概率， 实现光子传播全过程跟踪模拟， 统计分析光子信息， 从而完成三维模型特征提取。

图1 光传播过程示意图

1.1 模拟初始化

初始化即确定模拟的起始状态， 包括模型初始化与光子初始化两部分， 对于三维模型的初始化， 首先对三维模型进行尺度归一化、 旋转归一化处理[9]， 从而保证三维模型相对大小一致， 然后确定三维模型重心， 并以重心为原点建立三维空间坐标系， 将三维模型包围在标准约束空间中， 并使用方位角α以及偏转角β表示模型的轮廓， 图2给出了坐标系中方位角、 偏转角示意图。

图2 三维空间坐标系示意图

以光子束为最小模拟单位， 模拟光在约束空间中的传播历程， 光子初始化时光子束的入射位置(x,y,z)均匀分布在约束空间表面， 入射方向指向坐标原点， 并使用位置坐标点与(x,y,z)坐标轴的余弦值(Ux,Uy,Uz)表示入射方向。 此外， 为了实现光子束的能量信息统计， 初始化后的光子束均被设置w0能量。

1.2 计算移动步长

定义光子在传播方向未改变时的直线传播距离为移动步长， 移动步长受到包括各向异性、 折射率、 吸收系数、 反射率、 散射系数在内多种因素影响， 但主要由介质的散射系数、 吸收系数决定[10]， 式(1)表示了光子平均移动步长与散射系数、 吸收系数反比例关系， 式中Δs表示平均步长，σs和σa分别代表光散射系数、 吸收系数。

(1)

移动步长的大小具有随机性， 其概率密度函数遵守Beer定律[9]， 见式(2)

p(s)=(σa+σs)e-(σa+σs)s

(2)

(3)

结合概率密度以及式(1)反比例关系， 计算得出移动步长， 见式(3)， 式中随机变量ε∈[0, 1]。 此外， 光子移动单个步长的同时伴随着光子的能量转移， 体现光子束的能量损失， 单次移动后的光子束能量为

(4)

式中，w′表示移动后的光子束能量，w为移动前光子束能量。

1.3 确定散射方向

光子束在传播过程中因光散射而改变传播方向， 三维模型的形状越复杂， 结构变异连接处的光学性质差异越大， 光散射效果越明显， 因此在模拟过程中增加光散射过程有助于提高三维模型特征提取精度。 相函数是一种描述入射光单次散射偏转角度的概率分布函数， 其中Henyey-Greenstein相函数(H-G函数)计算方便[11]， 在蒙特卡罗法模拟中最为常用， H-G相函数表达式为

(5)

式中，β为偏转角；g∈[-1, 1]表示各向异性系数； 随机数ε∈[0, 1]。

各项向性系数g=-1时， 入射光全部背向散射； 反之g=1时， 入射光完全前向散射。 散射偏转角确定后， 在散射方位角取值范围内取随机数， 完成散射方向的改变。

1.4 模拟终止化

光子束出射后在被包围的约束空间中传播， 与模型表面、 空间约束界面碰撞， 产生折射、 反射等现象， 但光子束不可能无限传播， 为保证模拟随机性与模拟效率， 定义2种光子束终止化方式： (1)光子束能量低于临界阈值w″时终止； (2)光子束在约束空间表面和模型表面透射时终止。

光子束传播到三维模型和约束空间表面时， 可能发生透射或内反射现象， 取随机数值ε∈[0, 1]与Fresnel反射系数[12]R(ai)比较， 判断是否发生内反射， 若ε>R(ai)， 光子束透射界面， 反之光子束发生内反射。

(6)

式中， ni表示传播介质折射率； nt表示出射介质折射率。

1.5 分析并提取特征

光子束从初始化到终止化为止， 光子束全部传播信息均可被跟踪记录， 光子束特有的传播信息体现三维模型不同形态与内部结构。 单个光子束传播信息构成随机事件， 根据概率统计理论对大量随机事件统计分析， 即可完成有效信息的提取， 表1显示了三维模型特征提取过程中光子束传播信息统计量。

表1 光子束传播信息统计量

将不同统计特征组成特征矩阵， 对特征矩阵、 角度分布、 距离分布、 能量分布特征进行加权平均， 完成三维模型特征提取。

2 实验部分

选用美国普渡大学开发的ESB国际通用模型库作为三维模型实验数据集， ESB数据集[13]包括薄壁、 棱柱、 回转体3大类共计867个模型， 具有种类齐全、 覆盖范围广、 易于使用的特点， 部分三维模型如图3所示。

图3 ESB三维模型节选

实验中三维模型的光学特性参数为： 各向异性系数g=0.12， 光学参数σs=12 mm-1，σa=0.03 mm-1， 在Matlab2016中按照图4所示的流程对光子束传播过程进行模拟。

为实现三维模型特征提取效果的定性和定量评价， 实验以查准率、 查全率、 模型检索效率[14]等作为三维模型特征提取的评价标准。

3 结果与讨论

3.1 光子数量分析

光子数量是决定特征提取效果的关键影响因素， 模拟光子数量不足导致提取特征信息量少， 无法反应模型基本特征， 数量过多则特征提取时间长， 特征提取效率低。 图5显示了光子束数量对特征提取准确度、 提取效率的影响， 其中横坐标表示不同梯度光子束数量， 纵坐标以查准率表征提取效果， 以所需时间表征提取效率。

图4 实验模拟流程图

图5 模拟数量分析图

由图5可知， 特征提取时间随着光子束数量提高而增加， 同等变化幅度内呈现线性变化趋势， 表明光子模拟算法无记忆性； 查准率曲线先升高， 在光子束数量为5万左右时达到顶峰后呈现下降趋势， 这是因为随着光子束数量增加， 光子在模型内部传播密度增大， 表征模型特征点的统计量更加充分， 从而导致查准率升高， 但随着光子束数量逐步增加， 各统计量数据逐渐饱和且趋于正态分布， 模型特征的区分度减小， 查准率降低。 综合特征提取效果和提取效率两方面， 确定光子束数量范围为10 000～25 000之间， 此时， 特征提取时间需要30 s左右， 表明光子特征提取方法更适合离线三维模型检索与应用。

3.2 约束空间分析

约束空间是影响特征提取的重要因素， 一方面约束空间决定了光子束的发射位置， 另一方面约束空间内表面形状影响光子反射角度。 图6显示了不同形状约束空间下部分光子三维传播轨迹， 图6(a)为ESB模型库中三维模型， 图6(b, c, d)分别为圆球体、 正方体、 椭球体中光子传播轨迹， 其中红色代表被约束空间表面内反射的光子， 绿色表示被吸收的光子， 蓝色为透射出三维模型表面的光子。

图6 不同约束空间中光子分布示意图

由图6光子分布图可知， 在不同约束空间中， 光子传播轨迹有较大差异， 为量化约束空间对特征提取的影响效果， 选取结构差异较大的10种三维模型在不同约束空间中进行特征提取， 由表2统计结果可知， 在特征提取效果方面， 圆球体和正方体优于椭球体特征提取； 在稳定性方面， 圆球体性能最优。 此外， 使用圆球体约束空间， 三维模型旋转变化后与约束空间保持相对位置不变， 因此在初始化阶段三维模型不需要旋转归一化处理， 算法复杂度降低， 综合特征提取效果与流程两方面， 确定圆球体为最佳模拟约束空间。

表2 不同约束空间特征提取效果

3.3 算法比较

为验证光子模拟方法的特征提取效果， 使用不同特征提取方法进行对比实验， 3种比较算法分别为使用频率较高的D2分布算法、 距离-夹角算法以及频率域小波变换方法[15]。 使用球体约束空间并将光子模拟次数设置为1.2万次， 在ESB模型库中选取10种不同类别三维模型， 每种三维模型进行5次特征提取， 以5次模型检索的查准率-查全率平均值作为最终统计数据， 比较结果如图7所示。

由图7可知， 整体上本文特征提取效果与小波变换算法近似， 检索准确度高于D2分布与距离夹角算法， 这是因为基于形态的D2算法与距离-夹角算法只以模型表面信息点作为特征提取源， 小波变换和光子传播模拟方法能够深入模型内部， 体现表面与内部融合的双重特征； 为进一步比较本文算法与小波变换算法的精度， 图8显示了以内部形状复杂多腔体三维模型作为检索对象的对比效果， 此时小波算法检索准确度明显低于光子模拟算法， 这是由于在使用小波变换进行特征提取前三维模型必须要进过旋转归一化处理， 复杂三维模型在旋转归一化处理时会发生空间坐标系不重合现象， 而小波变换对模型位置信息变化较为敏感， 从而一定程度上降低了特征提取区分度。 以上两方面对比分析表明， 光子传播模拟算法的特征提取准确度较高， 能够满足三维模型特征提取与检索要求， 尤其适用于内部复杂三维模型特征提取。

图7 查全率查准率曲线

4 结 论

(1)蒙特卡罗法可实现光在不同介质中散射、 折射、 反射等现象的模拟， 将该方法应用到三维模型特征提取方面， 能够提取出融合模型表面特性与内部形态的有效特征， 是一种新的三维模型特征提取方法。

图8 复杂三维模型E测度分析结果

(2)采用蒙特卡罗光子模拟法并在ESB模型库中进行三维特征提取时， 最佳约束空间为球体， 10 000～25 000之间是光子束模拟数量的最佳取值范围。

(3)在球体约束空间且光子束数量为1.2万次的条件下， 光子模拟法的特征提取准确度高于小波变换、 距离夹角以及D2分布算法， 但特征提取效率有待提高， 该方法更适合复杂三维模型的离线特征提取与应用。