基于快速匹配算法的颈动脉超声图像的运动分析

1、基于快速匹配算法的颈动脉超声图像的运动分析                     作者：徐明才，陶进绪连娟，李昕，吴学咏，汪华【摘要】  颈动脉超声图像中的运动信息能够间接地反应颈动脉弹性等状况，结合颈动脉内中膜厚度(CIMT)能够为心脑血管疾病诊断提供定性与定量的依据。我们将改进后的金字塔快速匹配算法(modified block sum pyramid，MBSP)应用于颈动脉超声波

2、图像斑点跟踪获得运动信息。理论和实验结果都表明，改进后的金字塔块匹配算法能有效地减少运动跟踪的运算量，并且有着和改进前的金字塔块匹配算法相同的准确度。运动跟踪结果能够为医生诊断心脑血管疾病起到一定的辅助作用。 【关键词】  颈动脉;斑点跟踪;金字塔算法;运动分析;平均绝对差函数Abstract：The motion information in the carotid ultrasound image can reflect the condition of carotid whether it has hardened, expanded and contracted in the

3、 right way. Then it can help doctor to diagnose cardio cerebral vascular diseases together with the carotid intima-medial thickness(CIMT).We used the modified block sum pyramid (MBSP) algorithm to analyze the motion of carotid ultrasound images. The result shows that the modified algorithm effective

4、ly reduced the computation while speckle tracking. Besides that, the MBSP algorithm has the same accuracy with the block sum pyramid. The tracking result can assist the doctor to evaluate the disk of cardio cerebral vascular diseases to a certain degree.Key words：Carotid; Speckle tracking; Block sum

5、 pyramid; Motion analysis; Mean absolute difference1 引 言心脑血管疾病与颈动脉内中膜厚度及其运动特征等有着密切的相关性1。对于超声波图像而言，可以通过斑点跟踪的方法获得运动信息。斑点作为噪声不可避免地存在于超声波图像中。虽然斑点噪声影响超声波图像质量，但是根据斑点噪声的形成原理，可以通过跟踪斑点来获取图像中相关目标的运动信息。为研究组织的运动特征和弹性信息，进而辅助评估心脑血管疾病提供了一种可能。斑点跟踪算法一般包括两个方面：匹配和搜索。匹配的准则有许多种，常用的准则有归一化相关函数(NCCF),平均均方误差函数(MSD)以及平均绝对差函数

6、(MAD)等。因MAD不涉及乘法等复杂运算，所以较为常用，本研究也采用该匹配准则。搜索方法也有好多种，有全搜索法2，十字搜索法3，菱形搜索法4等。此外，在这些算法基础上出了许多快速算法，如金字塔算法5,连续淘汰法6,改进的自适应2BIT变换法7等。临床实践中医生除了用肉眼定性地去观察颈动脉的运动特征，有时需要定量测量出感兴趣的颈动脉某一目标的运动信息(如：短轴切面轮廓在一个心动周期内收缩和扩张的运动状况)。这样就能够更加可靠地评估心脑血管疾病和风险。由于辅助诊断对运动跟踪的实时性要求不是很高，而对运动跟踪的准确性要求较高。故我们选择能达到全局最优的全搜索法和金字塔快速匹配算法的组合进行斑点运动

7、跟踪，并且对金字塔快速匹配算法进行了改进，将其应用于颈动脉运动分析，高效准确地跟踪出目标兴趣点运动的相关信息。2 原理与方法2.1 金字塔块匹配算法设有块X和Y,块大小均为：2M×2M，对块X和Y均建立 M+1层的金字塔。从塔顶往下依次为第0，1，M层。X,Y的金字塔的第m 层分别记为Xm,Ym。第m-1层和第m层关系如方程(1)所示：Xm-1(i,j)=Xm(2i-1,2j-1)+Xm(2i-1，2j)+Xm(2i,2j-1)+Xm(2i,2j)(1)令：MADm(X,Y)=2mi=12mj=1|Xm(i,j)-Ym(i,j)|(2)可以得出结论(3)：MAD0(X，Y)MAD1(

8、X,Y)MAD2(X,Y)MADm(X,Y)MADm(X,Y)(3)两个块在生成了各自的金字塔之后，从顶层往下逐层MADm(X,Y)。如果发现其大于或者等于最新的MADmin时，说明本候选块一定不是最优块，舍弃该候选块继续搜索下一个候选块，如果比较到最底层,仍然比最新的MADmin还小时，将MADmin更新为当前候选块与块的MADm(X,Y)。直到所有的候选块均搜索完后，最小的MADmin的候选块相对于参考块的偏移量即为所要求解的运动向量。 金字塔匹配算法就是通过尽快地舍弃那些不是最优的候选块而减少运算量的。虽然建立金字塔会增加额外的计算量，但是由于大部分的块在比较时就被舍弃了，所以总的运算量

9、还是大幅地减少。 2.2 改进的金字塔匹配算法金字塔块匹配算法从顶层开始往下逐层比较时，层与层间的MAD值相差很大。当层数较大时，最后一层与倒数第二层的MAD值相差更大，因此，比较时仍有不少的计算冗余。为此，希望在现在的金字塔基础上不改变运算量，在层与层之间新增一些阈值，能够尽早地舍弃那些不是最优的块。为了叙述方便，以下讨论均以在最后一层(第M层)与倒数第二层(第M-1层)之间增加阈值为例。两个要比较的块分别记为X,Y,块X和Y的金字塔的第M-1层分别记为：XM-1，YM-1。设SM-1为第M-1层塔的所有像素点空间。将XM-1，YM-1按照相同的方式均分割为N个区域：S1M-1,S2M-1，

10、S3M-1，SkM-1，,SNM-1，1kN，见图1。具体的分割方法还需要根据实际情况进行合理的选择。分割出的N个区域要满足以下两个约束关系：Nk=1SM-1k=SM-1(4)SM-1kSM-1l=，1klN(5)图1 像素空间分割方法Fig 1 Pixels partition method图2 实验所采取的分割方法Fig 2 Experiments method同时对块X，Y的金字塔的第M层也划分为N个区域SM1,SM2,SM3,SMk，SMN，1kN。也要满足约束关系：Nk=1SMk=SM和SMkSMl=，1klN。SM-1k区域中的每一个点对应于SMk区域中的4个点，此对应关系和构建金

11、字塔时从M层到M-1层的4个点累加成1个点的对应关系一致。构造如下函数表达式：(0nN)MADM-1n(X，Y)=MADM-1(X，Y)-nk=1(i,j)SM-1k|XM-1(i,j)YM-1(i,j)|+nk=1(i,j)SMk |XM(i,j)YM(i,j)|(6)根据层间区域关系式(1)，上式又可以写成(7)式：MADM-1n(X，Y)=MADM-1(X，Y)-f(n)(7)其中f(n)可写成(8)式:f(n)=nk=1(i,j)SM-1k|XM-1(i,j)YM-1(i,j)|-1         |XM

12、(2i-1,2j-1)YM(2i-1,2j-1)|+|XM(2i-1,2j)YM(2i-1,2j)|+|XM(2i,2j-1)YM(2i,2j-1)|+|XM(2i,2j)YM(2i,2j)|(8)又据三角不等式及上下层点的对应关系知：f(n)0，且f(n)随着n的增大单调递减,据此得出(9)、(10)、(11)的3个结论：MADM-1(X，Y)=MADM-10(X,Y)(9)MADM-1N(X，Y)=MADM(X,Y)(10)MADM-10(X，Y)MADM-11(X，Y)MADM-1n(X，Y)MADM-1N(X，Y)(11)如果将SM-1像素空间一列一列依次从左往右分割时，所得到的有关结

13、论则和8是一致的。文献8中所提到的逐行更新算法可看作是本研究改进算法的特例。现将改进的金字塔匹配算法(以在两帧中跟踪一个兴趣点为例)归纳如下：记前后两帧图像分别为fi(x,y),fi+1(x,y)，被跟踪点坐标为(a,b)，前一帧图像点(a,b)所对应的块称为块，搜索范围为A。(1)相对于参考块运动位移为(0,0)的候选块与参考块两者的MAD值作为MADmin，并且确定后一帧图像中的起始候选块。MAD计算公式如(12)式所示：MADmin=2M-1-1x=-2M-1 2M-1-1y=-2M-1|fi(a+x,b+y)-fi+1(a+x,b+y)|(12)(2)依据金字塔算法构建当前候选块和参考

14、块各自的金字塔。金字塔总共假设有M+1层。沿金字塔至上而下逐层进行比较，因为这里只考虑在M-1层和M层间新增阈值，所以前M-1层的比较和原金字塔块匹配算法一样。当M-1层比较完后，发现MADM-1(X，Y)仍然小于最新的MADmin。此时则需要按照构造函数表达式生成一个新的中间阈值，再进行比较，其间发现新生成的阈值大于或者等于MADmin时，舍弃，直接进入步骤(3)进行下一个候选块的匹配。否则一直进行到最后一层为止。如果其结果仍然小于MADmin时，更新MADmin。(3)按照选取的搜索算法(如全搜索法)，确认下一点被候选点及相应的候选块。重复步骤(2)，直到A所有需要搜索的点搜索完为止。(4

15、)最后，对应于MADmin值最小的候选块即为最优候选块，其相对于参考块的偏移量，即为所要跟踪的运动矢量。2.3 利用MBSP算法分析颈动脉图像(1)取颈动脉短轴切面一个心动周期的图像序列。在第一帧图像中，沿颈动脉短轴轮廓，手动大致均匀地选取出若干个兴趣点，兴趣点应尽量选取在整个运动序列图像中亮度较为平稳且肉眼观察较为明显的斑点噪声颗粒。(2)利用改进的金字塔跟踪算法(MBSP)逐帧地跟踪所有兴趣点，得到各个兴趣点在各帧中的运动矢量，再将各帧上各点的运动矢量按照适合人眼观看的比例显示在各帧图像上。最后以动态的方式播放附有运动矢量的颈动脉图像序列，医生可以根据播放视频中的运动矢量定性地掌握兴趣点运

16、动大小、方向、扭转、是否出现异常等情况，做出一些如是否可能出现斑块硬化等初步诊断。(3)最后给出各个兴趣点的收缩扩张曲线，每条曲线给出每个点随时间的运动速度变化。设第j帧图像中跟踪到的N个点的分别为坐标(X1，Y1),(Xi,Yi),(XN,YN)。其均值为(X,Y)。记Vi=(Xi,Yi)-(X,Y),13 实验与结果3.1 实验结果实验数据由合肥105超声科提供，为一个心动周期的颈动脉短轴切片序列图像，共25帧。第1、11、22帧图像分别见图3、4、5。图3中小白框即为我们要跟踪的6个兴趣点，左上角为1号点，按顺时针依次排列。为了表述方便，这里只选取了6个兴趣点为例。分别用BSP和MBSP

17、算法对该颈动脉序列图像进行斑点运动跟踪，其中MBSP算法在第3层和第4层增加了3个新的阈值。将第3层的像素空间从左到右依次分割成4个相同的区域，见图2。实验中块大小取为16×16，搜索范围为8×8。图4附上了第11帧图像上各兴趣点的运动矢量图，两个算法跟踪出来的结果一样。图6给出了各个兴趣点(从左往右，从上至下依次为1到6)随时间变化的扩张收缩运动速度图。从图中我们可以看出，兴趣点2、3、4的运动表现出较为一致的运动，在一个心动周期内经历一个明显的扩张和之后的收缩过程。但是观察兴趣点1、5、6的运动速度曲线，发现其运动起伏相对来说没有这些规律，从而可以认为1、5、6附近区域

18、可能出现一些病变。这和医生的诊断结果：本试验者颈动脉左下部区域内膜形成了较为明显的斑块，而斑块区域收缩运动常出现紊乱是相互吻合的。3.2 性能分析(1)准确性。分析改进后的算法和改进前的金字塔匹配算法(BSP)，容易得出前后算法运动跟踪结果应该是一致的，差异只在于计算量上。实际跟踪的结果也证明了这一点。图7、8给出了前后两种算法的跟踪兴趣点1的运动轨迹图，图中横轴、纵轴分别对应于图像的X轴和Y轴坐标轴。(2)效率。表1给出了改进前后算法在上述实验过程中需要计算的各MAD值统计次数。计算第m层的MADm值所需要运算量(含加减法及绝对值)可以公式化9为：3.22m-1,计算MADmn运算量约为：5

19、n·(3·22m+2-1)/16。以第4层为例，如果直接采用表1 金字塔算法改进前后各MAD次数原来的金字塔算法，需要计算总次数为821 457;但采用改进后的算法需要的计算量仅为721 747。可见因为第3层与第4层3个新增阈值的引入，第4层判决过程的总体运算就减少12.11%，可以考虑将区域进行更细划分，但是要权衡因为阈值增加而造成的比较运算量。如果在第3层等层也引入同样机制的话，运算量将进一步减少。4 本研究将改进后的金字塔匹配算法应用于斑点运动跟踪，实验结果表明该算法能准确而更加有效地跟踪出运动矢量。并且将该方法应用于颈动脉运动分析，结合颈动脉的自身特点，给出有利于辅助医生诊断的扩张与收缩曲线等依据，为颈动脉的进一步分析奠定了一定基础。【参考文献】  1胡玮,张军,孔国喜，等. 超声检测颈内动脉内-中膜厚度及其形态结构与相关心脑血管疾病的关系J. 第四军医大学学报 ,2007,28(15):1431-1433.2Brünig M, Niehsen W. Fast full-search block matchingJ. IEEE Trans Circuits Syst Video