CN103815932B - 基于光流法和应变的超声准静态弹性成像方法 - Google Patents

Abstract

本发明提供一种基于光流法和应变计算的实时超声准静态弹性成像方法，首先采集经过波束合成的B型图像数据像素信息；接着对获得的B型图像数据像素信息进行多项式展开，使用多项式去逼近每个像素的邻域；然后通过对每个像素进行位移估计，得到整帧图像的位移场，尤其是轴向位移场；最后对获得轴向位移场进行卷积计算，得到B型图像数据的轴向应变场；再经过对应变场进行降噪声、可视化、彩色处理得到彩色弹性图像。本发明直接对B型图像数据进行计算，算法效率高，能做到实时计算，且鲁棒性较高。

Description

基于光流法和应变的超声准静态弹性成像方法

技术领域

本发明涉及一种一种超声成像方法，具体是基于光流法和应变的超声准静态弹性成像方法。

背景技术

超声波回波成像技术目前已经被广泛应用于军事、医疗等领域。超声波回波成像中关于弹性成像的概念最早由Ophir等人于1991年首次提出。之后，弹性成像技术在近二十年里得到了迅猛发展，它被称为继A型、B型、D型、M型超声之后的E型模式。超声弹性成像是通过超声成像系统进行生物组织弹性参数成像，超声弹性图能够提供传统B超图像所无法反映的生物组织弹性特征，对于肿瘤检测等临床应用有非常大的帮助。由于具有易于实现、适用于实时诊断和对组织无侵害性等优点，超声弹性成像技术受到了广泛的关注。

至今，弹性成像技术已有多种实现方式，但它们的基本原理都是测量生物组织对于某个激励(外部或内部)所做出的力学反映，一般都是通过测量生物组织的位移分布，分析组织位移场来重建组织的弹性参数，将生物组织的弹性信息(包括应变、杨氏模量、泊松比等)进行成像显示，以作为临床诊断的参照。

超声准静态弹性成像是指对待检测生物组织施加一个静态/准静态的外力，通常是用超声探头缓速挤压待测生物组织，并同时记录组织在受到挤压前后的超声信号，通过比对两帧超声信号，计算出组织的运动场，进而进行弹性参数的重建。Ophir最早提出的弹性成像概念即为超声准静态弹性成像，由于实现简单并易于实时成像，因此超声准静态弹性成像一直以来都是弹性成像领域的研发热点。

现有的弹性成像方法，计算位移场的方法大多基于互相关算法和改进的互相关算法，其使用的数据为RF射频数据，而不能直接使用B模式图像，由于处理RF射频数据需要很大的运算，因此造成了现有的方法计算效率低下的问题；光流法是指空间运动物体在观察成像平面上的像素运动的瞬时速度，利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。计算应变的方法大多基于最小二乘拟合算法和梯度法，使用这些方法计算出的应变，在精度方面还有提高的空间。

为了在得到准确的位移场的同时，提升算法的效率、减小运算开销，本发明提出了一种基于多项式展开的稠密光流算法。该基本原理是：第一步，使用多项式去逼近前后帧图像每个像素的邻域，这一步可以使用多项式展开转换高效地实现。第二步，通过观察一个确切的多项式如何经过转换，一个从多项式展开系数去确定位移场的方法，经过一系列的发明改进后得到一个鲁棒性好的实时算法。

发明内容

本发明的目的是提供一种基于光流法和应变的超声准静态弹性成像方法，克服了传统超声弹性成像过程中计算效率较低、不能实时成像、鲁棒性较差的问题。本发明采用的技术方案是：

一种基于光流法和应变计算的实时超声准静态弹性成像方法，包括以下步骤：

步骤1：获取经过波束合成的B型图像数据像素信息；

步骤2：对经过步骤1处理获得的B型图像数据像素信息进行多项式展开，使用多项式去逼近每个像素的邻域；

步骤3：通过对步骤2得到的每个像素的邻域进行位移估计，得到整帧图像的位移场；

步骤4：对步骤3获得的整帧图像的位移场进行卷积计算，得到B型图像数据整帧图像的应变场；

步骤5：对经过步骤4得到的应变场信息进行取对数法降噪声处理；

步骤6：对步骤5降噪后的应变场信息进行可视化、彩色处理得到彩色弹性图像。

本发明与已有技术相比具有的优点是：提供了一种基于光流法和应变的超声准静态弹性成像方法，位移场、应变计算过程简单，不但克服了现有技术计算过程效率低的问题；且能进行实时成像，具有较高的鲁棒性。

附图说明

图1为本发明超声准静态弹性成像的工作流程图。

图2为本发明的实施例滤波核梯度算子的示意图。

具体实施方式

下面结合具体附图和实施例对本发明作进一步说明。

为让本发明的上述特征和优点能更明显易懂，下文特举实施例，并配合附图作详细说明如下。

步骤1：获取经过波束合成的B型图像数据像素信息，通过整帧图像的每个像素的位置信息进行采集。

步骤2：对经过步骤1处理获得的B型图像数据像素信息进行多项式展开，使用多项式去逼近每个像素的邻域，得到整帧图像的每个像素的邻域信息。我们可以用多项式来对信号的轴向位置进行计算：f₁(x)＝x^TA₁x+b₁ ^Tx+c₁，其中x代表某个维度的位置信息，A₁是一个对称矩阵，b₁是一个向量，c₁是一个标量，T是矩阵的转置符号。

步骤3：通过对步骤2得到的每个像素的邻域进行位移估计，得到整帧图像的移位场，例如，第t帧的时候某一点的轴向位置是x₁，在t+1帧的时候某一点的轴向位置我们假设为x₂，则该点的运动为：d(x)＝x₂-x₁，这里我们仅关注轴向位移，忽略其他方向的位移，即其他方向的位移量为0。

具体的我们可以通过如下多项式进行计算：

f₂(x)＝f₁(x-d)＝(x-d)^TA₁(x-d)+b₁ ^T(x-d)+c₁＝x^TA₂x+b₂ ^Tx+c₂，

其中：A2＝A1，

b₂＝b₁-2A₁d，

c₂＝d^TA₁d-b₁ ^Td+c₁，d为步骤2中B型图像像素某个点在相邻图像中的轴向位移，其由公式d(x)＝(∑wA^TA)^-1∑wA^TΔb，推导而来，其中w(Δx)是某一点邻域范围内各像素点的权重函数，且Δb(x)＝-0.5(b₂(x)-b₁(x))。A＝A₁＝A₂，其为一个矩阵函数，简称矩阵。

当然我们也可以人为选取相邻间隔较大位移的两张整帧图像进行位移场的计算，这是需要人工手动筛选或自动设置一个预先阀值进行判断。在自动筛选设置一个预先阀值判断时，进行如下操作，取出中代表轴向变化平移量的Δx与预设的位移阈值较，若满足以下公式：

其中是预先设置的一个阀值。

这时我们通过如下公式来计算整帧图像的轴向位移，公式如下：

Δb(x)＝-0.5(b₂(x₃)-b₁(x₁))+A₁(x₁)Δ(x),其中x₃＝x₁+Δ(x)。

步骤4：设计一个滤波核对位移场进行卷积，用于计算位移场的偏导，得到整帧图像的应变场。如图2所示，得到应变场该梯度算子大小为N行*1列(N为正整数)。其计算方法如下：

$y(M+k)=\frac{25*(3M^4+6M^3-3M+1)*k-35*(3M^2+3M-1)*k^3}{(2M+3)*(2M+1)*(2M-1)*(M+2)*(M+1)*M*(M-1)}$

$y(M-k)=-\frac{25*(3M^4+6M^3-3M+1)*k-35*(3M^2+3M-1)*k^3}{(2M+3)*(2M+1)*(2M-1)*(M+2)*(M+1)*M*(M-1)}$

y(M)＝0，其中，y代表信号的轴向位移，N＝2*M+1，k为大于0且不大于M的整数(1≤k≤M)，经过与传统弹性成像对比这样处理的最大特点在于滤除造的同时可以确保信号的形状、宽度不变。

该滤波核是一个N行*1列的梯度算子；该滤波核是一个一维滤波核，N是滤波核的直径，M指代滤波核的半径，M为大于1的整数(在数字图像处理中，滤波核直径常等于滤波核半径乘以二再加一)；y(M-k)、y(M)、y(M+k)分别为在滤波核中第M-K行第一列、第M行第一列、第M+K行第一列的元素值。

步骤5：对经过步骤4得到的整帧图像的应变场信息进行降噪声处理。首先求得整帧图像的位移场偏导的绝对值dVR；再对绝对值使用取对数的方法(log(1+dVR))去除噪声；然后将其归一化到0～255区间的整数，得到降噪后的整帧图像的应变场。

步骤6：对步骤5去噪后的应变场信息进行可视化、彩色处理得到彩色弹性图像。例如，将由应变场构成的灰度图像转换成彩色图像，再将下系统设置的颜色系统投射到应变场，形成彩色的弹性图像。

对比以往弹性成像方法中求取位移场采用互相关法及只能针对RF射频数据的缺陷，通过上述本本发明实施例，可以直接通过对超声B型图像数据信息进行分析，其节约了对RF射频数据分析的过程，提高了效率，改善了实时性；且本发明步骤4过程中采用滤波核的梯度算子的方法，使应变场处理过程鲁棒性更好。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应该涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (3)

1.一种基于光流法和应变计算的实时超声准静态弹性成像方法，其特征在于，包括以下步骤：

步骤1：获取经过波束合成的B型图像数据像素信息；

步骤2：对经过步骤1处理获得的B型图像数据像素信息进行多项式展开，使用多项式去逼近每个像素的邻域；

步骤3：通过对步骤2得到的每个像素的邻域进行位移估计，得到整帧图像的位移场；

步骤4：对步骤3获得的整帧图像的位移场进行卷积计算，得到B型图像数据整帧图像的应变场；

步骤5：对经过步骤4得到的应变场信息进行取对数法降噪声处理；

步骤6：对步骤5降噪后的应变场信息进行可视化、彩色处理得到彩色弹性图像。

2.如权利要求1所述的基于光流法和应变计算的实时超声准静态弹性成像方法，其特征在于，所述步骤4的位移场卷积为一个预先设置的滤波核，用于计算位移场的偏导，得到整帧图像的应变场。

3.如权利要求2所述的基于光流法和应变计算的实时超声准静态弹性成像方法，其特征在于，所述滤波核是一个N行*1列的梯度算子，计算公式如下：

$y(M+k)=\frac{25*(3M^4+6M^3-3M+1)*k-35*(3M^2+3M-1)*k^3}{(2M+3)*(2M+1)*(2M-1)*(M+2)*(M+1)*M*(M-1)}$

$\begin{array}{c}y(M-k)=-\frac{25*(3M^4+6M^3-3M+1)*k-35*(3M^2+3M-1)*k^3}{(2M+3)*(2M+1)*(2M-1)*(M+2)*(M+1)*M*(M-1)}\\ y\left(M\right)=0;\end{array}$

其中，N＝2*M+1，k为大于0且不大于M的整数，1≤k≤M；该滤波核是一个一维滤波核，N是滤波核的直径，M指代滤波核的半径，M为大于1的整数；y(M-k)、y(M)、y(M+k)分别为在滤波核中第M-K行第一列、第M行第一列、第M+K行第一列的元素值。