CN101569543A - 弹性成像的二维位移估计方法 - Google Patents

Abstract

一种弹性成像的二维位移估计方法，包括步骤：对待处理的超声图像的第一射频帧与第二射频帧，用具有鲁棒性的匹配算法进行二维初始位移估计；用得到的二维初始位移估计对第一射频帧图像进行调整，使得相比于原始的第一帧，调整后的第一帧和原始的第二帧之间的相关性更高；将原始的第二射频帧与调整后的第一射频帧，通过亚像素位移估计算法进行位移估计；将所得的二维初始位移估计结果与亚像素位移估计算法所得位移估计结果进行合成处理，得到总体的二维位移估计。本发明的弹性成像的二维位移估计方法，采用粗略但鲁棒的初始估计步骤以及一个准确但需要图像相关性较高的精细搜索步骤，同时达到高速和准确。

Description

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

技术领域

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

背景技术

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

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

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

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

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

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

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

发明内容

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

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

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

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

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

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

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

附图说明

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

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

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

具体实施方式

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

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

如图1所示，本发明实施例的二维位移估计方法的包括下列步骤：

S1：对两连续超声图像的第一射频帧与第二射频帧，用具有鲁棒性的匹配算法进行二维初始位移估计；

S2：用步骤S1得到的二维初始位移估计对第一射频帧图像进行拉伸调整；

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

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

优选的，还可以包括：

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

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

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

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

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

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

$(u,v)=\arg \min \left(\underset{i=1}{\overset{b1}{\mathrm{\Σ}}}\underset{j=1}{\overset{b2}{\mathrm{\Σ}}}\right|B_1(i+u,j+v)-B_2(i,j)\left|\right),B_1\⋐f_n,B_2\⋐f_{n+1}---\left(1\right)$

其中b1和b2为块维度，通常b1＝b2。B1和B2为两射频帧中相应的两块。当B1为当前帧时，B2遍历搜索其2P*2P(P为搜索范围)的邻近区域从而进行图像匹配性能绝对误差和(SAD)的评估。

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

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

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

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

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

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

f(i+u*δt，j+v*δt，t+δt)＝f(i，j，t)                      (2)

其中，t较小同时(u，v)是像素(i，j)的水平和速度向量。当t足够小，计算等式(2)左端的泰勒估计：

$f(i+u*\mathrm{\δt},j+v*\mathrm{\δt},t+\mathrm{\δt})=f(i,j,t)+\frac{\∂f}{\∂i}\mathrm{\δi}+\frac{\∂f}{\∂j}\mathrm{\δj}+\frac{\∂f}{\∂t}\mathrm{\δt}+H---\left(3\right)$

忽略高阶量H并合并等式(2)和(3)，得到：

$0=\frac{\∂f}{\∂i}\mathrm{\δi}+\frac{\∂f}{\∂j}\mathrm{\δj}+\frac{\∂f}{\∂t}\mathrm{\δt}---\left(4\right)$

然后可得：

$\frac{0}{\mathrm{\δt}}=\frac{\∂f}{\∂i}\frac{\mathrm{\δi}}{\mathrm{\δt}}+\frac{\∂f}{\∂j}\frac{\mathrm{\δj}}{\mathrm{\δt}}+\frac{\∂f}{\∂t}\frac{\mathrm{\δt}}{\mathrm{\δt}}=\frac{\∂f}{\∂i}u+\frac{\∂f}{\∂j}v+\frac{\∂f}{\∂t}---\left(5\right)$

其中 $\frac{\∂f}{\∂i}=f_x,$ $\frac{\∂f}{\∂j}=f_y$ 和 $\frac{\∂f}{\∂t}=f_t$ 代表水平、垂直和时间方向上的导数。

将等式(5)写成矩阵形式

$\left[\begin{array}{cc}{f}_x& f_y\end{array}\right]\×\left[\begin{array}{c}u\\ v\end{array}\right]=-f_t---\left(6\right)$

为了进一步限定等式(6)，可做多个假设从而产生多个版本的光流算法。假设光流为本地常量，例如，一个3*3的块：

$\left[\begin{array}{ccc}p11& p12& p13\\ p21& p22& p23\\ p31& p32& p33\end{array}\right]$

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

$\left[\begin{array}{cc}{f}_x\left(p11\right)& f_y\left(p11\right)\\ f_x\left(p12\right)& f_y\left(p12\right)\\ \·\·\·& \·\·\·\\ f_x\left(p33\right)& f_y\left(p33\right)\end{array}\right]\×\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}{f}_t\left(p11\right)\\ f_t\left(p12\right)\\ \·\·\·\\ f_t\left(p33\right)\end{array}\right]---\left(7\right)$

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

$\left[\begin{array}{c}u\\ v\end{array}\right]=d=\arg \min {|\mathrm{Ad}-b|}^2---\left(8\right)$

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

A^TA＝A^Tb                          (9)

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

d＝inv(A^TA)*A^Tb                   (10)

此处A^T和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+Δt时段在指定点(i，j)的离散化处理。

$\frac{u_{i,j}^{t+\mathrm{\Δt}}-u_{i,j}^t}{\mathrm{\Δt}}=$

$\frac{1}{\mathrm{\Δx}}(c_{E_{i,j}}^t{\▿}_E{u}_{i,j}^t-c_{W_{i,j}}^t{\▿}_W{u}_{i,j}^t)+\frac{1}{\mathrm{\Δy}}(c_{N_{i,j}}^t{\▿}_N{u}_{i,j}^t-c_{S_{i,j}}^t{\▿}_S{u}_{i,j}^t)---\left(11\right)$

下标N、S、E和W分别表示像素周围的北南东西，符号

表示在其下标所示方向上最邻近梯度，其他系数由下式给出

$c_{N_{i,j}}^t=g_k\left({\left|{\▿}_N{u}_{i,j}^t\right|}^2\right),$ $c_{S_{i,j}}^t=g_k\left({\left|{\▿}_S{u}_{i,j}^t\right|}^2\right)$

$c_{E_{i,j}}^t=g_k\left({\left|{\▿}_E{u}_{i,j}^t\right|}^2\right),$ $c_{W_{i,j}}^t=g_k\left({\left|{\▿}_W{u}_{i,j}^t\right|}^2\right)$

其中选择g_k为 $g_k\left(\mathrm{\ξ}\right)=\frac{1}{1+\frac{\mathrm{\ξ}}{k}}$

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

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

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

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

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

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

Claims (10)

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

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

S2：用步骤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数字信号处理芯片或电脑的图形显示卡实现。

Images