基于SIFT特征关键点匹配的脑组织提取方法

江少锋 黄志鹏 杨素华 陈 震 张聪炫

(南昌航空大学测试与光电工程学院，南昌 330063)

引言

脑组织提取是指从脑部序列图像中将脑组织与颅骨、眼球、皮肤、脂肪等非脑组织分离出来的过程，是脑部磁共振图像处理的重要处理步骤之一。脑组织提取的研究可以简单地将其分为五大类：基于区域的、基于边界的、基于图谱的、基于学习的和基于混合的脑组织提取方法。

基于区域的脑组织提取方法采用灰度阈值分割或者分类的方法，将脑组织图像分成几个不同的区域，并且认为具有相同或相似特性的区域是相同的组织，然后采用区域合并或者形态学操作的方法对脑组织进行提取。基于区域的提取方法有聚类算法[1-2]、区域生长法[3]、图割(graph-cuts)法[4]、形态学方法[5-6]等。基于边界的脑组织提取方法通过找出一条封闭的曲线，将脑部划分为内部区域(脑组织)和外部区域(非脑组织)。典型的方法有BSE(brain surface extraction)算法[7]、可变点阵方法BET(brain extraction tool)算法[8]、活动轮廓模型方法[9-11]。基于图谱的脑组织提取方法先采用图像配准技术，将脑图像与标准图谱图像进行对齐，然后再进行脑组织的提取[12-13]。基于学习的脑组织提取方法，首先对训练样本提取灰度、空间坐标和梯度等信息作为特征，然后通过机器学习确定参数模型，进而得到实验测试样本的参数[14]。基于混合的脑组织提取方法试图将上述几种提取方法结合起来，从而得到更好的提取结果。混合的提取方法通常是一个由粗到精的过程。Segonne等结合了混合分水岭算法和可变点阵算法[15]，Jiang等结合了BET算法和混合水平集算法[16]，Galdames等结合了直方图阈值分割法和可变点阵法[17]，税午阳等结合了分水岭算法和水平集方法[18]，贾迪等结合了水平集算法和形态学方法[19]，Wang等结合了图谱法和机器学习方法[20]。

本研究将BET算法和图谱匹配相结合，提出了一个新的混合算法、实现脑组织稳定和精确提取。该方法是一个由粗到精的算法，在粗提取阶段通过特征匹配解决BET算法中的参数设置问题，从而得到稳定的脑组织边界；在精提取阶段以BET得到的脑组织边界作为初始轮廓进行graph cuts分割，由于初始轮廓非常接近脑组织，故解决了graph cuts方法对初始化敏感问题，实现了脑组织精确提取。

1 方法

1.1 BET简介

BET演化过程实质上就是通过3个力(平滑力u1、u2，扩张力u3)不断更新离散点阵位置，将初始轮廓上的顶点推到脑组织边界的过程。在图1中，A0、A1、A2分别为初始轮廓中3个连续顶点，O为初始轮廓的圆心，R是初始轮廓半径，Ac为A1和A2的中点，S为中线向量，Sn是垂直于A1和A2连线的垂线向量，At为垂点，St是Ac、At连线向量，由图中向量关系可知St=S-Sn，Sn的单位向量是eSn。

图1 BET演化力定义Fig.1 The forces defined in BET

本研究中3个力定义如下：

1)平滑力u1(见式1)。该力主要用于保持各顶点的等分，将A0拉向A1和A2中间，从而使得3个相邻的顶点等距离分布，并且该力平行于轮廓切线。

2)平滑力u2(见式2)。该力主要用于保持轮廓的平滑，作用于顶点A0，与u1垂直，当A0过于凸出时，u2会将其往内拉，这样u2在顶点演化过程中起到了平滑作用。θ是相邻顶点间的夹角，用θ表示轮廓的平滑度，即当θ<160°时，认为此处轮廓光滑度不够，这时简单地将θ推到180°；当θ≥160°时，则认为此处轮廓线足够光滑，无需进行处理。

3)轮廓演化力u3(见式3)。该力与u1垂直，主要在顶点演化过程中控制顶点扩张或收缩，将顶点推向脑组织边界。其中，Imin、Imax是以A0为起点，沿R方向寻找d1、d1/2距离内所有像素点的最小和最大灰度值，本研究d1取30。参数bt是0～1的一个常数，本研究取0.1；t2是一个灰度阈值，与脑组织边界部分低信号平均值有关。式(4)中的t4是用于调整灰度值的偏差对演化的影响，t4为脑白质和脑灰质的灰度平均值的粗略估算。这些参数中t2对BET提取影响最大，过大容易导致局部收敛，过小容易导致边界泄漏。

u1=St

(1)

t1=0.5(1+tanh(8(Imax-t4-t2)))bt+t2

(4)

这3个作用力的合力表达式为

u=α1u1+α2u2+α2lu3

(5)

1.2 模板的构建

本研究采用IBSR网站提供的标准MRI脑部图像作为模板(见图2(a))。直接对IBSR中MRI图像的手动脑组织提取结果(见图2(b))做BET处理，得到脑组织边缘的N个离散的点(见图2(c))，再将这N个离散点作为关键点，利用SIFT算法[21]对这些特征点提取描述子特征并进行存储，这样就集成了一个含有离散点阵位置、特征描述子等信息的模板特征。

图2 模板边界提取。(a)模板图像；(b)手动提取的脑组织；(c)BET提取的脑组织边界点Fig.2 The brain contour of template image.(a)Template image;(b)Brain tissue mask extracted by manual;(c)The brain contour obtained by BET

提取描述子特征方法如图3所示，左边是以关键点为中心的16×16的邻域窗口，将其划分为4×4的小窗口，计算每个小窗口的梯度方向直方图，并将梯度方向划分为8个方向，每个方向的大小是由该窗口内的像素点在该方向的梯度值经过高斯加权后叠加得到的，最终形成如图3右侧中4×4×8=128维的特征向量描述子。

图3 关键点的描述子Fig.3 The feature description of key point in SIFT

1.3 基于特征匹配的脑组织提取方法步骤

BET方法对参数t2较敏感，且提取结果常出现过平滑。为解决t2参数设置问题，本研究提出了基于特征匹配的脑组织提取方法，通过特征匹配找到已经位于脑组织边界的点，通过这些点自适应确定参数t2；为解决过平滑问题，在得到BET的结果后进行graph cuts精提取，从而得到精确的边界。

本研究所提出的基于特征匹配的脑组织提取方法本质上是一个由粗到精的过程，在读入脑部磁共振图像后，采用默认参数(t2=0.15，该参数较大，保证不会发生边界泄漏)进行BET脑组织迭代提取(见图4(a))，每次迭代后，以得到的轮廓点阵作为关键点提取SIFT特征，并和模板图像进行特征匹配，对于匹配上的点(图4(b)中的蓝色标记点)继续迭代演化，对于没有匹配上的点(图4(b)中的红色标记点)说明关键参数t2不合适，导致演化过慢，为此采用已匹配上的点去估计未匹配点的t2参数，其估计方法如式(6)所示。更新之后继续演化，直到达到最大迭代次数为止。

(6)

式中，Nm为已匹配的点的个数，Ii,min为第i个匹配点对应的Imin。

基于特征匹配的BET算法得到的脑组织轮廓已非常贴近脑组织(见图4(c))，但是存在过平滑现象。为得到更为精确的脑组织边界，本研究以该边界所含区域为初始值进行膨胀减腐蚀操作(形态学参数为菱形结构，半径为5)，定义一个环形感兴趣区域(ROI)，在该ROI中采用文献[4]的方法进行graph cuts精提取，并在式(7)～(11)中重新定义了graph cuts中边的权重定义，最终得到如图4(d)所示的脑组织边界。

图4 本研究算法脑组织提取过程。(a)BET提取中间结果；(b)匹配结果；(c)根据匹配结果继续BET演化结果：(d)graph cuts 精提取结果Fig.4 The processing steps of this brain extraction method. (a) The result of BET in the middle step; (b) The result of feature matching between template and target image; (c) The final result of BET; (d) The refined result of graph cuts

可以看出，本研究通过模板匹配，实现了BET演化过程中的重要参数t2的自动选择和脑组织的精确提取。该方法计算流程如图5所示。

图5 本研究算法流程Fig.5 The flowchart of this method

式中,Dpe为p点到脑组织边界的最近距离，M为图像的尺寸。Dt(p)使得位于脑组织边界附近的像素点权重值较小，而其余像素点权重值较大，这样对于远离边界的点，能够准确地被分割成非脑组织。对于边界附近的点须进一步通过IG(p)和IR(p)来区分。

(9)

式中：IG(p)是和ROI区域内的像素点灰度值相关的值；Ip代表像素点p的灰度值，μ代表ROI中图像的灰度平均值，Im代表ROI中图像的灰度最大值；k1为参数，可根据不同图像进行调整。IG(p)表明，Ip与μ差别越大，权重越高，一般脑组织边界上的点的灰度接近μ其权重低，非脑组织灰度远离μ其权重高。

式(10)、(11)是类似BET半全局搜索的约束项，表明若某点附近有很亮的点(如眼球位置)，则认为该点也不属于脑组织，给该点分配一个高权重，其中tm为脑组织均值估计，I(i)为沿中心方向距离该点为i的点的灰度。

1.4 特征匹配及匹配度计算

为尽可能消除误匹配和保持不变性，本研究提出一种基于空间距离弱约束投票的特征匹配计算方法，其表达式为

(12)

式中：Ti表示目标图像中脑组织边界中第i个顶点，Aj表示模板图像中脑组织边界中第j个顶点；FTi表示提取的Ti点的SIFT特征，FAj表示提取的Aj点的SIFT特征；Df表示两个顶点的特征距离，Dd表示两个顶点间的空间距离。VOTE表示对特征和空间距离先排序然后进行投票，只有特征和空间距离在前nA次投票均命中的才有效，但如果有多个匹配点，则取特征距离最近的那个点作为匹配点，而不是取空间距离最近的点作为匹配点，这样相当于采用距离作为一个弱约束，不仅能够在很大程度上消除误匹配，同时对图像的缩放和旋转具有较好的不变性。

一个简单采用基于空间距离弱约束投票特征匹配方法的二维脑组织边界匹配结果如图6所示，左图为脑组织模板图像，右图为目标图像。分别提取两幅图像脑组织边界的SIFT特征，然后采用基于空间距离弱约束投票的特征匹配方法进行特征匹配，可以看出多数边界点都匹配上，说明该方法对图像的缩放和旋转具有较好的不变性。同时可以看到，匹配上的点存在少量误匹配，本研究选取了19组数据进行匹配实验，这些数据在图像大小、灰度分布上存在较大差异，得到的误匹配率为10.1%。

图6 特征匹配结果(左为模板图像，右为目标图像，红色标记为未匹配点，蓝色标记为已匹配点)。(a)模板图像缩小0.8倍匹配结果；(b)模板图像旋转10°匹配结果Fig.6 The result of feature matching between templ-ate image (left) and target image(right)(Red points are unmatched points, blue points are matched). (a) The result of 0.8 down scaled template image; (b) The result of 10 degree rotated template image

1.5 实验方法

为了验证本方法的提取效果，采用哈佛大学医学院形态分析中心提供的IBSR测试数据中的20组标准T1加权磁共振图像作为实验数据。每组数据包含有约64个冠状切片，每片为256像素×256像素，切片厚度为3 mm。IBSR为每一幅脑组织图像都提供了相对应的标准的人工分割结果。采用本研究方法、BET算法[8]、ROBEX算法[14]和GCUT算法[4]等4种方法,对这20套数据进行脑组织提取对比实验。

1.6 评价方法

采用Jaccard 精度系数(JS)、Dice 精度系数(DS)、假阳性率(false positives rate，FP_Rate)、假阴性率 (false negatives rate, FN_Rate)、敏感性 (sensitivity, SE)和特异性 (specificity, SP) 进行评价比较，有

式中，IM代表提取结果中所有目标区域的像素点集合，IN代表人工分割标准的目标区像素点集合，I为元素都为1的图像集合。

2 结果

表1中列出了4种方法对20组实验数据进行提取，得到各个评价参数的平均值、标准差及分布范围，记录格式为平均值±标准差[范围]。从表1可以看出，6个评价参数中，本研究方法除FN_Rate和SE之外的参数都为最优，其中最大DS和JS的平均值表明精度最高，最小的方差表明数据分布紧密，算法对不同图像稳定性最好。FN_Rate最低的为GCUT方法，但其FP_Rate却为27.3%。通常来说，FN_Rate和FP_Rate越小，提取结果越准确。但只要提取结果包括所有大脑组织，FN_Rate就为0，而此时该结果中也可能包含大量的非脑组织。同时，如果提取结果总是包含在脑组织内，FP_Rate也等于0。所以说，一个算法仅仅具有较小的FN_Rate或FP_Rate并不一定代表其结果非常准确。一个准确、稳定的算法应该权衡考虑其FN_Rate或FP_Rate，本研究方法的FN_Rate单独较小不是最优，但综合考虑两者均低于5%，具有最好的均衡性。参数SE和SP也有相同的情况，SE和SP均衡性最好的是ROBEX算法，本研究方法次之。

由于本研究在定义ROI时，采用了形态学方法，而形态学元素的结构和大小对ROI有较大影响，为此本研究采用了3种形态学元素(圆形、菱形、正方形)，形态半径为5，对IBSR20组数据进行实验，得到平均DS系数分别为0.961、0.962、0.962，说明本方法对形态学元素形状不敏感。再采用菱形结构元素，形态半径分别为3、5、7各自进行试验，得到的平均DS分别为0.927、0.962、0.914。可见，形态学元素大小对提取结果有较大影响，形态学元素过大，使得ROI范围过大，容易产生边界泄漏；形态学元素过小，ROI范围小，容易陷入局部收敛；而中等大小的形态学参数，既能保证不产生边界泄漏，又不会陷入局部收敛。

图7定性列出了本方法、BET 算法、ROBEX算法和GCUT算法对某幅脑部MRI图像切片的处理结果。可以明显看出，本算法最接近手动提取结果，BET算法同时有边界泄露和局部收敛现象存在，GCUT算法处理结果中出现较多边界泄露，ROBEX算法处理结果中存在少量边界泄露且边界过于平滑。

表120组IBSR脑部MRI图像处理结果评价参数比较(均值±标准差)

Tab.1Comparisonofthismethodwithexistingbrainextractionapproaches,BET,GCUT,andROBEXusing20MRIvolumesfromIBSR(mean±SD)

方法DSJSSPSEFP_Rate%FN_Rate%BET0.849±0.0760.745±0.1100.976±0.0080.910±0.11222.87±7.959.0±11.3[0.666～0.928][0.499～0.866][0.963～0.991][0.698～0.999][11.8～39.8][0.1～30.2]GCUT0.880±0.0150 0.786±0.02400.973±0.0081.000±0.00027.34 ±3.910.01±0.02[0.844～0.907][0.729～0.829][0.957～0.985][0.999～1][20.5～37.1] [0.0～0.06]ROBEX0.940±0.0120.888±0.02100.987±0.0050.993±0.00511.9±2.750.67±0.46[0.912～0.954][0.839～0.911][0.976～0.993][0.981～0.993][8.31～17.00][0.21～1.92]本方法0.962±0.0080.926±0.01400.994±0.0040.972±0.0204.95±2.742.82±2.00[0.944～0.971][0.894～0.943][0.984～0.998][0.921～0.996][2.19～11.40][0.43～7.91]

图7 不同算法脑组织提取结果。(a)脑部MRI图像切片；(b)手动提取结果；(c) BET算法提取结果；(d) ROBEX算法提取结果；(e) GCUT算法提取结果；(f)本方法提取结果Fig.7 The extraction results with BET, ROBEX,GCUT and this method. (a) The original cerebral MRI slice; (b) The brain extracted by manual, (c) Result with BET; (d) Result with ROBEX; (e) Result with GCUT; (f) Result with this method

3 讨论

脑组织提取是一个二元分割问题，但是在脑部MRI图像中，脑组织和其他相邻组织的灰度缺乏足够的差异性，不同组织灰度又相互交叠，且常存在灰度分布不均、伪影现象，导致采用自动算法稳定精确提取脑组织难度很大[22]。为此，本研究尝试采用模板匹配和混合提取方法，以提高自动提取方法的稳定性和精度。

本研究首先在BET方法的初始化过程中采用了模板匹配方法，较好地控制了BET方法边界溢出(最低的FP_Rate)。本方法本质上是改进的BET方法，提高了BET方法的稳定性，但是过平滑现象仍然存在(见图4(c))，故需进一步采用后续处理，以得到精确的脑组织边界。本研究在后处理中采用的改进的Graph cuts方法，由于前处理提供了很好的初始边界，故很少出现较大的提取误差，保证了较高的稳定性。但是，算法中形态学大小参数对提取精度有一定影响，不能取太大或太小，经过大量实验，形态学半径取5较好。

对比实验中，本研究选取的两种提取对比算法(GCUT 和 ROBEX)和本研究方法均采用了graph cuts方法进行后续精提取，不同之处在于初始化方式和Graph cuts中的边权重的设置。其中,初始化主要影响graph cuts方法的稳定性，边权重影响graph cuts方法的精度。由于本方法初始化过程中采用了模板匹配方法，为后续精处理过程中的graph cuts方法提供了很好的初始边界，故很少出现较大的提取误差，保证了较高的稳定性。相比之下GCUT方法采用简单的阈值分割方法为graph cuts提供初始化，初始轮廓不够稳定，导致提取的结果存在较大的边界溢出(最高的FP_Rate)；而ROBEX方法采用学习方法提供的初始轮廓较为稳定，但是后续graph cuts方法没有针对脑组织特性进行优化，故提取结果存在少量边界泄露(较高的FP_Rate)，且边界过于平滑。相比之下，本方法在精提取的过程中，本方法根据脑组织提取特点重新定义了graph cuts的边权重，从而有着较高的准确性(最高的DS和JS)。可以看出本方法由于在初始化(稳定性)和边权重方面(精确性)均对graph cuts模型进行了优化，故能够保证较高的准确性和精度。对比实验中，本方法对所有处理图像在采用相同的参数下得到了最高的精度和最小的方差，经验证表明本方法的精确性和稳定性。当用户在用本方法提取脑组织时，不用考虑参数设定问题，大大方便了用户。

FN_Rate和FP_Rate两个参数具有一定的矛盾性，低的FN_Rate可能对应高的FP_Rate，反之亦然，为此需要根据实际应用来选择以哪个参数为主。比如，在fMRI图像配准应用中，不允许脑组织被误排除，就要求尽可能低的FN_Rate; 在脑组织测量应用中，则要求FN_Rate和FP_Rate都不能太高。

4 结论

本研究针对提取稳定性，提出了模板匹配BET方法。为graph cuts精提取提供稳定的初始轮廓；针对提取精度；改进了graph cuts精提取中的边权重定义，从而实现了精确稳定的脑组织提取。虽然本自动算法取得了较好的效果，但是由于在脑组织提取迭代过程中需多次进行匹配计算，导致运算量较大，故后续将研究采用并行方法以提高计算速度；同时进一步研究特征提取和特征匹配技术，进一步降低误配率，提高提取精度。