1 引　言

在我国每年数以万计的原发性脑肿瘤患者中, 神经胶质瘤的发病数占到了将近一半, 其中脑胶质瘤是最常见的颅内原发性肿瘤, 约占所有脑部肿瘤的30%–40%^[1-5]. 脑胶质瘤病人随着症状的发展有可能会出现昏迷或者死亡的情况. 当前脑胶质瘤的治疗方案中, 目标肿瘤区域的位置判断准确性对外科手术的术前规划至关重要, 通过手术术前规划确定脑胶质瘤位置、大小和形状, 选择最佳手术点进行开颅手术的程度以及手术器械的切入途径^[2], 为规划优化和个性化治疗提供重要信息^[3]. 传统没有三维图像指导的情况下, 大多数外科医生必须同时滚动浏览二维(2D)正交切片序列, 或者通过术前计算机断层扫描(computed tomography, CT)或磁共振成像(magnetic resonance imaging, MRI)在计算机上进行交互重建. 近年来, 混合现实(mixed reality, MR)技术可以减少外科手术从2D图像到3D图像的可视化转换障碍^[4], 引起了医学影像分析相关领域的广泛关注. MR技术通过将虚拟模型与真实场景进行叠加, 有效实现颅内结构的3D模型与患者头颅的虚实融合, 使用户与真实世界及虚拟模型实时交互, 为外科手术的术前计划和术中导航提供了直观、有效的技术支撑. 当前, 如何有效检测病灶的具体位置与大小, 实现病灶的有效分割与实时的三维重建, 是面向医学影像分析MR系统中的首要任务.

由于脑胶质瘤医学影像中存在灰度偏移场, 使得MRI肿瘤图像组织边缘模糊, 特别是脑胶质瘤位置随机性伴有肿瘤块效应, 导致很难精确观测肿瘤病变区域. 基于传统机器学习的检测与分割方法需要对图像特征进行人工分析和选择, 费时费力, 且算法泛化能力弱、鲁棒性差、分割效率低^[5]. 近年来, 深度学习方法广泛用于医学图像分析中. 本文基于深度学习方法, 进行脑胶质瘤的检测与分割, 有效提高病灶区域检测效率与分割准确度, 为脑胶质瘤医学影像分析混合现实系统提供位置、大小和形状精确的病灶掩码. 其次, 将MRI图像经过深度学习算法检测和分割得到的掩码数据, 在病灶模型与脑颅组织的三维重建算法中直接影响MR系统的性能与效果. 目前医学图像重建的方法主要有3D面绘制(surface rendering)^[6]与体绘制(volume rendering)^[7]. 体绘制可以绘制出包含丰富信息的三维实体, 但是计算量大、算法效率低. 面绘制可以有效地绘制三维实体的表面, 将感兴趣的部分提取出来, 因此面绘制成为医学图像三维重建的主流方法^[8]. 本文研究基于面绘制的脑胶质瘤医学影像多结构组织的三维绘制与优化方法, 通过运用不同颜色对分割结果的各种组织进行绘制, 使用不同蒙版进行可视化分割、提取与优化, 生成MR系统在同一坐标系下的3D虚拟模型. 此外, 为了实现物理世界与虚拟世界的有机融合, MR系统需要将三维重建好的虚拟模型准确定位到真实环境中的物理对象上, 并能根据真实场景中目标位置的变化, 实时地按照使用者的当前视角, 重新建立空间坐标系, 将虚拟场景准确地渲染到真实环境中. 这涉及到MR三维注册与跟踪关键技术, 目前主要包括基于标志物的和无标识的两类三维跟踪注册方法. 本文研究并实现了基于交互式无标识的三维注册方法, 同时为了增强脑胶质瘤影像分析MR系统的实用性, 进一步研究了MR多设备的第3视角空间实时共享方法, 实现了微软HoloLens设备与安卓手机的视角空间实时共享, 使用户可以直接通过安卓手机体验脑胶质瘤医学影像分析MR技术.

综上所述, 本文旨在将当前深度学习技术、混合现实技术与医学影像智能分析技术有机结合, 研究基于深度学习的脑胶质瘤精准检测与分割算法; 在此基础上, 研究脑胶质瘤医学影像面绘制算法与优化方法, 实现脑胶质瘤不同区域的三维重建与可视化展示, 研究并实现了基于交互式无标识的三维注册跟踪方法, 实现了基于标识图的MR移动端多设备的第3视角空间实时共享技术, 并通过实验与MR系统, 测试、验证了相关算法及系统的综合性能, 为医学影像分析、手术术前规划、医学教育与培训等提供了新方法与技术支撑.

2 相关工作

自1957年引入沉浸式环境以来, 人们开始对虚拟现实(virtual reality, VR)和增强现实(augmented reality, AR)环境在医疗领域的应用进行挖掘, 以帮助外科医生进行外科手术的训练、规划和导航^[9,10]. 在神经外科领域, 已经涌现出许多VR模拟器和工作台, 但用户必须戴上特定的VR设备, 使用控制器进行手术模拟, 这种沉浸式的体验只能运用于模拟, 并未与真实环境相叠加. Shamir等人^[11,12]描述的AR计划环境被证明在微创神经外科手术中也可促进术中导航, 随后AR系统的开发愈发普及. Lerotic等人^[13]利用一种新颖的非真实感技术可视化虚拟模型, 即在改变模型透明程度的基础上同时保持暴露表面的组织器官, 验证了虚实融合显示技术在外科手术导航的可行性. 此后, Elmi-Terander等人^[14]开发了一种增强现实手术导航(ARSN)系统, 该系统集成了术中成像系统和光学跟踪技术, 应用AR技术进行胸椎椎弓根螺钉置入.

但是, 交互性不强的增强现实技术还不能满足外科医生对术中导航位置进行实时调整的需求. 与传统的虚拟环境相比, 新兴的混合现实技术提供了更直观的可视化途径和交互方法. Wei等人^[15]运用混合现实指导进行了射频针植入手术的试验研究, 该系统提供了全息导航平台, 可指导外科医生执行肝肿瘤射频消融(RFA), 能够准确规划并精确靶向肿瘤, 展示了MR技术辅助手术的潜力.

Perkins等人^[16]提出了一个基于MR技术的乳房手术术前规划系统. 该系统通过Microsoft HoloLens将乳房MRI的3D虚拟模型投影到患者身上, 该技术旨在通过帮助外科医生确定肿瘤的范围来减少重复手术的次数. 这项初步研究在真实患者的肿瘤图像上进行实验, 其结果表明MR技术在指导外科手术时, 确实可以提供有关肿瘤位置信息, 起到令人鼓舞的效果. 但是, 由于计算时间长达30 min, 因此该系统的实时性较差, 在呈现重建的体积方面存在问题, 同时该系统互动性较差, 在手术过程中传达的信息量也比较有限.

Pratt等人^[17] 研发了混合现实系统, 重建腿骨、血管、皮肤和软组织的3D虚拟模型, 应用HoloLens对医学图像3D模型进行投影叠加, 允许外科医生在手术过程中, 可与虚拟3D对象进行实时交互, 证实了MR技术具有简化手术的潜力.

Bucioli等人^[18]提出了基于冠状动脉断层扫描(CT)图像重建的3D全息心脏的实时可视化模型. 虚拟心脏模型可以与许多HoloLens设备实时共享, 提供了基于MR环境的协作医疗诊断解决方案. 但是由于多台HoloLens的运用不仅使手术成本提高, 而且手术过程中的长期佩戴给外科医生带来了不适感.

本文研究了集脑胶质瘤智能检测与分割、三维重建及MR技术为一体的手术导航系统, 能自动检测和分割输入的MRI图像, 将三维重建的脑胶质瘤病灶区域叠加到病患区域, 并使用混合现实设备与Android端设备共享体验, 能够提高病灶定位精准度、增加互动体验并节约计算时间, 降低手术成本.

3 基于深度学习的脑胶质瘤分割算法

医学影像与自然图像不同, 医学影像作为多模态、深层次、边界模糊、梯度复杂的灰度图像, 包含较多的高分辨率信息, 同时也拥有人体内部结构固定、分布规律、语义简单明确的低分辩信息. 目前以全卷积经网络(fully convolutional networks, FCN)、UNet及其改进模型为主的深度学习框架虽然结合了物体类型识别与精准分割定位的高、低分辨信息, 能较好适用于医学影像分割, 但医学影像作为3D体数据图像, 使得有二维卷积和池化的结构框架更容易丢失深度信息. 所以为了更好地进行多尺度、多模态的预测与提取, 让分割边界更加精细以提高准确性, 本文运用了基于3D UNet^[19]的网络结构检测分割脑胶质瘤. 其过程主要包括MRI图像的预处理模块和检测分割模块, 如图1所示.

预处理模块对4个原始模态输入的MRI图像进行标准化、裁剪处理以及合并4个模态切片图像数据, 并保存在Numpy数组中. 由于MRI各个模态图像对比度的差异, 采用Z-score方式对图像零均值和单位标准差处理, 提高脑胶质瘤形状、大小等特征信息在各模态图像上的可比性, 加快参数优化速度收敛. 对4个模态图像(155×240×240)脑部周围的背景信息进行裁剪变为155×160×160, 避免肿瘤区域占比较小出现严重的数据不平衡, 提升网络模型分割性能.

检测与分割模块直接对预处理后的脑胶质瘤MRI图像的检测和分割, 并将检测与分割结果输出为3个嵌套子区域掩码(Whole Tumor、Core Tumor、Enhancing Tumor). 由于脑胶质瘤MRI图像是多模态的三维切片, 肿瘤在体积、形状、位置及精细度、内部结构等其他特征方面有较大变化, 加上大脑血管、脑浆等结构使得其具有丰富的形态特征. 使用二维的卷积池化或上采样将丢失部分空间信息^[20]. 所以使用三维结构将数据的空间信息利用起来, 提炼图层切片之间的隐含特征.

三维结构中的3D卷积层(conv)是2D卷积的第三维扩展, 输入F和卷积核W之间的操作定义为:

$ O_q^l(x, y, {\textit{z}}) = \partial \left(b_q^l + \sum\limits_k {\sum\limits_{r, s, t} {F_k^{l - 1}(x + r, y + s, {\textit{z}} + t)\times W_{kq}^l(r, s, t)} } \right) $

(1)

其中, $ \partial $ 激活函数使用的是经过校正的线性单位激活函数ReLU, 式中 $ F_k^{l - 1} $ 和 $ O_q^l $ 分别指的是第k层的第 $ l - 1 $ 个特征和第q层的第l个特征, $ W_{kq}^l $ 是指第l层的卷积核, $ b_q^l $ 为偏置项. 为了加快卷积的速度在每个卷积之后衔接批归一化(batch normalization, BN)操作, 可以使得学习快速进行, 不依赖初值, 抑制过拟合, 同时减少内部协变量偏移, 加快训练过程并提高准确性^[21]. 将其预测过程形式描述见算法1.

算法1. 归一化处理, 以mini-batch为学习单位, 按mini-batch进行正则化.

输入: mini-batch的 $ m $ 个输入数据集合: $B = \{ {x_1}, {x_2}, \ldots , {x_m}\}$ ; 待学习参数为: $ \gamma , \beta $ ;

输出: $ \left\{ {{y_i} = {{BN}_{\gamma , \beta }}({x_i})} \right\} $ .

begin

1. $\; {\mu _B} \leftarrow \dfrac{1}{m}\displaystyle\sum\limits_{i = 1}^m {{x_i}} $ /*对mini-batch的输入数据进行均值化*/

2. $ \sigma _B^2 \leftarrow \dfrac{1}{m}\displaystyle\sum\limits_{i = 1}^m {({x_i} - } {\mu _B}{)^2} $ /*计算mini-batch输入数据方差*/

3. $ \hat x{}_i \leftarrow \dfrac{{{x_i} - {\mu _B}}}{{\sqrt {\sigma _B^2 + \varepsilon } }} $ /*得到归一化结果, $ \varepsilon $ 为防止除零引入的极小量*/

4. ${y_i} \leftarrow \gamma \hat x{}_i + \beta \equiv {{B}}{{{N}}_{\gamma , \beta }}({x_i})$ /*对正则化后的数据进行缩放和平移变换*/

end

3D卷积层之后连接3D最大池化层, 利用不同通道之间的平移不变性对特征进行下采样, 取当前窗口的最大值作为新图的像素值, 然后深度上滑动, 得到多个特征图. 3D UNet的解码过程对应4次上采样, 将编码得到的高级语义特征图恢复到原图片的分辨率^[22]. 最后使用Softmax函数对预测结果进行分类, 得到病灶肿瘤的分割结果.

正由于该模型使用了该类三维卷积进行特征提取, 并通过三维最大池化层进行特征图下采样, 因此在编码过程中自动产生多尺寸特征图, 解决了传统全卷积分割网络中存在的输入尺寸受限、全局信息匮乏等问题^[23]. 此外, 该模型在解码过程使用了跳转连接机制(skip connection)对网络前后的多尺度特征图进行融合, 并通过上采样机制不断还原特征图, 因此全面地获取了图像的浅层特征和深层特征, 相比于传统全卷积分割网络具有更好的分割效果^[24].

4 混合现实三维重建算法 4.1 Marching cubes (MC)算法原理

MRI医学图像经过图像预处理和分割后, 图像序列将转换为一个规则的三维体数据场, 三维重建的任务就是要将分割提取的区域构建出对应组织或器官的三维模型, 实现三维可视化显示、操作及分析, 并通过几何单位拼接拟合物体表面来描述物体三维结构^[25].

Lorensen等人^[26]和Lai等人^[27]在1987年提出通过构造等值面来实现三维重建的经典方法, 以数据场中的相邻两层MRI切片上的各4个像素为顶点组成的立方体为最小等值面搜索单元(体素). 其算法的基本原理^[28]是: 首先提取一个体素, 将体素8个顶点的灰度值与阈值比较, 因此每个顶点会出现大于或小于两种情况, 整个体素就有28种构型组合, 排除立方体体素具有的对称、镜像等性质, 可简化为15种基本拓扑构型, 如图2所示.

传统面绘制方法中采用线性插值计算出与等值面相交的体元棱边上的交点计算坐标, 假设体素的棱边分别平行于x轴、y轴、z轴, 则三角面片的顶点坐标 $ V(x, y, {\textit{z}}) $ 见公式(2):

$ \left\{ {\begin{array}{*{20}{c}} {x = {x_1} + \dfrac{{(T - {I_1})\times({x_2} - {x_1})}}{{{I_2} - {I_1}}}} \\ {y = {y_1} + \dfrac{{(T - {I_1})\times({y_2} - {y_1})}}{{{I_2} - {I_1}}}} \\ {{\textit{z}} = {{\textit{z}}_1} + \dfrac{{(T - {I_1})\times({{\textit{z}}_2} - {{\textit{z}}_1})}}{{{I_2} - {I_1}}}} \end{array}} \right. $

(2)

其中, $ {V_1}({x_1}, {y_1}, {{\textit{z}}_1}) $ 和 $ {V_2}({x_2}, {y_2}, {{\textit{z}}_2}) $ 为体素棱边两个顶点坐标, $ {I_1} $ 和 $ {I_2} $ 为顶点在 $ {T_1} $ 和 $ {T_2} $ 处的灰度值, T是根据研究对象不同设置的阈值, 根据体元8个顶点与等值面的相对位置, 将这些交点按一定方式连接成等值面, 采用三角片组合曲面在体素内的逼近表示, 所有体素中的等值面构成整个数据场的等值面. 为解决三角面片之间的明暗度不连续变化的问题, 采用了中心差分法计算出梯度, 再通过采用线性插值法求解三角面片顶点法向量^[29], 体素棱边顶点处的梯度公式为:

$ \left\{ {\begin{array}{*{20}{c}} {{G_x} = \dfrac{{{f_{(x + r, y, {\textit{z}})}} - {f_{(x - r, y, {\textit{z}})}}}}{{2r}}} \\ {{G_y} = \dfrac{{{f_{(x, y + s, {\textit{z}})}} - {f_{(x, y - s, {\textit{z}})}}}}{{2s}}} \\ {{G_{\textit{z}}} = \dfrac{{{f_{(x, y, {\textit{z}} + t)}} - {f_{(x, y, {\textit{z}} - t)}}}}{{2t}}} \end{array}} \right. $

(3)

其中, r、s、t为体素棱边边长, $ {f_{(x, y, {\textit{z}})}} $ 为点 $ {V_{(x, y, {\textit{z}})}} $ 的灰度值, $ {G_x} $ 、 $ {G_y} $ 、 $ {G_{\textit{z}}} $ 为3个方向的梯度. 得到顶点梯度后, 通过索引值和索引表获取相交边, 再如公式(4)进行计算三角面顶点法向量 $ {N_{({N_x}, {N_y}, {N_{\textit{z}}})}} $ :

$ \left\{ {\begin{array}{*{20}{c}} {N{}_x = N{}_{x1} + \dfrac{{(T - {I_1})\times({N_{x2}} - {N_{x1}})}}{{{I_2} - {I_1}}}} \\ {N{}_y = N{}_{y1} + \dfrac{{(T - {I_1})\times({N_{y2}} - {N_{y1}})}}{{{I_2} - {I_1}}}} \\ {N{}_{\textit{z}} = N{}_{{\textit{z}}1} + \dfrac{{(T - {I_1})\times({N_{{\textit{z}}2}} - {N_{{\textit{z}}1}})}}{{{I_2} - {I_1}}}} \end{array}} \right. $

(4)

其中, $ {N_{({N_{x1}}, {N_{y1}}, {N_{{\textit{z}}1}})}} $ 和 $ {N_{({N_{x2}}, {N_{y2}}, {N_{{\textit{z}}2}})}} $ 为体素顶点的法向量. 在计算完体素后, 查询三角面片表, 确定组合方式得到等值面图像.

4.2 三角网格的简化

使用MC算法构造的组织表面模型一般都有数万数量的三角面形成的模型, 三角面的数量大小直接影响模型的存储容量和后续MR交互的实时性. 因此, 需要提高实时的交互能力就必须对模型进行网格简化^[29]. 在Hamann^[30]提出的基于三角形移去的模型简化方法中采用了几种简化机制^[31]: 抽样、自适应、细分、删除、顶点合并. 本文对模型使用删除方法, 重复地从网格中删除混淆顶点或者三角共面面片, 对体素之间产生的洞进行三角化. 这种算法运行简单、速度快, 多用于冗余的几何结构. 同时结合顶点合并^[32]的算法将三角面片中2个或3个顶点合并, 按材质合并网格或将网格合并为单个对象, 将网格分解为单独的元素, 以达到减少整个多边形面片数, 完成三角网格的简化.

归纳上述算法与过程, 将其用伪代码的形式描述该算法, 如算法2所示.

算法2. Marching cubes算法三维重建伪代码.

输入: MRI影像序列产生的三维离散体数据场顶点 $ V(x, y, {\textit{z}}) $ ;

输出: 三角面片顶点坐标法向量 $ ({N_x}, {N_y}, {N_{\textit{z}}}) $ 送入面绘制环节.

begin

1. for $k \leftarrow 1 \ldots {V_x}$ do /*分层读入三维离散数据场的数据*/

2. 　 for $j \leftarrow 1 \ldots {V_y}$ do

3. 　　 for $i \leftarrow 1 \ldots {V_{\textit{z}}}$ do

4.　　　 Read point value $(i, j, k), \ldots , (i + 1, j + 1, k + 1)$ /*读入体素顶点值*/

5.　　　 for ${{ii}} \leftarrow 1 \ldots 8$ do

6. 　　　　 if s[ii]>=value /*判断角点 ${V_0}, \ldots , {V_7}$ 标量值与等值面相对大小*/

7.　　　　　Index $ \leftarrow $ CASE_MASK[ii]

8. 　　　　 end if

9. 　　　 end for

10. 　　　Grayscale I, threshold T /*输入灰度值I与阈值T*/

11.　　　 $v \leftarrow {v_1} + \dfrac{{(T - {I_1})\times ({v_2} - {v_1})}}{{{I_2} - {I_1}}}$ , $ v \leftarrow x, y, {\textit{z}} $ /*计算体素棱边顶点坐标*/

12.　　　 $ G \leftarrow \dfrac{{{f_1} - {f_2}}}{{2r}} $ , $ r \leftarrow r, s, t $ /*中心差分法计算棱边顶点梯度向量*/

13. 　　　 $N \leftarrow N{}_1 + \dfrac{{(T - {I_1})\times ({N_2} - {N_1})}}{{{I_2} - {I_1}}}$ /*线性插值计算三角面片顶点法向量*/

14. 　　　Surface rendering $ \leftarrow G, N $

15. 　　 end for

16. 　 end for

17. end for

18. Optimize surface vertices

end

5 混合现实三维注册算法 5.1 基于交互式无标识的三维注册跟踪

混合现实利用三维注册技术实时检测出摄像头相对于真实场景的位姿状态, 确定所需要叠加的虚拟信息在投影平面中的位置, 并将这些虚拟信息实时显示在屏幕中的正确位置, 使得计算机生成的虚拟模型或信息完美融入到真实世界中, 以此来增强现实世界的视觉效果^[33-35]. 基于交互式无标识的注册方法是在HoloLens的视野范围内, 摄像机识别特定手势后与虚拟模型进行交互, 通过手势移动虚拟模型使其在视觉上实现与真实场景的融合^[36], 其原理如图3所示.

在图3显示的注册原理图^[37]中, ${C_{\rm world}}({x_w}, {y_w}, {{\textit{z}}_w})$ 为真实世界坐标系, ${C_{\rm camera}}({x_c}, {y_c}, {{\textit{z}}_c})$ 为HoloLens摄像机坐标, ${C_{\rm virtual}}({x_v}, {y_v}, {{\textit{z}}_v})$ 为虚拟模型坐标系, ${C_{\rm screen}}({x_s}, {y_s})$ 为显示屏坐标系.

如图4是基于HoloLens的三维注册坐标转换过程^[38], $T_{\rm virtual}^{\rm camera}$ 是虚拟模型坐标系到HoloLens摄像头坐标系的转换矩阵, $T_{\rm world}^{\rm camera}$ 是世界坐标系到HoloLens摄像头坐标系的转换矩阵, $T_{\rm camera}^{\rm screen}$ 是HoloLens摄像机坐标系到屏幕坐标系的转换矩阵. $T_{\rm world}^{\rm camera}$ 是三维空间中的刚性变换, 已知世界坐标系的一点为 $ {P_w}({x_w}, {y_w}, {{\textit{z}}_w}) $ , $ {P_c}({x_c}, {y_c}, {{\textit{z}}_c}) $ 为 $ {P_w} $ 在HoloLens摄像机坐标系下的三维坐标点, $ {P_c} $ 和 $ {P_w} $ 的转换关系如式(5)所示:

$ \left[ {\begin{array}{*{20}{c}} {{x_c}} \\ {{y_c}} \\ {{{\textit{z}}_c}} \\ 1 \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} R&{\overrightarrow T } \\ 0&1 \end{array}} \right]_{4 \times 4}}\left[ {\begin{array}{*{20}{c}} {{x_w}} \\ {{y_w}} \\ {{{\textit{z}}_w}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{R_{1x}}}&{{R_{1y}}}&{{R_{1{\textit{z}}}}}&{{T_x}} \\ {{R_{2x}}}&{{R_{2y}}}&{{R_{2{\textit{z}}}}}&{{T_y}} \\ {{R_{3x}}}&{{R_{3y}}}&{{R_{3{\textit{z}}}}}&{{T_{\textit{z}}}} \\ 0&0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_w}} \\ {{y_w}} \\ {{{\textit{z}}_w}} \\ 1 \end{array}} \right] $

(5)

其中, R为绕3个坐标轴的旋转矩阵, 包括变量 $ {R_x} $ 、 $ {R_y} $ 、 $ {R_{\textit{z}}} $ . T为三维平移向量, 二者是决定在空间中6个自由度的主要参数^[39]. 通过利用HoloLens摄像头通过对空间扫描, 确定R、T矩阵完成 $ T_{\rm world}^{\rm camera} $ 三维空间中的刚性变换.

$T_{\rm virtual}^{\rm camera}$ 的转换需要在注册中实时获得, 因此佩戴HoloLens进行注册时, 若突然加速或头部频繁晃动, 可能导致虚拟模型的位置发生变化, 增大注册误差^[39,40].

$T_{\rm camera}^{\rm screen}$ 是三维到二维平面的转换, 其过程为:

$ \left[ {\begin{array}{*{20}{c}} {{x_s}} \\ {{y_s}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} f&0&0 \\ 0&f&0 \\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\text{1}}&0&0&0 \end{array}} \\ {\begin{array}{*{20}{c}} 0&{\text{1}}&0&0 \end{array}} \\ {\begin{array}{*{20}{c}} 0&0&{\text{1}}&0 \end{array}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_c}} \\ {{y_c}} \\ {{{\textit{z}}_c}} \\ 1 \end{array}} \right] $

(6)

其中, f为焦距, 属于HoloLens摄像机内参, 整个坐标转换过程可以被表示为^[41]:

$ S = K[R|T]X $

(7)

其中, K为 $ f \cdot E $ , E为单位矩阵, R为绕3个坐标轴的旋转矩阵, T为三维平移向量, X为世界坐标点坐标. $ {P_s}({x_s}, {y_s}) $ 为点 $ {P_w}({x_w}, {y_w}, {{\textit{z}}_w}) $ 在HoloLens屏幕中得到的二维成像点. 点 $ {P_c}({x_c}, {y_c}, {{\textit{z}}_c}) $ 经过 $ T_{\rm world}^{\rm camera} $ 矩阵转化后为已知参数. $ T_{\rm camera}^{\rm screen} $ 的转换在三维注册过程中完成, 其作用是将摄像机中的虚拟模型呈现在二维屏幕上. 虚线部分表示世界坐标系和虚拟模型坐标系统一受 $ T_{\rm world}^{\rm camera} $ 和 $ T_{\rm virtual}^{\rm camera} $ 变换的影响而改变^[39-42].

将三维注册算法推导过程进行归纳, 如算法3所示.

算法3. 虚拟模型三维注册算法.

输入: 世界坐标系 $ {C_{\rm world}} $ 与虚拟模型坐标系 $ {C_{\rm virtual}} $ ; 设备参数: $ R, \vec t, f $ ;

输出: $ {C_{\rm screen}} \leftarrow K[R|t]X $ .

begin

1. $ {C_{\rm world}} \leftarrow {x_w}, {y_w}, {{\textit{z}}_w} $ /*通过Spatial mapping扫描得到世界坐标系数据*/

2. $ {C_{\rm virtual}} \leftarrow {x_v}, {y_v}, {{\textit{z}}_v} $

3. $ {C_{\rm camera}} \leftarrow T_{\rm world}^{\rm camera}{C_{\rm world}}, T_{\rm world}^{\rm camera} \leftarrow [R|\vec t] $ /*真实环境经 $ T_{\rm world}^{\rm camera} $ 转换为摄像头坐标系*/

3. $ {C_{\rm camera}} \leftarrow T_{\rm virtual}^{\rm camera}{C_{\rm virtual}} $ /*虚拟模型经实时矩阵 $ T_{\rm virtual}^{\rm camera} $ 投影在摄像头*/

4. $ ({x}_{s}, {y}_{s})\leftarrow K({x}_{c}, {y}_{c}, {{\textit{z}}}_{c}), K\leftarrow f\cdot E $ /*三维空间叠加投影在二维屏幕坐标系*/

5. $ S \leftarrow K[R|t]X $

end

5.2 移动MR端基于标识图的三维注册与位置转换

移动MR端的三维注册完成空间坐标转换, 需要Socket服务器处理两客户端并发送数据, 同时也使用基于识别图的注册方式, 因此需要HoloLens端导入Vuforia SDK来获取现实世界识别图的坐标信息, 并在一定时间内更新空间信息. 由于识别图显示在Android端屏幕, 所以当HoloLens端扫描识别图时, 虚拟模型在Android端的初始位置位于屏幕中央.

如图5所示为移动MR三维注册的空间坐标转换过程, 与混合现实设备HoloLens在屏幕上的二维坐标投影相似, ${C_{\rm Android}}({x_A}, {y_A})$ 为Android手机屏幕中央坐标系, ${C_{\rm PhoneCam}}({x_p}, {y_p}, {{\textit{z}}_p})$ 为Android端摄像头坐标系, Android手机上的摄像头位置的不同会让虚拟模型在手机屏幕上空间共享的精准度发生变化, $ {P_p}({x_p}, {y_p}, {{\textit{z}}_p}) $ 为点 $ {P_w}({x_w}, {y_w}, {{\textit{z}}_w}) $ 在Android手机摄像头中得到的三维坐标点, 转换关系见式:

$ \left[ {\begin{array}{*{20}{c}} {{x_p}} \\ {{y_p}} \\ {{{\textit{z}}_p}} \\ 1 \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {{R'}}&{{T'}} \\ 0&1 \end{array}} \right]_{4 \times 4}}\left[ {\begin{array}{*{20}{c}} {{x_w}} \\ {{y_w}} \\ {{{\textit{z}}_w}} \\ 1 \end{array}} \right] $

(8)

基于ARCore的Android端会跟踪虚拟模型空间的锚点位置, 通过Android摄像头扫描定位后确定 $ {R'} $ 、 $ {T'} $ 矩阵参数, 完成真实世界到Android摄像头坐标系的转换. 同时 $ T_{\rm virtual}^{\rm PhoneCam} $ 矩阵转换也需要在注册中实时获得.

而 $T_{\rm PhoneCam}^{\rm Android}$ 是三维到二维平面的转换:

$ \left[ {\begin{array}{*{20}{c}} {x_A^{'}} \\ {y_A^{'}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{f'}}&0&0&0 \end{array}} \\ {\begin{array}{*{20}{c}} 0&{{f'}}&0&0 \end{array}} \\ {\begin{array}{*{20}{c}} 0&0&1&0 \end{array}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_p}} \\ {{y_p}} \\ {{{\textit{z}}_p}} \\ 1 \end{array}} \right] $

(9)

$ \left[ {\begin{array}{*{20}{c}} {{x_A}} \\ {{y_A}} \\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{ - \sin \theta }&0 \\ {\sin \theta }&{\cos \theta }&0 \\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {x_A^{'}} \\ {y_A^{'}} \\ 1 \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {pos.{x_p}} \\ {pos.{y_p}} \\ 1 \end{array}} \right] $

(10)

其中, $ {f'} $ 为焦距, 是Android端摄像机内参. $ {P_A}({x_A}, {y_A}) $ 为点 $ {P_p}({x_p}, {y_p}, {{\textit{z}}_p}) $ 在Android手机端得到的二维成像点, 在经过二维逆时针旋转( $ \theta = {90^ \circ } $ )后, 使虚拟模型与手机横屏状态下一致, 得到成像点后减去手机摄像头的位置坐标( $ pos.{x_p} $ , $ pos.{y_p} $ ), 最终使得虚拟模型在Android端位置与在HoloLens中看到的位置一致, 即可完成HoloLens与ARCore空间坐标的统一.

归纳上面的算法和过程, 将移动端空间共享过程表达为伪代码形式, 如算法4.

算法4. 移动混合现实三维注册空间共享算法伪代码.

输入: 头部坐标位置 $ {C_{\rm camera}} $ , 移动设备空间位置 $ {C_{\rm Android}} $ , 移动设备摄像头位置 ${C_{\rm PhoneCam}}$ ;

输出: 共享目标物体大小、位置、角度及状态数据.

begin

/*初始化位置、大小、角度及状态, 设置变化阈值Threshold $ \leftarrow $ 0.01f及刷新时间updateTime $ \leftarrow $ 0.1f */

1. if Time.time-lastTime > updateTime /*判断位置、大小、角度及状态更新*/

2. 　 if Mathf.Abs(position.gap) > Threshold || Mathf.Abs(rotate.gap) > Threshold || Mathf.Abs(scale.gap) > Threshold || Mathf.Abs(state.gap) > Threshold

3. 　　New Vector3 (px, py, pz) $ \leftarrow $ transform.localPosition

4. 　　New Vector3 (rx, ry, rz) $ \leftarrow $ transform.localEulerAngle

5.　　 New Vector3 (sx, sy, sz) $ \leftarrow $ transform.localScale

6. 　　 foreach var itemin sliders /*遍历滑条值控制目标物体状态*/

7. 　　　Item.onValueChanged.Addlistener(OnSlide)

8. 　　　 for $i \leftarrow 1 \ldots clarity.Length$ do

9. 　　　　New Value(color.a) $ \leftarrow $ material.color.a

10. 　　　 end for

11.　　 end for

12. 　　SendTransformInfo() /*向房间发送目标物体位置、大小、角度及状态信息*/

13.　　 Object[] parameters $ \leftarrow $ New object[px, py, pz, rx, ry, rz, sx, sy, sz, a]

14. 　　SendTransformPhone /*向移动设备发送目标物体位置、大小、角度及状态信息*/

15. 　　ShowPhoneTransform /*同步移动混合现实设备及虚拟模型状态空间信息*/

16.　　 New Transform.localPosition $ \leftarrow $ ARcoreCam.position

17. 　　New Transform.localEulerAngles $ \leftarrow $ ARcoreCam.EulerAngles(rx-90, ry, rz)

18.　　 New Transform.localScale $ \leftarrow $ ARcoreCam.Scale

19. 　　New parent.localPosition $ \leftarrow $ Vector3(px, py, pz)-CamCenterVector3(pos.x, pos.y, 0)

20.　　 New parent.localEulerAngles $ \leftarrow $ Vector3(rx, ry, rz)

21. 　　New parent.localScale $ \leftarrow $ Vector3(sx, sy, sz)

22. 　　New parent. color $ \leftarrow $ Value(color.a)

23.　 end if

24. end if

end

6 实验与结果分析 6.1 3D UNet分割算法实验 6.1.1 数据集

本文采用国际医学图像计算和计算机辅助干预协会(MICCAI)的BraTS挑战赛数据集. MICCAI自2012年开始举办BraTS挑战赛^[43], 在该比赛中, 主办方提供了一个同时具有高(WHO3~4)、低级别(WHO1~2)胶质瘤的脑部MRI数据集, 该数据集中提供每个病例的4个模态的MRI图像序列, 同时提供了每个病例的多区域脑胶质瘤分割区域掩码, 所有数据都是通过医疗专家手工标注得到的真实掩码(the ground truth, GT).

BraTS挑战赛数据集是脑胶质瘤分割领域目前最权威的公开数据集. 在该数据集中, 病人的脑部MRI主要包括四种序列: T1序列(T1-weighted)、T1ce序列(T1 contrast enhanced)、T2序列(T2-weighted)和FLAIR序列(FLAIR). 同时也包含了脑胶质瘤各个子区域的分割掩码——Seg序列^[5]. 在病理上将其分为4个主要肿瘤区域, 即: 肿瘤周围水肿(peritumoral edema, ED), 增强肿瘤区域(enhancing tumor, ET), 非增强肿瘤核心(non-enhancing tumor, NET)以及坏死核心(NCR). 为了更好地在临床上得到应用, BraTS比赛将上述的4种肿瘤区域再合并为3种互相包容的肿瘤区域, 作为分割的真实掩码, 即: 肿瘤周围水肿指定为“全肿瘤(whole tumor, WT)”区域, 掩码值为2 (如图6中绿色区域); 非增强肿瘤核心和坏死核心合并指定为“肿瘤核心(tumor core, TC)”区域, 掩码值为1(如图6中黄色区域); 增强肿瘤区域指定为“增强肿瘤核心(enhancing tumor core, ET)”区域, 掩码值为4 (如图6中红色区域); 对于正常组织和非脑部区域的掩码值为0.

6.1.2 评价指标

在脑胶质瘤分割任务中, 分割结果评估一般基于4种度量方式进行评估, 包括Dice score, 敏感性(sensitivity), 特异性(specificity), 以及Hausdorff距离. 其中Dice score与传统的分类评估指标F值近似相同, 用于衡量脑胶质瘤分割算法的精确程度, 如公式(11):

$ Dice(P, T) = \frac{{2 \times \left| {{P_1} \wedge {T_1}} \right|}}{{\left| {{P_1}} \right| + \left| {{T_1}} \right|}} = \frac{{2 \times TP}}{{2 \times TP + FP + FN}} $

(11)

其中, 将预测结果分为4类, 包括 $ TP $ (true positive)为真阳性, $ TN $ (true negative)为真阴性, $ FP $ (false positive)为假阳性, $ FN $ (false positive)为假阴性. 而预测结果区域 $ {P_1} = TP + FP $ , 真实掩码区域 $ {T_1} = TP + FN $ , 预测背景区域 $ {P_0} = TN + FN $ 和真实的背景区域 $ {T_0} = TN + FP $ .

敏感性(真阳性率, 公式(12))和特异性(真阴性率, 公式(13))表示模型预测正确的背景区域与真实背景区域的比例, 衡量了模型对背景区域的敏感程度.

$ \textit{Sens}(P, T) = \frac{{\left| {{P_1} \wedge {T_1}} \right|}}{{\left| {{T_1}} \right|}} = \frac{{TP}}{{TP + FN}} $

(12)

$ \textit{Spec}(P, T) = \frac{{\left| {{P_0} \wedge {T_0}} \right|}}{{\left| {{T_0}} \right|}} = \frac{{TN}}{{TN + FP}} $

(13)

在脑胶质瘤分割中, Hausdorff距离是用来描述模型预测掩码与真实标签之间相似程度的一种量度, 如公式(14):

$ Haus(P, T) = \max {\text{\{ }}\mathop {{\text{sup}}}\limits_{p \in \partial {P_1}} \mathop {{\text{inf}}}\limits_{t \in \partial {T_1}} {{d(p, t), }}\mathop {{\text{sup}}}\limits_{t \in \partial {T_1}} \mathop {{\text{inf}}}\limits_{p \in \partial {P_1}} {{d(t, p)\} }} $

(14)

在上述4种脑胶质瘤的分割评估方法中, 比较重要的是Dice score和Hausdorff_95距离, 它们分别用于衡量模型是否能够最大程度地分割脑胶质瘤以及是否能精确分割脑胶质瘤. 实验结果对比

6.1.3 实验结果对比

本实验平台的配置主要参数为: Intel(R)Xeon(R)Silver 4210 CPU, 主频2.20 GHz, GeForce RTX 2080 Ti GPU, 64 GB内存, 11 GB显存, 64位Windows 10操作系统以及深度学习算法库PyTorch、图像算法库OpenCV、SimpleITK等.

在实验过程中, 为了体现三维重建3D模型的准确性和验证本文使用的3D UNet的有效性, 进行了不同分割模型结构的对比实验, 使用了FCN-8s、2D UNet、2D Hybrid UNet、2D DeepResUNet与3D UNet进行对比, 首先, 对输入数据的155层二维图像进行预测, 在所有层都经过2D模型预测之后, 再通过创建全黑的三维矩阵将所有层的预测结果进行叠加和WT、TC、ET区域的拼接合并, 以及真实GT标注的WT、TC、ET拼接.

实验的数据集共335份病例样本, 其中训练集285份(包括HGG病例210例、LGG病例75例、无症状病例1例), 每份病例包含4个模态类型和1份医疗专家手工标注得到的真实掩码, 包含测试集50例. 模型训练的输入通道为4, batch size设置为2, 使用Adam进行优化, momentum为0.9, 初始学习率和权值衰减分别为0.03和0.0001, 并使用BCEDiceLoss作为模型损失函数. 在训练过程中使用早停法不仅能够解决epoch数量需要手动设置的问题, 还能有效避免网络发生过拟合, 从而获得较好的实验结果. 使用相同实验参数得到训练各个模型, 在测试集上的预测结果对比如图7所示.

从图7中可以看出, 对于全肿瘤区域、肿瘤核心及增强肿瘤核心区域的分割, 3D UNet都展现了其优越性, 正由于脑胶质瘤MRI图像是多模态的三维切片, 使用三维卷积、三维池化和跳转机制使其Dice Score指标分别能达到0.897、0.849、0.791, 比FCN、2D UNet、2D Hybrid UNet、2D DeepResUNet 这些二维卷积深度学习模型的效果更好如图7(a)所示, 柱形图中蓝色为全肿瘤区域, 橙色为肿瘤核心区域, 灰色为增强肿瘤核心区域, 同时在图7(b)和图7(c)也表明, 3D UNet在Sensitivity指标(WT: 0.897、TC: 0.913、ET: 0.838)及Specificity指标(WT: 0.882、TC: 0.875、ET: 0.828)上也具有较高的精确度. Hausdorff距离是描述该模型预测掩码与标准标签掩码之间的相识距离评价指标, 如图7(d)所示3D UNet在3个区域的Hausdorff距离分别为1.125、0.855、1.257. 图7的实验结果充分表明3D UNet更有效地获取了图像的浅层特征和深层特征, 比FCN-8S、2D UNet等二维网络有更好的检测性能.

同时, 本文验证了不同batch size对脑胶质瘤检测分割结果的影响, 分别完成了batch size=2, 4, 16时, 不同模型对全肿瘤区域(WT)分割效果, 如图8所示. 从图中可以看出, 增大batch size能一定程度上提高分割效果, 但总体上2D网络的分割效果还是差于3D UNet的分割效果. 此外, batch size参数的选择影响模型的优化程度和速度, 因此, 综合上考虑, 本文设置batch size为2, 对脑胶质瘤进行检测分割.

本文进一步将3D UNet模型与近年来一些典型文献在BraTS挑战赛数据集下的实验结果进行对比, 主要考虑最重要的Dice Score指标和Hausdorff_95距离指标, 其结果如表1所示. 从表1可以看出, 本文3D UNet模型在Dice Score指标上对于肿瘤核心(TC)的分割效果好于目前典型文献中的方法, 而在对全肿瘤区域(WT)与增强肿瘤核心区域(ET)的分割效果基本接近于典型文献中的方法; 但是, 在Hausdorff_95距离指标上, 本文3D UNet模型对于全肿瘤区域(WT)、肿瘤核心(TC)和增强肿瘤核心区域(ET)的分割效果明显优于目前典型文献中的方法. 因而, 整体上表明本文3D UNet模型具有较好分割性能.

为了进一步从三维可视化角度直观展示分割效果, 本文从测试集中选取几例病例在FCN、2D UNet、2D Hybrid UNet、2D DeepResUNet与3D UNet网络结构下的分割结果分别进行掩码可视化, 与专家标注的真实掩码(GT)进行对比, 如表2所示, 其中第1列为真实掩码可视化结果, 最后一列为3D UNet分割的可视化结果. 从表2较清晰地表明, 在肿瘤核心区域(TC, 黄色区域)和增强肿瘤核心区域(ET, 红色区域)上, 3D UNet模型的预测掩码更接近各病例的真实掩码GT标准, 在全肿瘤区域(WT, 绿色区域)的分割效果也比其他模型稍好些. 总体上, 3D UNet在3个区域上均有较好的分割效果, 并有效地减少各个区域错分和多分的情况, 与真实的脑胶质瘤掩码分布基本一致, 为后续MR系统提供了更加精确的模型来源.

6.2 三维重建实验

基于面绘制算法, 将3D UNet网络模型检测和分割得到脑胶质瘤掩码, 同时将脑胶质肿瘤结构的阈值随掩码输出, 利用不同颜色蒙版材质来表示不同子区域, 使用“Threshold”“Flood filling”“Smoothing”等操作通过设定好的阈值完成掩码到3D模型的面绘制, 当阈值不同其模型重建效果也不同, 如图9(a)当阈值范围为100–200时, 重建结果为大脑表层; 当阈值范围为200–300时, 大脑组织结构基本重建如图9(b); 随着阈值范围的不断增大, 大脑重建过程由外到里展现不同的组织区域, 直到阈值范围为500~600时, 可清晰看到全肿瘤结构组织如图9(e). 本文将阈值范围设置为20~700, 完整地重建整个脑组织模型, 如图9(f)所示, 以便在MR系统中对大脑及脑瘤的组织结构进行全面的观察.

由于3D UNet输出掩码结果为像素级特征图, 同一肿瘤区域的掩码值相同, 意味着其阈值也相同. 肿瘤重建结果与大脑组织结构重建模型如图10(a)所示, 绿色区域为阈值是1时的全肿瘤区域模型(掩码值为1), 同样地, 黄色区域模型为阈值是2时肿瘤核心区域模型(掩码值为2), 红色区域为阈值是4时的增强肿瘤核心区域模型(掩码值为4), 通过图右下角的脑胶质瘤肿瘤在水平面(左上)、冠状面(左下)及矢状面(右下)切片, 同时通过滑条可观察脑胶质瘤各个子区域在脑部切片区域的具体位置.

为了对比, 本文选用医学体数据的光线投射法^[7]重建脑胶质瘤的三维模型, 其基本原理是: 从图像的每一个像素, 沿着视线方向发射出一条光线, 光线穿越整个图像序列, 沿着光线方向等距离采样, 利用插值计算出图像序列的采样点颜色值与不透明度, 同时根据光线吸收模型进行累加, 最终获得渲染图像的颜色和不透明度. 当观察方向发生变化时, 采样点的前后关系也将变化, 需要重新进行采样, 因而其计算量极为庞大. 本实验三维重建结果如图10(b)所示.

面绘制渲染重建后的大脑和肿瘤模型实际上是一个导出时子区域完全分离的3D模型, 并且属于同一坐标系下. 体绘制在重建过程中不产生等值面, 本质仍为三维体数据场在二维面的图像, 导致得到的三维模型无法导出到混合现实框架中作为虚拟对象交互使用^[52,53], 因此只能采用面绘制的方式获取MR系统实验的三维模型.

由于面绘制结果是将整个三维模型以三角面的形式表现. 虽然场景中的3D模型的面数越多, 在视觉效果上越精细动人, 但也意味着会消耗更多的计算和存储资源, 对设备的负担越大, 导致帧速率大幅降低, 会出现画面延迟卡顿现象. 为了提高渲染效率和实时性, 需要对模型进行优化, 即在保证逼真效果的情况下, 减少三角面数, 达到与模型最逼近的效果, 同时减小模型内存大小, 使导入MR设备中的模型视觉效果最佳. 因此, 通过Blender将面绘制STL格式转换为FBX格式, 进一步对3D模型进行优化.

本文使用几何元素删除法中的顶点删除及区域合并方法模型进行优化, 在循环遍历模型全部顶点后, 删除那些影响较大的顶点, 然后对删除区域进行重新网格化, 并对网格孔洞进行修护, 优化结果如表3所示.

表3中最后一列为本实验三维重建模型的原始顶点数和面数, 其3D显示的每秒传输帧数(fps)仅为12. 当模型三角顶点数量减少约70002、三角面的数量减少约140000、优化程度为39.48%时, 模型的纹理与细节都有很好的保留, 并且实时画面稳定, 帧率为34, 满足系统的实时性要求; 当模型三角顶点数量约为35002、三角面数量约为70000、优化程度为69.74%时, 帧率为55, 此时实时性较高, 但模型表面多处突起慢慢消失, 趋于平整; 当模型三角顶点数量减少约1002、三角面的数量减少约2000、优化程度为99.14%, 此时模型纹理基本消失, 圆润程度有了很大改变并产生了许多突刺; 模型三角顶点数量减少约118、三角面的数量减少约232、优化程度为99.90%时, 模型可清晰的看到面面之间的共同棱边, 突起明显. 因此, 在同时兼顾模型渲染质量和系统运行性能情况下, 实验使用优化程度为39.48%的模型.

6.3 脑胶质瘤医学影像分析MR系统实验 6.3.1 医学影像混合现实系统设计

医学影像混合现实系统可用于脑胶质瘤患者病情分析与辅助诊断, 外科手术术前讨论与手术规划; 可在手术过程中实时指导肿瘤位置, 也能用于医学辅助教学, 帮助医生掌握手术过程中的关键技术. 系统设计流程框图如图11所示. 首先, 从患者头部MRI影像获取脑胶质瘤空间信息, 经过基于深度学习的智能检测与分割算法, 获得肿瘤主要区域的分割掩码, 然后将大脑及脑胶质瘤掩码进行3D虚拟模型绘制重建, 同时创建混合现实环境, 导入所需SDK工具包, 将重建大脑模型及肿瘤子区域进行全息可视化显示. 在此基础上, 开发MR系统人机交互界面, 主要包括交互手势、3D运动、方向拖拽及透明度可视化等, 将大脑及脑胶质瘤病灶区域融合到真实患者的头部位置, 进行全方位检查与诊断, 并从不同视角空间同步共享.

脑胶质瘤医学影像分析MR系统硬件平台主要包括Microsoft HoloLens头盔和Android手机. HoloLens设备的硬件参数为Intel 32位全息处理单元(HPU1.0)CPU, HoloLens Graphics GPU/HPU、64 GB闪存和2 GB RAM、1 GB显存、Windows 10.0.17763.1879 ×32位. Android手机的主要参数为: 小米8、高通骁龙845 8核2.8 GHz处理器、6 GB运行内存、64 GB存储器、2248×1080分辨率、Android 10操作系统. 软件开发工具主要包括高通AR开发工具包Vuforia SDK、谷歌AR工具包ARCore SDK、Unity 3D等

6.3.2 第1视角虚实融合实验

用户戴上HoloLens设备, 使用脑胶质瘤MR系统, 获得第1视角场景. 当启动MR系统后, 大脑和脑胶质瘤的3D虚拟模型就会出现在用户第1视角中, 如图12(a)所示, 用户可通过滑条对重建的3D虚拟模型实现整体透明及肿瘤全透明, 在不遮挡真实场景情况下进行虚实融合; 图12(b)展示了不同脑胶质瘤子区域在脑部区域的位置及大小. 如图12(c)所示, 在融合阶段用户可以使用“轻击并按住”手势获取目标框, 进行平移和旋转运动获取虚拟模型坐标及真实世界坐标系, 并与大脑3D虚拟模型进行互动, 随后借助混合现实设备的空间扫描能力提取稀疏特征进行目标识别, 将目标渲染到真实世界并实时跟踪进行虚实融合. 基于无标识的注册渲染将虚拟模型固定在患者头部上, 并保持静止, 通过调整模型的透明度可选择最佳的透明放置, 以手势旋转缩放匹配不同位姿的融合. 在整个建模和实施步骤中都会保留3D虚拟肿瘤的比例及位置信息, 因此可以更好地使用模型投影进行外科手术的术前规划及术中指导.

6.3.3 移动手机端第3视角空间同步设计

(1) 人机交互界面搭建

为了实现HoloLens端与Android手机端第3视角空间实时共享, 在Unity 3D中通过Canvas容器设计人机交互界面, 添加HoloLens端与Android手机端的通信连接与视觉共享的界面按钮. 如图13(a)为HoloLens端系统的UI主界面, 包括创建Socket服务器房间和上传、下载空间锚点, 为Vuforia扫描Android端识别图界面提供了识别按钮, 同样也可对虚拟模型进行交互手势、3D运动、方向拖拽及透明度可视化等功能. 图13(b)为Android手机的子界面, 主要实现通信连接和断开功能. 同时为了Android手机端更好的共享体验, 设置了隐藏UI的最小化按钮, 扩大可用屏幕宽度并突出显示效果等功能.

(2) 同步信息处理

混合现实设备HoloLens和支持ARCore的手机各自具有空间连续定位功能, 将两种设备统一空间坐标并通过建立Socket网络通信方法同步大小、位置、角度和状态数据. 同时混合现实系统通过创建空间锚点来确定虚拟模型位置, 一个空间锚点代表世界中一个必须时间追踪的重要点. 当设备获得更多的空间数据时, 空间锚点在需要时会根据另一个锚点来调整位置, 以保证每一个锚点在真实世界中保持精确的位置. 将3D模型的其他数据与空间锚点关联, MR系统会将3D虚拟模型相关信息保存在其锚点处, 用户就可在后续的使用中找到锚点处的虚拟模型的相关信息. 将空间锚点及其相关环境信息从HoloLens端传输到其他移动设备端, HoloLens里的虚拟模型便能分享到移动设备端(Android手机), 两个设备就能识别到同一个位置, 实现空间位置共享.

(3)空间位置校准

基于外科手术在手术过程中对模型精度和模型位置叠加契合度的要求, 本研究的第3视角测试采用Android端在识别图注册时的位置信息作为系统的评价标准, 旨在为手术精准导航. 位置校准实验过程为: 佩戴HoloLens至舒适状态后打开Android手机程序准备第3视角静态误差测试, 在Android端开发过程中事先准备了与手机尺寸大小(154.9 mm×74.8 mm×7.6 mm)相同的虚拟模型, 在更改Android端摄像头相对手机中心位置( $ pos.{x_p} $ , $ pos.{y_p} $ )后, 点击Android端UI界面上的手机坐标按钮, 将手机虚拟模型固定放置(记为(0 mm, 0 mm)), 记录每次放置位置与真实手机放置位置坐标, 相同相机位置放置5组数据, 保存虚拟模型偏离真实手机的测试结果平均值于表4. 当手机相机位置距离手机屏幕中心坐标为(−61 mm, 25 mm)时模型偏离较远向上偏离7.4 mm、向右偏离1.5 mm, 随着参数设置越接近手机摄像头真实位置, 模型偏离越小, 三维注册效果越好.

6.3.4 脑胶质瘤MR分析实例

本实验以脑胶质瘤医学影像的诊断、分析为案例, 显示外科主治医生利用MR系统与医生团队、患者、患者家属等分析脑胶质瘤的状况, 有效实现病情沟通、手术术前规划等. 主治医生戴上HoloLens混合现实设备, 启动HoloLens端的MR程序, 在UI界面使用“轻击”手势, 点击连接Socket服务器并创建房间, 等待其他用户的登录. 其他医生及病人家属通过Android手机启动MR应用程序, 在弹出的客户端UI界面, 点击连接服务器并加入, 此时手机位置为第3视角同步的初始位置. 在HoloLens端通过UI界面的功能按钮控制脑胶质瘤患者的3D虚拟模型, 同时手机端可从不同的位置观测各子区域在患者头部情况, 完成视角共享, 典型场景如图14所示.

如图14(a)所示, 主治医生创建房间并为患者家属或协助医生分析MRI图像上脑胶质瘤的大小、区域位置等情况, 方便家属和协助医生理解, 有效地解决医生与患者或医学信息量不足的医生之间的沟通不畅问题. 在图14(b)中, 主治医生通过操作手势以MR系统第1视角为患者及协助医生分别展示了增强肿瘤区域、肿瘤核心区域和全肿瘤区域结构及脑部肿瘤整体区域的三维模型, 通过第1视角中的人机交互界面可分别调节肿瘤区域透明度及脑部整体透明. 同时, 用户通过手机端直观地了解患者肿瘤子区域大小及在脑部位置, 如图14(c), 该方式便于制定可靠的术前计划, 同时基于手机端第3视角的移动MR设备极大地节约了手术设备成本, 为外科手术术中导航提供便利.

为了分析MR系统性能, 本文应用Windows Performance Analyzer对上述MR实验系统CPU占用率和系统实时性进行测试. 如图15所示, MR系统将在HoloLens设备上消耗大量的计算资源. 图15(a)是系统运行所有进程的CPU使用情况, 其中MR系统及其服务的CPU使用率总体保持平稳. 图15(b)是MR系统应用程序和按钮操作进程单独进行分析的结果. 在前28 s系统进行房间创建、透明度调整、模型缩放和旋转操作, MR系统在CPU使用率上有明显的抖动, 但整体保持平稳. 在40 s时, Android端用户登录加入房间, CPU占用率为17.07%. 在56~62 s时, MR系统的HoloLens端完成扫描识别图, 同步空间信息, CPU使用率由平均19%上升到30%左右. 运行过程中, 系统没有发生明显的卡顿现象, 保证了系统交互过程的流畅性和虚实融合的实时性.

7 结　论

本文研究的医学影像智能诊断与混合现实技术系统, 以辅助诊断、三维重建、MR三维注册及实时同步共享为重点, 协助神经外科医生进行脑胶质瘤的诊断、检测与切除手术. 本文提出基于深度学习模型3D UNet的MR脑胶质瘤定位与区域分割算法, 采用基于面绘制的方法对脑胶质瘤影像进行多结构组织的三维绘制与优化, 提出了交互式无标识的移动混合现实三维注册跟踪与视觉空间共享算法, 设计、实现了脑胶质瘤医学影像分析MR系统, 有效支撑脑胶质瘤的辅助诊断与治疗过程, 极大地减轻了外科医生审视MRI图像的视觉疲劳程度, 提高医生的审片效率, 同时也提高了病灶区域定位的准确性. 随后通过重建大脑和肿瘤子区域后, 将其实现到HoloLens混合现实设备和Android设备中, 有助于神经外科医生们直观了解肿瘤范围, 从而可以对肿瘤进行几乎完全的切除, 提高手术的成功率, 降低手术风险.

本文最后实现了脑胶质瘤医学影像分析MR原型系统, 在实验室环境下展示了其工作过程与原理, 目前的实验显示了良好的应用前景. 然而, 目前的三维注册虚实融合技术存在一定的技术制约, 例如空间扫描精度注册配准上需要进一步的改善, 高性能的硬件设备也能够有效提高准确率. 同时MR系统的实时性、术中指导的即时性等问题, 需要更多实际的临床测试与评估. 下一步将进一步完善系统功能与性能, 优化检测分割模型, 研究检测分割、三维重建与MR系统的快速融合接口技术, 加强远程辅助诊断功能, 为精准微创治疗提供有效的可视化平台.