CN114129188B - 一种超声造影灌注流向估计方法与成像系统 - Google Patents

一种超声造影灌注流向估计方法与成像系统 Download PDF

Info

Publication number
CN114129188B
CN114129188B CN202111372351.1A CN202111372351A CN114129188B CN 114129188 B CN114129188 B CN 114129188B CN 202111372351 A CN202111372351 A CN 202111372351A CN 114129188 B CN114129188 B CN 114129188B
Authority
CN
China
Prior art keywords
point
points
streamline
perfusion
flow direction
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
CN202111372351.1A
Other languages
English (en)
Other versions
CN114129188A (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.)
Saset Chengdu Technology Ltd
Original Assignee
Saset Chengdu Technology Ltd
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 Saset Chengdu Technology Ltd filed Critical Saset Chengdu Technology Ltd
Priority to CN202111372351.1A priority Critical patent/CN114129188B/zh
Publication of CN114129188A publication Critical patent/CN114129188A/zh
Application granted granted Critical
Publication of CN114129188B publication Critical patent/CN114129188B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/38—Registration of image sequences
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/06—Measuring blood flow
    • 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/481—Diagnostic techniques involving the use of contrast agent, e.g. microbubbles introduced into the bloodstream
    • 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/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10016—Video; Image sequence
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10132—Ultrasound image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20092—Interactive image processing based on input by user
    • G06T2207/20104—Interactive definition of region of interest [ROI]

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Pathology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • Biophysics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Animal Behavior & Ethology (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Hematology (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明涉及超声成像技术领域，特别涉及一种新的超声造影灌注流向估计方法与成像系统。方法具体包括以下步骤：S1，基于配准后的超声造影图像序列，生成TIC曲线数据；S2，对所述TIC曲线数据进行拟合，从拟合结果中提取灌注参数，并生成二维的灌注参数矩阵，所述灌注参数矩阵包含基于达峰时间的参数矩阵；S3，在所述基于达峰时间的参数矩阵中，依据造影剂从达峰时间取值小的区域到达峰时间取值大的区域估计超声造影灌注流向。本发明基于TIC曲线提出了一种新的超声造影成像模式，实现了超声造影灌注流向的识别和可视化。

Description

一种超声造影灌注流向估计方法与成像系统
技术领域
本发明涉及超声成像技术领域，特别是涉及一种超声造影灌注流向估计方法与成像系统。
背景技术
目前国内常用的超声造影剂是意大利博莱科公司生产的第二代造影剂声诺维，声诺维采用六氟化硫微气泡，并在微气泡表面包以稳定剂，使微气泡在低机械指数下超声共振而不破裂，从而产生较强的谐波信号。超声造影可以显示直径﹤100μm的小血管，而常规超声只能显示直径﹥200μm的血管，因此，超声造影能够更精细地显示组织的血流灌注情况。
具体的，颈部淋巴结相关病理如下：正常淋巴结含有淋巴门结构，血管及输出淋巴管从淋巴门穿行，进入淋巴门的动脉呈分支状分布至皮质，形成毛细血管环。经肘正中静脉团注射声诺维造影剂后，正常淋巴结的超声造影显示造影剂自淋巴门开始增强，呈离心性增强，并迅速均匀分布于整个淋巴结。而发生肿瘤转移的淋巴结，由于肿瘤细胞通过淋巴管首先聚集于边缘窦，继而逐渐侵入淋巴结实质；同时，肿瘤细胞可破坏淋巴结正常的动静脉结构，诱导新生血管形成，所以转移淋巴结超声造影多表现为从周边开始向中心扩散的向心性强化，同时可见局部高灌注、低灌注区或无灌注区；一般高灌注或低灌注区为肿瘤组织，而无灌注区为坏死组织；若病变累及整个淋巴结，常表现为弥漫分布的强弱不等的灌注区。
不同的病例，会在灌注特征(造影成像)和形态(普通二维灰阶超声)等各方面表现出不同。因此，需要在灌注时，了解灌注的情况，并实时通过图像显示出来，以用于病例的分析和判断。
发明内容
本发明在采用由TIC曲线得到的参数进行二维的可视化成像的基础上，对淋巴结造影成像中的灌注流向进行估计和可视化，提出了一种超声造影灌注流向估计方法与成像系统。
为了实现上述发明目的，本发明提供了以下技术方案：
一种超声造影灌注流向估计方法，具体包括以下步骤：
S1，基于配准后的超声造影图像序列，生成TIC曲线数据；
S2，对所述TIC曲线数据进行拟合，从拟合结果中提取灌注参数，并生成二维的灌注参数矩阵，所述灌注参数矩阵包含基于达峰时间的参数矩阵；
S3，在所述基于达峰时间的参数矩阵中，依据造影剂从达峰时间取值小的区域到达峰时间取值大的区域估计超声造影灌注流向。
作为本发明的优选方案，步骤S3具体包括：
S31，在所述达峰时间的参数矩阵中，依据造影剂从达峰时间取值小的区域到达峰时间取值大的区域，得到流线图；
S32，将所述流线图转换为向量场，得到向量场矩阵，所述向量场矩阵中，流线经过的点有非零的向量，其他点都为零向量；
S33，通过所述向量场矩阵预判超声造影灌注的流向并进行成像显示。
作为本发明的优选方案，步骤S31具体包括：
S311，在所述达峰时间的参数矩阵中，选取达峰时间小的区域中的点为起点，在所述起点临近区域内按照搜索条件搜索候选点，并基于搜索出的预定标准候选点和搜索条件，进行迭代搜索，获取一系列的预定标准候选点；
S312，将所述一系列的预定标准候选点依次连接起来，形成估计灌注流向的流线图；
所述搜索条件为：当前候选点的达峰时间值大于前一候选点的达峰时间值，并且所述当前候选点与所述前一候选点的时间差大于时间阈值。
作为本发明的优选方案，步骤S311中，选取起点的步骤包括：
S3111，基于所述达峰时间的参数矩阵中的最大值和最小值计算判定阈值；
S3112，将所述达峰时间的参数矩阵中小于判定阈值的坐标点设置为起点。
作为本发明的优选方案，步骤S311中基于搜索出的预定标准候选点和搜索条件，进行迭代搜索，获取一系列的预定标准候选点具体包括以下步骤：
Q3111，获取以当前候选点为中心点的8邻域点，若从所述8邻域点中得到预定标准候选点，则搜索结束，否则将当前候选点和以当前候选点为中心点的8邻域点加入群pointGroup中；
Q3112，以群pointGroup中的各点为第二当前候选点，获取第二当前候选点的4邻域点，若从所述4邻域点中得到预定标准候选点，则搜索结束，否则将第二当前候选点和以第二当前候选点为中心点的4邻域点加入群pointGroup中；
Q3113，重复步骤S3112，依次得到多个预定标准候选点，直到群pointGroup的长度大于限制值。
作为本发明的优选方案，步骤S312中，所述流线图中流线的方向用所述达峰时间的参数矩阵成像图中的HSV中色调“H”分量来表示。
作为本发明的优选方案，所述步骤S32中具体包括以下步骤：
遍历流线图中每一条流线，使流线中的点的向量方向指向流线中的下一点，生成流线对应的向量矩阵。
作为本发明的优选方案，当流线中的点有多个向量方向时，利用流线计数矩阵选择流线计数最大的下一点为选中指向；
所述流线计数矩阵标记了每个点被流线经过的次数。
作为本发明的优选方案，步骤S2中还包括以下步骤：
对所述TIC曲线数据进行拟合前，对所述TIC曲线数据进行粗筛选，对所述TIC曲线数据进行拟合后，对所述TIC曲线数据进行细筛选；
所述粗筛选是指：按照预设的筛选范围去掉TIC曲线数据中的异常数据；或者将所述TIC曲线数据进行缩放，按照预先设置的缩放筛选条件去掉TIC曲线数据中的异常数据；
所述细筛选是基于NRMSE和修正的余弦相似度的计算值进行筛选得到的，其中，
RMSE为均方根误差，max(fi)是TIC曲线数据拟合后的最大值。
基于相同的构思，还提出了一种超声造影灌注流向估计成像系统，包括成像设备和至少一个处理器，以及与所述至少一个处理器通信连接的存储器；
所述成像设备根据所述处理器输出的结果显示图像；
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够执行上述任一项所述的方法。
与现有技术相比，本发明的有益效果：
本发明基于TIC曲线提出了一种新的超声造影成像模式，实现了超声造影灌注流向的识别和可视化，便于快速分析病灶处的血流信息协助医生做出判断，对在造影视频中利用计算机技术重建体内细血管的研究有一定的借鉴意义。
附图说明：
图1为实施例1中一种超声造影灌注流向估计方法的流程图；
图2为实施例1中视频ROI的人工标定示意图；
图3为实施例1中视频帧运动配准的方法流程图；
图4为实施例1中搜索区域示意图；
图5为实施例1中利用Gamma模型进行曲线拟合的仿真图；
图6为实施例1中基于原始真实值进行判断筛选的仿真图；
图7为实施例1中基于缩放数据进行判断筛选的仿真图；
图8为实施例1中TTP参数成像图；
图9为实施例1中AUC参数成像图；
图10为实施例1中搜索中对中心点的一层一层包围点示意图；
图11为实施例1中依次连接各点生成流线示意图；
图12为实施例1中流线图显示示意图；
图13为实施例1中流线图及流线计数成像图对比图；
图14为实施例1中流线计数成像图的log增强图；
图15为实施例1中LIC成像方式示意图；
图16为实施例1中LIC成像图。
具体实施方式
下面结合试验例及具体实施方式对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。
实施例1
一种超声造影灌注流向估计方法，流程图如图1所示，包括以下步骤：
S1，基于配准后的超声造影图像序列，生成TIC曲线数据；
S2，对TIC曲线数据进行拟合，从拟合结果中提取灌注参数，并生成二维的灌注参数二维矩阵，可成像显示，灌注参数二维矩阵中包含达峰时间的参数二维矩阵；
S3，在达峰时间的参数二维矩阵中，依据造影剂从达峰时间取值小的区域到达峰时间取值大的区域估计超声造影灌注流向。
其中，步骤S1具体包括以下步骤：
S101，视频ROI的人工标定。超声造影视频内容包含超声灰阶B模式和造影模式，超声医生会在B模式的参考帧上面绘制ROI区域，如图2左侧图中的轮廓线，该标记过程即为视频ROI的人工标定。标记的ROI区域的作用在于，在B模式的参考帧上面绘制ROI区域后，视频图像序列中的每一帧根据ROI区域并通过B模式图像序列进行运动配准，人工标定的ROI区域可以根据轮廓线生成一幅Mask二值图像来标记该区域，如图2右侧图中的白色区域。
S102，视频帧运动配准。
视频帧运动配准概述为：在超声造影成像的采集时间内，淋巴结不会发生显著的形变。虽然超声探头的位置相对固定，但病人还是会因为呼吸等身体牵引的运动造成颈部淋巴结的微小移动，造成ROI感兴趣区域的微小运动位移。因此，选择了基于“MAD(平均绝对误差)最小化”和菱形搜索法的块匹配方法对获取的视频帧进行配准。
平均绝对误差(Mean of Absolute Difference,MAD)的计算公式为：
其中，C代表从待配准图像上截取的图像块，R代表参考帧的参考图像块，M和N分别代表模板图像块宽高，i和j代表图像块内像素相对于左上角的坐标。
视频帧运动配准的方法流程图如图3所示，具体包括以下步骤：
B1021，选择K帧参考帧；对一个视频会随机选择了K帧作为候选固定参考帧，都进行配准操作得出配准结果，再选出最优结果。这里选择K为10。
B1022，采用“固定参考帧”法配准；
进一步的，B1022，采用“固定参考帧”法配准的理由为：
因为在给定的视频数据里，某一帧或连续的几帧由于探头位置的移动或病人身体的大幅度移动而造成图像不连续，这就会导致下一帧根本没有淋巴结，用异常结果进行后续的配准将会导致“连锁反应”，最终匹配失败。所以选择了“固定的参考帧”配准法，即利用一帧的ROI区域在其他所有帧上进行配准。淋巴结的位置相对固定，基本不会发生大的位移和形变，所以可以确定一个固定的配准搜索范围。此外，人工标定的ROI区域太小会不能包含组织结构部分，会造成图像块的配准失败，所以搜索模板应该相比人工标定的ROI区域做相应扩展。
如图4所示，若人工标定的ROI区域是中间的椭圆形“1”的话，则将搜索模板变大，扩展为中间的矩形“2”。而搜索范围为以其扩展的外层矩形“3”。这里两次扩展都是上下左右各增加40个像素(触及图像边缘即停止)。
B1023，选择最佳参考帧及其对应的配准结果。
先将所有候选的参考帧都执行一遍运动配准操作，然后选择最佳的参考帧，以及最佳参考帧相应的配准结果。选择策略如下：
(1)配准后会算出配准后的MAD曲线数据均值μ和标准差σ；
(2)对10组(由于选择了K帧参考帧，K取10，因此此处相应选择10组)均值μ和方差σ进行归一化处理(除以其对应的最小值)，避免由于它们取值差异较大造成筛选的误差，
μ＝μ/μmin
σ＝σ/σmin
(3)以“0.6*μ+0.4*σ”作为筛选标准，取其最小值所对应的参考帧，即得到了最佳参考帧及其运动配准结果。
其中，步骤S2概述为：超声造影视频中每一帧的ROI区域提取出一个像素值，n帧就会形成一个长度为n的TIC曲线数据。参数成像就是将视频的ROI区域每一个像素得到一条TIC曲线，针对其中每一个灌注参数(如达峰时间TTP)均可以生成一个二维参数矩阵，再经过颜色编码就可以变为一幅参数成像图像。
步骤S2具体包括以下步骤：
1)逐像素点曲线拟合
A、二维平滑窗口的选择。
实验会遇到一个问题。TIC曲线会出现很多不符合灌注趋势的曲线。
一种原因是信号本身有噪声(通过曲线的拟合处理可以减弱)；一种原因组织内部本来就是并非固定不动的，会随着呼吸发生轻微的运动变化(很细微)。
第一种原因引起的情况能够通过公式拟合进行减弱。对于第二种原因引起的情况，这里考虑利用一个二维平滑窗口，如果将呼吸引起的细微运动的轨迹笼罩住，那窗口内的平均值的变化就能符合灌注趋势了。平滑窗口的大小选择的是“29*29”。
B、TIC曲线数据的提取
时间强度曲线TIC通常是造影模式观察灌注强度变化的常用方法。每一个像素的取值首先会从RGB颜色空间映射到HSV颜色空间，然后利用V分量得到像素的强度值。这样便得到了一条组成曲线的数据。
实验中HSV颜色空间三分量值域分别为，H分量为[0,360]，S分量为[0,1]，V分量为[0,1]。这里提取每一个像素的强度值时，提取颜色V分量后统一缩放至[0,255]。
这里首先采用S-G平滑滤波(Savitsky-Golay滤波器)，多项式阶数是2，平滑窗口大小是81(原视频帧率是10fps，所以这代表8.1秒)。滤波后的数据再进行Gamma曲线拟合。
利用HSV方式提取的数据再经过S-G平滑滤波后被称为“TIC曲线原数据”。
C、基于Gamma函数的曲线拟合
平滑滤波可以消除一部分异常扰动数据，优化曲线拟合的结果。滤波之后利用Gamma函数进行曲线拟合。
在得到的1-D TIC曲线的基础上，利用Gamma模型进行曲线拟合，使得TIC曲线更平滑，便于提取相关的成像参数。图5是gamma拟合的仿真图，其中，较平滑的是拟合曲线，另外的是拟合前数据。
Gamma模型公式为：
Y(t)＝Atce-kt+B
其中，A>＝0，C>0，k>0，B是初始值即t＝0时刻的y值。B为固定值，不属于可变参数。
作为本发明的优选方案，在TIC曲线拟合前后，还对TIC曲线进行了筛选。
提取的TIC曲线原数据，需要进行筛选，因为并不是ROI区域内的所有像素均有血流经过。对于没有很好灌注趋势的像素点应该去除。粗筛选是曲线拟合之前，利用数据原本的特点进行一次快速地筛选；细筛选就是利用曲线拟合结果的好坏进行筛选，必须粗细筛选都通过才能继续拟合，否则是失败的点。
第一步，粗筛选。
粗筛选是在曲线拟合前对直接的TIC曲线原数据先进行简单快速的判断，去掉存在异常的曲线原数据。大致分两个思路，一个是基于真实值进行筛除，一个是将TIC曲线原数据缩放到[0，m]的值，再进行筛除，m属于一个正整数，这里m为8。
(1)基于真实值判断，就是挑选TIC曲线原数据的峰值、起始值判断其大小。而其中时间强度曲线中的强度值是HSV颜色空间的“V”分量得到的，范围统一为[0，255]。
比如实验中，基于真实值筛选，留下的数据为：峰值不能小于30，曲线起始值不大于100，如图6所示，峰值不能太小，因为太小造影剂没灌注效果；起始值不能太大，太大就是起始亮度太大，在造影里就是组织。
(2)基于缩放值判断，将曲线原数据值缩放到[0，8]，就可以对数据中的第二峰值、后面(1/2时间之后)最大值、起伏过程中的上升幅度、起伏过程中的异常上升总次数等衡量数据进行定量化。这是非缩放到固定范围的数据做不到的。
比如实验中，第二峰值在数据后面不能大于5，在前面不能相比之前的波谷上升超过3。后端1/5时间段最大值不能是6。上升幅度大于3为异常上升，异常上升次数不能多于3个。基于缩放值进行筛选的仿真示意图如图7所示。
当基于真实值判断和基于缩放值判断后均保留下来的TIC曲线数据，才通过了粗筛选。
第二步，细筛选。
这一步就是用Gamma公式去拟合粗筛选成功的曲线原数据，然后根据拟合情况再次对数据进行筛选。
拟合度准则如下：
A、NRMSE
即归一化的RMSE,使RMSE除以拟合值最大值。
RMSE即均方根误差，即
其中，ri是原始数据，fi是拟合后的数据。
B、adjusted cosine
即修正的余弦相似度，公式：
u即原数据r与拟合数据f合在一起共同的平均值。
使用修正的余弦相似度需要两个时间序列长度一样，时间点一一对应。
最终，只有当NRMSE≤0.15，并且adjust cosine≥0.75两个条件都满足的拟合数据才能保留下来。
当曲线成功通过粗筛选，又通过细筛选，则细筛选公式的结果作为最终的TIC曲线数据用于提取灌注参数特征。
2)提取2-D灌注参数特征
经过双重筛选后的TIC曲线，可以提取灌注参数，而参数大致分两类：(1)基于时间的，如TTP(达峰时间)、MTT(平均渡越时间)；(2)基于数值的，如AUC(曲线下面积)、PI(峰值)。
每条曲线可得一组灌注参数，而所有曲线对于一个灌注参数可形成二维的数据矩阵，即“灌注参数二维矩阵”。
其中，对于拟合失败的曲线，采取人工设置值的方式给一个默认值。
3)灌注参数成像的可视化
得到灌注参数二维矩阵后，根据二维矩阵中的数据进行成像，生成参数矩阵成像图并着色。
首先，将数据的取值范围调整到0～240。因为HSV里面的H分量，[0，240]刚好是由红到蓝。
假设有一组灌注参数二维矩阵数据，最大值为maxValue,最小值为minValue。对于灌注参数二维矩阵中的任一坐标对应求得映射值inputValue；
基于映射公式计算出映射值，映射值计算公式为：
然后，基于该映射值进行颜色编码，对于一个坐标点，将转换后的映射值作为HSV的“h”分量，另外两个分量S和V设置成相同的默认值“1”，则实现了灌注参数成像的可视化。在灌注参数成像的可视化过程中有两种成像着色方法，其中一个是采用“大到小，蓝到红”的方法，另一个是采用“小到大，蓝到红”的方法。
对于TTP(达峰时间)来说，达峰时间越早医生越关注，采用偏红色的颜色编码。这一类数据采用“大到小，蓝到红”的方法。求得某坐标点的映射值outputValue后，则该坐标点的HSV颜色值为(floor(outputValue)，1，1)，例如，求得某一个坐标点的映射值outputValue为220.2，则该坐标点的HSV颜色值为(220，1，1)。对灌注参数二维矩阵中每个坐标点都按照该方法进行转换，遍历参数矩阵，就得到了参数矩阵成像图。图8给出了TTP参数成像图，ROI区域内为拟合好的点，其中两个黑点是拟合失败的点。其中，floor()是向下取整函数。
对于AUC(曲线下面积)来说，则值越大表示越可能灌注总量大就越是受关注的。这一类数据采用“小到大，蓝到红”的方法。求得某坐标点的映射值outputValue后，则该坐标点的HSV颜色值为(240-floor(outputValue)，1，1)，例如，求得某一个坐标点的映射值outputValue为220.4，则该坐标点的HSV颜色值为(20，1，1)。对灌注参数二维矩阵中每个坐标点都按照该方法进行转换，遍历参数矩阵，就得到了参数矩阵成像图。图9给出了AUC参数成像图，ROI区域内为拟合好的点，其中两个黑点是拟合失败的点。
其中，步骤S3概述为：淋巴结良性病例一般会先从淋巴门开始灌注，有些恶性病例可能会因为内部组织病变而从外围开始灌注。
灌注方向的估计是本发明的关键创新点，这里考虑利用TTP(达峰时间)参数矩阵进行判断。估计方法的规则是依据造影剂从TTP小的区域流向TTP大的区域。
步骤S3具体包括以下步骤：
1)流线种子点的选取
在TTP二维数据矩阵中查找最大值maxValue和最小值minValue，求得灌注比较早的区域的判定阈值Th，判定阈值Th的计算公式如下：
Th＝minValue+(maxValue-minValue)*frac
其中Th即灌注比较早的区域的判定阈值，frac为“0”到“1”之间的小数。这里设置frac＝0.3。
将小于判定阈值Th的TTP二维数据点全部设置为种子点。
2)每条流线的生成
利用TTP参数二维数据矩阵计算流线需要利用当前点与其邻域点之间的关系。以达峰时间TTP得到的2-D参数矩阵为例，以秒为单位，一般小数点后一位数。通过设置一个起始点，在该点的邻域内搜寻TTP大于该点的值且时间差大于“0.2”或“0.3”秒的最小值，因为这些坐标的曲线达峰时间接近当前点又不会被认为同时达到峰值。一步一步进行迭代，每一个点用相同的方法寻找下一个候选点，可以估计出造影剂逐渐增强的大致灌注轨迹。
我们的搜索范围并不是采用固定的8邻域，而是以当前点为中心，一层一层(1、2到3的顺序)地搜寻外围的邻域点，从而找到一个TTP大于该点且超过一定数值的最小值。图10中给出了搜索过程中对中心点的一层一层包围点示意图，可称为“寻点二维矩阵”，通过“寻点二维矩阵”实现了灌注方向估计。采用图10的邻域设置，从起点坐标依次找到下一个点的坐标，然后把这些坐标点依次连接起来，就得到了流线。以图10的邻域设置为例说明以当前点P搜索下一个候选点的过程如下：
A、设当前点为P，并设置预定标准为下一点是该点值是比当前点大0.2的点的最小值。
B、搜索P点的8邻域点，如图中标记为“1”的那些点。如果找到预定标准的候选点就结束；否则将当前点和其8邻域点加入一个点群pointGroup，然后进行第C步；
C、继续搜索pointGroup点群的4邻域点，如图中标记为“2”的那些点。参考点依旧为中心点P(即第A步的点)。如果找到预定标准的点就结束；否则就将其周围点加入pointGroup点群，重复本步骤，搜索图中标记为“3”的那些点，直到pointGroup点群的长度大于限制值即终止。
搜索后生成流线的示意图如图11所示，一条流线的例子，从A点开始，用上述过程一次可以找下一个点是B。点B可以找到点C，直到F点结束(下一个点找不到终止)就形成了一条流线。
3)流线的显示
通常一条流线可以利用链表的形式给出，即一条流线的点组成一个列表list。
对找到的流线list中的下标相邻的两个点，假如在图像物理空间不相邻，利用Bresenham算法可以找到一条近似直线的路径连接起来。
但是，方向怎么表示？原图本来就小，再画箭头几乎看不清。于是，利用TTP参数成像图原本的HSV中色调“H”分量来表示先后关系。流线图示意如图12所示，左侧为原始TTP参数成像图像，中间为提取的TTP较小区域作为流线搜索起始区域，右侧为流线结果图。
另外，对该ROI区域内的每一个点统计经过的流线，即每条流线经过一个点就将该坐标的计数加1，所有流线叠加可以生成一个流线计数二维矩阵。其意义是，流线经过比较多的区域有很大可能是有血管的区域，通过流线计数二维矩阵可以找到流线经过比较多的区域。流线图与流线计数成像图如图13所示，左为流线图，右为流线计数成像图，由于流线计数图各线条太暗，于是利用log增强技术将灰度较暗的细节进行了增强，并显示，流线计数成像图的log增强图如图14所示。
有关流线计数图的计算方式，大致可以描述为：
第一步，遍历所有利用bresenham算法扩充后的流线列表list；
第二步，遍历每条流线list中的每一个点，当前点被流线流过一次，则相应流线计数矩阵中对应的坐标点数字加1。
4)流线图转向量场和LIC显示
通过找种子点可以生成很多流线，利用这些流线可以生成一个二维向量场，并利用LIC(线积分卷积)的方式进行显示，LIC成像方式如图15所示。
一条流线从前到后依次从一个点指向下一个点，走过的点可以赋值成向量形式，在生成向量场过程中，只有流线经过的点有非零的向量，其他点向量都为零向量。遍历所有流线的所有点，就可以得出一个向量场矩阵。
设置向量场矩阵大小为(h，w，2)，即高宽与原图像相同而每个元素两个值{dx,dy}表示向量大小。同时生成一张mask标记各坐标是否设置了向量(即同时生成了流线计数二维矩阵)。则向量场矩阵生成过程为：
(1)当遍历某条流线到某一点，向量场矩阵中该坐标未设置向量，则对该坐标设置向量指向流线下一个点(没有下一点就设置为零向量)；
(2)当向量矩阵中某个点有多个向量指向的时候，利用流线计数矩阵选择流线计数最大的下一点为最优指向(流线计数矩阵多个坐标的流线计数相等时，选择最先遍历到的流线的下一点)。
向量场生成的LIC图如图16所示，向量场生成LIC图时要对零向量区域进行一步抑制，缩小其亮度范围，并且利用参数成像图的颜色进行着色。而在实验中图像已经很小，不能用箭头标记向量场。只能将计算结果显示出来。结果中也进行了着色来表示向量方向。
以上显示和描述了本发明的基本原理和主要特征及本发明的优点，对于本领域技术人员而言，显然在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
此外，应当理解，虽然本说明书按照实施方式加以描述，但并非实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (8)

1.一种超声造影灌注流向估计方法，其特征在于，具体包括以下步骤：
S1，基于配准后的超声造影图像序列，生成TIC曲线数据；
S2，对所述TIC曲线数据进行拟合，从拟合结果中提取灌注参数，并生成二维的灌注参数矩阵，所述灌注参数矩阵包含基于达峰时间的参数矩阵；
S3，在所述基于达峰时间的参数矩阵中，依据造影剂从达峰时间取值小的区域到达峰时间取值大的区域估计超声造影灌注流向；
步骤S3具体包括：
S31，在所述达峰时间的参数矩阵中，依据造影剂从达峰时间取值小的区域到达峰时间取值大的区域，得到流线图；
S32，将所述流线图转换为向量场，得到向量场矩阵，所述向量场矩阵中，流线经过的点有非零的向量，其他点都为零向量；
S33，通过所述向量场矩阵预判超声造影灌注的流向并进行成像显示；步骤S31具体包括：
S311，在所述达峰时间的参数矩阵中，选取达峰时间小的区域中的点为起点，在所述起点临近区域内按照搜索条件搜索候选点，并基于搜索出的预定标准候选点和搜索条件，进行迭代搜索，获取一系列的预定标准候选点；
S312，将所述一系列的预定标准候选点依次连接起来，形成估计灌注流向的流线图；
所述搜索条件为：当前候选点的达峰时间值大于前一候选点的达峰时间值，并且所述当前候选点与所述前一候选点的时间差大于时间阈值。
2.如权利要求1所述的一种超声造影灌注流向估计方法，其特征在于，步骤S311中，选取起点的步骤包括：
S3111，基于所述达峰时间的参数矩阵中的最大值和最小值计算判定阈值；
S3112，将所述达峰时间的参数矩阵中小于判定阈值的坐标点设置为起点。
3.如权利要求1所述的一种超声造影灌注流向估计方法，其特征在于，步骤S311中基于搜索出的预定标准候选点和搜索条件，进行迭代搜索，获取一系列的预定标准候选点具体包括以下步骤：
Q3111，获取以当前候选点为中心点的8邻域点，若从所述8邻域点中得到预定标准候选点，则搜索结束，否则将当前候选点和以当前候选点为中心点的8邻域点加入群pointGroup中；
Q3112，以群pointGroup中的各点为第二当前候选点，获取第二当前候选点的4邻域点，若从所述4邻域点中得到预定标准候选点，则搜索结束，否则将第二当前候选点和以第二当前候选点为中心点的4邻域点加入群pointGroup中；
Q3113，重复步骤S3112，依次得到多个预定标准候选点，直到群pointGroup的长度大于限制值。
4.如权利要求1所述的一种超声造影灌注流向估计方法，其特征在于，步骤S312中，所述流线图中流线的方向用所述达峰时间的参数矩阵成像图中的HSV中色调“H”分量来表示。
5.如权利要求1所述的一种超声造影灌注流向估计方法，其特征在于，所述步骤S32中具体包括以下步骤：
遍历流线图中每一条流线，使流线中的点的向量方向指向流线中的下一点，生成流线对应的向量矩阵。
6.如权利要求5所述的一种超声造影灌注流向估计方法，其特征在于，
当流线中的点有多个向量方向时，利用流线计数矩阵选择流线计数最大的下一点为选中指向；
所述流线计数矩阵标记了每个点被流线经过的次数。
7.如权利要求1-6任一所述的一种超声造影灌注流向估计方法，其特征在于，步骤S2中还包括以下步骤：
对所述TIC曲线数据进行拟合前，对所述TIC曲线数据进行粗筛选，对所述TIC曲线数据进行拟合后，对所述TIC曲线数据进行细筛选；
所述粗筛选是指：按照预设的筛选范围去掉TIC曲线数据中的异常数据；并且将所述TIC曲线数据进行缩放，按照预先设置的缩放筛选条件去掉TIC曲线数据中的异常数据；
所述细筛选是基于NRMSE和修正的余弦相似度的计算值进行筛选得到的，其中，
，
RMSE为均方根误差，是TIC曲线数据拟合后的最大值；
所述修正的余弦相似度计算公式为
其中，是余弦相似度，/>是原始数据，/>是拟合后的数据，u是原始数据与拟合后的数据合在一起共同的平均值。
8.一种超声造影灌注流向估计成像系统，其特征在于，包括成像设备和至少一个处理器，以及与所述至少一个处理器通信连接的存储器；
所述成像设备根据所述处理器输出的结果显示图像；
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够执行权利要求1至7中任一项所述的方法。
CN202111372351.1A 2021-11-18 2021-11-18 一种超声造影灌注流向估计方法与成像系统 Active CN114129188B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111372351.1A CN114129188B (zh) 2021-11-18 2021-11-18 一种超声造影灌注流向估计方法与成像系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111372351.1A CN114129188B (zh) 2021-11-18 2021-11-18 一种超声造影灌注流向估计方法与成像系统

Publications (2)

Publication Number Publication Date
CN114129188A CN114129188A (zh) 2022-03-04
CN114129188B true CN114129188B (zh) 2023-09-15

Family

ID=80390459

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111372351.1A Active CN114129188B (zh) 2021-11-18 2021-11-18 一种超声造影灌注流向估计方法与成像系统

Country Status (1)

Country Link
CN (1) CN114129188B (zh)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101779968A (zh) * 2009-01-07 2010-07-21 株式会社东芝 医用图像处理装置和超声波图像取得装置
JP2014200616A (ja) * 2013-04-10 2014-10-27 株式会社東芝 医用画像処理装置、ｘ線診断装置及び医用画像処理プログラム
CN104363834A (zh) * 2012-06-07 2015-02-18 株式会社东芝 图像处理装置以及x射线诊断装置
CN104688272A (zh) * 2015-03-06 2015-06-10 西安交通大学 一种基于单像素tic源的超声造影灌注参量成像方法
CN104688269A (zh) * 2015-03-06 2015-06-10 西安交通大学 一种时间强度曲线的呼吸运动补偿及双峰拟合的方法
WO2018130601A2 (en) * 2017-01-12 2018-07-19 Stichting Katholieke Universiteit Extracting flow information from a dynamic angiography dataset
CN109952063A (zh) * 2016-11-14 2019-06-28 皇家飞利浦有限公司 用于表征造影剂流的肝脏灌注的系统和方法
CN110458836A (zh) * 2019-08-16 2019-11-15 深圳开立生物医疗科技股份有限公司 一种超声造影成像方法、装置和设备及可读存储介质
CN111110277A (zh) * 2019-12-27 2020-05-08 深圳开立生物医疗科技股份有限公司 超声成像方法、超声设备及存储介质

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110150309A1 (en) * 2009-11-27 2011-06-23 University Health Network Method and system for managing imaging data, and associated devices and compounds
US9192357B2 (en) * 2013-02-19 2015-11-24 Kabushiki Kaisha Toshiba Method and system for quantitative vectorial perfusion based upon blood flow direction using 4D medical imaging

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101779968A (zh) * 2009-01-07 2010-07-21 株式会社东芝 医用图像处理装置和超声波图像取得装置
CN104363834A (zh) * 2012-06-07 2015-02-18 株式会社东芝 图像处理装置以及x射线诊断装置
JP2014200616A (ja) * 2013-04-10 2014-10-27 株式会社東芝 医用画像処理装置、ｘ線診断装置及び医用画像処理プログラム
CN104688272A (zh) * 2015-03-06 2015-06-10 西安交通大学 一种基于单像素tic源的超声造影灌注参量成像方法
CN104688269A (zh) * 2015-03-06 2015-06-10 西安交通大学 一种时间强度曲线的呼吸运动补偿及双峰拟合的方法
CN109952063A (zh) * 2016-11-14 2019-06-28 皇家飞利浦有限公司 用于表征造影剂流的肝脏灌注的系统和方法
WO2018130601A2 (en) * 2017-01-12 2018-07-19 Stichting Katholieke Universiteit Extracting flow information from a dynamic angiography dataset
CN110458836A (zh) * 2019-08-16 2019-11-15 深圳开立生物医疗科技股份有限公司 一种超声造影成像方法、装置和设备及可读存储介质
CN111110277A (zh) * 2019-12-27 2020-05-08 深圳开立生物医疗科技股份有限公司 超声成像方法、超声设备及存储介质

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于LIC的海流动态模拟与可视化;刘振东等;《海洋通报》;第38卷(第02期);179-185 *

Also Published As

Publication number Publication date
CN114129188A (zh) 2022-03-04

Similar Documents

Publication Publication Date Title
Luo et al. Micro-vessel image segmentation based on the AD-UNet model
CN110503649B (zh) 一种基于空间多尺度U-net与超像素修正的肝脏分割方法
Garnavi et al. Automatic segmentation of dermoscopy images using histogram thresholding on optimal color channels
CN109431531B (zh) 基于灌注成像的血管分割方法及装置和计算机装置
JP2004032684A (ja) 医療画像において腫瘤や実質組織変形をコンピュータを用いて検出する自動化した方法と装置
CN108615239B (zh) 基于阈值技术和灰度投影的舌图像分割方法
Qureshi et al. Detection of glaucoma based on cup-to-disc ratio using fundus images
JP2008520345A (ja) 超音波画像における病変の検出及び分類方法、及びそのシステム
WO2011066689A1 (en) Method and device for detecting bright brain regions from computed tomography images
Goyal et al. Multi-modality image fusion for medical assistive technology management based on hybrid domain filtering
CN116071337A (zh) 一种基于超像素分割的内镜图像质量评价方法
Schenk et al. Automatic high-speed video glottis segmentation using salient regions and 3D geodesic active contours
CN114066911A (zh) 一种基于乳腺和肿瘤分割的bpe自动化提取方法和系统
CN112712540B (zh) 一种基于ct影像的肺部支气管提取方法
CN117274216B (zh) 基于水平集分割的超声下颈动脉斑块检测方法及系统
Zhang et al. Retinal spot lesion detection using adaptive multiscale morphological processing
CN117522862A (zh) 一种基于ct影像肺炎识别的图像处理方法及处理系统
CN111091578B (zh) 一种针对血管医学影像的分割方法
CN114129188B (zh) 一种超声造影灌注流向估计方法与成像系统
CN108629780B (zh) 基于颜色分解和阈值技术的舌图像分割方法
CN116523924A (zh) 一种医学实验用数据处理方法及系统
CN113298754B (zh) 一种前列腺组织轮廓线控制点的检测方法
CN113658193B (zh) 基于信息融合的肝脏ct影像肿瘤分割方法
Bouzidi et al. BrainSeg3D to detect Multiple sclerosis lesions using magnetic resonance imaging
Almi'ani et al. A modified region growing based algorithm to vessel segmentation in magnetic resonance angiography

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
GR01 Patent grant
GR01 Patent grant