CN105266849B - 实时超声弹性成像方法和*** - Google Patents
实时超声弹性成像方法和*** Download PDFInfo
- Publication number
- CN105266849B CN105266849B CN201410326830.3A CN201410326830A CN105266849B CN 105266849 B CN105266849 B CN 105266849B CN 201410326830 A CN201410326830 A CN 201410326830A CN 105266849 B CN105266849 B CN 105266849B
- Authority
- CN
- China
- Prior art keywords
- gpu
- image
- biological tissue
- msup
- mrow
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明提供一种实时超声弹性成像方法和***,成像方法主要包含下述步骤:获取生物组织受外力缓速挤压后产生形变的不同状态下的B模式图像;选定两帧图像的ROI区域;使用光流法求取位移场,在计算过程中将需加速运算的含重复计算的子步骤的参数信息传输给至第一GPU数据缓存模块;由GPU工作组进行需加速运算的各子步骤的运算处理;将各子步骤的运算结果信息带入光流法运算过程,最终求出ROI区域的位移场即光流场;获取图像ROI区域的轴向应变场;对轴向应变场信息进行降噪声处理;对获取的信息进行彩色化处理;对ROI区域中的对象区域进行弹性分级。该实时超声弹性成像方法运算速度快,精度较高,同时具备较强的鲁棒性。
Description
技术领域
本发明涉及超声波回波成像领域,尤其是一种实时超声弹性成像方法和***。
背景技术
超声波回波成像技术目前已经被广泛应用于军事、医疗等领域。超声波回波成像中关于弹性成像的概念最早由Ophir等人于1991年首次提出。之后,弹性成像技术在近二十年里得到了迅猛发展,它被称为继A型、B型、D型、M型超声之后的E型模式。超声弹性成像是通过超声成像***进行生物组织弹性参数成像,超声弹性图能够提供传统B超图像所无法反映的生物组织弹性特征,对于肿瘤检测等临床应用有非常大的帮助。由于具有易于实现、适用于实时诊断和对组织无侵害性等优点,超声弹性成像技术受到了广泛的关注。
现有的弹性成像方法中,计算位移场的方法大多基于互相关算法和改进的互相关算法,其使用的数据为RF射频数据,而不能直接使用B模式图像,由于处理RF射频数据需要很大的运算,因此造成了现有的方法计算效率低下的问题;光流法是指空间运动物体在观察成像平面上的像素运动的瞬时速度,利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系,从而计算出相邻帧之间物体的运动信息的一种方法。计算应变的方法大多基于最小二乘拟合算法和梯度法,使用这些方法计算出的应变,在精度方面还有提高的空间。
1998年,Negahdaripour将光流(Optical Flow)定义为动态图像的几何变化和辐射度变化的全面表示。光流法利用图像序列中的像素强度数据的时域变化和相关性来确定各自像素位置的运动,表达了图像灰度在时间上的变化与景象中物体结构及其运动的关系。用于评价光流法精度的一个经典标准来自于Middlebury,其用于光流法评价的数据集有12张图像,除此之外还有Sintel和KITTI。
光流法的基本原理如下:
令I(x,y,t)表示在t时刻图像上像素(x,y)处的亮度(或颜色),那么光流法的目的就是求出在t+1时刻的图像上,该像素相对于原来(x,y)的位移量(u,v),用方程表示即:
I(x+u,y+v,t+1)=I(x,y,t) (1)
其中u和v是待求解的位移量。该方程被称为Brightness Constancy Model(亮度恒定模型)。
在经典的光流法中,一般利用一阶泰勒展开作为工具来建立图像梯度和位移之间的关系,这一步骤通常被称为线性化。具体原理如下:
假设图像的亮度是连续的,如图2的一维示例,曲线1表示frame1(帧1)中的图像,曲线2表示frame2(帧2)中的图像,待求的位移是向右的箭头对曲线进行一阶泰勒展开,其实就是假设曲线的局部是线性的,这样可以考察如图2的向右箭头、向上箭头、粗线段组成的三角形。并非是粗线段的长度,而是其斜率。这样可以得到图中所示的关系,注意负号是因为斜率其实表示的是钝角的tan值。这样一来就建立了图像的导数和位移之间的关系,注意是图像在空间上的导数,是图像在时间上的导数。
将图2中的方程用在二维的图像上,对于每个像素,可以写出以下方程:
Ixu+Iyv+It=0 (2)
其中,Ix和Iy是图像沿空间的x和y两个方向的导数,即图像梯度的两个分量;It是图像沿时间变化的导数,可以用两帧图像的差异来近似;u和v是该像素沿x和y方向的位移,也就是待求的光流未知量。该方程是Brightness Constancy Model(亮度恒定模型)线性化的结果,被称为Gradient Constraint Equation(梯度约束方程)。
需要注意的是,上面这个模型是建立在两个假设基础上:第一,图像变化是连续的;第二,位移不是很大。如果这两点假设不成立,即图像不连续且位移很大,则无法将位移和图像梯度联系起来,泰勒展开得到的近似结果将很差。在实际操作中,可以很容易使上述两个假设其成立。关于第一点假设,一般可以预先对图像进行高斯平滑,使其变化较为平缓;关于第二点假设,一般对图像降分辨率建立金字塔,通过Coarse-To-Fine(由粗到精)的方式去求解。
基于上面的方程,产生了最经典的两个光流法:Lucas-Kanade方法和Horn-Schunck方法,他们分别从不同的角度增加了求解该方程的稳定性。Lucas-Kanade方法是将每个像素周围的一些像素考虑进来,每个像素的未知量单独求解,是一种局部光流方法;而Horn-Schunck方法是将上面的方程纳入到一个规则化的框架中,将局部平滑的优先级加入光流计算,所有像素的未知位移量之间相互依赖,需要用全局优化的方法求解。
当前如何优化现有弹性成像***,是技术人员需要考虑的研发方向,现有超声***大都使用CPU进行全局的运算,本发明运用GPU进行弹性成像***中图像处理运算,能够提高现有超声弹性成像***的运算速度,减少CPU的负载从而提高了原有***的稳定性。同时,随着日益增长的诊断需求,如何快速准确进行辅助医生进行疾病的诊断,也是医疗领域的发展方向。
发明内容
本发明的第一目的在于提供一种实时超声弹性成像方法,运算速度快,精度较高,同时具备较强的鲁棒性。本发明还同时提出了一种实时超声弹性成像***。本发明采用的技术方案是:
一种实时超声弹性成像方法,包括下述步骤:
S10.对生物组织进行缓速挤压时对生物组织目标区域进行超声波扫查并接收回波信号;
S20.对步骤S10中接收的回波信号进行处理,形成线数据;
S30.对步骤S20得到的线数据进行处理,形成生物组织受外力缓速挤压后产生形变的不同状态下的B模式图像;
S40.对生物组织受外力缓速挤压后产生形变的不同状态下的两帧B模式图像选定两个在不同帧图像的同一位置区域,即选定两帧图像的ROI区域;
S50.通过第一CPU数据处理模块对该两帧图像的ROI区域使用光流法求取位移场;并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数据缓存模块,进行缓存;
S60.第一GPU数据缓存模块将接收到的子步骤的参数信息传输给各自的GPU数据处理模块进行数据处理,S50中的含重复计算的各子步骤的运算分配给至少一个GPU核工作单元,由至少一个GPU核工作单元组成GPU工作组;
S70.将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存中,等待所有GPU工作组的参数信息映射完成,然后各GPU工作组进行步骤S50中的各子步骤的运算处理;
S80.将GPU工作组的运算结果信息传输至第二GPU数据缓存模块;
S90.第二GPU数据缓存模块将运算结果信息传输至第二CPU数据处理模块,第二CPU数据处理模块读出运算结果信息并带入光流法运算过程,最终求出ROI区域的位移场即光流场;
S100.使用低通滤波器对轴向光流场求取应变,得到图像ROI区域的轴向应变场;
S105.对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理;
S110.对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物组织彩色弹性图像。
所述步骤S20中,通过波束合成模块对回波信号具体进行A/D转换、同相位叠加处理;
所述步骤S30中,通过B模式图像处理模块对步骤S20得到的线数据具体进行降噪、滤波、提高信噪比处理。
进一步地,所述步骤S50中,含重复计算的子步骤的参数信息是指:构建图像金字塔、中心梯度计算、计算目标图像及其导数的形变、估计双变量、估计位移这些计算的参数信息。
进一步地,所述步骤S100中,低通滤波器的核如下列公式所示:
该滤波核的大小为N行*1列,N为正奇数,其中,M=(N-1)/2,k为大于0且不大于M的整数,1≤k≤M,y代表滤波核不同位置的值。
更进一步地,滤波核的N取值不大于31。
进一步地,所述步骤S110具体为:首先将由轴向应变场构成的灰度图像经波段合成转换成彩色图像,再将该彩色图像映射到自然颜色***、孟塞尔颜色***、PANTONE色卡颜色***、TILO管理颜色***,或者自定义的颜色***的至少一种,形成彩色化的生物组织弹性图像。
进一步地,步骤S110之后还包括步骤S120,对ROI区域中的对象区域进行弹性分级,当生物组织弹性低于预设门限时显示相应等级的报警图标。
步骤S120中的弹性分级具体包括:
首先,提取出ROI区域中的对象区域;
然后,统计对象区域所有像素值分别位于各弹性值区间的个数;各弹性值区间依弹性值由低到高分布;
再计算出每段弹性值区间的像素个数占该对象区域总像素的比例;
将占对象区域总像素的比例最高的弹性值区间的弹性级别判为该对象区域生物组织的弹性级别。
本发明提出的一种实时超声弹性成像***,包括:
换能器,用于对生物组织进行缓速挤压时对生物组织目标区域进行超声波扫查并接收回波信号;
波束合成模块,用于对接收的回波电信号进行A/D转换、同相位叠加处理,形成线数据;
B模式图像处理模块,用于对上述线数据进行处理,形成生物组织受外力缓速挤压后产生形变的不同状态下的B模式图像;
第一CPU数据处理模块,用于对生物组织受外力缓速挤压后产生形变的不同状态下的两帧B模式图像的同一位置区域即ROI区域使用光流法求取位移场;并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数据缓存模块;
第一GPU数据缓存模块,用于缓存接收到的含重复计算的子步骤的参数信息并传输给各自的GPU数据处理模块进行数据处理;
一个或多个GPU数据处理模块,用于将上述各子步骤的运算分配给至少一个GPU核工作单元,将至少一个GPU核工作单元组成GPU工作组,将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存中,等待所有GPU工作组的参数信息映射完成,然后各GPU工作组进行光流法中含重复计算的各子步骤的运算处理;
第二GPU数据缓存模块,用于缓存GPU数据处理模块的运算结果信息并传输至第二CPU数据处理模块;
第二CPU数据处理模块,用于读出运算结果信息并带入光流法运算过程,最终求出ROI区域的位移场即光流场;使用低通滤波器对轴向光流场求取应变,得到图像ROI区域的轴向应变场;对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理;
彩色定量显示模块,用于对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物组织彩色弹性图像;
图像辅助预诊断模块,用于对ROI区域中的对象区域采用统计的方法进行弹性分级,当生物组织弹性低于预设门限时控制显示单元显示相应等级的报警图标;
显示单元,用于显示彩色化的生物组织弹性图像,以及弹性等级对应的报警图标。
本发明的优点在于:
1.采用本发明的光流法求取的位移场精度高,能够显著提高图像处理的精度。
2.采用GPU加速的方法解决了光流法效率低下的问题,因此本发明既能求得准确的生物组织弹性,又满足了实际使用中弹性计算的实时性要求。
3.使用对生物组织的弹性进行分级的方法,相对于直接使用阈值判定法更科学合理。
4.采用GPU加速运算大量重复子步骤信息,能够提高***的鲁棒性,减少***卡顿的现象,提高***的流畅性;同时降低了CPU的工作负荷,有利于低配置的超声诊断设备进行超声弹性成像处理或版本的升级。
附图说明
图1为本发明的结构组成示意图。
图2为求取位移场的光流法的原理示意图。
图3为本发明的滤波核示意图。
图4为本发明的流程图。
具体实施方式
下面结合具体附图和实施例对本发明作进一步说明。
本发明所提出的实时超声弹性成像方法,由CPU与多核GPU(Graphic ProcessingUnit)联合实现,该方法步骤如下:
S10.通过超声换能器对生物组织进行缓速挤压时对生物组织目标区域进行超声波扫查并接收回波信号;
S20.通过波束合成模块对步骤S10中的接收的回波电信号进行A/D转换、同相位叠加等处理,形成线数据;
S30.通过B模式图像处理模块对步骤S20得到的线数据进行降噪、滤波、提高信噪比等处理,形成生物组织受外力缓速挤压后产生形变的不同状态下的B模式图像;
S40.通过使用ROI框(感兴趣区域)对生物组织受外力缓速挤压后产生形变的不同状态下的两帧B模式图像选定两个在不同帧图像的同一位置区域;
ROI(Region Of Interest),即感兴趣区域,在图像处理中是指从待处理图像中提取出的要处理的区域。ROI的使用一方面可以突出框内区域的特征,另一方面又可以提高图像处理的速度。本***预设了一个ROI框,用户可以调整ROI框的大小和位置以满足不同的要求。
S50.通过第一CPU数据处理模块对该两帧图像的ROI区域使用光流法求取位移场;并在计算过程中将需要加速运算的含大量重复计算的子步骤的参数信息传输至第一GPU数据缓存模块,进行缓存;其中需加速运算的含大量重复计算的子步骤的参数信息是指:构建图像金字塔、中心梯度计算、计算目标图像及其导数的形变、估计双变量、估计位移等计算的参数信息。
光流法求取位移场公式如下:
其中小写的u和v分别表示每个像素沿x和y方向的位移;大写的U和V分别表示由所有像素的u和v组成的位移场;t代表时间参数;λ代表一个常量,代表权重系数,λ≤6;10-6为此函数的一个优化系数。
S60.第一GPU数据缓存模块将接收到的需加速运算的子步骤的参数信息传输给各自的GPU数据处理模块进行数据处理,S50中需加速运算的各子步骤的运算分配给至少一个GPU核工作单元,由至少一个GPU核工作单元组成GPU工作组。需加速运算的含大量重复计算的子步骤的参数信息是指:构建图像金字塔、中心梯度计算、计算目标图像及其导数的形变、估计双变量、估计位移这些计算的参数信息。例如:
构建图像金字塔的参数信息,包括:上下文、上一层图像的宽与高、下一层图像的宽与高、上一层图像数据、下一层图像数据、图像通道数、图像位深、图像步长等。
中心梯度计算的参数信息,包括:图像金字塔该层的图像数据、上下文、图像宽与高、图像步长、x方向导数、y方向导数、x方向导数步长等。
计算目标图像及导数的形变的参数信息,包括:图像金字塔该层的图像数据、上下文、图像宽与高、图像纹理、图像步长、x与y方向位移量等。
估计双变量的参数信息,包括:初步估计的位移数据、上下文、全局线程、本地线程、tau值、x与y方向位移量等。
估计位移的参数信息,包括:上下文、初步估计的位移数据、全局线程、本地线程、theta值、图像宽与高、图像步长、容错数据、x与y方向位移量等。
S70.将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存中,等待所有GPU工作组的参数信息映射完成,然后各GPU工作组进行步骤S50中需加速运算的各子步骤的运算处理;GPU数据处理模块单元由至少一个的GPU工作组组成,且GPU工作组包含至少一个GPU工作单元。
S80.将GPU工作组的运算结果信息传输至第二GPU数据缓存模块;
此步骤将需要GPU加速的各子步骤经过S60、S70运算的运算结果信息传输至第二GPU数据缓存模块;
S90.第二GPU数据缓存模块将运算结果信息传输至第二CPU数据处理模块,第二CPU数据处理模块读出运算结果信息并带入光流法运算过程,最终求出ROI区域的位移场即光流场;(用光流法计算出的位移场也叫光流场)。
将第二GPU数据缓存模块中的信息由第二CPU数据处理模块读出,并将处理结果带入光流法的运算过程。待所有需要GPU加速的处理子步骤均完成并将各自的信息带入光流法,且完成计算后,图像ROI区域的光流场即求得。
在实际处理中,上述方法由第一CPU数据处理模块、第一GPU数据缓存模块、GPU数据处理模块、第二GPU数据缓存模块、第二CPU数据处理模块共同进行处理。这样,将光流法中构建图像金字塔、中心梯度计算、计算目标图像及其导数的形变、估计双变量、估计位移等巨大繁重而又重复的工作,交由GPU工作组的至少一个的GPU核工作单元来处理,当然一般而言GPU核工作单元的数量大于2个,通过这样的方法大大提高了运行速度,对于提高算法的实时性、扩展算法的应用领域具有非常重要和积极的作用,减少了***卡顿现象,从而提高了超声仪器的整体***的稳定性,提高了使用人员操作感观的舒适度。
S100.使用低通滤波器对轴向光流场求取应变,得到图像ROI区域的轴向应变场;此步骤在第二CPU数据处理模块中进行。
将光流法得到的位移场分解成X方向位移场和Y方向位移场,即横向位移场与轴向位移场。由于轴向位移场对应着探头按压、组织受力的方向,因此本发明更关注轴向位移场和轴向应变场。使用一个低通滤波器对轴向位移场进行卷积,计算轴向位移场的偏导,得到图像ROI区域的轴向应变场,从而反映ROI区域的组织弹性。
本发明采用的低通滤波器的核如下列公式所示:
如图3所示,该滤波核的大小为N行*1列(N为正奇数),其中,M=(N-1)/2,k为大于0且不大于M的整数(1≤k≤M),y代表滤波核不同位置的值。
由于本发明所采用的光流法的精度高,离群值较少,因此该滤波核中的N可以取较小的值,即滤波核可以设计得较小,这样既不会降低应变计算的效果,又能提高效率,还能保存一些图像细节。本发明采用的滤波核的N取值不大于31。
S105.对经过滤波得到的图像ROI区域的轴向应变场信息进行取对数法降噪声处理;此步骤在第二CPU数据处理模块中进行。
对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理。首先对上个步骤得到的图像ROI区域轴向应变场取绝对值abs(dVR);再对绝对值使用取对数的方法log(1+abs(dVR))去除噪声;然后将其归一化到0~255区间的整数,得到降噪后的图像ROI区域的轴向应变场。
S110.对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物组织彩色弹性图像。此彩色化定量显示的步骤在彩色定量显示模块中完成。
首先将由轴向应变场构成的灰度图像经波段合成转换成彩色图像,再将该彩色图像映射到自然颜色***、孟塞尔颜色***、PANTONE色卡颜色***、TILO管理颜色***,或者自定义的颜色***中的至少一种,形成彩色化的生物组织弹性图像。
S120.对ROI区域中的对象区域采用统计的方法进行弹性分级,当生物组织弹性低于预设门限时显示相应等级的报警图标。弹性分级在图像辅助预诊断模块中进行,显示单元可显示相应的报警图标。
生物组织弹性低于预设门限时,分为4个等级,即严重、中等、轻微、警告,测量值达到相应等级即显示相应的图标。
首先,使用阈值分割方法或其他图像分割方法提取出ROI区域中组织弹性较小的区域,该组织弹性较小的区域为疑似病变区域,也就是关注的对象区域。该疑似病变区域通常是ROI区域的一部分子区域,也有可能充满整个ROI区域(这时病变部位非常大)。然后,根据***预设值Value_verylow、Value_low、Value_middle、Value_high,统计该对象区域所有像素弹性值分别位于[0,Value_verylow)、[Value_verylow,Value_low)、[Value_low,Value_middle)、[Value_middle,Value_high)四段区间的个数,再计算出每段区间的像素个数占该组织弹性较小的区域总像素的比例。若[0,Value_verylow)区间的比例最高,则判定该ROI区域中疑似病变区域的组织弹性非常低,***上的图标显示严重;若[Value_verylow,Value_low)区间的比例最高,则判定该ROI区域中疑似病变区域的组织弹性很低,***上的图标显示中等;若[Value_low,Value_middle)区间的比例最高,则判定该ROI区域中疑似病变区域的组织弹性较低,***上的图标显示轻微;若[Value_middle,Value_high)区间的比例最高,则判定该ROI区域中疑似病变区域的组织弹性有些低,***上的图标显示警告。使用该统计的方法对生物组织的弹性进行分级,相对于直接使用阈值判定法更科学合理。
本发明提出的实时超声弹性成像***如图1所示,包括:
换能器,用于对生物组织进行缓速挤压时对生物组织目标区域进行超声波扫查并接收回波信号;
波束合成模块,用于对接收的回波电信号进行A/D转换、同相位叠加处理,形成线数据;
B模式图像处理模块,用于对上述线数据进行处理,形成生物组织受外力缓速挤压后产生形变的不同状态下的B模式图像;
第一CPU数据处理模块,用于对生物组织受外力缓速挤压后产生形变的不同状态下的两帧B模式图像的同一位置区域即ROI区域使用光流法求取位移场;并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数据缓存模块;
第一GPU数据缓存模块,用于缓存接收到的含重复计算的子步骤的参数信息并传输给各自的GPU数据处理模块进行数据处理;
一个或多个GPU数据处理模块,用于将上述各子步骤的运算分配给至少一个GPU核工作单元,将至少一个GPU核工作单元组成GPU工作组,将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存中,等待所有GPU工作组的参数信息映射完成,然后各GPU工作组进行光流法中含重复计算的各子步骤的运算处理;
第二GPU数据缓存模块,用于缓存GPU数据处理模块的运算结果信息并传输至第二CPU数据处理模块;
第二CPU数据处理模块,用于读出运算结果信息并带入光流法运算过程,最终求出ROI区域的位移场即光流场;使用低通滤波器对轴向光流场求取应变,得到图像ROI区域的轴向应变场;对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理;
彩色定量显示模块,用于对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物组织彩色弹性图像;
图像辅助预诊断模块,用于对ROI区域中的对象区域采用统计的方法进行弹性分级,当生物组织弹性低于预设门限时控制显示单元显示相应等级的报警图标;
显示单元,用于显示彩色化的生物组织弹性图像,以及弹性等级对应的报警图标。
显示单元可以是常规的台式机显示器,也可以是触摸屏显示器或手机等终端接收显示单元。
Claims (9)
1.一种实时超声弹性成像方法,其特征在于,包括下述步骤:
步骤S10.对生物组织进行缓速挤压时对生物组织目标区域进行超声波扫查并接收回波信号;
步骤S20.对步骤S10中接收的回波信号进行处理,形成线数据;
步骤S30.对步骤S20得到的线数据进行处理,形成生物组织受外力缓速挤压后产生形变的不同状态下的B模式图像;
步骤S40.对生物组织受外力缓速挤压后产生形变的不同状态下的两帧B模式图像选定两个在不同帧图像的同一位置区域,即选定两帧图像的ROI区域;
步骤S50.通过第一CPU数据处理模块对该两帧图像的ROI区域使用光流法求取位移场;并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数据缓存模块,进行缓存;
步骤S60.第一GPU数据缓存模块将接收到的子步骤的参数信息传输给各自的GPU数据处理模块进行数据处理,S50中的含重复计算的各子步骤的运算分配给至少一个GPU核工作单元,由至少一个GPU核工作单元组成GPU工作组;
步骤S70.将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存中,等待所有GPU工作组的参数信息映射完成,然后各GPU工作组进行步骤S50中的各子步骤的运算处理;
步骤S80.将GPU工作组的运算结果信息传输至第二GPU数据缓存模块;
步骤S90.第二GPU数据缓存模块将运算结果信息传输至第二CPU数据处理模块,第二CPU数据处理模块读出运算结果信息并带入光流法运算过程,最终求出ROI区域的位移场即光流场;
步骤S100.使用低通滤波器对轴向光流场求取应变,得到图像ROI区域的轴向应变场;
步骤S105.对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理;
步骤S110.对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物组织彩色弹性图像。
2.如权利要求1所述的实时超声弹性成像方法,其特征在于:
所述步骤S20中,通过波束合成模块对回波信号具体进行A/D转换、同相位叠加处理;
所述步骤S30中,通过B模式图像处理模块对步骤S20得到的线数据具体进行降噪、滤波、提高信噪比处理。
3.如权利要求1所述的实时超声弹性成像方法,其特征在于:
所述步骤S50中,含重复计算的子步骤的参数信息是指:构建图像金字塔、中心梯度计算、计算目标图像及其导数的形变、估计双变量、估计位移这些计算的参数信息。
4.如权利要求1所述的实时超声弹性成像方法,其特征在于:
所述步骤S100中,低通滤波器的核如下列公式所示:
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mi>y</mi>
<mo>(</mo>
<mi>M</mi>
<mo>+</mo>
<mi>k</mi>
<mo>)</mo>
<mo>=</mo>
<mfrac>
<mrow>
<mn>25</mn>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<msup>
<mi>M</mi>
<mn>4</mn>
</msup>
<mo>+</mo>
<mn>6</mn>
<msup>
<mi>M</mi>
<mn>3</mn>
</msup>
<mo>-</mo>
<mn>3</mn>
<mi>M</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>*</mo>
<mi>k</mi>
<mo>-</mo>
<mn>35</mn>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<msup>
<mi>M</mi>
<mn>2</mn>
</msup>
<mo>+</mo>
<mn>3</mn>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>*</mo>
<msup>
<mi>k</mi>
<mn>3</mn>
</msup>
</mrow>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mi>M</mi>
<mo>+</mo>
<mn>3</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mn>2</mn>
<mi>M</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mn>2</mn>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mi>M</mi>
<mo>+</mo>
<mn>2</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mi>M</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
<mo>*</mo>
<mi>M</mi>
<mo>*</mo>
<mo>(</mo>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mfrac>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>y</mi>
<mo>(</mo>
<mi>M</mi>
<mo>-</mo>
<mi>k</mi>
<mo>)</mo>
<mo>=</mo>
<mo>-</mo>
<mfrac>
<mrow>
<mn>25</mn>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<msup>
<mi>M</mi>
<mn>4</mn>
</msup>
<mo>+</mo>
<mn>6</mn>
<msup>
<mi>M</mi>
<mn>3</mn>
</msup>
<mo>-</mo>
<mn>3</mn>
<mi>M</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>*</mo>
<mi>k</mi>
<mo>-</mo>
<mn>35</mn>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<msup>
<mi>M</mi>
<mn>2</mn>
</msup>
<mo>+</mo>
<mn>3</mn>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>*</mo>
<msup>
<mi>k</mi>
<mn>3</mn>
</msup>
</mrow>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mi>M</mi>
<mo>+</mo>
<mn>3</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mn>2</mn>
<mi>M</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mn>2</mn>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mi>M</mi>
<mo>+</mo>
<mn>2</mn>
<mo>)</mo>
<mo>*</mo>
<mo>(</mo>
<mi>M</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
<mo>*</mo>
<mi>M</mi>
<mo>*</mo>
<mo>(</mo>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mfrac>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>y</mi>
<mo>(</mo>
<mi>M</mi>
<mo>)</mo>
<mo>=</mo>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
该低通滤波器的核的大小为N行*1列,N为正奇数,其中,M=(N-1)/2,k为大于0且不大于M的整数,1≤k≤M,y(M+k)、y(M-k)、y(M)分别代表低通滤波器的核不同位置的值。
5.如权利要求4所述的实时超声弹性成像方法,其特征在于:
低通滤波器的核的N取值不大于31。
6.如权利要求1所述的实时超声弹性成像方法,其特征在于:
所述步骤S110具体为:首先将由轴向应变场构成的灰度图像经波段合成转换成彩色图像,再将该彩色图像映射到自然颜色***、孟塞尔颜色***、PANTONE色卡颜色***、TILO管理颜色***,或者自定义的颜色***的至少一种,形成彩色化的生物组织弹性图像。
7.如权利要求1~6中任一项所述的实时超声弹性成像方法,其特征在于:
步骤S110之后还包括步骤S120,对ROI区域中的对象区域进行弹性分级,当生物组织弹性低于预设门限时显示相应等级的报警图标。
8.如权利要求7所述的实时超声弹性成像方法,其特征在于:步骤S120中的弹性分级具体包括:
首先,提取出ROI区域中的对象区域;
然后,统计对象区域所有像素值分别位于各弹性值区间的个数;各弹性值区间依弹性值由低到高分布;
再计算出每段弹性值区间的像素个数占该对象区域总像素的比例;
将占对象区域总像素的比例最高的弹性值区间的弹性级别判为该对象区域生物组织的弹性级别。
9.一种实时超声弹性成像***,其特征在于,包括:
换能器,用于对生物组织进行缓速挤压时对生物组织目标区域进行超声波扫查并接收回波信号;
波束合成模块,用于对接收的回波电信号进行A/D转换、同相位叠加处理,形成线数据;
B模式图像处理模块,用于对上述线数据进行处理,形成生物组织受外力缓速挤压后产生形变的不同状态下的B模式图像;
第一CPU数据处理模块,用于对生物组织受外力缓速挤压后产生形变的不同状态下的两帧B模式图像的同一位置区域即ROI区域使用光流法求取位移场;并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数据缓存模块;
第一GPU数据缓存模块,用于缓存接收到的含重复计算的子步骤的参数信息并传输给各自的GPU数据处理模块进行数据处理;
一个或多个GPU数据处理模块,用于将上述各子步骤的运算分配给至少一个GPU核工作单元,将至少一个GPU核工作单元组成GPU工作组,将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存中,等待所有GPU工作组的参数信息映射完成,然后各GPU工作组进行光流法中含重复计算的各子步骤的运算处理;
第二GPU数据缓存模块,用于缓存GPU数据处理模块的运算结果信息并传输至第二CPU数据处理模块;
第二CPU数据处理模块,用于读出运算结果信息并带入光流法运算过程,最终求出ROI区域的位移场即光流场;使用低通滤波器对轴向光流场求取应变,得到图像ROI区域的轴向应变场;对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理;
彩色定量显示模块,用于对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物组织彩色弹性图像;
图像辅助预诊断模块,用于对ROI区域中的对象区域采用统计的方法进行弹性分级,当生物组织弹性低于预设门限时控制显示单元显示相应等级的报警图标;
显示单元,用于显示彩色化的生物组织弹性图像,以及弹性等级对应的报警图标。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410326830.3A CN105266849B (zh) | 2014-07-09 | 2014-07-09 | 实时超声弹性成像方法和*** |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410326830.3A CN105266849B (zh) | 2014-07-09 | 2014-07-09 | 实时超声弹性成像方法和*** |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105266849A CN105266849A (zh) | 2016-01-27 |
CN105266849B true CN105266849B (zh) | 2017-10-17 |
Family
ID=55137283
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410326830.3A Active CN105266849B (zh) | 2014-07-09 | 2014-07-09 | 实时超声弹性成像方法和*** |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105266849B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107280703B (zh) * | 2016-07-22 | 2019-12-24 | 珠海医凯电子科技有限公司 | 基于gpu平台的实时3d超声扫描变换方法 |
CN106725609A (zh) * | 2016-11-18 | 2017-05-31 | 乐普(北京)医疗器械股份有限公司 | 一种弹性检测方法和装置 |
CN106805997B (zh) * | 2016-12-26 | 2020-08-07 | 乐普(北京)医疗器械股份有限公司 | 一种弹性成像方法和装置 |
CN106901776B (zh) * | 2017-01-11 | 2019-07-26 | 中国人民解放军第三军医大学第三附属医院 | 基于可变滤波器长度的超声弹性成像方法 |
CN107958214A (zh) * | 2017-11-21 | 2018-04-24 | 中国科学院深圳先进技术研究院 | Ecg信号的并行分析装置、方法和移动终端 |
CN110332987B (zh) * | 2019-08-22 | 2020-09-01 | 广东电网有限责任公司 | 一种声纹信号成像方法及麦克风阵列信号的成像方法 |
CN114881923A (zh) * | 2022-03-30 | 2022-08-09 | 什维新智医疗科技(上海)有限公司 | 一种基于超声图像弹性信号的结节回声分析装置 |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6520915B1 (en) * | 2000-01-28 | 2003-02-18 | U-Systems, Inc. | Ultrasound imaging system with intrinsic doppler capability |
CN1298290C (zh) * | 2004-07-23 | 2007-02-07 | 清华大学 | 超声弹性成像的平衡测压装置 |
CN101569543B (zh) * | 2008-04-29 | 2011-05-11 | 香港理工大学 | 弹性成像的二维位移估计方法 |
WO2013026141A1 (en) * | 2011-08-19 | 2013-02-28 | The University Of British Columbia | Elastography using ultrasound imaging of a thin volume |
CN102626327B (zh) * | 2012-04-26 | 2014-02-19 | 声泰特(成都)科技有限公司 | 基于接收端空间复合的超声弹性成像及压力反馈方法 |
CN104349817B (zh) * | 2012-05-29 | 2017-12-15 | 皇家飞利浦有限公司 | 用于辐射治疗中的经改进的门控效率和动态裕度调节的基于弹性成像的方法 |
CN103040488B (zh) * | 2012-12-21 | 2014-06-04 | 深圳大学 | 一种实时超声弹性成像位移估计方法和*** |
CN103815932B (zh) * | 2014-02-17 | 2016-06-22 | 无锡祥生医学影像有限责任公司 | 基于光流法和应变的超声准静态弹性成像方法 |
-
2014
- 2014-07-09 CN CN201410326830.3A patent/CN105266849B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN105266849A (zh) | 2016-01-27 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105266849B (zh) | 实时超声弹性成像方法和*** | |
US11229419B2 (en) | Method for processing 3D image data and 3D ultrasonic imaging method and system | |
US5795296A (en) | Pipeline process for automatically measuring object boundary from ultrasound image samples | |
US9113826B2 (en) | Ultrasonic diagnosis apparatus, image processing apparatus, control method for ultrasonic diagnosis apparatus, and image processing method | |
CN106127711A (zh) | shearlet变换和快速双边滤波器图像去噪方法 | |
CN101887581B (zh) | 图像融合方法及设备 | |
JP7296171B2 (ja) | データ処理方法、装置、機器及び記憶媒体 | |
DE102012108121A1 (de) | Verfahren und System für ultraschallgestützte automatische Erkennung, Quantifizierung und Nachverfolgung von Pathologien | |
JP2013542046A5 (zh) | ||
CN101909165B (zh) | 一种基于混合测度的视频数据宽景成像方法 | |
CN111820948B (zh) | 胎儿生长参数测量方法、***及超声设备 | |
Chen et al. | Real-time freehand 3D ultrasound imaging | |
CN113994367A (zh) | 用于生成合成弹性成像图像的方法和*** | |
CN114565572A (zh) | 一种基于图像序列分析的脑出血ct图像分类方法 | |
KR100760251B1 (ko) | 초음파 영상 처리 시스템 및 방법 | |
CN107146211A (zh) | 基于线扩散函数和双边滤波的视网膜血管图像降噪方法 | |
EP2425784B1 (en) | Providing a color Doppler mode image in an ultrasound system | |
TW202027090A (zh) | 醫療矢面影像的取得方法、神經網路的訓練方法及計算機裝置 | |
US20220249060A1 (en) | Method for processing 3d image data and 3d ultrasonic imaging method and system | |
CN107076820A (zh) | 用于评估和提高在精细结构分析数据中的数据质量的方法 | |
CN110930394B (zh) | 测量肌肉肌纤维束线斜率和羽状角的方法及终端设备 | |
CN111062935B (zh) | 一种乳腺肿瘤检测方法、存储介质及终端设备 | |
CN111242853B (zh) | 基于光流处理的医学ct图像去噪方法 | |
Cong et al. | Fast and automatic ultrasound simulation from ct images | |
US9307948B2 (en) | Image processing apparatus, method, and program |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CP01 | Change in the name or title of a patent holder |
Address after: 214028 Changjiang Road in Jiangsu province Wuxi new area, industrial park five, 51 No. 53 No. 228 block Patentee after: Wuxi CHISON medical Polytron Technologies Inc Address before: 214028 Changjiang Road in Jiangsu province Wuxi new area, industrial park five, 51 No. 53 No. 228 block Patentee before: Xiangsheng Medical Image Co., Ltd., Wuxi |
|
CP01 | Change in the name or title of a patent holder |