WO2009132516A1

WO2009132516A1 - 弹性成像的二维位移估计方法

Info

Publication number: WO2009132516A1
Application number: PCT/CN2009/000454
Authority: WO
Inventors: 郑永平; 周永进; 黄铮铭
Original assignee: 香港理工大学
Priority date: 2008-04-29
Filing date: 2009-04-28
Publication date: 2009-11-05
Also published as: CN101569543A; CN101569543B

Description

弹性成像的二维位移估计方法技术领域

本发明涉及一种超声弹性成像技术，尤其与弹性成像中的二维位移估计方法有关。背景技术

弹性成像是通过比较来自在某种应力作用下的介质的图像从而获取组织或材料（介质）弹性比对的技术。所述图像可通过任意已知模式获取，包括但不限于超声波、光学、 X光以及 MRI图像等等。这些图像可在表面沿着纵向或横向的截面或通过投影获取。同时，应力作用也可以有多种形式存在，包括准静态压缩、拉伸和剪切、多向震动、渗透压等等。这些应力作用可以来自外部源或内部源，后者例如通过心跳或呼吸导致的组织移动。在以下的描述中，将釆用在准静态压缩下获取的超声波图像，但是本发明可以应用于其它各种类型的图像以及其它类型的应力作用方式。其中超声波图像可以是视频图像或射频图像。以下的描述以超声波射频图像为例。

超声弹性成像的核心在于利用在压缩前后得到的超声射频（RF)信号来估计超声图像中的局部位移。该估计可为一维或二维。在超声弹性成像中，通常采用超声波换能器来轻微压缩感兴趣区域从而在使得沿着压缩方向的不同材料之间产生形变。形变可以发生在轴向和横向。通过估计位移场方向图，可进一步计算沿压缩方向的形变分布。假设压缩是弹性的，形变场可指示材料弹性参数在沿压缩方向的分布。

在现有的弹性成像技术中，主要存在下列的四个不足。

第一，弹性成像对比度与和压缩前后之间的信号连续性之间的折中问题；通常，为保持压缩前后信号高度相关（即信号连续性）并控制噪音水平，形变不可过大，但是，形变过小可导致所生成的图像欠缺对比度；其中，该对比度指示在所感兴趣的区域内，不同硬度组织的在形变图像中的亮度区别。因此，由于需要更加优化的不同组织部分的视觉效果，最好得到较高的对比度，此为弹性成像中待解决的主要问题。

第二，较之于脱机操作，实际操作中的实时性要求更高。诸多报告显示在超声场中压缩前后对象位移的实时估计很有必要性的却其实现具有挑战性是，其中有题目为 "A modified block matching method for real-time freehand strain imaging"的作者为 Y. Zhu T. Hall 的发表于 Ultrasonic Imaging, 24： 161-176, 2002年的报告，以及题目为 "Clinical assessment of real-time freehand elasticity imaging system based on the combined autocorrelation method "β勺作者为 T. Shiina, Μ. Yamakawa, Ν. Nitta, Ε. Ueno, Τ. Matsumara, S. Tamano和 T. Mital e的报告，涉及了这个问题。但上述文献中釆用较复杂的算法来达到高准确性，这与实现所需要的实时处理之间总是存在矛盾的，因此在弹性成像领域发明具有高速处理能力且同时具有高准确性的算法是非常必要的。

第三、压缩所诱发的组织位移的模型简化问题。传统上，仅有轴向位移可细化估计，横向位移只是通过简单内插得到。因此，形变的横向组成部分会通过去相关噪音扰乱轴向位移估计。对于实时应用的场合，用复杂的位移模型并基于二维超声来完全求解不同方向上的形变是不现实的，因此必须进行简化假设。这些假设因而产生出不同的模型以及进一步不同的位移估计方法。对此问题的解决方案之一是 E. Konofagou和 J. Ophir的美国专利 6270459B1 , "横向位移，横向形变和形变比率的成像和估计方法" 。但是，其只将其散斑跟踪方法局限在互相关算法范畴，而互相关算法的计算是非常耗时的。

第四、从位移数据中获取优化的形变图像的问题。在弹性成像中，无论内部 /外部或震动 /准静态组织压缩的情况下，轴向形变都被视为压缩时产生的轴向分布的梯度。但是，梯度操作自身将在位移数据中会放大噪音。虽然可采用滤波技术能够减轻这种问题（F. Kallel和 J. Ophir 1997; W. Luo, J. Bai, P. He 2004) 。但大体 ±说，使用现有方法不可能同时达到兼顾平滑性和边界清晰度的要求。发明内容

针对现有技术中存在的不足，本发明目的在于提供一种弹性成像的二维位移估计方法。

一种弹性成像的二维位移估计方法，包括步骤-

S 1 :对待处理的两帧图像的第一帧与第二帧，用具有鲁棒性的匹配算法进行二维初始位移估计；

S2: 用步骤 S1得到的二维初始位移估计对第一帧图像进行调整，使得相比于原始的第一帧，调整后的第一帧和原始的第二帧之间的相关性更高；

S3：将原始的第二帧与调整后的第一帧，通过亚像素位移估计算法进行位移估计；

S4:将步骤 S1所得的二维初始位移估计结果与步骤 S3所得位移估计结果进行合成处理，得到总体的二维位移估计。

本发明优点在于，本发明的二维位移估计方法，其中包括粗略但可靠的初始估计步骤以及准确但需要数据相关性较高的细化检测步骤，从而通过使用获取自经过机械扰动前后的软材料的两幅超声图像生成并显示二维弹性对比图像。本发明在组织的弹性成像中快速进行鲁棒二维位移估计，克服了现有技术中存在的诸多不足，很好地考虑了弹性成像对比度与和压缩前后之间的信号连续性之间的折中；且在一定程度上符合实际操作中的实时性要求；简化了模型；而且能够从位移数据中获取优化的形变图像。附图说明

所包括的附图提供了对本发明的进一步理解，其被并入到本说明书中并构成为本说明书的一部分，所述附图示出了本发明的实施例并与说明书一起用于解释本发明的原理。在附图中相同的附图标记表示相同的部件。在附图中：

图 1是本发明实施例的弹性成像的二维位移估计方法流程图；

图 2A-图 2K为本发明实施例的弹性成像的二维位移估计方法处理一组射频帧的过程及其对应的结果；

图 3A-图 3K为本发明实施例的弹性成像的二维位移估计方法处理一组在较大压缩下获取自胸腔弹性图像的射频帧的典型结果。

具体实施方式在以下的描述中，将采用在准静态压缩下获取的超声波图像，但是本发明可以应用于其它各种类型的图像以及其它类型的扰动。

本发明尤其涉及采用超声射频信号的位移估计新方法以及用于在该位移估计之后获取形变的新方法。以两连续的射频帧（图像）为例，首先采用具有鲁棒性的匹配算法计算粗略位移分布，如块匹配算法（Block Matching Algorithm, 简称 BMA) ，然后根据位移估计得到的结果来对其中一帧（假设选用第一射频帧）进行调整生成新的帧，该调整例如通常为拉伸。由于调整后的第一射频帧较之原始的第一射频帧与第二射频帧具有更佳的相关性，可使用密集光流方法，如 Lucas-Kanade光流方法（LKOF)对调整后的第一射频帧和原始的第二射频帧进行计算得到轴向和横向上的具有亚像素精度的位移估计。

如图 1所示，本发明实施例的二维位移估计方法的包括下列步骤- S1 :对两连续超声图像的第一射频帧与第二射频帧，用具有鲁棒性的匹配算法进行二维初始位移估计； S2: 用步骤 SI得到的二维初始位移估计对第一射频帧图像进行拉伸调整； S3：将原始的第二射频帧与调整后的第一射频帧，通过亚像素位移估计算法进行二维位移估计；

S4:将步骤 S1所得位移估计结果与步骤 S3所得的位移估计结果进行合成处理，得到较二维初始位移估计精确的总体二维位移估计。

优选的，还可以包括：

步骤 S5 : 重复进行下列过程直至符合设定条件：用步骤 S3得到的亚像素位移估计对步骤 S2得到的调整后的第一射频帧图像进行进一步调整，例如拉伸调整；并用第二射频帧与进一步调整后的第一射频帧通过亚像素位移估计算法进行二维位移估计；再将所得到的亚像素二维位移估计累加到步骤 S4得到的二维位移估计。

步骤 S5之后还可以进行步骤 S6: 对步骤 S4或步骤 S5得到的图像，通过各向异性扩散方法进行微分处理。

本发明使用的 BMA，首先由 Jain引入图像处理中。该技术将图像分为多个矩形块，针对各个块，在另一图像（例如射频帧）中检索最佳估计，从而估计位移向量。由于各个块足够小，所以旋转和缩放可通过待匹配的两图像中的相应区域分段转换来进行大致估计。简单说，各种版本的 BMA, 包括：（1 )可能使用不同匹配标准，例如最小化方差和，均方误差和绝对平均误差，或最大化互相关函数以及弹性极限差分分类；（2)较之于穷尽搜索，使用不同的检索策略虽然损失准确性但能够减少计算量，例如使用三步搜索、菱形搜索、二维对数算法、正交搜索和分级搜索等等；（3 )使用拉普拉斯图像金字塔或小波变换进行多解析度块匹配；其他变换包括浮点精度块匹配，块面积变化的块匹配和可变形块匹配等等。

该块匹配首次由 Levinson等人采用来提取组织弹性信息。 Zhu等人提供了块匹配方法的修正版本，包括适应性搜索窗调节和实时建立弹性成像的新质量控制标准。

为了保证块匹配结果，本发明的二维位移估计方法，可使用亚像素精度遍历搜索 BMA。遍历搜索 BMA也是块匹配最为直接的方法。

具体的说，对于两帧射频帧 fn和 fn+l , 如果采用最小化的绝对误差作为匹配标准，则二维位移向量为：

(u , ν ) ( 1 ) 其中

当 B1为当前帧时， B2遍历搜索其 2P*2P (P为搜索范围）的邻近区域从而进行图像匹配性能绝对误差和（SAD) 的评估。

为了达到更高的初始位移估计，还可以用 BMA对进行了像素插值的两帧图像进行计算，从而得到亚像素的位移分辨率。比如在每两个像素点之间插一个新的像素点，就可以得到半像素的位移分辨率。

根据 BMA移动估计结果调整后的第一射频帧与原始的第二射频帧图像具有更佳的相关性，因此保证了后续可用 LK-OF进行亚像素移动估计的效果。这样，在实际应用中就可以在保证准确的基础上使用大于 1%的形变来获得相对高的形变对比度。本发明中，其他块匹配方法也可在步骤 S1中使用。

歩骤 S2中，用步骤 S1得到的二维初始位移估计对第一射频帧图像进行拉伸调整，使得相比于原始的第一帧，调整后的第一帧和原始的第二帧之间的相关性更高；

然后进行步骤 S3: 将原始的第二射频帧与经步骤 S2调整后的第一射频帧，通过亚像素位移估计算法进行位移估计。

密集光流方法可以实现以亚像素精度计算位移场，但是却缺乏鲁棒性或需要图像相关性较髙。本发明步骤 S1使用 BMA完成了具有鲁棒性的初始二维位移估计后，可以进行具有鲁棒性的亚像素精度的位移场计算。

众所周知，光流方法已成功用于血管弹性成像，例如【Wanetal.2001 and Maurice et al. 2005】， Pellot-Barakat等人在取自类似组织的图像的射频数据上以及活体胸腔组织中测试了光流方法，该结果优于取自基于相关技术的比较对象。首次提出于 1981 年的 Lucas-Kanade方法仍是最为广为应用的双帧微分光流方法。光流通常指相对于在图像平面上图像强度的可见移动的瞬间移动场。如果将 fn写为 f(i,j，t)，其中，函数 (i，j)为像素位置， t为时间参数，在下一射频帧中，光流的强度常数假设为：

f(i + u *St,j + v^t,t + St) = f(i, j,t) (2) 其中， t较小同时 (ιι,ν)是像素 (i,j)的水平和速度向量。当 t足够小，计算等式（2)左端的泰勒估计- f(i + u*dt,j + v^St,t + St) = f(i,j,t) +—Si + ^-Sj + ^-dt + H (3) di dj dt

忽略高阶量 H并合并等式（2)和（3) ，得到- f _s. df

0 = ^~δί + St 4 di

然后可得：

df St ― df

+ df

St 一 di dt dj St dt dt d dj d 其中 ' dj- - - /,和 d t ■ 代表水平、垂直和时间方向上的导数

将等

为了

假设光流为本地常量，例如，一个 3*3的块：

ρΠ p\2 pl3

ρ2\ ρ22 ρ23

ρ3 \ ρ32 ρ33

从等式 (6)中得到 9个等式。

在等式（7) 中的三个向量或矩阵通常分别记为 A， d和 b，由此

最后的平方问题引出 Lucas-Kanade等式

A^TA = A^Tb (9) 当 A^TA 得到充分调节时，也就是说，矩阵 ^²¹^ 的特征值 λ ΐ或 λ 2足够大，等式（9) 就有个简单解。

d = inv ( A ^T A) * A ^Tb (10) 此处 A^r和 inv(*)分别代表 A的转置矩阵以及矩阵求逆函数。

美国专利 6277074使用类似粗到细的策略作为其发明，但是并未解决上述所有弊病。在其细化步骤中使用互相关算法导致该方法效率低于本发明。本发明可实现具有更多运算能力的实时弹性成像。

本发明的二维位移估计方法，同时估计水平形变以及轴向形变；同时通过在 BMA和 LK-OF模型中调整窗口参数，可结合局部信号连续性来获得更准确的二维位移估计，同时无须在射频 A线之间插值。

当压缩前后的图像去相关太过严重时，可重复进行 LK-OF和图像调整过程，来在所述框架中达到更佳结果。这意味着 LK-OF算法的结果首先将用于对第一射频帧进行调整从而可在第二轮使用再次调整过的第一射频帧和原始的第二射频帧进行 LK-OF计算得到更加准确的结果。本发明步骤 S5的过程，可反复进行直至得到设定的条件达到时，例如达到 LK-OF的最终优化结果时。将来自 BMA和 LK-OF的二维位移估计结果合成之后，就得到了通过第一射频帧和第二射频帧估计的总体位移估计图。合成算法可以是简单的加和或选用加权加和算法，比如可以给 BMA的结果较大的权重，又比如给重复迭代的前几次较大的权重。现有技术的文献中有许多不同的方法来从位移估计图得到代表弹性的形变图。本发明步骤 S6，采用各向异性扩散方法来平滑进行微分处理并保持主要边界清晰。扩散方法不仅可以起到平滑作用，同时保持图像中的边缘，该边缘对探测图像中不同物体的边界。可采用不同的扩散方法达到此目的。因此，本发明解决了如何平滑形变图像同时保持组织边界。

对于扩散方法，例如，经典 Perona-Malik扩散方法的实质在于在 t到 t+ A t时段在指定点（i, j ) 的离散化处理。

u i ,j ^u ij

A t

(11)

1 1

Ax '·' '^{J i ,J} Ay

下标 N、 S、 E和 W分别表示像素周围的北南东西，符号 V表示在其下标所示方向上最邻近梯度，其他系数由下式给出

2 2

c_E' _iJ = gk i ^ s )， ^c = S_k ( ^ w ^u ' )

其中选择 gk为 ) - i

1 + ί- k

在迭代进行过程中，高梯度项权重将大大高于低梯度项的输出权重，因此通过调整参数 k，可保留强边界区域同时平滑其他部分。

需要说明的是，各算法均可由 FPGA芯片硬件实现，可由 DSP数字信号处理芯片实现，也可由电脑的图形显示卡实现。

下面结合图 2A-图 2K和图 3A-图 3K，显示两组使用本发明二维位移估计方法获得的典型超声图像结果。图 2Α-图 2Κ显示一个具有很小变形的情况，而图 2Α-图 2Κ显示一个具有较大变形的图像组。在两组图像中，图 2Α-图 2Κ和图 3Α-图 3Κ分别表示如下的内容。图 2Α和图 3Α为第一射频帧，图 2Β和图 3Β为第二射频帧，图 2C和图 3C为采用 ΒΜΑ 半像素精度产生的垂直位移场估计，图 2D和图 3D为采用 ΒΜΑ半像素精度产生的水平位移场估计（全黑色表示几乎没有水平位移），图 2Ε和图 3Ε为根据位移估计通过将第一射频帧向第二射频帧调整得到的新 "帧"，图 2F和图 3F为由调整后的第一射频帧和原始的第二射频帧釆用 LK-OF估计得到的垂直位移，图 2G和图 3G由调整后的第一射频帧和原始的第二射频帧采用 LK-OF估计得到的水平位移，图 2H和图 3H最终的垂直位移场，图 21和图 31最终的水平位移场，图 2J和图 3J为垂直形变图，图 2K和图 3K为二十次各向异性扩散后的垂直形变图。

如说明书所述，本发明公开的方法即可实现高速计算的同时保持计算准确性。此外，可同时在垂直和水平方向得到位移和形变图像。

本发明采用粗略但鲁棒的初始估计步骤以及一个准确但需要图像相关性较高的的精细搜索步骤。这个创新举措使该方法同时达到高速和准确。

以上，是为了本领域技术人员理解本发明，而对本发明所进行的详细描述，但可以想到，在不脱离本发明的权利要求所涵盖的范围内还可以做出其它的变化和修改，这些变化和修改均在本发明的保护范围内。

Claims

权利要求

1、一种弹性成像的二维位移估计方法，包括步骤-

S1 : 对待处理的两帧图像的第一帧与第二帧，用具有鲁棒性的匹配算法进行二维初始位移估计；

S3：将原始的第二帧与调整后的第一帧，通过亚像素位移估计算法进行位移估计； S4: 将步骤 S1所得的二维初始位移估计结果与步骤 S3所得位移估计结果进行合成处理，得到总体的二维位移估计。

2、如权利要求 1所述的位移估计方法，其特征在于，步骤 S1中，所述二维初始位移估计是对像素插值后的图像进行从而得到亚像素分辨率的二维初始位移估计。

3、如权利要求 1所述的位移估计方法，其特征在于，步骤 S2中，所述调整为图像拉伸调整或是基于平滑后的二维初始位移估计。

4、如权利要求 1所述的位移估计方法，其特征在于，步骤 S4中，所述合成处理是加和算法或加权加和算法。

5、如权利要求 1所述的位移估计方法，其特征在于，所述图像是超声波射频图像、超声波视频图像、光学图像、光学层析图、 X光图像或核磁共振图像。

6、如权利要求 1-5所述的位移估计方法，其特征在于，还包括步骤 S5，重复进行下列过程直至符合设定条件：用步骤 S3得到的亚像素位移估计对步骤 S2得到的调整后的第一射频帧图像进行进一步调整;并用第二射频帧与进一步调整后的第一射频帧通过亚像素位移估计算法进行位移估计；再将所得到的亚像素位移估计累加到步骤 S4得到的位移估计。

7、如权利要求 1-5所述的位移估计方法，其特征在于，还包括：

步骤 S6: 通过各向异性扩散方法进行微分处理。

8、如权利要求 1所述的位移估计方法，其特征在于，所述鲁棒性的匹配算法为 BMA 块匹配算法。

9、如权利要求 1 所述的位移估计方法，其特征在于，所述亚像素位移估计算法为 Lucas-Kanade光流方法。

10、如权利要求 1所述的位移估计方法，其特征在于，所述各算法是由 FPGA芯片硬件、 DSP数字信号处理芯片或电脑的图形显示卡实现。