Freehand三维超声体数据重建的最新研究进展

范阳阳 倪 倩 温铁祥 辜 嘉 王 磊 谢耀钦 刘金根

1(中国科学院深圳先进技术研究院，广东 深圳 518055)2(武汉理工大学信息工程学院，武汉 430070)3(深圳市福田区中医院(广州中医药大学深圳医院)，广东 深圳 518034)

Freehand三维超声体数据重建的最新研究进展

范阳阳1,2倪 倩3温铁祥1#*辜 嘉1#王 磊1#谢耀钦1#刘金根2

1(中国科学院深圳先进技术研究院，广东 深圳 518055)2(武汉理工大学信息工程学院，武汉 430070)3(深圳市福田区中医院(广州中医药大学深圳医院)，广东 深圳 518034)

Freehand三维超声成像系统因其扫描方式自由、能提供更大的成像视角和更高的图像分辨率等优点，比较符合医生习惯和手术室环境，成为近年来超声影像引导介入手术的主要研究方向。Freehand三维超声从一系列不规则的二维B超图像入手，重构器脏结构的三维体数据，并进行三维渲染显示。体数据重建是Freehand三维超声成像系统的关键技术，对提高重建图像质量有着重要的作用。对Freehand三维超声体数据的重建算法进行归纳和分类，比较分析其中的典型算法，最后对Freehand三维超声体数据重建算法的研究状况和未来的发展方向进行展望。

Freehand三维超声；体数据重建；重建算法

引言

三维超声成像技术应用于临床已有50余年时间，三维超声成像的概念最初由Buan和 Greewood在1961年提出，他们采用一系列平行的人体器官二维超声截面图像，用叠加的方式得到这个器官的三维图像。这个概念被提出之后，越来越多的研究人员投入到相关领域的研究工作中。目前，三维超声成像方式主要有两种[1]：一种是基于二维面阵探头(也称“三维容器探头”)的容积成像方法，如图1(a)所示。该方法能直接得到扫描部位的三维容积数据，成像实时性高。与传统二维超声探头相比，容积探头的主要缺点是成本高、扫描深度、角度以及图像分辨率有限。另一种是利用传统的二维超声诊断设备，结合某种空间定位机构(主要包括机械臂定位、光学定位以及电磁定位)来获取一系列二维超声图像以及相应的空间位置信息，然后通过算法重建出三维超声体数据。后者是目前实现三维超声图像引导介入手术的主要方法，主要包含3个阶段：二维超声图像的获取、三维体数据的重建、三维渲染显示。二维超声图像的获取有两类方法：一是预定轨迹采样，也叫“机械运动扫描法”，主要是将超声探头固定在机械臂上，由计算机控制步进电机，以驱动探头按照事先规划好的数据采集轨迹和采样频率，以一定的位移或角度增量步进，对器官进行采样。用这种方式所得到的图像间的位置或角度关系都是事先确定的，具体的采集轨迹有线性扫描、扇形扫描和旋转扫描3种形式，如图1(b)所示。二是随机采样，也叫自由扫查(Freehand)，即由操作者手持带空间定位装置的超声探头，按一定的顺序、以相对自由的方式对器脏进行采集，采集到的图像位置和角度都是相对任意的，如图1(c)所示。从临床应用角度来看，Freehand扫查方式更易于与手术器械相结合，在超声图像引导介入手术方面，Freehand 三维超声有着广泛的临床应用前景。

图1 三维超声成像扫描方式。(a)容积探头方式[9]；(b)电机驱动方式[10]；(c)自由式方式[1]Fig.1 The scanning of 3D ultrasound imaging: (a) 2D matrix probe-based scanning[9]; (b) Mechanical scanning[10]; (c) Freehand scanning[1]

近年来，Freehand三维超声成像系统被广泛地应用于医学成像领域。2003年，葡萄牙的Combia大学研制了VOLUS系统[2]，用于人体心脏左心室诊断和研究；2011年，Inglis等研制了一种超声内镜系统[3]，用于对人体食管进行自由式的三维超声重建。在国内，2005年，海军总医院研究了三维超声导航系统在神经外科手术中的应用[4]；2007年，清华大学研制了用于肝癌介入治疗的三维超声导航系统[5]；2011年，哈尔滨工业大学在基于三维超声图像的穿刺手术机器人辅助系统方面进行了研究[6]；2013年，吉林大学对经阴道三维超声在妇科疾病诊断中的应用进行了深入研究[7]。2014年，Jiang等对半自动化三维注释的乳腺超声成像系统的开发和可行性进行了研究[8]。因此，对Freehand三维超声成像系统的研究具有重要的科研价值，引起人们的广泛关注。

1 Freehand三维超声重建

1.1 三维超声成像原理

Freehand三维超声成像，是利用传统的二维超声扫描设备，结合光学或磁定位系统，通过Freehand扫描方式，获得一系列二维超声图像以及相应的空间位置信息，然后进行三维超声体数据的重建，最后对重建的三维体数据进行渲染显示。其中，三维体数据重建是实现Freehand三维超声高精度成像的关键技术环节之一。根据重建目的的不同，超声图像三维重建可分为两大类：一类是基于表面重建的方法，另一类是基于体数据重建的方法。第1类方法需要对脏器做精确的轮廓分割，无法显示脏器内部的详细解剖内容，故很少被临床采用。第2类方法是用一定数目的体素，按相应的空间位置排列，构成三维立体图像，可对人体结构的所有组织信息进行重建，是目前使用较多的三维重建方法。

目前，用于临床医学B超扫描的超声探头有许多种，它们有不同的分类方式。例如，按波束控制方式分为线阵型、扇型、凸阵型和相控阵等，按诊断部位的不同可分为眼科探头、心脏探头、腹部探头、颅脑探头和儿童探头等。在Freehand三维超声成像中，根据扫查部位的不同，选择相应的超声探头。例如，浅表结构(眼、甲状腺、颈部血管、睾丸)一般采用线阵探头进行平行扫查，腹部(肝、胆、胃)、胸腔、妇科一般采用凸阵探头进行扇形扫查。对器脏进行扫描，得到二维超声图像序列，然后利用体数据重建算法重构器脏的三维图像。例如：Laurence等利用相控阵探头扫描人体脑部，得到脑肿瘤的二维B超图像，然后利用基于像素的方法重建脑肿瘤的三维图像[11]；Huang等用线阵探头对人体前臂进行自由式扫描，之后利用基于体素的自适应均方距离加权方法对其进行三维重建[12]；郭强等人用扇形探头对人体心脏进行旋转扫描，得到心脏的二维超声图像，之后利用基于函数的马尔科夫随机场插值方法重建心脏的三维结构[13]；Huang等用凸阵探头对胎儿模型进行自由式扫描，得到二维超声图像，之后利用基于函数的贝塞尔插值法对它们进行三维体数据重建[14]等。

三维超声成像因其图像显示直观、精确测量结构参数、准确定位病变组织和缩短医生诊断所需时间等优点，被广泛地用于人体器官病变的诊断，如乳腺癌、肝脏肿瘤、胎儿、颈动脉疾病、肾脏肿瘤等。近年来，研究人员对三维超声成像提出了一些新的方法。例如，目前对于乳腺癌的初步筛查和诊断主要有3种方式[15-16]：基于二维超声(Ultrasound，US)、三维磁共振(Magnetic Resonance Imaging, MRI)以及二维压缩投影的乳腺钼靶(Mammography，MAM)设备。最近研究人员提出三维超声断层成像技术(three dimensional ultrasound tomography，3D-UST)，并基于乳房的外部曲面结构和内部组织特征融合3D-UST和3D-MRI影像，为乳腺癌的筛查和诊断提供更准确的资料。3D-UST是通过上下移动环型超声探头阵列来逐层扫描乳房切面，基于数据重建得到最终的三维乳腺影像[17-18]。相比于US、MRI和MAM，3D-UST的相对优势有：一是成像时间短，影像分辨率高，且不需要挤压乳房；二是3D-UST生成3组数据，分别反映乳腺组织的超声波反射信息、传播速度和衰减信息，能够用于区分肿瘤；三是3D-UST在肿瘤良恶性区分和乳腺密度测量上与MRI和MAM具有高度的一致性。因此，UST有逐步替代US、MRI和MAM，成为一种常规的低廉的乳腺成像设备。

1.2 体数据重建过程

Freehand三维超声体数据重建的过程如图2所示，主要包括3个步骤：体数据结构构建、采样像素重新分配和体数据中体素值计算。

图2 自由式三维超声重建过程[19]Fig.2 The process of Freehand 3D ultrasound reconstruction[19]

三维超声重建的第1步，是根据二维超声图像信息，确定重建体数据的尺寸规格，具体包括体数据的坐标原点、维度大小和体素间的物理间隔。文献[20-21]曾分别提出，用图像关键帧或更复杂的主成分分析等方法，确定重建体数据结构的大小。这里采用一种基于包围盒(bounding box)的技术快速确定重建区域的大小，且无需预先确定或限制重建扫描的区域。包围盒只由它的最小点(Xmin,Ymin,Zmin)和最大点(Xmax,Ymax,Zmax)完全确定。

三维超声重建的第二步，是对二维平面上的像素进行重新分配，即遍历二维超声平面上的每一个像素点，并根据它们位置信息的变换矩阵，把像素映射到三维体数据中。如果不止一个像素同时落入到同一个体素中，则需要根据一定的规则，选取合适的值(如平均值、最大值、最先(最后)到达的值等)。在Freehand三维超声系统中，称从二维超声图像到三维体数据的映射为“前向映射”，如图3所示。

图3 前向映射和后向映射[19]Fig.3 Forward mapping and backward mapping[19]

前向映射可表示为

ui=MXi

(1)

式中，Xi表示在超声平面中的像素坐标，M表示前向变换矩阵，ui表示重建体数据中对应的体素坐标。

根据Freehand三维超声扫描坐标系统设置，前向变换矩阵可进一步分解为

(2)

式中，MRP表示从超声图像扫描平面坐标系(P)到超声探头定位器坐标系(R)的变换矩阵，MTR表示从坐标系(R)到跟踪系统自身世界坐标系(T)的变换，MCT表示从坐标系(T)到重建体数据坐标系(C)的变换。

变换矩阵MRP是未知的，但可通过空间校准方法获得，并且一旦校准之后就固定不变。变换矩阵MTR可从实时跟踪系统直接获得，这是一个实时动态变化的矩阵。矩阵MCT的作用是将坐标系统T和坐标系统C对齐，以方便重建后体数据的操作和显示。为便于矩阵的统一操作运算，本研究统一使用4×4的齐次变换矩阵。

前向映射的优点是算法运行效率高。除了前向映射，还可以采用后向映射来完成像素的重分配。后向映射定义为前向映射的逆，可表示为

(3)

后向映射将遍历重建体数据中的每一个体素ui，找到一组离该体素最近的像素集合S={(di,Xi)}(表示像素Xi离当前体素距离为di的集合)，然后通过一定的法则(距离加权、最邻近采样、加权中值等)，对当前体素进行赋值。

除了可以使用前向映射和后向映射方法之外，也可以使用向量运算的方法来完成二维超声平面上的点到三维立体网格中的映射，如图4所示。

图4 用向量方法实现二维数据到三维空间的映射Fig.4 Vector-based method of the mapping from 2D to 3D

在图4中，O(0,0,0)是世界坐标系的原点，C是重建三维体数据的坐标原点，P、P1、P2三点确定一个二维超声图像平面，P为坐标原点,向量PP1、PP2分别为二维平面的两个相互垂直的坐标轴，Xi(x,y)为二维超声平面上的任一点，ui为二维坐标Xi在重建三维体数据中的对应点。首先，要获得向量PXi的坐标，这可以根据Xi在二维平面两坐标轴上的投影来计算，有

(4)

从而，可以得到向量OXi的坐标，有

(5)

然后，根据向量关系可得到ui的坐标，有

(6)

计算出ui的坐标之后，将Xi的像素值赋值给ui，从而完成二维平面的点到三维立体的映射。

由于Freehand三维超声的采样数据是稀疏的，像素分配之后不可避免地会在重建体数据中留下空白区域。因此，重建过程的第3步就是对空白体数据区域进行填补(hole-filling)。针对已知数据进行插值的方法有很多，但基本原理都是利用周围已知的像素值去插值出体素网格中未知的体素值。

2 Freehand三维超声重建算法的分类

目前，用于Freehand三维超声重建的算法有很多，这些算法所完成的工作主要是像素的映射和体素的估算。像素的映射，是指根据空间变换关系，将二维超声图像中的像素值映射到重建三维体数据中的过程；体素的估算，是指在重建的三维体数据中，根据灰度值已知的体素，完成对重建体数据中空体素的计算。根据这些算法的具体实现，三维超声体数据重建方法可以分为3类[22-23]：基于体素的方法(voxelbasedmethod,VBM)、基于像素的方法(pixelbasedmethod,PBM)和基于函数的方法(functionbasedmethod,FBM)。下面将对这些算法进行详细的分析和描述。

2.1VBM

2.1.1VBM变种算法

VBM遍历重建体素网格中的每一个体素，根据后向映射方法找到二维平面图像中与该位置相符的像素(一个或多个)，并把这些像素值按一定的算法赋给相应的体素。根据像素与体素的对应情况，VBM又有许多变种算法。

1)仅用一个像素值来确定一个体素的体素值，这种方法也称为最邻近体素法[24](voxelnearestneighbor,VNN)。VNN将遍历重建体数据中的每一个体素，并计算得到离该体素点最近邻的采样像素点，这主要是通过计算点到面的距离确定，假设超声图像所在的平面为p·n=d，其中n为该平面法向量，p为该平面上的任一点，d为该平面到坐标系统原点的距离，那么重建体数据中任一个体素点q到该平面的距离为d′=q·n-d。根据点到面的距离计算公式，VNN的算法流程可表示为：算法1，VNN重建算法→输入，重建体数据中的每个体素点→输出，体素的灰度值→操作。

具体操作：对每一个体素点q，寻找最邻近的像素点p。首先，计算q到每个二维超声平面图像的法向量an，an=(q·n-d)n，最小的法向量即为最邻近平面；找到最邻近平面之后，寻找距离该法向量最邻近的像素,即为最邻近像素点p,p=q-an。将最邻近像素点p的值赋值给对应的体素点q。

2)由多个像素值确定一个体素的体素值，这类方法称为VBM插值算法(voxelbasedmethodswithinterpolation,VBMI)。VBMI方法将遍历每一个体素点，首先记录距离当前体素点小于某一邻域半径的一组像素点集合，然后根据一定计算法则，利用该集合内的像素灰度值给当前体素点进行赋值，具体流程可表示为：算法2，VBMI重建算法→输入，重建体数据的每一个体素点→输出，体素的灰度值→操作。

具体操作：遍历输入的重建体数据的每一个体素点ui；根据点到面的距离，寻找最邻近二维超声平面上的一组像素集合{Xi}；根据距离加权法则，计算该像素集合的加权值，f(ui)=∑i∈NiωiXi；用计算获得的加权值给对应的体素点赋值。

2.1.2VBMI插值算法

图5给出了VBMI重建方法的原理，其中白色的空心点表示当前待重建的体素点，黑色的实心点表示在采样超声图像平面上的像素点，设定的邻域是半径大小为radius的球体。当前VBMI用到的插值方法主要有距离加权(distanceweighted,DW)插值方法、探头轨迹(probetrajectory,PT)插值方法和中值滤波(medianfilter，MF)插值方法。

图5 VBMI重建原理[19]Fig.5 The reconstruction principle diagram of VBMI[19]

1)DW插值方法[25-27]：距离加权算法通过对局部邻域内的所有像素值的距离加权平均来计算体素值。计算公式如下：

(7)

式中，I(Vc)是当前体素的灰度值，I(Vk)是体素局部邻域内第k个像素的像素值，Wk为第k个像素的权值，等于像素到体素距离dk的倒数。

DW插值算法有多种改进版本[12,28]，包括距离平方加权(squaredistanceweighted,SDW)方法、自适应距离加权(adaptivedistanceweighted,ADW)插值方法和自适应距离平方加权(adaptivesquaredistanceweighted,ASDW)。它们的基本原理与DW算法相同，不同之处在于加权系数Wk的计算。

2)PT插值方法[29]：该方法先计算通过体素中心的虚拟平面，然后找到距离体素最近的两个超声平面(位于体素的两侧)，并记录体素到两个平面的距离，之后由这两个超声平面估计超声探头的运动轨迹，记录探头运动轨迹与这两个超声平面的交点，在每个超声平面上分别找到距离这两个交点最近的4个像素，先对这4个像素分别使用双线性插值，得到每个交点的像素值，最后再对这两个交点用距离加权算法计算当前体素的灰度值。

3)MF插值方法[30-31]：它利用中值滤波器原理，将体素局部邻域内的所有像素灰度值的中值作为该体素的灰度值。常用的MF算法有标准中值滤波(standardmedianfilter，SMF)插值、距离加权中值滤波(distanceweightedmedianfilter，DWMF)插值和自适应加权中值滤波(adaptiveweightedmedianfilter,AWMF)插值。

这里，用{Xi}表示体素点ui局部邻域内像素的集合，di表示像素到体素的距离，ωi表示像素权重，I(ui)表示体素的体素值，I(Xi)表示像素的灰度值。

对于SMF，其原理可表示为

(8)

式中，median{…}表示中值滤波。

而DWMF和AWMF，先在SMF的每一项上乘以对应权值ωi，然后再取其中值。对于DWMF，其权值函数ωi=1/(di+α)2；而对于AWMF，其权值函数ωi=c-σ2di/μ·r。其中，c为正常数，r为空体素的邻域半径，σ为邻域内已知体素值的标准差，μ为邻域内已知体素值的均值。

2.2PBM

PBM方法通常包括两个步骤：像素重新分配(distributionstep，DS)和孔洞填补(hole-fillingStep，HFS)。在DS步骤，将遍历输入二维超声图像上的每一个像素，根据前向映射方法，找到重建的体数据网格中对应的体素，并将该像素按一定的算法分配给对应的体素。在HFS步骤，将遍历每个体素，找到没有赋值的空缺体素，并对其进行插值计算。对DS步骤，最常采用的方法是最邻近像素赋值方法(pixelnearestneighbor,PNN)[32-35]，它将最邻近的像素值赋值给对应的体素。如果多个像素对应于同一个体素，则通常取这些像素值的平均值赋值给对应的体素。通常，由于扫描数据的不规则性和稀疏性，完成DS阶段后，重建的三维标准网格中仍存在空缺的体素，所以需要进行HFS步骤，对空缺体素进行填补。对HFS步骤，根据权重函数的不同可有多种算法，主要包括距离加权算法和快速步进方法(fastmarchingmethod，FMM)[36]。传统的PNN算法在HFS阶段采用距离加权进行计算，其算法流程可表示为：算法3，PNN重建算法→输入，二维超声图像→输出，无空缺的重建体数据→操作。

步骤1：像素重分配阶段(DS)。遍历输入的二维超声图像上的每一个像素点，寻找最邻近的体素点。如果该点体素值为空，则直接赋值voxelvalue=pixelvalue；如果该点体素值不为空，则将多个像素点对该体素的贡献值进行平均，voxelvalue=(voxelvalue+pixelvalue)/2。

步骤2：空体素填充阶段(HFS)。遍历待填充的重建体数据的每一个体素点ui；如果该体素点为空，则以该体素为圆心，确定一个半径为r的圆形邻域，寻找该邻域内体素值不为空的体素集合{Xi}；根据距离加权法则，对邻域内已知体素值进行加权平均，f(ui)=∑i∈NiωiXi；将计算获得的加权均值赋值给对应的空体素点。

对于FMM算法，首先确定待填充的孔洞区域，然后利用快速步进方法遍历孔洞区域的每一个体素进行填补，其基本原理如图6所示，其中，Ω表示重建体素值中的待填充区域，∂Ω为待填充区域的边界。FMM算法的初始位置从位于体素未知区域的初始边界∂Ω0开始，朝着∂Ω的法向量方向移动，并在移动的过程中计算当前边界点上的体素值，如此反复，直至整个未知区域被填充。对处于FMM移动边界上的点p，可采用泰勒展开的方式，对以p为中心、半径大小为ε的球形区域中已知体素进行插值，有

(9)

(10)

2.3FBM

FBM是根据二维超声图像中的像素点与其空间位置建立某种函数关系，再根据该函数计算出重建网格中体素的灰度值。FBM选择一个特殊函数(如多项式) 并确定其系数， 使一个或多个函数曲线通过所输入的像素点。然后，根据这些函数对规则间隔点上的体素进行估值。所选插值函数不同，FBM可以有多种实现，主要包括径向基函数(radialbasisfunction，RBF)[37-38]插值、马尔科夫随机场(Markovrandomfield，MRF)[13,39-40]插值、贝塞尔(Bezier)[14]插值以及基于核回归(kernelregression，KR)[41]方法的插值。

2.3.1RBF插值

该方法的基本思想是利用采样图像上的像素点构建分段的径向基函数，然后利用该函数对重建体数据中的规则网格点进行拟合求值。RBF插值一般分为两步进行：第1步，将体数据分割为小段区域(尺寸不一定要相同，但要确保每个矩形区域有足够的有效采样点)，根据该窗口内的数据拟合该段的RBF函数；第2步，利用拟合出的RBF函数，对重建体数据中的规则网格点进行赋值。假设每个矩形区域内的非零体素数目为N，其对应的体素值为Yi(i=1,2,…,N),它们在体素阵列中的位置为Li=(xi,yi,zi)，则RBF算法将寻找样条函数S(L)，使其尽可能地接近这些数据点，并使此样条函数尽可能平滑。因此，S(L)要满足下列要求：

(11)

式中，等式右边第1项为数据点的样条偏差，第2项是数据平滑项，变量λ决定这两部分的比例关系。S(L)的解为

(12)

式中，T(L)为趋向函数，ai为标量系数，R(L,Li)为径向基函数。

通常T(L)应为一个常数，且T(L)=a0，此时径向基函数为

(13)

样条函数的系数可以通过求解下列线性方程组来获得，有

(14)

(15)

通过N+1个方程求解未知系数ai，由系数ai值和样条函数的表达式，可以得出待插值体素的体素值。

2.3.2MRF插值

该方法运用统计方法，结合贝叶斯框架，估算重建三维立体中的体素值，可表示为

(16)

式中，f表示重建数据的估计值，B(L)=[b1(L),b2(L),…,bN(L)]T是由三维空间基函数bi构成的向量，U=[u1,u2,…,uN]T是一个N×1的系数向量。每一个基函数bi定义为由一个已知函数h位移而得到， h是一个三线性插值函数。

因此，要估算出f，第一步是根据观察数据V={Vi}，Vi=(Yi,Li))对系数U做估算。根据贝叶斯理论，采用最大后验概率方法来估算U，可得出优化问题如下：

(17)

式中，Y为已知三维体数据灰度值Yi的集合，U是未知系数。

假设Y中所有元素都独立同分布，则对似然函数有

(18)

研究发现，超声采样数据服从Rayleigh分布。根据Rayleigh分布的概率密度函数，可以得出似然函数p(Y|U)的表达式为

(19)

然后，假设系数U的先验概率服从Gibbs[42]分布，根据Gibbs与MRF的等价关系[43-44]，可得先验概率p(U)的表达式为

(20)

最后，根据获得的系数U计算三维网格的体素灰度值，完成三维超声体数据重建。

2.3.3Bezier插值

该方法利用Bezier曲线，完成对三维超声体数据的重建。n阶Bezier曲线的表达式为

(21)

从式(21)可以看出，P0,P1,…,Pn一共n+1个控制点决定这条Bezier曲线。其中，Pi为曲线第i+1个控制点，P0为曲线的起始点，Pn为曲线的终点。

通常采用三阶Bezier曲线，即可完成超声三维重建。在实现Bezier插值时，采用一个可以包含4帧二维超声图像的控制窗口来控制生成Bezier曲线的控制点个数。假设P1～P4分别表示控制窗口内的4帧二维超声图像上的像素点，且它们在各自图像中对应的位置相同，将P1～P4转换到三维网格的体元中，它们在三维空间对应的体素为V1～V4，这4个体素在三维空间的坐标分别为VP1(x1,y1,z1)、VP2(x2,y2,z2)、VP3(x3,y3,z3)、VP4(x4,y4,z4)。根据Bezier曲线定义，由4个体素的坐标,可以构造一条基于位置的Bezier曲线，其表达式如下：

(22)

同时，由4个体素的灰度值，可以构造一条基于体素灰度值的Bezier曲线，有

(23)

得到这个三次Bezier曲线的表达式后，就可以插值计算曲线上其他体素的灰度值。假设曲线上某一体素VP(t)的坐标为(x，y，z)，变量t则等于待插值体素V(t)与第一个控制点V1的距离对整条曲线起点V1与终点V4之间距离做归一化所得结果，有

(24)

将t代入曲线公式，就可得出对应的体素灰度值。而对于重叠区域，则采用距离加权来计算体素值。

2.3.4KR插值方法

该方法采用如下的数据来观察模型，有

(25)

式中，r(·)是核回归函数，Li=(Li0,Li1,Li2)是体素的三维坐标，εi是一个均值为0、方差为σ2的高斯噪声，Yi表示超声数据采样值。

假设待求体素点L离已知数据采样点Li很近，根据N阶泰勒展开式，可得：

r(Li)≈r(L)+r′(L)(Li-L)+

β0+β1(Li-L)+β2(Li-L)2+…+βN(Li-L)N

(26)

根据最小二乘法，可得最优问题如下：

(27)

式中，M是邻域窗口内已知体素的个数，KH(·)是核函数，H是平滑参数。核函数KH(·)通常为指数函数或高斯函数[45]，只需满足下列条件

∫tK(t)dt=0 ∫t2K(t)dt=c

(28)

式(27)可用矩阵[46]表述为

(29)

式中

(30)

(31)

(32)

式中，Y是所有体素值组成的向量集，b是所有βi参数组成的向量集，Wl是对角矩阵，对角线上的元素为KH(·)，其他元素为零。

根据最小二乘法，可求得式(29)的解为

(33)

3 Freehand三维超声重建算法分析

对于Freehand三维超声重建算法性能的分析，可以从以下方面考虑：重建图像效果、重建速度和重建精确度等。以算法的重建时间和重建后的均方误差(meansquareerror,MSE)[47]为指标，分析重建算法的性能。

根据文献[23,41]，可以对Freehand三维超声体数据重建算法做出如下分析。

1)从重建精度来看，因为VNN方法可以直接将体素最近的像素灰度值作为体数据值，所以在图像边缘、纹理保持等方面表现较好，但该方法的重建误差比较大，尤其是在采用间隔更大的采样数据集上。ADW算法类似于均值滤波，能够有效地降低超声图像的斑点噪声，但同时也模糊了图像的细节信息和小组织器官的轮廓边缘。PNN算法容易在DS分配区域和HFS区域的边缘产生不连续，会产生图像伪影，影响重建质量。FMM算法因引入图像边缘的结构信息进行空洞区域的填补，所以较好地克服PNN算法存在的不连续问题。MRF算法能对超声图像中的斑点噪声进行滤除，但图像的边缘信息也容易被模糊。与VNN方法类似，KR算法在超声图像的边缘、纹理保持方面表现较好，且重建出的图像清晰。

2)从重建速度来看，文献[41]指出，KR算法在重建精度上最高，但在运行时间上却最慢，时间长达6h之久。比较而言，PBM(如PNN)算法的运行时间是比较少的，一般在十几分钟，VBM(如VNN、ADW)算法的运行时间居中，将近1h，而FBM(如MRF、KR)算法的运行时间最长，长达几小时。因此，现有的Freehand三维超声重建算法的实时性不高，造成这种现象的主要原因是用于重建的二维B超图像序列较多，且生成的三维体数据模型较大。例如，腹部超声体模的B超图像序列有102张，生成的三维体数据大小为131×425×421；人体肝脏的B超图像序列有167张，生成的三维体数据大小为350×338×387。生成的体数据范围越大，算法的计算量就越大，从而运行时间越慢。另外，这些算法的计算过程是由CPU完成的，而CPU多采用串行计算方式，所以速度较慢。

4 未来的研究方向

综上分析可知，在重建图像效果、重建精度和重建速度上，每种重建算法各有利弊。例如：KR算法的重建精度最好，但其重建速度却是最慢的；VNN算法可以较好地保持图像边缘，但其重建误差比较大；PNN算法重建时间居中，但它容易在重建图像上产生伪影，图像精度不高。因此，目前存在的Freehand三维超声体数据重建算法在重建速度和重建精度方面很难兼容，在抑制斑点噪声和保持图像边缘信息两方面也存在相互矛盾的问题。所以，在今后的研究工作中，研究人员应该重点在以下几个方面展开研究和改进：

1)要提供算法的自适应调整能力。目前的算法并不能保证适用于所有的应用场景或图像，或者需要手动调整许多参数。因此，未来的三维超声重建算法应能自动地分析采集超声图像数据，并对不同的场景做出自适应的调整，满足不同场景对高质量重建图像的要求。

2)体数据重建图像质量有待提高。目前，超声图像三维重建结果或多或少地存在抑制斑点噪声和保持图像边缘信息相矛盾的问题。理想的重建算法应该能够有效地抑制斑点噪声，同时尽可能地保持图像中的边缘信息。因此，重建算法的性能还有待研究人员进一步改进。

3)体数据重建算法的复杂度有待降低。在三维超声重建的临床应用领域，特别是超声图像引导的计算机辅助介入手术系统，往往要求三维重建算法具有比较高的实时性。但是，现有的重建算法，尤其是重建图像精度比较好的算法，都普遍存在时空复杂度过高的问题。因此，研究人员要采取有效措施，使重建算法在保证图像质量的同时，让时间和空间复杂度有大幅的降低，如采用图形处理单元(GPU)[48]，提高算法的运行速度。

4)体数据重建算法的客观评测方法需要研究人员建立和完善。目前，研究人员对重建算法的评测还偏重于主观评测，并没有结果一致、得到研究人员普遍认可的针对三维超声体数据重建算法的客观评价体系。因此，建立一个针对体数据重建领域的标准化客观评测体系，将会是以后体数据重建研究的一个方向。

5 结语

到目前为止，研究人员主要从3个方向来研究Freehand三维超声体数据重建这一问题：VBM重建、PBM重建和FBM重建。现有的Freehand三维超声重建算法在重建速度和重建精度上很难兼容，在抑制斑点噪声和保持边缘信息方面也存在着矛盾。因此，探索重建效果好、重建精度高且实时性高的Freehand三维超声重建算法，是未来Freehand三维超声用于计算机辅助介入手术中的重点研究方向之一，值得专业人员深入研究。

[1] 胡伟.自由臂超声三维图像快速重建算法研究 [D]. 广州:华南理工大学, 2013.

[2]VarandasJ,BaptistaP,SantosJ,etal.VOLUS-avisualizationsystemfor3Dultrasounddata[J].Ultrasonics, 2004, 42(1): 689-694.

[3]InglisS,DerekC,JohnPL.Anovelthree-dimensionalendoscopicultrasoundtechniqueforthefreehandexaminationoftheoesophagus[J].UltrasoundinMedicine&Biology, 2011, 37(11): 1779-1790.

[4] 田增民,孙君昭,杜吉祥,等. 3D超声导航系统在神经外科手术中的应用 [J]. 第二军区大学报, 2015, 26(12): 1384-1387.

[5] 徐静,杨向东,朱森强,等. 用于肝癌介入治疗的术中三维超声导航系统 [J]. 中国生物医学工程学报, 2007, 26(5): 719-723.

[6] 孙银山. 基于三维超声图像的穿刺手术机器人辅助系统研究 [D]. 哈尔滨:哈尔滨工业大学, 2011.

[7] 何宇. 经阴道三维超声在妇科疾病诊断中的应用 [D]. 吉林:吉林大学, 2013.

[8]JiangWeiwei,LiAnhua,ZhengYongping.Asemi-automated3Dannotationmethodforbreastultrasoundimaging:systemdevelopmentandfeasibilitystudyonphantoms[J].UltrasoundinMedicine&Biology, 2014, 40(2): 434-446.

[9] 李曦，陈光杰. 三维超声成像原理及其发展 [J]. 中国医学物理学杂志，2001，18(3)：136-140.

[10]FensterA,BaxJ,NeshatH,etal.Advancementsandbreakthroughsinultrasoundimaging: 3Dultrasoundimaginginimage-guidedintervention[M].Rijeka:InTechOpenAccessPublisher, 2013.

[11]MercierL,PetreccaK,AraujoD,etal.Onlinedatabaseofbraintumorimages[EB/OL].http://dx.doi.org/10.1118/1.4709600/ 2012- 05-22/2016-05-11.

[12]HuangQinghua,ZhengYongping.Anadaptivesquared-distanceweightedinterpolationforvolumereconstructionin3Dfreehandultrasound[J].Ultrasonics, 2006, 44(1): 73-77.

[13] 郭强,许健,杨新,等. 基于MAP-MRF的旋转3维超声心动图断层重建及降噪 [J]. 中国图象图形学报, 2005, 10(10): 1281-1288.

[14]HuangQinghua,HuangYanping,HuWei,etal.Bezierinterpolationfor3Dfreehandultrasound[J].IEEETransactionsonHuman-MachineSystems, 2015, 45(3): 385-392.

[15]KuhlCK,SchradingS,LeutnerCC,etal.Mammography,breastultrasound,andmagneticresonanceimagingforsurveillanceofwomenathighfamilialriskforbreastcancer[J].JournalofClinicalOncology, 2005, 23(33): 8469-8476.

[16]WarnerE,PlewesDB,ShumakRS,etal.Comparisonofbreastmagneticresonanceimaging,mammography,andultrasoundforsurveillanceofwomenathighriskforhereditarybreastcancer[J].JournalofClinicalOncology, 2001, 19(15): 3524-3531.

[17]DuricN,LittrupP,LiC,etal.Breastultrasoundtomography:bridgingthegaptoclinicalpractice[C] //BoschJG.MedicalImaging2012:UltrasonicImaging,Tomography,andTherapy.SanDiego:SPIE, 2012, 8320: 83200O1-83200O9.

[18]Delphinusmedicaltechnologies[EB/OL].http://www.delphinusmt.com/2014-09-19/2016-07-20.

[19] 温铁祥. 自由式三维医学超声重建关键技术研究 [D]. 深圳:中国科学院大学, 2014.

[20]GeeAH,PragerRW,TreeceGM,etal.Engineeringafreehand3Dultrasoundsystem[J].PatternRecognitionLetters, 2003, 24(4): 757-777.

[21]RaulSJE,MarcosMF,PabloCM,etal.Atheoreticalframeworktothree-dimensionalultrasoundreconstructionfromirregularlysampleddata[J].UltrasoundinMedicine&Biology, 2005, 29(2): 255-269.

[22]SolbergOV,LindsethF,TorpH,etal.Freehand3Dultrasoundreconstructionalgorithms—areview[J].UltrasoundinMedicine&Biology, 2007, 33(7): 991-1009．

[23] 孙银山,吴冬梅,杜志江,等.Freehand3维超声重建算法综述 [J]. 中国图象图形学报, 2010, 15(11): 1561-1568.

[24]SherebrinS,FensterA,RankinRN,etal.Freehandthreedimensionalultrasound:implementationandapplications[C] //ProceedingsoftheMedicalImaging1996:PhysicsofMedicalImaging.NewportBeach:SPIE, 1996 : 296-303.

[25]BergS,TorpH,MartensD,etal.Dynamicthree-dimensionalfreehandechocardiographyusingrawdigitalultrasounddata[J].UltrasoundinMedicine&Biology, 1999, 25(5): 745-753.

[26]MartensD,GiljaOH.TheEchoPAC-3DSoftwareforImageAnalysis[M].SingaporeCity:WorldScientificPublishing, 2005: 305-336.

[27]TrobaughJW,TrobaughDJ,RichardWD.Three-dimensionalimagingwithstereotacticultrasonography[J].ComputerizedMedicalImagingandGraphics, 1994, 18(5): 315-323.

[28]HuangQinghua,LuMinhua,ZhengYongping,etal.Specklesuppressionandconstrastenhancementinreconstructionoffreehand3Dultrasoundimagesusinganadaptivedistanceweightedmethod[J].AppliedAcoustics, 2009, 70(1): 21-30.

[29]PierrickC,PierreH,XavierM,etal.Probetrajectoryinterpolationfor3Dreconstructionoffreehandultrasound[J].PreprintSubmittedtoMedicalImageAnalysis, 2007, 11(16): 604-615.

[30]HuangQinghua,ZhengYongping.Medianfiltersusedforvolumereconstructioninfreehand3Dultrasound[C]//Proceedingsofthe27thAnnualConferenceoftheEngineeringinMedicineandBiologySociety.Shanghai:IEEE2005 : 1826-1829.

[31]HuangQinghua,ZhengYongping.Volumereconstructionoffreehandthree-dimensionalultrasoundusingmedianfilters[J].Ultrasonics, 2008, 48(3): 182-192.

[32]CarrJC,StallkampJL,FynesMM,etal.Designofaclinicalfreehand3Dultrasoundsystem[C]//ProceedingsintheMedicalImaging2000:UltrasonicImagingandSignalProcessing.SanDiego:SPIE, 2000 : 14-25.

[33]NelsonTR,PretoriusDH.Interactiveacquisitionanalysisandvisualizationofsonographicvolumedata[J].InternationalJournalofImagingSystemsTechnology, 1997, 8(1): 26-37.

[34]SanR,MartínM,AlberolaC,etal.FreehandultrasoundreconstructionbasedonROIpriormodelingandnormalizedconvolution[C] //Proceedingsofthe6thInternationalConferenceonMedicalImageComputingandComputer-AssistedIntervention.Berlin:Springer, 2003 : 382-390.

[35]SanR,MartinM,CaballeroPP,etal.Atheoreticalframeworktothree-dimensionalultrasoundreconstructionfromirregularlysampleddata[J].UltrasoundinMedicine&Biology, 2003, 29(2): 255-269.

[36]WenTiexiang,ZhuQingsong,QinWenjian,etal.AnaccurateandeffectiveFMM-basedapproachforfreehand3Dultrasoundreconstruction[J].BiomedicalSignalProcessingandControl, 2013, 8 (6): 645-656.

[37]RohlingR,GeeAH,BermanL.Radialbasisfunctioninterpolationforfreehand3Dultrasound[J]. 1999, 1613 : 478-483.

[38] 何琪,孙丰荣,耿俊卿,等. 基于径向基函数的超声体数据重建仿真研究 [J]. 系统仿真学报, 2010, 22(8): 1992-1996.

[39]SanchesJM,MarquesJS.Arayleighreconstruction/interpolationalgorithmfor3Dultrasound[J].PatternRecognitionLetters, 2000, 21(10): 917-926.

[40]SanchesJM,MarquesJS.Amultiscalealgorithmforthree-dimensionalfreehandultrasound[J].UltrasoundinMedicine&Biology, 2002, 28(8): 1029-1040.

[41]ChenXiankang,WenTiexiang,LiXingmin,etal.Reconstructionoffreehand3Dultrasoundbasedonkernelregression[J].BioMedicalEngineeringOnLine, 2014, 13: 124-138.

[42]GemanS,GemanD.Stochasticrelaxationgibbsdistributionsandthebayesianrestorationofimages[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 1984, 6(6):721-741.

[43]KindermannR,SnellJL.Markovrandomfieldsandtheirapplications[EB/OL].http://www.ams.org/books/conm/001/conm001-endmatter.pdf/1980/2016-08-19.

[44]LiSZ.Markovrandomfieldmodelingincomputervision[M].Berlin:Springer-Verlag, 1995.

[45]SilvermanBW.Densityestimationforstatisticsanddataanalysis[M].NewYork:Chapman&Hall, 1986.

[46]TakedaH,FarsiuS,MilanfarP.Kernelregressionforimageprocessingandreconstruction[J].IEEETransactionsonImageProcessing, 2007, 16(2): 349-366.

[47]RohlingRN,GeeAH,BermanL.Acomparationoffreehandthree-dimensionalultrasoundreconstructiontechniques[J].MedicalImageAnalysis, 1999, 3(4): 339-359.

[48] 马鸣飞. 基于GPU加速的三维超声成像系统 [D]. 哈尔滨：哈尔滨工业大学， 2011.

The Latest Research Progress of Freehand 3D US Volume Reconstruction

Fan Yangyang1,2Ni Qian3Wen Tiexiang1#*Gu Jia1#Wang Lei1#Xie Yaoqin1#Liu Jingen2

1(ShenzhenInstitutesofAdvancedTechnologyChineseAcademyofSciences,Shenzhen518055,Guangdong,China)2(SchoolofInformationEngineering,WuhanUniversityofTechnology,Wuhan430070,China)3(ShenzhenFutianHospitalofChineseMedicine(ShenzhenHospitalofGuangzhouUniversityofChineseMedicine),Shenzhen518034,Guangdong,China)

Freehand 3D ultrasound imaging system is more in line with doctor′s habits and the surrounding of operating room for its simple, useful, wider field of view and higher image resolution, which makes it one of the research focuses in the field of guidance of interventional surgery based on ultrasound imaging. Freehand 3D ultrasound produces 3D volume data of anatomical objects from a sequence of irregularly located 2D B-mode ultrasound images, and renders the reconstructed 3D volume data. Volume reconstruction is one of the key procedures in freehand 3D ultrasound imaging system and plays an important role in improving the reconstructed image quality. This paper summed up three kinds of volume reconstruction methods for freehand 3D ultrasound image, incuding voxel-based method, pixel-based method and function-based method. After that, some algorithms are presented with their merits and drawbacks. Finally, the research progress is summarized and future research directions are suggested.

freehand 3D ultrasound; volume data reconstruction; reconstruction algorithms

10.3969/j.issn.0258-8021. 2017. 01.012

2016-05-10， 录用日期:2016-07-29

国家自然科学基金(61401451)；中国博士后科学基金(2016M590827)；广东省创新团队项目(2011S013)

R318

A

0258-8021(2017) 01-0092-11

# 中国生物医学工程学会高级会员(Senior member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)， E-mail: tx.wen@siat.ac.cn