基于磁导航的徒手超声三维重建

杨凡帆,张勤俭,范巨峰,李海源,晁明辉,褚浩杰

(1.北京信息科技大学 仪器科学与光电工程学院,北京100192;2.北京信息科技大学 机电工程学院,北京100192;3.首都医科大学附属北京朝阳医院,北京100020;4.北京邮电大学 自动化学院,北京100876)

0 引言

超声成像是目前临床诊断和手术中最流行的医疗成像方式之一[1],将超声影像结合医学图像处理技术和三维重建技术可以获得更丰富的病灶信息。近年来,徒手三维超声成像因其扫描方式灵活、能提供更大的视野以及便于医生操作,而成为超声影像引导手术的主要研究方向,在图像引导类以及导航类的微创手术中具有实际意义。

常用的三种获取超声影像的方式为机械式扫描[2]、徒手扫描(带位置传感器[3-4]或不带传感器[5-6])和容积探头方式[7]。其中,带位置传感器的徒手超声三维重建方式鲁棒性更强、更方便、成本更低,目前使用更为广泛[8]。这种方式通常是由传统的超声设备结合某种空间定位装置来获取一系列二维超声图像及其空间方位信息,再结合标定技术和三维重建技术,实现超声影像的三维重建。市面上已经有一些较为成熟的空间定位方法,包括机械(臂)定位、声学定位、电磁定位和光学定位[9]。其中,电磁定位体积小、无遮挡、精度和价格适中,比较适用于术中手术设备需要移动空间且移动范围不大的穿刺手术。

徒手扫描方式获取到的是空间中非规则分布的二维序列图像,与计算机断层扫描(computed tomography,CT)、磁共振成像(magnetic resonance imaging,MRI)等获取的规则分布的断层图像不同,超声三维重建需要先计算二维图像到超声探头的转换信息和一系列空间坐标转换,这一过程即为超声探头标定。探头标定的基本原理是通过扫描一个已知几何特征及其空间位置的模型,根据这些特征在超声图像中和在实际模型中的位置关系,来求解未知标定矩阵中的参数。标定流程包括标定模型制作、图像特征提取和标定方程求解3个步骤[10]。根据标定模型内特征的形状可将标定模型分为点式、线式和面式三大类,经典的标定模型有交叉线形、三线形、N线形等[11-14]。N线形模型制作简单,便于扫描,且应用广泛,符合本文对重建的要求。

徒手超声重建是将超声探头扫描图中的二维像素点重建到三维体数据中,重建流程包括:构建体数据、像素的重新分配和计算体数据值。构建体数据是预先确定一个三维重建空间,可以通过主成分分析法和包围盒等方式。分配像素值是根据求解的坐标变换,将二维平面上的像素插入三维空间中,可以利用不同的插值方法将像素分配到重建区域。计算体数据值是对插值后的三维体进行填充,将空值补充完整,得到一个完整的三维数据。如何在图像质量和计算时间之间取得平衡是目前三维成像重建算法面临的挑战[15]。

1 设备和坐标系统

图1展示了超声重建过程,先扫描标定模型完成探头标定,再扫描待重建区域进行三维重建。

图1 磁导航徒手超声三维重建流程Fig.1 Process of freehand ultrasound 3D reconstruction based on magnetic navigation

通过视频采集卡采集B超探头画面,将磁导航位置传感器固定在B超探头上,在模型表面模拟医生扫体,每一幅超声图像都对应一个空间方位数据(x,y,z,α,β,γ)。其中,(x,y,z)代表位置传感器相对于电磁发射器的坐标,单位为mm;(α,β,γ)为位置传感器相对于电磁发射器的方位角、俯仰角和滚转角,单位为°。利用方向和位置信息构建三维重建空间,再对该区域进行插值重建,最后将三维数据可视化得到超声的三维影像。

超声影像重建中涉及的坐标系如图2所示。空间导航设备的位置传感器被固定在一个传统的超声探头上。重建需要将二维超声图像中的像素点转化到三维模型空间。具体过程是先将图像(像素)坐标系(用p来表示)转化到传感器坐标系(用s来表示),然后转化到电磁发射器坐标系(默认为世界坐标系,用w来表示),再转换到模型坐标系(用m来表示),由此可将像素点放入三维模型空间。得到扫描序列图之间的实际物理关系,重建结果才有实际意义。

图2 坐标转换过程Fig.2 Coordinate conversion process

将超声帧中的一个像素(u,v)转换到模型中的一个点(x,y,z)的过程如式(1)所示(默认模型坐标系为重建坐标系)。

M=Tm←w·Tw←s·Ts←p·P

(1)

式中:Ta←b代表从坐标系b向坐标系a的空间变换矩阵;P为像素点在二维超声图中的坐标;M为像素点在模型中的坐标。

2 探头标定实验

2.1 标定模型

标定模型选择使用最广泛且更适用于临床的N线形模型,本文基于N线形模型思路提出一种“三层简易N线”标定模型,以少量的点计算出图像中多组横纵方向上的数据,以此降低标定环节的复杂度。标定模型的壳体是使用树脂材料通过三维打印的方式制成,4条鱼线分为3层固定在壳体内部,第1层是两条平行线,两线间距|AC|=|BD|=25 mm,其余两层是方向相反的斜线,每层之间距离相等,间距为10 mm。以模型中心位置为例进行说明,扫描时超声探头方向如图3(a)所示,探头与AB垂直且与AC平行,E、F、G、H为超声平面在模型中心位置与鱼线的交点,对应的俯视示意图如图3(b)。扫描结果如图3(c),提取标定点后得图3(d)。为便于后续的说明,连接GH,点F投影到线段GH上得到点F′,可知|EF|=|FF′|=1/2|EF′|。

图3 标定模型和标定点Fig.3 Calibration model and calibration point

2.2 标定计算

本文将点的连通域中心作为其坐标,标定的目的是计算不同坐标系之间的变换关系。

变换矩阵Tw←s可由位置传感器的输出数据(x,y,z,α,β,γ)得到,如式(2)所示:

(2)

将式(1)展开为如下形式:

(3)

式中:sx和sy分别为横向像素比例和纵向像素比例,两者均未知;Tm←w和Ts←p未知。未知参数和矩阵的求解顺序为:首先,通过标定点之间的关系求出超声图像的像素比例sx和sy;然后,通过模型中的点在传感器坐标系下的对应关系求出Tm←w;最后,结合之前的求解结果,通过标定点在图像和模型中的对应关系求出Ts←p。

2.2.1 求解像素比例

超声图像的像素比例单位为mm/像素,即图中的一个像素代表多少毫米。不同种类探头的像素比例不同,且一幅超声图中横向和纵向的比例也不同,计算式为

(4)

式中:X为超声图像中像素的横向距离,x为模型中与之对应的距离;Y为超声图像中像素的纵向距离,y为模型中与之对应的距离。

从超声图像中计算出的距离和在标定模型中测量得到的距离均存在一定误差,在实际求解像素比例时可选取n组距离进行计算,求其平均值作为最后的结果,如式(5)所示:

(5)

2.2.2 求解Tm←w

首先通过测量得到标定模型上的点在模型坐标系下的坐标M和在传感器坐标系下的坐标S,然后通过转换矩阵Tw←s将S转换为电磁发射器坐标系下的坐标W,最后利用最小二乘法求解式(6)中的坐标转换矩阵Tm←w。

M=Tm←w·W

(6)

2.2.3 求解Ts←p

对式(1)移项得式(7),可用多幅图中的F点在图像和模型中的坐标关系求出Ts←p的最小二乘解。

(7)

式中:PF为F点在超声图像中的坐标(sxu,sxv),可直接通过计算得出。MF为F点在模型中的坐标(xF,yF,zF),可利用相似三角形原理得到式(8),则(xF,yF,zF)如式(9)所示,其中z坐标值固定为15.5 mm。

(8)

(9)

至此,标定环节结束,得到了图像像素点到模型坐标系下的变换。由此空间坐标变换可以将像素点分配到三维重建空间。

3 三维重建及可视化

包围盒(bounding box)技术是构建重建空间的简易方法,它通过填充一系列二维超声帧的像素,再找到最大点和最小点来计算体积,构成一个由空体素组成的三维体。基于此,本文提出简易快速包围盒方法:在一次扫描得到的序列中直接利用所有坐标中的最小值和最大值快速确定重建区域范围。如图4所示。

图4 简易快速包围盒算法流程Fig.4 Simple and fast bounding box algorithm flow

根据插值方式的不同可以将超声影像三维重建分为基于体素的方法(voxel based method,VBM)、基于像素的方法(pixel based method,PBM)和基于函数的方法(function based method,FBM)。VBM遍历每个体素,选择其周围的像素值填入,常用的最近邻域体素法(voxel nearest neighbor,VNN)将距当前体素最近的像素的值插入。FBM对二维超声图像和其对应的三维重建空间建立函数关系,将像素值映射到目标体素中去。PBM遍历图像中每个像素,将像素值分配给体素。最邻近像素算法(pixel nearest neighbor,PNN)是一种经典的PBM,也是本文使用的重建方式。图5展示了PNN的两个基本步骤:像素重新分配和孔洞填补。步骤①遍历二维像素值,通常与权值一起存储,计算最邻近的一组像素集合的加权平均值赋值给三维体素。步骤②遍历三维体素,对邻域内非空体素进行加权平均后赋值给相应的空体素点。

图5 PBM流程示意图Fig.5 Diagram of the PBM process

重建的结果是一个三维格式,常见的可视化方法分为两大类:基于面绘制的方法和基于体绘制的方法。本文采用经典的光线投射(ray casting,RC)算法进行可视化。当光线穿过体数据时,RC方法会沿着光线方向对经过的每个体素的颜色和阻光度进行积累,最后在平面上绘制出投影。

4 实验

实验采用NDI的3D Guidance trakSTAR 电磁导航定位仪,深圳迈瑞公司DP-7 超声设备,使用epiphan(DVI2USB 3.0)的视频采集卡连接电脑采集B超探头画面。

本文使用图3(d)计算像素比例,将相应的数值代入式(5)计算像素比例得式(10):

(10)

式中:XGF、XFH、XGE、XEH、XGH表示超声图像中两点之间的横向距离;YEF、YFG、YFH、YEG、YEH表示超声图像中两点之间的纵向距离。由图3(a)可知,上述距离对应其在标定模型中的实际距离为12.5 mm、12.5 mm、12.5 mm、12.5 mm、25 mm和10 mm、10 mm、10 mm、20 mm、20 mm。

表1记录了求解坐标转换矩阵Tm←w时所选取的标定模型上的8个点分别在两个不同坐标系下的空间位置。在计算时需要统一单位。先利用Tw←s将S下的点转换到W下,再将表1中的点带入式(6)中求出Tm←w。

表1 点在模型坐标系下的位置和用传感器测得的位置

通过式(8)和式(9)计算出F点在模型中的坐标MF,在求解Ts←p时需计算多幅超声图像的F点在模型中的坐标并带入式(7)求最小二乘解。至此得到标定关系Tm←wTw←sTs←p。

标定探头时可将探头固定在滑轨上,保持探头平行于AC和BD且垂直于模型底面,由AC边扫描至BD边。选取扫描后较为清晰的31幅超声图像,其中10幅图像用于上述标定计算,剩余21幅图像用于测试标定误差。比较F点在标定模型中的实际坐标值与计算坐标值,以此来测试探头标定误差。图6为误差曲线,x轴的平均误差为3.28 mm,y轴的平均误差为0.53 mm,z轴误差忽略不计(探头固定在滑轨上平行移动)。以较小的妇科肿瘤为例,肿瘤直径一般是2～3 cm,穿刺的误差要求通常在毫米级别内,医生通常会努力将穿刺误差控制在2 mm或更小的范围内。本文探头标定误差(不含z轴)平均为1.90 mm,初步达到了手术需求。

图6 误差分析Fig.6 Error analysis

因为不同病灶区和影像设备的清晰度不一样,像核磁或CT影像清晰且易分割,但超声画面伪影大、噪点多,往往很模糊。若提前将目标区域分割出来会减少重建压力。本文通过点选出目标区域,再利用样条曲线拟合轮廓后进行阈值分割。RC方法和最大密度投影(maximum intensity projection,MIP)方法是基于体绘制的两种常见方法。本文用球体模型进行实验,用两种可视化方法对比分割前后的重建效果。结果如图7所示。

图7 目标模型的可视化结果Fig.7 Visualization of the target model

根据结果可以明显看出:模型内部有干扰介质,未经分割的模型不能很好地展示目标区域,分割后重建效果更好,且RC方法的可视化效果更好。

5 结束语

图像引导的微创介入手术逐渐成为当今接受度更高的手术方式之一,也一直是国际前沿的研究热点之一。本文对徒手超声三维重建的关键技术进行了研究,主要包括超声标定、超声重建及三维可视化等内容。研究改进了标定模型,简化了标定环节,并优化了构建重建区域的方法,在一定程度上简化了整个徒手超声重建过程。

为了提高重建效率,可以考虑从以下几个方面进行优化,如:在数据采集完毕后对图像进行增强或分割来优化重建范围;优化标定模型,降低标定环境的复杂度;对重建和可视化的算法公式进行优化。

在未来的工作中,将会采集更多超声数据,并把这种简易的标定思路运用到自动标定的研究中。此外,还会基于超声影像研究自适应的分割方式,以此来提高重建的效率。