吕超凡, 言颖杰, 林 力,3, 柴 岗, 鲍劲松
(1. 东华大学 机械工程学院,上海 201620; 2. 上海交通大学医学院附属第九人民医院 整复外科,上海 200011; 3. 上海交通大学 塑性成形技术与装备研究院,上海 200030)
下颌角截骨手术用于解决患者的下颌角肥大、左右两侧不对称和面部上下宽度比例不协调等问题,是近年来较为热门的颅面整形手术.医生获取患者的下颌骨三维模型后,在考虑手术安全性的前提下,结合面部协调性和患者的美观需求,确认下颌角截骨平面,并基于此设计个性化截骨定位导板.在手术过程中,医生在截骨定位导板的辅助下,沿着截骨平面切除患者的下颌骨突出肥大区域,改善患者的下颌角形态,以达到美观效果.若截骨平面定位不理想, 术中可能会发生下颌骨升支意外骨折、下牙槽神经受损和双侧不对称等问题[1].截骨平面的设计需要医师有丰富的临床经验,并且耗时较长,平均每个病人需要花费26 min.因此,快速准确的截骨平面自动设计对于下颌角整形手术具有重要意义.
不同患者的下颌骨形态差异较大,难以通过数学建模的方式得到特定患者的截骨平面.而深度学习方法能够学习样本数据的内在规律,解决复杂的模式识别难题,因此,本文利用深度学习实现截骨平面的自动设计.截骨平面设计可以视为一个回归问题,即输入一个三维模型,输出截骨平面参数.然而该任务难以训练,因此本文先对输入的三维模型进行语义分割,即从病人的下颌角三维模型中找到突出肥大区域,然后由语义分割结果得到截骨平面.目前,基于深度学习的三维模型语义分割方法可分为基于体素、基于多视图和基于点云的方法.基于体素的方法[2-3]主要瓶颈是低分辨率体素会丢失较多信息,而高分辨率体素会导致计算和存储成本急剧增加.基于多视图的方法[4-5]主要局限是三维到二维的信息丢失问题.而基于点云的方法[6-12]能够最大程度地保留三维模型的空间位置信息和几何结构信息.因此,本文采用基于点云的语义分割算法实现下颌角截骨平面设计.
基于点云的算法是直接在点云上进行三维处理任务.文献[6]开创了点云深度学习的先河,提出了利用共享多层感知机(MLP)提取逐点特征,并通过最大池化聚合全局特征,从而解决点云的无序性问题, 但是其缺乏提取局部特征的能力.为了提取潜在的几何形状,文献[7]利用最远点采样和分组策略来提取本地特征; 文献[8]通过边卷积来聚合局部特征;文献[9]对三维点云之间的几何关系编码方式进行研究.尽管这些方法能够捕获局部上下文信息,但是在不同类别点的交界处仍不能产生细粒度的分割,原因为这些特征聚合模块通常采用最大池化或平均池化来聚合特征,使得区域内的点具有语义一致性.为了提升边界的分割性能,文献[10-11]在聚合局部特征时引入了注意力机制,但是非本地特征的提取能力较差.
本文结合下颌骨点云语义分割任务的特点,即要求截骨区与非截骨区过渡区域的精细化分割,提出一种基于Transformer的点云语义分割网络.利用基于注意力机制的本地特征提取层,网络能够提取下颌骨点云中的局部精细结构,从而对截骨平面附近点云进行细粒度分割;利用基于 Transformer的非本地特征提取层,网络能够提取点云的全局上下文信息,并将逐点特征映射到一个高维的语义分割特征空间,进一步提高了下颌骨语义分割精度.试验结果表明,与其他算法相比,本文算法能够较好地预测下颌骨点云的截骨区域,截骨区与非截骨区的过渡区域最为平整,分割结果更好.
针对下颌角整形手术,医生在获取颅骨电子计算机断层扫描(CT)数据后,利用CT三维重建技术得到患者下颌骨的三维模型;其次,手工标记下颌角中的重要解剖点、线、面和下牙槽神经位置,作为截骨手术规划的参考;再次,在考虑手术安全性的前提下,结合面部协调性和患者的美观需求,确认下颌角截骨位置;最后,根据患者的下颌角切除区域设计相应的手术导板[13].整个流程如图1(橙色部分)所示.
在下颌角整形手术术前规划中,截骨平面的设计是最为关键且耗时的环节.本文利用深度神经网络自动预测截骨平面,截骨平面的智能规划流程如图1(蓝色部分)所示.首先,本文利用表面点采样获取下颌骨和下牙槽神经的点云数据;然后,利用点云语义分割网络预测下颌骨点云中的截骨区域,并计算截骨平面;最后,以下牙槽神经点云作为参考,对截骨平面进行调整,以确保截骨平面不会切到神经,从而保证手术的安全性.
图1 下颌角截骨术前规划的总体流程Fig.1 Overall process of preoperative planning for mandibular angle osteotomy
图2 下颌骨点云语义分割网络Fig.2 Semantic segmentation network for mandibular point clouds
在获取患者下颌骨和下牙槽神经的三维模型后,首先利用表面点采样将其转换为固定点数的点云数据,然后对下颌骨点云进行归一化处理:将点云重心平移至坐标原点,并对其进行缩放,使点云中x、y和z方向上的最值在[-1,1]之间.点云分割网络的输入为归一化后的下颌骨点云,在得到预测结果后将点云缩放至原比例,用于后续截骨平面的调整.
1.2.1网络结构 下颌骨点云语义分割算法的总体框架如图2所示.该网络旨在将输入点云(点的数量为N)转换为新的高维特征向量,以获取语义丰富、互相关联的逐点特征,并基于此进行点云语义分割任务.首先,通过局部特征提取模块对各个点的局部几何结构信息进行编码;其次,将这些特征输入到Transformer层中,学习不同点特征之间的依赖关系,丰富点特征的语义信息,并通过共享MLP将点特征转换为1 024维;再次,将最大池化与平均池化分别应用于逐点特征,并将两个池化操作的输出结果连接在一起,以得到一个有效表示输入点云的全局特征;最后,对全局特征进行复制后,将其与逐点特征连接,并通过共享MLP(512,256, 3)(数值表示各层神经元个数)和Softmax激活函数得到各个点的预测类别.
1.2.2局部特征提取模块 给定点云P={Pi|i=1, 2, …,N},各个点Pi的属性为其三维坐标pi∈R1×3(文中R的上标为矢量和张量的维度),也可以包含一些额外属性(例如颜色和法向量,本文试验中只输入点的坐标),记为qi∈R1×d.
局部特征提取模块作用于点云P中的每个点如图3所示,该模块将逐点特征升维至de.对于点Pi,首先通过K最近邻(KNN)算法获取该点的邻域A(i)=Pj,j=1, 2, …,k;然后对点Pi及其邻近点进行空间关系编码.文献[9]的研究表明,将中心点与邻近点的关系编码为如下形式,其对于分类和分割等任务的效果最好.
ei,j=concat(pi,pj,pi-pj,di,j)
(1)
式中:pj为点的坐标;di,j为点Pi和点Pj的欧氏距离;concat为向量的连接.
编码邻域的空间关系后得到十维的向量ei,j,将其与点Pj的额外信息qj连接,然后通过映射函数ρ对其进行特征提取,最后聚集该邻域的特征,得到指定维度的特征向量.局部特征提取函数表示为
Fi=R{fi,1,fi,2, …,fi,k}=
R{ρ(concat(ei,j,qj)) |Pj∈A(i)}
(2)
式中:Fi∈R1×de为聚合后的局部特征;ρ由共享MLP、批标准化和ReLU激活函数组成;R{·}为特征的聚合函数,用于汇总邻近点的特征;fi,k∈R1×de为边特征.
图3 局部特征提取模块Fig.3 Module of local feature aggregation
以往研究多用最大池化[7]或平均池化[14]聚合相邻特征,但这样会丢失许多信息.因此,有研究尝试将注意力机制应用于相邻特征的聚合,自动学习邻域中各个点对于当前点的重要程度[10-11].
边缘点是不同类别点云之间的过渡区域.在下颌骨语义分割任务中,边缘点预测标签是否正确比较重要,原因为下颌骨分割面是利用边缘点得到的.然而边界上的特征提取通常不明确,原因为其混合了不同类别点的特征,不同类别的特征跨边界传播将会导致边缘点的分割结果不佳.语义分割算法的本质是利用深度神经网络将原始数据映射到一个高维的特征空间.在这个特征空间中,同类别点的特征向量应尽可能接近,不同类别点的特征向量应尽可能远离.因此,可以利用注意力池化来聚合邻域A(i) 的特征集合{fi,1,fi,2, …,fi,k}, 如图4所示.
利用注意力池化可以学习邻域中各个点对于中心点的注意力得分,从而减缓特征的跨区域传播,使过渡区域的分割结果更好.注意力池化分为计算注意力得分和加权求和两个步骤.
在点Pi的邻域A(i)中,任意邻近点特征fi,j的注意力权重αi,j是通过MLP学习到的逐通道注意力得分,表示为
αi,j=μ(concat(fi,j,pi-pj)),Pj∈A(i)
(3)
式中:μ由共享MLP、批标准化和激活函数组成;αi,j与fi,j的维度相同.
在获得各个邻近点特征的注意力权重后,通过加权求和汇总邻近点特征为
(4)
1.2.3Transformer层 在下颌骨点云语义分割任务中,为了完成细粒度的分割,需要提取各个点的本地特征和非本地特征.本地特征能够提供局部区域的几何形状信息和位置信息;非本地特征能够提供点的全局信息,建立点云中任意点与当前点的依赖关系,使同类别的点在特征空间中尽可能接近.
Transformer具有全局信息建模能力,且对输入数据或特征具有顺序无关性,因此被广泛应用于自然语言处理和时序数据处理等任务[15-17].点云作为三维空间中的一组无序点集可用Transformer处理.使用局部特征提取模块提取各个点的本地特征后,通过Transformer层提取点的非本地特征,具体如下所述.
Transformer由自注意力机制、前馈层和残差连接组成,如图5所示.通过局部特征提取模块提取的点特征记为Fe∈RN×de.Transformer层的输入为点云P中各个点的坐标和特征.
图5 Transformer层Fig.5 Transformer layer
自注意力机制是注意力机制的一种形式,用于计算输入数据内部的语义关联.自注意力机制主要由查询向量Q、键向量K和值向量V组成,计算方式如下:
(Q,K,V)=Fe(θq,θk,θv)
(5)
Q,K∈RN×da,V∈RN×de
式中:θq,θk和θv是共享MLP的可学习参数,为提高计算效率,在试验中设置da=de/4.
注意力机制不会保留输入数据的位置信息,因此需要通过位置编码使注意力层适应输入数据的空间位置关系.在自然语言处理和图像领域中,常用的位置编码方式为绝对位置编码[15]和相对位置编码[18].在三维点云的处理任务中,点的坐标本身包含空间位置信息,但是单独一个点的坐标并不能体现该点在整个点云中的位置,因此本文把点云P的重心pc作为参考点,建立各个点与整个点云的相对位置关系,即
(6)
β=δ(pi-pc), ∀pi∈P
(7)
式中:β为位置编码;δ为一个R1×3→R1×da的共享MLP.经过试验发现位置编码对于注意力得分的生成比较重要,因此本文将位置编码β分别与矩阵Q和K相加,即
Q′=Q+β
(8)
K′=K+β
(9)
Fout=ρ(Fe-Fa)+Fe
(10)
下颌神经管位于术区附近,神经管内有下牙槽神经和伴行血管.下牙槽神经损伤是严重的术后并发症之一,若截骨过程中不慎伤及神经,将导致下颌、下唇部位,甚至牙龈感觉的丧失.因此,一个合格的术前设计方案,其截骨面与神经的最近距离必须控制在安全范围内.
利用深度学习算法预测的结果可能不符合截骨手术规划要求,尤其是截骨方案的安全性问题.因此,利用所述算法预测截骨区域后,本文以下颌骨神经点云为参考,判断截骨平面是否需要进行调整.左右两边的计算方式相同,因此以左截骨平面为例.对于截骨区中的每个点,计算其到非截骨区的最小距离,若该距离小于安全距离阈值ds,则将该点视为截骨平面附近的点.在得到截骨平面附近的点集PB后,利用平面拟合得到截骨平面L.最后计算L到神经管的距离dmin,若该距离大于ds,则平移L.算法的具体描述如下.
输入:左截骨区点云Pα∈RN1×3,非截骨区点云Pε∈RN2×3,左神经点云Pe,截骨面到神经的安全距离阈值ds.
输出:截骨平面L.
1.计算Pα和Pε的逐点距离,得到距离矩阵D∈RN1×N2.
2.建立空集PB=∅,来保存截骨平面附近的点.
3.fori←1 toN1do
4.di=100
5. forj←1 toN2do
6. ifDi,j 7.di=Di,j 8. ifdi<0.1 then 9. 将点Pα,i加入到集合PB中. 10.对点集PB进行平面拟合,得到截骨平面L={ax+by+cz+m=0}. 11.计算点集Pe到L的最小距离dmin. 12.ifdmin 本文数据集为上海交通大学医学院附属第九人民医院整复外科提供的240个病例,每个病例的数据包含术前下颌骨、左右两侧的截骨区域和下牙槽神经,其数据格式为立体光刻(STL)模型.制作带标签数据集的过程为:① 对输入数据进行表面点采样;② 对截骨区点云进行平面拟合,以获取两侧的截骨平面;③ 利用截骨平面将术前下颌骨点云分为3个部分,并对各个区域的点赋予不同的标签,术前下颌骨区域标记为0,左、右截骨区域分别标记为1和2;④ 对点云进行归一化处理. 当点数较多时,训练时占用的显存和训练时间都会急剧增加;当点数较少时,会给截骨平面的拟合带来误差.因此,本文将下颌骨点云的点数设为 4 096,下牙槽神经的点数设为500,处理后的点云如图6所示.下牙槽神经点云是训练时的安全性参考指标,不作为点云语义分割网络的输入.按照比例4∶1∶1将原始数据分为训练集、验证集以及测试集. 图6 带标签的下颌骨点云Fig.6 Mandible point cloud with labels 本文试验在Ubuntu16.04操作系统、Tesla V100显卡(32 GB显存)的服务器上进行.采用 Pytorch 实现所述框架,使用交叉熵损失函数.训练批次设为8,迭代次数为200次,采用Adam优化器,初始学习率为0.001,学习率衰减方式为余弦退火,权重衰减为0.000 1. 试验中对分割网络的评价指标为准确率、精确率、召回率和F1值.本文将截骨区点云作为正类,非截骨区点云作为负类,在下颌骨点云数据中,正负类样本比例约为1∶10,因此相比于精确率,F1值对于预测结果的评价更为准确.选取PointNet++、DGCNN和PCT等点云语义分割网络与本文算法进行对比.在进行算法的对比试验时,试验环境保持一致. 下颌角截骨方案的主要评价指标为截骨对称性、安全性和美观性[19].测量过程中涉及的几个解剖标志点如图7所示,各评价指标的定义和测量方式如下. 截骨对称性:下颌角点是评价面部轮廓外观的关键点,在评价截骨对称性时,可以在下颌骨双侧测量点C与点G的距离,两侧dC G的差值为ΔdC G,若ΔdC G<3 mm, 则左右两边对称,否则不对称. 截骨安全性:在截骨手术时要保证截骨区域离下颌角神经有一定距离,避免手术过程中对神经造成损伤.在评价截骨安全性时,可以测量截骨平面到神经的最小距离ds,若ds>4 mm,则视为安全,否则不安全. 截骨美观性:下颌角美观性的评价依据是下颌角角度θ(髁顶点、下颌角点和颏下点的夹角),θ在115° 左右为正常,下颌角肥大患者的θ在110° 左右,θ在120° 左右为美观.因此,当θ∈(115°, 125°)时,视为美观. 图7 下颌骨的解剖标志点Fig.7 Anatomical landmarks of mandible 不同算法的各项评估指标在表1中列出.试验结果表明,所提方法的各项评估指标均达到最佳性能,在下颌骨语义分割试验中,本文所述方法优于先前的网络框架. 表1 不同方法在下颌骨语义分割测试集上的评估结果 本文根据测试集中40个病例的截骨方案预测结果,在医学影像控制系统中对下颌角模型进行虚拟截骨后, 对相关数据进行测量,结果如表2所示. 表2 预测截骨方案的评价指标Tab.2 Evaluation index of predicted osteotomy plan 在40个病例的预测结果中,只有2个不符合要求,试验结果表明本文算法的截骨手术规划效果良好. 在测试集中选取1个病例,不同方法对该病例的预测结果如图8所示(未做截骨平面微调),图中红点为预测错误的点.可知,在DGCNN和PCT算法的预测结果中,预测标签错误的点大多集中在截骨区与下颌骨主体的过渡区域,原因是这些点在局部特征聚合时混合了不同类别点的特征,而本文算法通过注意力池化减缓了特征的跨区域传播.此外,加入相对位置编码的Transformer层也能更好地建立点与点之间的依赖关系,使得同类点在特征空间中尽可能接近,提高其语义一致性,从而提高下颌骨点云的分割精度.试验结果表明,与其他两种模型相比,本文算法的语义分割精度最高,非截骨区与下颌骨主体的过渡区域也最平整. 图8 不同点云语义分割算法的预测结果可视化Fig.8 Visualization of prediction results of different point cloud semantic segmentation algorithms (1) 邻近点的数量.邻近点的数量k用于确定提取局部特征时局部区域中点的数量,结果如表3所示.当k=16时,模型性能最佳.当邻域较小(k=4)时,模型可能没有足够的上下文信息以准确预测各个点的标签;当邻域较大(k=64)时,每个点特征中包含较多邻近点的信息,可能会引入过多的噪声,导致模型的准确性下降.因此,k选取16最为合适. 表3 不同邻近点数量的比较 (2) 特征聚合函数.特征聚合函数用于聚合局部区域中邻近点的特征.本文将注意力池化与平均池化和最大池化进行比较,试验结果如表4所示.试验结果表明,注意力池化能够提高点云语义分割的性能. (3) 位置编码.进行3组对比试验,即没有位置编码、使用绝对位置编码和使用相对位置编码.相对位置编码如前文所述,绝对位置编码是将式(7)中的δ(pi-pc)改为δ(pi),结果如表5所示.试验结果表明,如果没有位置编码,性能将大大下降;使用绝对位置编码,精度比不使用绝对位置编码高;相对位置编码产生最佳结果. 表4 不同特征聚合函数的比较Tab.4 Comparison of different feature aggregation functions % 表5 不同位置编码方式的比较Tab.5 Comparison of different position encoding methods % 提出一个结合本地特征提取模块和非本地特征提取模块的下颌骨点云语义分割网络,并基于下颌角整形手术的诊断病例,构造了下颌骨语义分割数据集,训练出一个能有效预测截骨区域的语义分割模型.引入注意力机制的局部特征提取模块能够获取精细的局部特征,减缓特征的跨区域传播,使得截骨区与非截骨区的过渡区域更为平滑.加入相对位置编码的Transformer层能够获取下颌骨点云的全局上下文信息,自动学习点特征间的依赖关系,进一步提高分割精度.试验结果表明,本文算法预测的截骨区域十分接近医师手动标记的截骨区域,可以作为医师进行截骨手术规划的参考.2 试验验证
2.1 数据集构建
2.2 试验环境设置
2.3 网络性能评估指标
2.4 临床验证评估指标
3 试验结果与分析
3.1 网络性能分析
3.2 临床验证结果分析
3.3 预测结果可视化
3.4 消融研究
4 结语