CN103735287B - 一种血管内超声弹性成像二维多级混合位移估计方法 - Google Patents

一种血管内超声弹性成像二维多级混合位移估计方法 Download PDF

Info

Publication number
CN103735287B
CN103735287B CN201310645210.1A CN201310645210A CN103735287B CN 103735287 B CN103735287 B CN 103735287B CN 201310645210 A CN201310645210 A CN 201310645210A CN 103735287 B CN103735287 B CN 103735287B
Authority
CN
China
Prior art keywords
point
level
compression
window
displacement
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
Application number
CN201310645210.1A
Other languages
English (en)
Other versions
CN103735287A (zh
Inventor
钱馨然
崔崤峣
顾天明
焦阳
吕铁军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Suzhou Institute of Biomedical Engineering and Technology of CAS
Original Assignee
Suzhou Institute of Biomedical Engineering and Technology of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Suzhou Institute of Biomedical Engineering and Technology of CAS filed Critical Suzhou Institute of Biomedical Engineering and Technology of CAS
Priority to CN201310645210.1A priority Critical patent/CN103735287B/zh
Publication of CN103735287A publication Critical patent/CN103735287A/zh
Application granted granted Critical
Publication of CN103735287B publication Critical patent/CN103735287B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Abstract

本发明公开了一种血管内超声弹性成像二维混合位移估计方法,本发明提出的二维混合位移估计方法,结合了数据点的多级搜索和单级跟踪的方法,体现了合理平衡计算量和计算精度的优点。在算法初级阶段,利用二维多级搜索的方法来获取轴向和横向偏移量,而避免在高级阶段对横向方向进行搜索。在高级阶段,为最大程度减小计算量,则利用质量指导和零点相位差跟踪方法来估计轴向位移。这种算法可以很好地适用于血管内超声成像的实时成像要求,以及成像精度。

Description

一种血管内超声弹性成像二维多级混合位移估计方法
技术领域
[0001] 本发明涉及超声成像技术领域,具体涉及的是一种血管内超声弹性成像二维多级 混合位移估计方法。
背景技术
[0002] 血管内超声成像(intravascularultrasound,IVUS)是20世纪80年代末迅速发 展起来的将无创性的超声技术和有创性的心导管技术相结合的一种新的介入式超声成像 技术,是目前唯一商业化用于临床检测的可以实时提供患者冠状动脉血管横截面图像的检 查手段。具体方法是在心导管检查过程中,经导丝将超声探头送至心血管腔内进行探测,继 而回撤超声导管,可观测到管腔大小、管壁结构的厚度和形态等,观测其收缩和舒张的功能 变化,还能对病变进行定性分析甚至是精确测量分析。不仅能反映血管内腔的变化,同时也 能反映含有斑块的血管横断面结构以及斑块的性状等等,提供了冠状动脉造影技术所不能 提供的重要信息。目前IVUS成像技术已可对冠状动脉甚至更细小血管进行血管内成像,可 定性、定量地提供动脉壁微结构灰度图像,对于冠状动脉粥样硬化与狭窄等心血管疾病的 诊断与治疗具有重要意义。临床应用经验表明该技术具有直观、准确等优点,被认为是诊断 冠心病新的"金标准"。
[0003] 目前应用于IVUS系统的主要技术有虚拟组织学成像(VirtualHistology IVUS,VH-IVUS)和复合散斑成像(IntegratedBackscatterIVUS,IB-IVUS)等,均是基于超 声回波强度差别来区分富含脂质的低回声斑块即软斑块和富含纤维成分的高回声斑块即 硬斑块两种,但对于脂类或者半脂类的斑块不能够很好的进行定性,而且存在导管位置偏 离轴心和运动不均匀等伪影影响。
[0004] 1991年,Ophir和其研究小组初次提出了超声成像技术--弹性成像 (Elastography)以来,逐步成为当前国际医学超声领域最具价值的热点课题之一。弹性成 像方法不仅可以提供组织硬度的图像,也就是关于病变的组织特征的信息,以此来定位病 变,而且可以根据不同组织间弹性系数不同,在受到外力压迫后组织发生变形的程度不同, 将受压前后回声信号移动幅度的变化转化为实时彩色图像。将超声弹性成像方法运用到血 管内超声成像技术上,可以分辨高低应变区域,从而识别冠状动脉血管内不同斑块的相对 生物力学特征,辅助判断斑块破裂难易程度,进一步的提高心血管疾病的检测手段。
[0005] 自弹性成像提出以来,就有许多方法不断发展起来,总体来说可以分为两大类型: 一是一维方法,主要包括两种:1.基于梯度运算的应变估计算法,即先对组织的位移进行 估计,然后对位移估计做差分处理,得到组织的应变分布,如时域互相关算法、相位检测法、 过零点跟踪法等;2.直接的应变估计算法,不以组织位移为中间估计值,直接得到组织的 应变分布,如自适应伸展算法、谱相关算法等。具体方法是利用超声探头,采集组织压缩前 和压缩后的探头每一条扫描线的射频信号,再分别利用上述提到的算法对每一条扫描线进 行对应组织的应变分布。除分别具有自身的优缺点外,一维方法仅考虑换能器每一条扫描 线接收到的组织轴向方向运动的数据,而没有考虑组织在横向的运动,因此并没有包括组 织运动的全部信息。而现有的二维方法,优点是可以同时测量运动物体的轴向和横向运动 分量,不像一维运动估计算法那样受限于测量角度等问题。但是现有的二维算法中,如图像 相减算法(Heinetal.,1993),就是用一帧图像与另一帧图像作相减运算,如果两幅图像 中相同位置处的像素的亮度没有变化,则该像素在相减运动结果图像中亮度为零,表示组 织没有发生运动;而亮度有变化的区域,在相减运算结果图像中具有不同亮度,表示组织在 这个区域有运动。这种方法对存在较大差别的两幅图像的检测效果还可以,但亮度级数太 少,对比度不高。优点是具有实时性,缺点是不能得到运动的定量信息,检测不出运动方向。 或者二维区域匹配算法,以及基本相关算法,本质上是利用图像序列之间的一种对应关系, 找到目标区域的变化,但这种算法计算量非常大,十分耗时,不能很好的应用于实时的血管 内超声成像系统,不能满足其实时性要求。
发明内容
[0006] 本发明的目的在于克服现有技术存在的以上问题,提供一种血管内超声弹性成像 二维多级混合位移估计方法,很好地解决平衡计算时效和精度的问题。
[0007] 为实现上述技术目的,达到上述技术效果,本发明通过以下技术方案实现:
[0008] -种血管内超声弹性成像二维混合位移估计方法,包括以下步骤:
[0009] 步骤1)利用商用超声换能器得到待测生物组织的压缩前的一帧二维射频信号;
[0010] 步骤2)手持所述超声换能器沿着探头纵向对所述生物组织施加一个微小的挤压, 得到压缩后的一帧二维射频信号;
[0011] 步骤3)分别对从步骤1和步骤2中得到的组织压缩前、后的射频信号进行计算, 得到基带解析数据,具体步骤为先对射频数据进行下采样减少数据点,为避免下采样时组 织位移信息丢失,压缩前、后的信号采样频率应尽可能高,而下采样时,在不违背麦奎斯特 采样定律的前提下,下采样频率也应该尽可能高,再进行希尔伯特变换,乘以exp( ]Y\其中 w。为超声换能器的中心频率,则得到基带解析信号pre(X,y)和post(X,y),其中x表示沿 探头轴向方向采样点,y表示沿探头横向方向采样点;
[0012] 步骤4)对步骤3中得到的压缩前、后基带解析信号进行信号处理,分别得到压缩 前、后信号的幅值8 1(1,7)和82(1,7);
[0013] 步骤5)确定第一级窗口的数量及大小,可取9个,每行每列各3个;
[0014] (1)压缩前的目标区域第一级窗口大小为轴向或横向采样点的五分之一;
[0015] (2)压缩后的搜索区域第一级窗口大小略大于轴向或横向采样点的五分之一,可 取四分之一;
[0016] 步骤6)利用步骤5中确定的参数,对步骤4中得到的压缩前、后信号的幅值 ai(x,y)和&2(1,5〇进行第一级窗口搜索,搜索时对压缩前、后帧进行相关系数最大值匹配, 相关系数公式如式(1)所示;
[0017] (1)在最低分辨率下,压缩前和压缩后帧中的每一个窗口中,每四个米样点作为一 个对象点,图2中的上一图中的方形点为压缩前窗口的中心,对齐到压缩后的搜索区域进 行相关系数计算,找到最佳匹配的矩阵形区域并做记录;
[0018] (2)将分辨率提高一级,每两个采样点作为处理对象点,上一步中找到的外圆方形 点(图2中的下一图中的点)周围的x方向三个点作为这一步的匹配候选点,计算相关系数, 找到最佳匹配点,并做记录;
[0019] (3)分辨率最高的一级,不间隔点进行搜索,将上一步方形点周围的三个点作为候 选点,且计算相关系数时不再间隔任何采样点,找到最佳匹配点,完成第一级搜索,只有在 最低分辨率时,候选匹配位置尺寸与搜索区域的大小相关;在下一步分辨率提高的情况下, 对相关系数的要求高于计算量,候选匹配区域定为三个,得到9个独立的轴向-横向位移估 计值;
[0020]
Figure CN103735287BD00071
[0021] 其中I辦I#分别是横向和轴向位移,%分别是压缩前和压缩后的图像幅值, %.,:. %分别是%,_ %的内插窗口平均值,T是窗口大小;
[0022] 步骤7)确定第二级的窗口数量及大小,其数量与算法参数无关,其中有九个为第 一级9个窗口的子集,因此第一级搜索得到的轴向和横向位移估计就可以传递到第二级, 第二级其余的窗口则均匀分布整个帧,第二级运动的方法和第一级相同,但间隔因子初始 即取每两个采样点作为一个对象点,压缩前的目标区域窗口大小为轴向或横向采样点的 十五分之一,而压缩后的搜索区域窗口大小则略大于目标区域,可取轴向或横向采样点的 十分之一;
[0023] 步骤8)重复步骤6中的(1)、(2)、(3)步,但第二级的匹配过程由动态跟随来完 成,首先第二级的9个窗口由第一级的输出确定一个初始点pi(x0,y0),其余的窗口初始值 则由第一级的结果大致定向,定位于每一个窗口中心,将pl存储于点集S中,以及一个队列 L集中,第二步,计算pi的四个邻点,并将其加入到S集中,52 =命7,/?2, />4, /?5人 j〇2 =(M+1,jO),j〇3 =(i-1,jO),= (i,jO+1),j〇5 = (i,jO-1);此时S集中有四个点待处理,根据式(1),选出下一个处理点P2:/2=max C(/4), ),其中C即为互相关系数最大,此时用P2将队列L集更新,紧接着搜索当前点 P2的邻点,如此经过连续不断的递进搜索,L集队列中永远是相关系数最大值,直到判断条 件S集不再有更新,则该区域搜索完毕,得到更普遍分布的第二级轴向-横向位移估计值;
[0024] 步骤9)步骤8完成后,利用Savitzky-Golay滤波方法,对得到的横向估计进行平 滑;
[0025] 步骤10)第三级搜索之前,对第二级得到的横向位移估计值利用双线性插值法,得 到第三级的横向位移估计,因为横向位移对噪声敏感,窗口越小,这种影响越大,因此,第三 级不再对横向位移估计做更进一步的精确;
[0026] 步骤11)确定第三级窗口数量及大小,第三级需要完成整幅帧的搜索,因此数量由 窗口大小确定,窗口大小在轴向方向可取四十分之一,第二级获得的轴向位移估计作为第 三级窗口的子集,进行结果传递;
[0027] 步骤12)第三级因为窗口较小,此时默认每个窗口的弹性模量为常数,位移可以近 似看做深度的线性函数,窗口的位移可以假定参考窗口中心点的位移,但由于散斑和背向 散射可能会引起小窗口中心位移估计局部产生巨大变差,因此利用对数压缩方法,对基带 信号进行处理,以此来减小幅值的变化,如式2和式3所示;
[0028] 窃x從软pa " (2)
[0029] 0 ™s〇g^ircx ^arg(p^t{^ (3)
[0030] 其中,c表示压缩因子。步骤13)因为第三级窗口的数量较多,不再利用互相关系 数作为匹配准则,而是改用压缩前、后信号相关性最大时,其相位差趋近于零点的特点,进 行牛顿迭代快速收敛的方法,如式4所示,其中arg(prelcig(x,y))和arg(postlcig(x,y))为步 骤12中得到的压缩前、后对数压缩后的相值,W(x,y)为采样点对位移估计的贡献值,即权 函数,表示前一步得到的位移值,&^41:为当前所求的位移值,w。是换能器的中心频率, 而且因为信号的采样点为离散点,且经过下采样处理,因此实际中很难完全趋近于零点,可 以通过限制迭代次数来达到精确位移估计的目的,同时不会进入死循环,迭代次数可以由 具体的信号采样点数量及所要求的精度最终确定;
[0031]
Figure CN103735287BD00081
[0034] 步骤14)步骤13中的W(x,y),可取1,也可取步骤4得到的压缩前、后信号的幅值 a: (X,y)和a2 (X,y)对应窗口的和,如式5所示;
[0035]
Figure CN103735287BD00082
(5)
[0036] 步骤15)利用步骤8中的跟踪方法结合步骤13中提到的相位差方法,将轴向位移 估计传递到剩下的每一个第三级窗口,完成算法的最后一步搜索,即可得到整幅帧图的位 移估计值。
[0037] 进一步的,所述步骤2中的微小的挤压控制在1%的数量级。
[0038] 本发明的有益效果是:
[0039] 本发明提出的二维混合位移估计方法,结合了数据点的多级搜索和单级跟踪的方 法,体现了合理平衡计算量和计算精度的优点。在算法初级阶段,利用二维多级搜索的方法 来获取轴向和横向偏移量,而避免在高级阶段对横向方向进行搜索。在高级阶段,为最大程 度减小计算量,则利用质量指导和零点相位差跟踪方法来估计轴向位移。这种算法可以很 好地适用于血管内超声成像的实时成像要求,以及成像精度。
附图说明
[0040] 图1是算法总体流程不意图;
[0041] 图2为第一级的搜索策略示意图图;
[0042] 图3为第二级、第三级的搜索策略流程图。
具体实施方式
[0043]下面将参考附图并结合实施例,来详细说明本发明。
[0044]一种血管内超声弹性成像二维混合位移估计方法,包括以下步骤:
[0045]步骤1)利用商用超声换能器得到待测生物组织的压缩前的一帧二维射频信号;
[0046]步骤2)手持该超声换能器沿着探头纵向对该组织施加一个微小的挤压(一般控制 在1%的数量级),得到压缩后的一帧二维射频信号;
[0047]步骤3)分别对从步骤1和步骤2中得到的组织压缩前、后的射频信号进行计算,得 到基带解析数据,具体步骤为先对射频数据进行下采样减少数据点,为避免下采样时组织 位移信息丢失,压缩前、后的信号采样频率应尽可能高,而下采样时,在不违背麦奎斯特采 样定律的前提下,下采样频率也应该尽可能高,再进行希尔伯特变换,再乘以exp( ]Y\其中 w。为超声换能器的中心频率,则得到基带解析信号pre(x,y)和post(X,y),其中x表示沿 探头轴向方向采样点,y表示沿探头横向方向采样点;
[0048] 步骤4)对步骤3中得到的压缩前、后基带解析信号进行信号处理,分别得到压缩 前、后信号的幅值8 1(1,7)和82(1,7);
[0049]步骤5)确定第一级窗口的数量及大小,可取9个,每行每列各3个,(1)压缩前的 目标区域第一级窗口大小为轴向或横向采样点的五分之一;(2)压缩后的搜索区域第一级 窗口大小略大于轴向或横向采样点的五分之一,可取四分之一;
[0050] 步骤6)利用步骤5中确定的参数,对步骤4中得到的压缩前、后信号的幅值 ai(x,y)和&2(1,5〇进行第一级窗口搜索,搜索时对压缩前、后帧进行相关系数最大值匹配, 相关系数公式如式(1)所示,搜索策略可参照图1,( 1)在最低分辨率下,压缩前和压缩后帧 中的每一个窗口中,每四个采样点(间隔3个点)作为一个对象点,图2的上一图中的方形点 为压缩前窗口的中心,对齐到压缩后的搜索区域(图2的下一图)进行相关系数计算,找到最 佳匹配的矩阵形区域(图2的下一图中外圆方形点为其中心点)并做记录;(2)将分辨率提 高一级,每两个采样点作为处理对象点,上一步中找到的外圆方形点周围的x方向(轴向)三 个点作为这一步的匹配候选点,计算相关系数,找到最佳匹配点(图2中的下二图中外圆方 形点),并做记录;(3)分辨率最高的一级,不间隔点进行搜索,将上一步方形点周围的三个 点作为候选点,且计算相关系数时不再间隔任何采样点,找到最佳匹配点(图2中的下三图 中外圆方形点),完成第一级搜索。只有在最低分辨率时,候选匹配位置尺寸与搜索区域的 大小相关;在下一步分辨率提高的情况下,对相关系数的要求高于计算量,候选匹配区域定 为三个,得到9个独立的轴向-横向位移估计值;
[0051]
Figure CN103735287BD00091
[0052]其中分别是横向和轴向位移,%分别是压缩前和压缩后的图像幅值, %,%分别是%的内插窗口平均值,T是窗口大小;
[0053]步骤7)确定第二级的窗口数量及大小,其数量与算法参数无关,其中有九个为第 一级9个窗口的子集,因此第一级搜索得到的轴向和横向位移估计就可以传递到第二级, 第二级其余的窗口则均匀分布整个帧,第二级运动的方法和第一级相同,但间隔因子初始 即取每两个采样点(间隔一个采样点)作为一个对象点,压缩前的目标区域窗口大小为轴向 或横向采样点的十五分之一,而压缩后的搜索区域窗口大小则略大于目标区域,可取轴向 或横向米样点的十分之一;
[0054] 步骤8)重复步骤6中的(1)、(2)、(3)步,但第二级的匹配过程由动态跟随来完 成,首先第二级的9个窗口由第一级的输出确定一个初始点pi(x0,y0),其余的窗口初始值 则由第一级的结果大致定向,定位于每一个窗口中心,将pl存储于点集S中,以及一个队列 L集中,第二步,计算pi的四个邻点,并将其加入到S集中,52 =命7,/?2, />4, /?5人 j〇2 =(M+1,jO),j〇3 =(i-1,jO),= (i,jO+1),j〇5 = (i,jO-1);此时S集中有四个点待处理,根据式(1),选出下一个处理点P2:/2=max C(/4),C〇d5)),其中C即为互相关系数最大,此时用P2将队列L集更新,紧接着搜索当前点P2的邻点,如此经过连续不断的递进搜索,L集队列中永远是相关系数最大值,直到判断条 件S集不再有更新,则该区域搜索完毕,得到更普遍分布的第二级轴向-横向位移估计值, 具体流程可参考图3 ;
[0055] 步骤9)步骤8完成后,利用Savitzky-Golay滤波方法,对得到的横向估计进行平 滑;
[0056] 步骤10)第三级搜索之前,对第二级得到的横向位移估计值利用双线性插值法,得 到第三级的横向位移估计,因为横向位移对噪声敏感,窗口越小,这种影响越大,因此,第三 级不再对横向位移估计做更进一步的精确;
[0057] 步骤11)确定第三级窗口数量及大小,第三级需要完成整幅帧的搜索,因此数量由 窗口大小确定,窗口大小在轴向方向可取四十分之一,第二级获得的轴向位移估计作为第 三级窗口的子集,进行结果传递;
[0058] 步骤12)第三级因为窗口较小,此时默认每个窗口的弹性模量为常数,位移可以近 似看做深度的线性函数,窗口的位移可以假定参考窗口中心点的位移,但由于散斑和背向 散射可能会引起小窗口中心位移估计局部产生巨大变差,因此利用对数压缩方法,对基带 信号进行处理,以此来减小幅值的变化,如式2和式3所示;
[0059] »lilCl -t c x ^re(^y)i)^p〇x arg(pr^^y))) ⑵
[o060]pas!;^(5: ^log(! 4-c^ ^}')m(3)
[0061] 其中,c表不压缩因子。
[0062] 步骤13)因为第三级窗口的数量较多,不再利用互相关系数作为匹配准则,而是改 用压缩前、后信号相关性最大时,其相位差趋近于零点的特点,进行牛顿迭代快速收敛的方 法,如式4所示,其中arg(prelcig(x,y))和arg(postlcig(x,y))为步骤12中得到的压缩前、后 对数压缩后的相值,W(x,y)为采样点对位移估计的贡献值,即权函数表示前一步得到 的位移值,4%0为当前所求的位移值,w。是换能器的中心频率,而且因为信号的采样点为 离散点,且经过下采样处理,因此实际中很难完全趋近于零点,可以通过限制迭代次数来达 到精确位移估计的目的,同时不会进入死循环,迭代次数可以由具体的信号采样点数量及 所要求的精度最终确定;
[0063]
Figure CN103735287BD00101
[0064]
[0065]
Figure CN103735287BD00111
[0066] 步骤14)步骤13中的W(x,y),可取1,也可取步骤4得到的压缩前、后信号的幅值a: (X,y)和a2 (X,y)对应窗口的和,如式5所示;
[0067]
Figure CN103735287BD00112
(5)
[0068] 步骤15)利用步骤8中的跟踪方法结合步骤13中提到的相位差方法,将轴向位移 估计传递到剩下的每一个第三级窗口,完成算法的最后一步搜索,即可得到整幅帧图的位 移估计值。

Claims (2)

1. 一种血管内超声弹性成像二维混合位移估计方法,其特征在于,包括以下步骤: 步骤1)利用商用超声换能器得到待测生物组织的压缩前的一帧二维射频信号; 步骤2)手持所述超声换能器沿着探头纵向对所述生物组织施加一个微小的挤压,得 到压缩后的一帧二维射频信号; 步骤3)分别对从步骤1)和步骤2)中得到的组织压缩前、后的射频信号进行计算,得 到基带解析数据,具体步骤为先对射频数据进行下采样减少数据点,为避免下采样时组织 位移信息丢失,压缩前、后的信号采样频率应尽可能高,而下采样时,在不违背麦奎斯特采 样定律的前提下,下采样频率也应该尽可能高,再进行希尔伯特变换,乘以exp ( ·]Υ\其中w。 为超声换能器的中心频率,则得到基带解析信号pre(x, y)和post (X,y),其中X表示沿探头 轴向方向采样点,y表示沿探头横向方向采样点; 步骤4)对步骤3)中得到的压缩前、后基带解析信号进行信号处理,分别得到压缩前、 后信号的幅值S1 (X,y)和a2(x, y); 步骤5)确定第一级窗口的数量及大小,取9个,每行每列各3个; (1) 压缩前的目标区域第一级窗口大小为轴向或横向采样点的五分之一; (2) 压缩后的搜索区域第一级窗口大小略大于轴向或横向米样点的五分之一; 步骤6)利用步骤5)中确定的参数,对步骤4)中得到的压缩前、后信号的幅值ai (X,y) 和a2 (X,y)进行第一级窗口搜索,搜索时对压缩前、后帧进行相关系数最大值匹配,互相关 系数公式如式①所示; (1) 在最低分辨率下,压缩前和压缩后帧中的每一个窗口中,每四个米样点作为一个对 象点,方形点为压缩前窗口的中心,对齐到压缩后的搜索区域进行相关系数计算,找到最佳 匹配的矩阵形区域并做记录,标记为外圆方形点; (2) 将分辨率提高一级,每两个采样点作为处理对象点,上一步中找到的外圆方形点周 围的X方向三个点作为这一步的匹配候选点,计算相关系数,找到最佳匹配点,标记为外圆 方形点; (3) 分辨率最高的一级,不间隔点进行搜索,将上一步外圆方形点周围的三个点作为候 选点,且计算相关系数时不再间隔任何采样点,找到最佳匹配点,完成第一级搜索,只有在 最低分辨率时,候选匹配位置尺寸与搜索区域的大小相关;在下一步分辨率提高的情况下, 对相关系数的要求高于计算量,候选匹配区域定为三个,得到9个独立的轴向-横向位移估 计倌:
Figure CN103735287BC00021
其中分别是横向和轴向位移,1%分别是压缩前和压缩后的图像幅值, 纹p 分别是◎i'» 的内插窗口平均值,T是窗口大小; 步骤7)确定第二级的窗口数量及大小,其数量与算法参数无关,其中有九个为第一级 9个窗口的子集,因此第一级搜索得到的轴向和横向位移估计就传递到第二级,第二级其余 的窗口则均匀分布整个帧,第二级运动的方法和第一级相同,但间隔因子初始即取每两个 采样点作为一个对象点,压缩前的目标区域窗口大小为轴向或横向采样点的十五分之一, 而压缩后的搜索区域窗口大小则略大于目标区域; 步骤8)重复步骤6)中的(1)、(2)、(3)步,但第二级的匹配过程由动态跟随来完成, 首先第二级的9个窗口由第一级的输出确定一个初始点pi (xO, yO),其余的窗口初始值则 由第一级的结果大致定向,定位于每一个窗口中心,将pi存储于点集S中,以及一个队列 L集中,第二步,计算pi的四个邻点,并将其加入到S集中,吳命7,/¾ />4, /?5人/?2 =UO +1, /0),/?3 = UO -1, /0),/>4 =(成 jO +1),/?5 =(成 -1);此时 S 集 中有四个点待处理,根据式①,选出下一个处理点P2:/2 = max C(/4),(成)),其中C即为求取互相关系数,此时用P2将队列L集更新,紧接着搜索当前点 P2的邻点,如此经过连续不断的递进搜索,L集队列中永远是相关系数最大值,直到判断条 件S集不再有更新,则该区域搜索完毕,得到更普遍分布的第二级轴向-横向位移估计值; 步骤9)步骤8)完成后,利用Savitzky-Golay滤波方法,对得到的横向估计进行平滑; 步骤10)第三级搜索之前,对第二级得到的横向位移估计值利用双线性插值法,得到第 三级的横向位移估计,因为横向位移对噪声敏感,窗口越小,这种影响越大,因此,第三级不 再对横向位移估计做更进一步的精确; 步骤11)确定第三级窗口数量及大小,第三级需要完成整幅帧的搜索,因此数量由窗口 大小确定,窗口大小在轴向方向取四十分之一,第二级获得的轴向位移估计作为第三级窗 口的子集,进行结果传递; 步骤12 )第三级因为窗口较小,此时默认每个窗口的弹性模量为常数,位移近似看做深 度的线性函数,窗口的位移假定参考窗口中心点的位移,但由于散斑和背向散射可能会引 起小窗口中心位移估计局部产生巨大变差,因此利用对数压缩方法,对基带信号进行处理, 以此来减小幅值的变化,如式②和式③所示;
Figure CN103735287BC00031
其中,c表示压缩因子; 步骤13)因为第三级窗口的数量较多,不再利用互相关系数作为匹配准则,而是改用压 缩前、后信号相关性最大时,其相位差趋近于零点的特点,进行牛顿迭代快速收敛的方法, 如式④所示,其中arg(prelcig(x, y))和arg(postlcig(x, y))为步骤12)中得到的压缩前、后对 数压缩后的相值,W(x,y)为采样点对位移估计的贡献值,即权函数,ΔΧ&表示前一步得到 的位移值,当前所求的位移值,w。是换能器的中心频率,而且因为信号的采样点 为离散点,且经过下采样处理,因此实际中很难完全趋近于零点,通过限制迭代次数来达到 精确位移估计的目的,同时不会进入死循环,迭代次数由具体的信号采样点数量及所要求 的精度最终确定;
Figure CN103735287BC00041
步骤14)步骤13)中的W(x,y)取步骤4)得到的压缩前、后信号的幅值&1(1,7)和 a2(x,y)对应窗口的和,如式⑤所示;
Figure CN103735287BC00042
步骤15)利用步骤8)中的动态跟随方法结合步骤13)中提到的相位差方法,将轴向位 移估计传递到剩下的每一个第三级窗口,完成算法的最后一步搜索,即可得到整幅帧图的 位移估计值。
2.根据权利要求1所述的血管内超声弹性成像二维混合位移估计方法,其特征在于: 所述步骤2)中的微小的挤压控制在1%的数量级。
CN201310645210.1A 2013-12-05 2013-12-05 一种血管内超声弹性成像二维多级混合位移估计方法 Active CN103735287B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310645210.1A CN103735287B (zh) 2013-12-05 2013-12-05 一种血管内超声弹性成像二维多级混合位移估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310645210.1A CN103735287B (zh) 2013-12-05 2013-12-05 一种血管内超声弹性成像二维多级混合位移估计方法

Publications (2)

Publication Number Publication Date
CN103735287A CN103735287A (zh) 2014-04-23
CN103735287B true CN103735287B (zh) 2015-11-18

Family

ID=50492413

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310645210.1A Active CN103735287B (zh) 2013-12-05 2013-12-05 一种血管内超声弹性成像二维多级混合位移估计方法

Country Status (1)

Country Link
CN (1) CN103735287B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105748100B (zh) * 2014-12-19 2019-04-16 深圳开立生物医疗科技股份有限公司 准静态超声弹性成像位移计算方法和装置
EP3067091B1 (de) * 2015-03-13 2020-07-29 BIOTRONIK SE & Co. KG Dislokationssensor
CN106651868A (zh) * 2016-08-31 2017-05-10 沈阳东软医疗系统有限公司 一种位移测量的方法及装置
CN107198545B (zh) * 2017-06-06 2020-05-05 苏州国科昂卓医疗科技有限公司 生物组织的弹性位移及应变估计方法、装置
CN108784736B (zh) * 2018-05-23 2020-02-14 成都信息工程大学 一种二维迭代的超声弹性成像应变估计方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1586408A (zh) * 2004-08-20 2005-03-02 清华大学 一种多尺度的生物组织位移估计方法
CN1586411A (zh) * 2004-08-06 2005-03-02 清华大学 一种二维综合互相关的生物组织位移估计方法
CN101569543A (zh) * 2008-04-29 2009-11-04 香港理工大学 弹性成像的二维位移估计方法
CN102048560A (zh) * 2010-12-14 2011-05-11 深圳市蓝韵实业有限公司 一种采用双尺度的生物组织位移估计方法
CN102764141A (zh) * 2012-07-20 2012-11-07 中国科学院深圳先进技术研究院 弹性成像方法和系统及其中的生物组织位移估计方法和系统
CN102824194A (zh) * 2011-06-14 2012-12-19 深圳迈瑞生物医疗电子股份有限公司 一种弹性成像中的位移检测方法及装置
CN102920480A (zh) * 2012-11-26 2013-02-13 重庆理工大学 一种超声弹性成像性能增强方法
CN102973296A (zh) * 2012-11-16 2013-03-20 清华大学 一种血管组织位移估算方法
CN103040488A (zh) * 2012-12-21 2013-04-17 深圳大学 一种实时超声弹性成像位移估计方法和系统

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080144902A1 (en) * 2006-10-25 2008-06-19 Aloka Co., Ltd. Optimal block searching algorithm for tissue displacement estimation in elasticity imaging
US8920323B2 (en) * 2010-10-13 2014-12-30 Wisconsin Alumni Research Foundation Coupled axial and lateral displacement estimation for elasticity imaging

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1586411A (zh) * 2004-08-06 2005-03-02 清华大学 一种二维综合互相关的生物组织位移估计方法
CN1586408A (zh) * 2004-08-20 2005-03-02 清华大学 一种多尺度的生物组织位移估计方法
CN101569543A (zh) * 2008-04-29 2009-11-04 香港理工大学 弹性成像的二维位移估计方法
CN102048560A (zh) * 2010-12-14 2011-05-11 深圳市蓝韵实业有限公司 一种采用双尺度的生物组织位移估计方法
CN102824194A (zh) * 2011-06-14 2012-12-19 深圳迈瑞生物医疗电子股份有限公司 一种弹性成像中的位移检测方法及装置
CN102764141A (zh) * 2012-07-20 2012-11-07 中国科学院深圳先进技术研究院 弹性成像方法和系统及其中的生物组织位移估计方法和系统
CN102973296A (zh) * 2012-11-16 2013-03-20 清华大学 一种血管组织位移估算方法
CN102920480A (zh) * 2012-11-26 2013-02-13 重庆理工大学 一种超声弹性成像性能增强方法
CN103040488A (zh) * 2012-12-21 2013-04-17 深圳大学 一种实时超声弹性成像位移估计方法和系统

Also Published As

Publication number Publication date
CN103735287A (zh) 2014-04-23

Similar Documents

Publication Publication Date Title
CN103735287B (zh) 一种血管内超声弹性成像二维多级混合位移估计方法
KR101868381B1 (ko) 의료용 초음파 이미징에서의 전단파 정보의 해석
CN101066211B (zh) 用于在超声波系统中显示信息的方法
CN101317774B (zh) 超声波诊断装置
US6884216B2 (en) Ultrasound diagnosis apparatus and ultrasound image display method and apparatus
US8538103B2 (en) Medical image processing device, medical image processing method, medical image diagnostic apparatus, operation method of medical image diagnostic apparatus, and medical image display method
US9330461B2 (en) Image-based method for measuring elasticity of biological tissues and system thereof
CN103239258B (zh) 采用超声波的同轴切变波表征
CN102695457B (zh) 超声波诊断装置及测量内中膜的厚度的方法
Rappaport et al. Assessment of myocardial regional strain and strain rate by tissue tracking in B-mode echocardiograms
EP2529666B1 (en) Ultrasonic diagnosis device and method used therefor to track measurement point
Liebgott et al. Transverse oscillations for tissue motion estimation
CN103327904B (zh) 超声波摄像装置和超声波摄像方法
EP2189117B1 (en) Region setting for intima media thickness measurement in an ultrasound system
JP5726081B2 (ja) 超音波診断装置及び弾性画像の分類プログラム
US8300909B2 (en) Ultrasonographic device and ultrasonographic method
CN100383554C (zh) 心壁应变成像
CN105748100B (zh) 准静态超声弹性成像位移计算方法和装置
US20190365354A1 (en) Ultrasonic blood flow parameter displaying method, and ultrasonic imaging system therefor
Nilsson et al. A method to measure shear strain with high spatial resolution in the arterial wall non-invasively in vivo by tracking zero-crossings of B-mode intensity gradients
Tournoux et al. Estimation of radial strain and rotation using a new algorithm based on speckle tracking
Sade et al. Second-generation tissue Doppler with angle-corrected color-coded wall displacement for quantitative assessment of regional left ventricular function
US20130158403A1 (en) Method for Obtaining a Three-Dimensional Velocity Measurement of a Tissue
US8777860B2 (en) Method for evaluation of renal vascular perfusion using power doppler ultrasonography
Park et al. Arterial elasticity imaging: comparison of finite-element analysis models with high-resolution ultrasound speckle tracking

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant