CN111528912A - 一种超声弹性成像方法、装置及*** - Google Patents
一种超声弹性成像方法、装置及*** Download PDFInfo
- Publication number
- CN111528912A CN111528912A CN202010447279.3A CN202010447279A CN111528912A CN 111528912 A CN111528912 A CN 111528912A CN 202010447279 A CN202010447279 A CN 202010447279A CN 111528912 A CN111528912 A CN 111528912A
- Authority
- CN
- China
- Prior art keywords
- data
- frame
- elastography
- matching
- calculating
- 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.)
- Pending
Links
- 238000002091 elastography Methods 0.000 title claims abstract description 85
- 238000000034 method Methods 0.000 title claims abstract description 42
- 238000003384 imaging method Methods 0.000 claims abstract description 61
- 238000012216 screening Methods 0.000 claims abstract description 49
- 238000011156 evaluation Methods 0.000 claims abstract description 21
- 238000007781 pre-processing Methods 0.000 claims abstract description 9
- 238000006073 displacement reaction Methods 0.000 claims description 32
- 238000003825 pressing Methods 0.000 claims description 29
- 230000010363 phase shift Effects 0.000 claims description 20
- 238000005070 sampling Methods 0.000 claims description 18
- 238000002113 ultrasound elastography Methods 0.000 claims description 17
- 239000000523 sample Substances 0.000 claims description 15
- 238000004590 computer program Methods 0.000 claims description 9
- 238000002604 ultrasonography Methods 0.000 claims description 8
- 230000006835 compression Effects 0.000 claims description 5
- 238000007906 compression Methods 0.000 claims description 5
- 238000012285 ultrasound imaging Methods 0.000 claims description 5
- 230000004323 axial length Effects 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 9
- 230000002159 abnormal effect Effects 0.000 abstract description 8
- 238000004364 calculation method Methods 0.000 description 32
- 238000001914 filtration Methods 0.000 description 17
- 230000004927 fusion Effects 0.000 description 9
- 238000004422 calculation algorithm Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000007726 management method Methods 0.000 description 4
- 230000007704 transition Effects 0.000 description 4
- 239000000969 carrier Substances 0.000 description 3
- 238000001514 detection method Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000004393 prognosis Methods 0.000 description 2
- 208000003028 Stuttering Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000003139 buffering effect Effects 0.000 description 1
- 229910010293 ceramic material Inorganic materials 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/48—Diagnostic techniques
- A61B8/485—Diagnostic techniques involving measuring strain or elastic properties
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5215—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5269—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving detection or reduction of artifacts
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Heart & Thoracic Surgery (AREA)
- Surgery (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Veterinary Medicine (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明涉及超声成像技术领域,公开了一种超声弹性成像方法,包括以下步骤:获取弹性成像原始数据,对弹性成像原始数据进行数据预处理得到IQ数据,并对IQ数据进行缓存;从缓存的多帧IQ数据中提取部分IQ数据作为待筛选IQ数据集;对待筛选IQ数据集中IQ数据进行两两组合,得到多个IQ数据组;计算每一IQ数据组的结构相似性评估系数,取结构相似性评估系数在设定阈值范围内的IQ数据组作为初筛IQ数据组;计算初筛IQ数据组中各IQ数据组的成像质量参数,并根据成像质量参数从初筛IQ数据组中筛选出最佳组合;根据最佳组合进行弹性成像,得到弹性图像。本发明提供的超声弹性成像方法具有弹性成像稳定、无异常帧,且资源占用率低的技术效果。
Description
技术领域
本发明涉及超声成像技术领域,具体涉及一种超声弹性成像方法、装置、***以及计算机存储介质。
背景技术
超声弹性成像(ultrasonicelastography)是近年来发展比较快的一种新的超声成像技术,该方法由J.Ophir在1991年提出。超声弹性成像为临床疾病诊断提供了更加丰富的功能信息,因此该方法受到业界和相关研究机构的广泛关注,并体现出了较高的商用价值。
超声按压弹性成像算法,是基于相邻两帧数据之间在相同应力作用下存在位移,组织由于弹性大小的不同从而产生的应变不同来实现的。通过不同颜色反映出组织的弹性大小(应变大小),实现弹性成像。超声按压弹性成像技术,由于需要手动操作,成像质量与操作手法的相关性较大。高质量的成像以按压力的连续性及稳定性为前提,但是在实际操作过程中很难做到按压力度的绝对稳定,更难保证采集到的每帧图像应力是固定不变的。采集到数据中的相邻帧数据之间的存在应力差异性,就可能导致计算出来的弹性图像(应变图像)颜色变化较大。此外,操作过程中应力大小的变化可能导致采集数据的相邻帧之间出现去相关的现象,导致弹性图像异常。这些较大的颜色变化、弹性图像的异常结果在实时显示过程中会导致弹性图像不稳定的现象,影响临床判断。
当前针对弹性成像不稳定的问题,一般的做法是,弹性成像完成后,获取弹性图像的应变均值与匹配参数的融合结果,应变均值与匹配参数融合结果为弹性图像的成像质量参数,根据应变均值与匹配参数的融合结果来选择是否输出当前弹性图像,即:筛除成像异常的弹性图像,输出符合质量要求的弹性图像,从而提升弹性成像质量。然而,根据弹性成像原始数据计算出弹性图像,需要通过弹性成像算法实现,弹性成像算法是一个计算量较大的算法,计算出弹性图像之后,再判断是否启用该弹性图像的方法存在效率低、硬件资源占用损耗大的问题。尤其是当出现较多弃用数据时,容易导致实时弹性图像出现卡顿的假象,弃用数据也是通过弹性成像算法获取的弹性结果,这一过程占用了运算资源,但是得到的却是无用的异常成像结果,一定程度上减缓了弹性成像的帧率。
发明内容
本发明的目的在于克服上述技术不足,提供一种超声弹性成像方法、装置、***以及计算机存储介质,解决现有技术中针对应力不稳定所导致的弹性成像异常问题进行质量提升时,效率低、资源占用率高的技术问题。
为达到上述技术目的,本发明的技术方案提供一种超声弹性成像方法,包括以下步骤:
获取弹性成像原始数据,对所述弹性成像原始数据进行数据预处理得到IQ数据,并对所述IQ数据进行缓存;
从缓存的多帧IQ数据中提取部分IQ数据作为待筛选IQ数据集;对所述待筛选IQ数据集中IQ数据进行两两组合,得到多个IQ数据组;
计算每一所述IQ数据组的结构相似性评估系数,取结构相似性评估系数在设定阈值范围内的IQ数据组作为初筛IQ数据组;
计算所述初筛IQ数据组中各IQ数据组的成像质量参数,并根据所述成像质量参数从所述初筛IQ数据组中筛选出最佳组合;
根据所述最佳组合进行弹性成像,得到弹性图像。
本发明还提供一种超声弹性成像装置,包括处理器以及存储器,所述存储器上存储有计算机程序,所述计算机程序被所述处理器执行时,实现所述超声弹性成像方法。
本发明还提供一种超声弹性成像***,包括所述超声弹性成像装置,还包括超声探头、B超成像装置以及显示装置;
所述超声探头用于分时探测所述超声弹性原始数据以及B超原始数据;
所述B超成像装置用于根据所述B超原始数据进行B超成像,得到B超图像;
所述显示装置用于对所述B超图像以及所述弹性图像进行显示。
本发明还提供一种计算机存储介质,其上存储有计算机程序,所述计算机该程序被处理器执行时,实现所述超声弹性成像方法。
与现有技术相比,本发明的有益效果包括:本发明首先通过弹性成像原始数据获取IQ数据;然后,并未直接通过IQ数据进行弹性成像计算,而是先对IQ数据进行组合筛选,筛选出结构相似性评估系数较大的初筛IQ数据组,筛除掉具有显著差异的IQ数据组,保证了采集数据的稳定性,进而筛选出成像质量较好的最佳组合,对弹性成像效果进行预后判断,保证筛选出的最佳组合的成像效果;最后基于最佳组合进行弹性成像计算,得到稳定的弹性图像。本发明。本发明在计算弹性图像之前,对采集的数据进行结构相似性分析判断,保证了输出的弹性图像稳定显示,减少操作手法差异或不当造成弹性图像显著差异的现象,不做冗余计算,节省运算资源;在计算弹性成像前确定用于计算弹性成像的有效帧数据的最佳组合,保证了输出的弹性图像无异常帧,提升弹性成像帧率,不做冗余计算,节省运算资源。
附图说明
图1是本发明提供的超声弹性成像方法一实施方式的流程图;
图2是本发明提供的生成实时按压曲线一实施方式的取样块取样位置图;
图3是本发明提供的超声弹性成像***一实施方式的结构示意图;
附图标记:
101、B超成像装置;102、超声弹性成像装置;103、显示装置;104、超声探头。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
实施例1
如图1所示,本发明的实施例1提供了超声弹性成像方法,包括以下步骤:
S1、获取弹性成像原始数据,对所述弹性成像原始数据进行数据预处理得到IQ数据,并对所述IQ数据进行缓存;
S2、从缓存的多帧IQ数据中提取部分IQ数据作为待筛选IQ数据集;对所述待筛选IQ数据集中IQ数据进行两两组合,得到多个IQ数据组;
S3、计算每一所述IQ数据组的结构相似性评估系数,取结构相似性评估系数在设定阈值范围内的IQ数据组作为初筛IQ数据组;
S4、计算所述初筛IQ数据组中各IQ数据组的成像质量参数,并根据所述成像质量参数从所述初筛IQ数据组中筛选出最佳组合;
S5、根据所述最佳组合进行弹性成像,得到弹性图像。
本实施例首先通过弹性成像原始数据获取IQ数据;具体的,IQ数据指对弹性成像原始数据进行IQ调制后得到的数据,将弹性成像原始数据分为两路:I路和Q路,分别进行载波调制,两路载波相互正交,相互正交的两路载波即相位相差90度的载波,I路和Q路分别进行载波调制后一起发射,得到IQ数据,有利于提高频谱利用率。然后,在计算弹性图像之前,先对缓存的一定量的IQ数据进行组合筛选,筛选出结构相似性评估系数较大的初筛IQ数据组,筛除掉具有显著差异的IQ数据组,保证用于计算弹性图像的IQ数据组存在较好应变位移结果,保证了采集数据的稳定性;在初次筛选的基础上,进行进一步的筛选,筛选出成像质量较好的最佳组合,对弹性成像效果进行预后判断,保证筛选出成像效果最佳的最佳组合;最后基于最佳组合进行弹性成像计算,得到稳定、高质量的弹性图像。
本实施例提供的超声弹性成像方法,避免了计算出异常的弹性图像结果,节省了计算异常结果的运算时间;得到的弹性图像之间过渡平稳,差异性小,减少操作手法差异或不当造成弹性图像显著差异的现象;同时,不做冗余计算,节省运算资源。
具体的,在进行弹性成像之前,首先需要对弹性成像原始数据进行信号处理以及缓存,具体阐述如下。
优选的,对所述弹性成像原始数据进行数据预处理得到IQ数据,并对所述IQ数据进行缓存,具体为:
对所述弹性成像原始数据依次进行带通滤波、正交解调、低通滤波、二次采样以及低通滤波,得到IQ数据;
实时缓存最新的第一设定帧数的IQ数据,得到缓存IQ数据集。
弹性成像原始数据经过带通滤波、正交解调、低通滤波、二次采样以及低通滤波获取到连续的IQ数据,本实施例中IQ数据的大小为384*800,其中384代表线束,800代表深度,成像深度为5cm。探头中心频率为7.5MHz。
本实施例中第一设定帧数设置为100帧,实时缓存最新的100帧IQ数据,包括I数据m_pSrc_I_data_int[100]以及Q数据m_pSrc_Q_data_int[100];使用帧号管理方法对最新的100帧IQ数据在第一缓冲区中的存储地址进行管理。本实施例中:最新帧在第一缓冲区中的地址为m_nFrameNO[0],m_nFrameNO[0]的值85,则对应的最新帧IQ数据保存在m_pSrc_I_data_int[85]与m_pSrc_Q_data_int[85]中。IQ数据在第一缓冲区内的存储方式为递增存储,即:成像刚启动时第一帧存放在m_pSrc_I_data_int[0]、m_pSrc_Q_data_int[0]中,第二帧存放在m_pSrc_I_data_int[1]、m_pSrc_Q_data_int[1]中,…,第100帧存放在m_pSrc_I_data_int[99]、m_pSrc_Q_data_int[99]帧中;随后,第101帧替换了第1帧存放在m_pSrc_I_data_int[0]、m_pSrc_Q_data_int[0]中,依次循环存储。
对IQ数据进行缓存后,即可进行组合筛选以及弹性成像了,然而,本实施例在进行弹性成像的同时,利用缓存的IQ数据生成实时按压曲线,实现了弹性图像和实时按压曲线的同时获取。具体阐述如下。
优选的,采用并行线程,在根据缓存的IQ数据进行弹性成像同时,根据缓存的IQ数据生成实时按压曲线。
本实施例在进行弹性成像计算的同时,生成实时按压曲线,实时按压曲线的计算与弹性成像计算采用并行线程实现。通过并行线程实现实时按压曲线与弹性成像的数据筛选以及成像计算的并行工作,保证弹性图像连续、稳定的同时也保证了按压曲线的实时性。实时按压曲线的计算生成具体阐述如下。
优选的,根据缓存的IQ数据生成实时按压曲线,具体为:
从缓存的IQ数据中提取最新帧IQ数据和次新帧IQ数据;
设定多个同一深度、同一大小的取样块;
通过所述取样块分别从所述最新帧IQ数据块和所述次新帧IQ数据块中取样多对相应的数据块,得到多组数据块组;
计算每一组所述数据块组的匹配系数;
根据所述匹配系数提取最新帧IQ数据与次新帧IQ数据之间的相对相移;
根据所述相对相移获取最新帧IQ数据与次新帧IQ数据之间的相对位移;
对各组数据块组所对应的相对位移进行加权求和,得到实时相对位移;
根据所述实时相对位移生成所述实时按压曲线。
本实施例在计算实时按压曲线时,首先判断第一缓冲区是否更新数据,如果是,则实时从第一缓冲区中提取最新帧与次新帧,m_nFrameNO[0]与m_nFrameNO[1]的值为最新帧与次新帧在第一缓冲区的地址。通过固定的取样块提取最新帧与次新帧的匹配系数,通过匹配系数提取最新帧与次新帧之间的相移,通过相移获取相对位移,相对位移反映应力的大小与方向,这里将相对位移存储到按压曲线存储变量m_pDdisplace_line[100]的m_pDdisplace[0]中,按压曲线存储变量中的对应位移依次后移一个值,生成实时按压曲线。
对最新帧与次新帧中相应数据块组的匹配系数求取,具体如下。
优选的,计算每一组所述数据块组的匹配系数,具体为:
采用相关窗在所述数据块组中进行互相关匹配,得到互相关系数:
R12=∫∫DS1(t+v,x+w)*S2(t+v,x+w)*dvdw
其中,R12为所述数据块组的互相关系数,t为待匹配的数据块的中心的轴向坐标,x为待匹配的数据块的中心的横向坐标,v为相关窗的横向长度,w为相关窗的轴向长度,S1(t+v,x+w)为数据块组中第一个数据块在(t+v,x+w)像素点处的IQ数据值,S2(t+v,x+w)为数据块组中第二个数据块在(t+v,x+w)像素点处的IQ数据值,D=v*w;
根据所述匹配系数提取最新帧IQ数据与次新帧IQ数据之间的相对相移,具体为:
其中,为最新帧IQ数据与次新帧IQ数据的相对相移,I1(i,j)为第一个数据块的正交解调数据的虚部,Q1(i,j)为第一个数据块的正交解调数据的实部,I2(i,j)为第二个数据块的正交解调数据的虚部,Q2(i,j)为第二个数据块的正交解调数据的实部;
根据所述相对相移获取最新帧IQ数据与次新帧IQ数据之间的相对位移:
其中,D为最新帧IQ数据与次新帧IQ数据的相对位移,Δt为相移需要的时间,c为超声波在组织中的传播速度,ω为超声波的中心频率对应的角速度,f0为中心频率。
具体的,本实施中使用三个同一深度、同大小的取样,取样块的位置如图2所示,,图中分别以P1、P2、P3表示,图中标识的尺寸的单位为像素点。分别计算着三个取样块代表的匹配点的匹配系数,然后根据匹配系数提取相移,然后相移提取对应的位移,最后对这三个匹配点对应的位移进行加权求和,对应的权值可以但不限于0.3*D1+0.4*D2+0.3*D3,其中D1、D2、D3分别对应三个取样块获取的位移。
本实施例在进行上述实时按压曲线计算的同时,进行弹性成像计算,弹性成像之前,需要对缓存的IQ数据进行组合筛选,以保证成像质量,具体阐述如下。
优选的,从缓存的多帧IQ数据中提取部分IQ数据作为待筛选IQ数据集,对所述待筛选IQ数据集中IQ数据进行两两组合,得到多个IQ数据组,具体为:
从所述缓存IQ数据集中提出取第二设定帧数的IQ数据作为待筛选IQ数据集;
IQ数据缓存完成后,就需要对缓存IQ数据集中IQ数据进行组合筛选。首先从缓存IQ数据集中提取第二设定帧数的IQ数据,本实施例中第二设定帧数为5,因此可以提取出10组IQ数据组。
IQ数据组合完成后,需要对IQ数据组进行初次筛选,分别对这10个组合数据进行结构相似性评估,筛选出结构相似性评估系数在设定阈值范围之间的IQ数据组,具体阐述如下。
优选的,计算每一所述IQ数据组的结构相似性评估系数,具体为:
计算所述IQ数据组中每一帧IQ数据的包络数据:
其中,data为IQ数据的包络数据,I为IQ数据的正交解调数据的虚部,Q为IQ数据的正交解调数据的实部;
根据所述包络数据计算所述IQ数据组中两帧IQ数据之间的结构相似性评估系数:
其中,x和y分别为同一组IQ数据组中的两帧IQ数据的包络数据,SSIM(x,y)为同一组IQ数据组中的两帧IQ数据之间的结构相似性评估系数,Ux为包络数据x的均值,Uy为包络数据y的均值,σx 2为包络数据x的方差,σy 2为包络数据y的方差,σxy为包络数据x和为包络数据y的协方差,C1、C2均为常数。
本实施例中,C1为6.5025,C2为58.5225。本实施例所求取的SSIM(x,y)的范围为0~1,当x与y相同时,该系数为1,x与y的差异性越大,SSIM(x,y)越小。本实施例中设定阈值范围为0.8~0.95,即筛选出结构相似性评估系数在0.8~0.95范围内的IQ数据组进行下一筛选环节。
初次筛选完成后,在初次筛选的基础上进行二次筛选,具体阐述如下。
优选的,计算所述初筛IQ数据组中各IQ数据组的成像质量参数,并根据所述成像质量参数从所述初筛IQ数据组中筛选出最佳组合,具体为:
所述成像质量参数包括匹配系数以及应变均值,计算所述初筛IQ数据组中各IQ数据组的匹配系数以及应变均值;
筛选出匹配系数大于相关阈值的IQ数据组,然后从筛选出的IQ数据组中选取应变均值最大的设定数量的IQ数据组作为所述最佳组合。
二次筛选包括两个筛选环节:互相关匹配环节以及应变均值匹配环节,对初筛IQ数据组中每一IQ数据组进行取样检测,得到对应的应变均值与匹配系数,取匹配系数大于0.8的IQ数据组中应变均值最大的前N3个IQ数据组作为最佳组合,进行弹性成像计算。下面分别对两个筛选环节进行具体阐述。
优选的,计算所述初筛IQ数据组中各IQ数据组的匹配系数,具体为:
将所述IQ数据组中两帧IQ数据分别设置为基准帧和匹配帧;
将所述基准帧均分为设定数量的基准数据块,在所述匹配帧中划分出与所述数据块位置一一对应的检索框;
计算所述基准帧中每一基准数据块与相应检索框中多个匹配数据块的互相关系数:
RA12(t,x,n,m)=∫∫DS1(t+v,x+w)*S2(t+n+v,x+m+w)*dvdw
其中,RA12(t,x,n,m)为互相关系数,v为相关窗的横向长度,w为相关窗的轴向长度,检索框中多个匹配数据块按设定步长依次移动选取,n为横向步进长度,m为轴向步进长度,t为匹配数据块的中心的轴向坐标,x匹配数据块的中心的横向坐标,S1(t+v,x+w)为基准数据块在(t+v,x+w)像素点处的IQ数据值,S2(t+n+v,x+m+w)为匹配数据块在(t+n+v,x+m+w)像素点处的IQ数据值,D=v*w;
计算基准帧的自相关系数:
RA11(t,x,0,0)=∫∫D|S1(t+v,x+w)|2dvdw
其中,RA11(t,x,0,0)为基准帧的自相关系数;
计算匹配帧的自相关系数:
RA22(t,x,n,m)=∫∫D|S2(t+n+v,x+m+w)|2dvdw
其中,RA22(t,x,n,m)为匹配帧的自相关系数;
根据每一组进行匹配的数据块对应的所述互相关系数、所述基准帧的自相关系数、所述匹配帧的自相关系数计算所述基准帧与所述匹配帧的匹配系数:
其中,C(t,x,n,m)为所述基准帧与所述匹配帧的匹配系数;
计算多组数据块对应的匹配系数的平均值作为最终的匹配系数。
互相关匹配环节,对每一个IQ数据组,在匹配帧中使用互相关匹配的方法匹配最佳帧,匹配系数最大的位置对应最佳帧;本实施例中,为了简化运算将IQ数据组划分为4*12个数据块进行匹配计算;4*12个数据块均匀的分布在整帧IQ数据中,其中4为横向分布,12为轴向分布;在匹配帧中检索匹配基准帧中的4*12个数据块相应位置的检索框;基准帧中的块P的大小为5*73,其中5为横向大小,73为轴向大小;在匹配帧中同一位置的检索框S的大小为9*81,其中9为横向大小,91为轴向大小;匹配步进为1个像素点,数据块P与检索框S的大小根据需求设定;基准数据块与检索框中匹配数据块的匹配系数,采用上述公式计算即可。
匹配系数的最大值对应的匹配数据块即为最佳匹配数据块,最佳匹配数据块的中心点即为最佳匹配点,将最佳匹配数据块的n、m转换为以基准数据块中心位置为原点的位移坐标。
匹配结束后,针对每一基准数据块的最佳匹配数据块获取到4*12个最佳匹配点,从而获取到的4*12个匹配系数,最后对4*12个匹配系数求均值,筛选出匹配系数大于0.8的IQ数据组进行下一步应变均值检测。
为了进一步简化运算,应变均值计算同样基于将IQ数据组均分为多个数据块进行,同时,由于匹配系数计算时已经进行了数据块划分,因此应变均值计算时,采用同样的划分方式,同样将IQ数据组均分为4*12个数据块提取对应的应变均值,避免重复划分增加运算量。应变均值的计算具体阐述如下。
优选的,计算所述初筛IQ数据组中各IQ数据组的应变均值,具体为:
将所述IQ数据组中两帧IQ数据分别均分为设定数量的数据块,得到多组数据块组;
计算每一组所述数据块组的匹配系数;
根据所述匹配系数提取所述IQ数据组中两帧IQ数据之间的相对相移;
根据所述相对相移获取所述IQ数据组中两帧IQ数据之间的相对位移;
计算所述相对位移的轴向梯度,得到应变分布值;
计算各数据块所对应的应变分布值的平均值得到所述应变均值。
具体的,本实施例中,对IQ数据组中数据块组进行相对位移进行计算时,采用计算实时按压曲线时,最新帧与次新帧的数据块组的相对位移计算方法实现。计算得到相对位移后,使用Sobel梯度算子计算其轴向梯度,得到应变分布值。本实施例中梯度算子为:[1 0-1;2 0 -2;1 0 -1]。
针对每一数据块组均可得到一个应变分布值,对同一IQ数据组的多个数据块的应变分布值求均值,得到应变均值。每个IQ数据组都可以获取到一个应变均值,最后对IQ数据组按应变均值从大到小的顺序进行排序,取排在前5个的应变均值最大的IQ数据组作为最佳组合,进行弹性计算。若进入这环节的IQ数据组数量不足5帧,则全部输出,或返回为空。
获取到最佳组合后,即可根据最佳组合进行弹性成像,具体阐述如下。
优选的,根据所述最佳组合进行弹性成像,得到弹性图像,具体为:
计算所述最佳组合的应变分布,将应变分布映射到对应的颜色表,得到弹性图像。
对最佳组合进行计算得到弹性图像时,需要计算应变分布值。本实施例中最佳组合的应变分布值计算,采用二次筛选时,计算初筛IQ数据组中IQ数据组的应变分布值的方法进行计算。同样是将最佳组合划分为多个的网格块(即数据块)进行计算,首先获取网格块的相移分布,然后获取位移分布,最后通过梯度计算模块获取到应变分布,将应变分布结果映射到对应的颜色表,输出为弹性图像。过程类似,在此也就不再重复描述。
求取相移分布前,优选对最佳组合进行高斯滤波,具体的,本实施例中对最佳组合中两帧IQ数据进行了一个窗为3*5的高斯滤波,然后使用窗大小为5*78的块在9*81的块中共进行互相关匹配后,依次得到获取相移分布、位移分布、应变分布。提取相移分布后,同样优选进行滤波,本实施例设置一个大小为5*5的高斯滤波窗,一个5*3的均值滤波窗。获取应变分布后,同样优选进行滤波处理,本实施例设置一个5*5的均值滤波窗以及与图像输出结果与平滑程度相关的高斯滤波窗,大小为5*9,滤波系数为1.5,最后得到弹性图像。
为了实现输出的弹性图像的良好过滤,本实施例还对弹性图像进行余晖处理,具体如下。
优选的,根据所述最佳组合进行弹性成像,得到弹性图像,还包括:
对当前弹性图像与历史弹性图像进行加权融合,得到最终的弹性成像。
将当前弹性图像与历史弹性图像进行加权融合得到最终的弹性成像结果,完成对弹性图像的实时过渡,输出过渡平缓的实时弹性图像。具体的,加权融合的权值根据用户选择的档位设置。
实施例2
本发明的实施例2提供了超声弹性成像装置,包括处理器以及存储器,所述存储器上存储有计算机程序,所述计算机程序被所述处理器执行时,实现实施例1提供的超声弹性成像方法。
本发明实施例提供的超声弹性成像装置,用于实现基于超声弹性成像方法,因此,超声弹性成像方法所具备的技术效果,超声弹性成像装置同样具备,在此不再赘述。
实施例3
如图3所示,本发明的实施例3提供了超声弹性成像***,以下简称本***,包括实施例2提供的超声弹性成像装置102,还包括超声探头104、B超成像装置101以及显示装置103;
所述超声探头104用于分时探测所述超声弹性原始数据以及B超原始数据;
所述B超成像装置101用于根据所述B超原始数据进行B超成像,得到B超图像;
所述显示装置103用于对所述B超图像以及所述弹性图像进行显示。
由于按压弹性成像,需要B超图像进行引导,因此本***对B超图像和弹性图像同时成像,B超成像与弹性成像的原始数据是共用一个超声探头104在不同模式下采集的。由于B超成像与弹性成像需要的发射模式与接收模式中需要调节的参数不同,本***交替控制超声探头104工作在B超成像模式与弹性成像模式两个模式下,本实施例中获取B超原始数据与超声弹性原始数据的帧数比为1:1,采集到超声弹性原始数据以及B超原始数据进入两个并行装置,即B超成像装置101和超声弹性成像装置102,最后进入显示装置103进行显示,具体阐述如下。
超声探头104内置10:前端控制装置。控制超声探头104的发射与接收,控制超声探头104工作在B超原始数据与弹性成像原始数据之间逐帧切换,弹性成像与B超原始数据的帧数为1:1,使用同一探头交替采集的。
B超成像装置101用于将前端获取的B超原始数据处理为B超图像。B超成像装置101采用现有技术实现即可,这里不做过多赘述。
超声弹性成像装置102,包含信号预处理模块、实时按压曲线生成模块、第一缓存模块、第二缓存模块、组合筛选模块、弹性成像模块以及弹性余晖模块共七个模块。
信号预处理模块用于对弹性成像原始数据进行预处理,包括带通滤波、正交解调、低通滤波、二次采样、低通滤波一系列模块,获取到可用于弹性成像计算的IQ数据。
第一缓存模块为第一缓冲区的缓存模块,用于缓存信号预处理模块输出的IQ数据,以便用于后续的实时按压曲线计算、组合筛选、弹性成像计算。该模块用于缓存第一设定帧数N1的IQ数据,使用帧号管理模块管理存入数据,通过帧号管理指示数据流入的顺序,用于循环存储输入的数据,存储满时最新帧数据替换第一缓冲区里时间最久数据的位置,缓存区内永远只保存当前帧数据与历史N1-1帧。这里取第一缓冲区的大小为100帧IQ数据,按最新数据与历史数据的顺序(例如:0-99)将对应的数据存放地址存储在帧号管理(a[100])中,其中a[0]始终存放最新帧数据的地址,a[1]、a[2]、…、a[99]依次存放之前输入数据。
按压曲线生成模块用于实时的根据第一缓存模块中最新帧数据与次新帧数据提取相对位移,并提取最新帧的按压位移,将获取的按压位移值保存到按压曲线缓存变量中。
第二缓存模块通过判断第一缓存模块是否更新,如果更新,从第一缓存模块更新数据中提取第二设定帧数N2的IQ数据,对这N2帧IQ数据进行两两组合;
组合筛选模块为本发明的关键技术点,在组合筛选模块对多个IQ数据组进行结构相似性评估(SSIM),取结构相似性评估系数在设定阈值范围内的IQ数据组进入下一环节,在该环节可以剔除由于按压力度过小造成信噪比不够的IQ数据组,也可以剔除由于按压力度过大导致出现去相关的IQ数据组;然后对IQ数据组进行取样检测,每一帧均匀取样4*12个数据块,其中4为横向取样,12为轴向取样;分别计算每个数据块组的应变均值以及匹配系数;与已有方法不同的是,这里取匹配系数大于0.85的组合;从相关匹配系数大于0.85的组合里,依次取设定数量N3个应变均值最大的IQ数据组,作为最佳组合输出给弹性成像模块;N3可以为N2或小于N2大于0的自然整数;
弹性成像模块为弹性图像计算装置,将组合筛选模块输出的最佳组合进行计算获取最佳的弹性图像。
弹性余晖模块完成对弹性成像模块输出的弹性图像的缓存以及一定比列的加权融合,加权融合的参数以及参与的帧数由余晖档位设置;余晖处理可以实现弹性图像的良好过渡,不会给用户造成突兀的影像现象。最后输出弹性图像到显示装置103。
显示装置103,用于对B超图像、弹性图像以及实时按压曲线进行融合显示。
实施例4
本发明的实施例4提供了计算机存储介质,其上存储有计算机程序,所述计算机程序被处理器执行时,实现实施例1提供的超声弹性成像方法。
本发明实施例提供的计算机存储介质,用于超声弹性成像方法,因此,超声弹性成像方法所具备的技术效果,计算机存储介质同样具备,在此不再赘述。
以上所述本发明的具体实施方式,并不构成对本发明保护范围的限定。任何根据本发明的技术构思所做出的各种其他相应的改变与变形,均应包含在本发明权利要求的保护范围内。
Claims (10)
1.一种超声弹性成像方法,其特征在于,包括以下步骤:
获取弹性成像原始数据,对所述弹性成像原始数据进行数据预处理得到IQ数据,并对所述IQ数据进行缓存;
从缓存的多帧IQ数据中提取部分IQ数据作为待筛选IQ数据集;对所述待筛选IQ数据集中IQ数据进行两两组合,得到多个IQ数据组;
计算每一所述IQ数据组的结构相似性评估系数,取结构相似性评估系数在设定阈值范围内的IQ数据组作为初筛IQ数据组;
计算所述初筛IQ数据组中各IQ数据组的成像质量参数,并根据所述成像质量参数从所述初筛IQ数据组中筛选出最佳组合;
根据所述最佳组合进行弹性成像,得到弹性图像。
2.根据权利要求1所述的超声弹性成像方法,其特征在于,还包括:采用并行线程,在根据缓存的IQ数据进行弹性成像同时,根据缓存的IQ数据生成实时按压曲线。
3.根据权利要求2所述的超声弹性成像方法,其特征在于,根据缓存的IQ数据生成实时按压曲线,具体为:
从缓存的IQ数据中提取最新帧IQ数据和次新帧IQ数据;
设定多个同一深度、同一大小的取样块;
通过所述取样块分别从所述最新帧IQ数据块和所述次新帧IQ数据块中取样多对相应的数据块,得到多组数据块组;
计算每一组所述数据块组的匹配系数;
根据所述匹配系数提取最新帧IQ数据与次新帧IQ数据之间的相对相移;
根据所述相对相移获取最新帧IQ数据与次新帧IQ数据之间的相对位移;
对各组数据块组所对应的相对位移进行加权求和,得到实时相对位移;
根据所述实时相对位移生成所述实时按压曲线。
4.根据权利要求1所述的超声弹性成像方法,其特征在于,计算每一所述IQ数据组的结构相似性评估系数,具体为:
计算所述IQ数据组中每一帧IQ数据的包络数据:
其中,data为IQ数据的包络数据,I为IQ数据的正交解调数据的虚部,Q为IQ数据的正交解调数据的实部;
根据所述包络数据计算所述IQ数据组中两帧IQ数据之间的结构相似性评估系数:
其中,x和y分别为同一组IQ数据组中的两帧IQ数据的包络数据,SSIM(x,y)为同一组IQ数据组中的两帧IQ数据之间的结构相似性评估系数,Ux为包络数据x的均值,Uy为包络数据y的均值,σx 2为包络数据x的方差,σy 2为包络数据y的方差,σxy为包络数据x和为包络数据y的协方差,C1、C2均为常数。
5.根据权利要求1所述的超声弹性成像方法,其特征在于,计算所述初筛IQ数据组中各IQ数据组的成像质量参数,并根据所述成像质量参数从所述初筛IQ数据组中筛选出最佳组合,具体为:
所述成像质量参数包括匹配系数以及应变均值,计算所述初筛IQ数据组中各IQ数据组的匹配系数以及应变均值;
筛选出匹配系数大于匹配阈值的IQ数据组,然后从匹配系数大于匹配阈值的IQ数据组中选取应变均值最大的设定数量的IQ数据组作为所述最佳组合。
6.根据权利要求5所述的超声弹性成像方法,其特征在于,计算所述初筛IQ数据组中各IQ数据组的匹配系数,具体为:
将所述IQ数据组中两帧IQ数据分别设置为基准帧和匹配帧;
将所述基准帧均分为设定数量的基准数据块,在所述匹配帧中划分出与所述数据块位置一一对应的检索框;
计算所述基准帧中每一基准数据块与相应检索框中多个匹配数据块的互相关系数:
RA12(t,x,n,m)=∫∫DS1(t+v,x+w)*S2(t+n+v,x+m+w)*dvdw
其中,RA12(t,x,n,m)为互相关系数,v为相关窗的横向长度,w为相关窗的轴向长度,检索框中多个匹配数据块按设定步长依次选取,n为横向步进长度,m为轴向步进长度,t为匹配数据块的中心的轴向坐标,x匹配数据块的中心的横向坐标,S1(t+v,x+w)为基准数据块在(t+v,x+w)像素点处的IQ数据值,S2(t+n+v,x+m+w)为匹配数据块在(t+n+v,x+m+w)像素点处的IQ数据值,D=v*w;
计算基准帧的自相关系数:
RA11(t,x,0,0)=∫∫D|S1(t+v,x+w)|2dvdw
其中,RA11(t,x,0,0)为基准帧的自相关系数;
计算匹配帧的自相关系数:
RA22(t,x,n,m)=∫∫D|S2(t+n+v,x+m+w)|2dvdw
其中,RA22(t,x,n,m)为匹配帧的自相关系数;
根据每一组进行匹配的数据块对应的所述互相关系数、所述基准帧的自相关系数、所述匹配帧的自相关系数计算所述基准帧与所述匹配帧的匹配系数:
其中,C(t,x,n,m)为所述基准帧与所述匹配帧的匹配系数;
计算多组数据块对应的匹配系数的平均值作为最终的匹配系数。
7.根据权利要求5所述的超声弹性成像方法,其特征在于,计算所述初筛IQ数据组中各IQ数据组的应变均值,具体为:
将所述IQ数据组中两帧IQ数据分别均分为设定数量的数据块,得到多组数据块组;
计算每一组所述数据块组的匹配系数;
根据所述匹配系数提取所述IQ数据组中两帧IQ数据之间的相对相移;
根据所述相对相移获取所述IQ数据组中两帧IQ数据之间的相对位移;
计算所述相对位移的轴向梯度,得到应变分布值;
计算各数据块所对应的应变分布值的平均值得到所述应变均值。
8.一种超声弹性成像装置,其特征在于,包括处理器以及存储器,所述存储器上存储有计算机程序,所述计算机程序被所述处理器执行时,实现如权利要求1-7任一所述的超声弹性成像方法。
9.根据权利要求1所述的超声弹性成像***,其特征在于,包括如权利要求8所述的超声弹性成像装置,还包括超声探头、B超成像装置以及显示装置;
所述超声探头用于分时探测所述超声弹性原始数据以及B超原始数据;
所述B超成像装置用于根据所述B超原始数据进行B超成像,得到B超图像;
所述显示装置用于对所述B超图像以及所述弹性图像进行显示。
10.一种计算机存储介质,其上存储有计算机程序,其特征在于,所述计算机该程序被处理器执行时,实现如权利要求1-7任一所述的超声弹性成像方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010447279.3A CN111528912A (zh) | 2020-05-25 | 2020-05-25 | 一种超声弹性成像方法、装置及*** |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010447279.3A CN111528912A (zh) | 2020-05-25 | 2020-05-25 | 一种超声弹性成像方法、装置及*** |
Publications (1)
Publication Number | Publication Date |
---|---|
CN111528912A true CN111528912A (zh) | 2020-08-14 |
Family
ID=71969537
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010447279.3A Pending CN111528912A (zh) | 2020-05-25 | 2020-05-25 | 一种超声弹性成像方法、装置及*** |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111528912A (zh) |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103845081A (zh) * | 2012-11-28 | 2014-06-11 | 深圳迈瑞生物医疗电子股份有限公司 | 超声弹性成像***和方法、实时动态帧间处理方法 |
CN104739442A (zh) * | 2013-12-25 | 2015-07-01 | 深圳迈瑞生物医疗电子股份有限公司 | 压力弹性成像位移检测方法、装置和超声成像设备 |
CN109745073A (zh) * | 2019-01-10 | 2019-05-14 | 武汉中旗生物医疗电子有限公司 | 弹性成像位移的二维匹配方法及设备 |
-
2020
- 2020-05-25 CN CN202010447279.3A patent/CN111528912A/zh active Pending
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103845081A (zh) * | 2012-11-28 | 2014-06-11 | 深圳迈瑞生物医疗电子股份有限公司 | 超声弹性成像***和方法、实时动态帧间处理方法 |
CN104739442A (zh) * | 2013-12-25 | 2015-07-01 | 深圳迈瑞生物医疗电子股份有限公司 | 压力弹性成像位移检测方法、装置和超声成像设备 |
CN109745073A (zh) * | 2019-01-10 | 2019-05-14 | 武汉中旗生物医疗电子有限公司 | 弹性成像位移的二维匹配方法及设备 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
EP4343675A1 (en) | Visual image enhancement generation method and system, device, and storage medium | |
US20210236095A1 (en) | Ultrasound ct image reconstruction method and system based on ray theory | |
US20210407371A1 (en) | Method, apparatus and device for acquiring luminance compensation data | |
CN1134808A (zh) | 用于超声波诊断的显示图像的方法和装置 | |
CN109191416A (zh) | 基于稀疏字典学习和剪切波的图像融合方法 | |
CN109557540A (zh) | 基于目标散射系数非负约束的全变差正则化关联成像方法 | |
CN111208509A (zh) | 一种超宽带雷达人体目标姿态可视化增强方法 | |
CN116385436B (zh) | 一种基于ct图像的胆结石辅助检测*** | |
CN112529870A (zh) | 基于源域与频域结合的多尺度CNNs肺结节假阳性剔除方法 | |
CN111862038A (zh) | 一种斑块检测方法、装置、设备及介质 | |
CN109982090A (zh) | 一种结合灰度熵和盲反卷积的采样率自适应分块压缩感知方法 | |
CN104739442B (zh) | 压力弹性成像位移检测方法、装置和超声成像设备 | |
CN113095358A (zh) | 一种图像融合方法及*** | |
Wang et al. | Lightweight multiple scale-patch dehazing network for real-world hazy image | |
CN1936926A (zh) | 一种基于稀疏变换的图像盲源分离方法 | |
CN112102217A (zh) | 一种可见光图像与红外图像快速融合方法及*** | |
CN102760293B (zh) | 基于距离矩阵的图像质量评价方法 | |
CN111528912A (zh) | 一种超声弹性成像方法、装置及*** | |
Li | (Retracted) Infrared image filtering and enhancement processing method based upon image processing technology | |
CN114419074B (zh) | 4k医疗影像处理方法 | |
CN109615591A (zh) | 一种基于gpu并行加速的三维块匹配降噪方法 | |
Ma et al. | Glomerulus extraction by using genetic algorithm for edge patching | |
CN106251304B (zh) | 动态图像分段***和方法 | |
CN114821493A (zh) | 基于计算机视觉、ais和雷达的船舶信息显示方法及*** | |
CN101739684B (zh) | 基于彩色分割和像素显著性估计的视差估计方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
RJ01 | Rejection of invention patent application after publication | ||
RJ01 | Rejection of invention patent application after publication |
Application publication date: 20200814 |