CN104463844B - 基于brox光流法的超声造影图像血管灌注区提取方法 - Google Patents

Abstract

本发明公开了一种基于brox光流法的超声造影图像血管灌注区提取方法，按照以下步骤进行：1)获取图像帧，形成超声造影图像序列；2)图像预处理，减少每帧超声造影图像的斑点噪声；3)基于brox光流法估计相邻两帧超声造影图像的运动场，并根据估计出的偏差量对相邻两帧超声造影图像进行校正；4)利用阈值分割法提取出血管灌注区域的点。其效果是：本方法采用了对噪声不敏感、鲁棒性强的brox光流法估计出相邻图像帧的运动位移并进行校正，基于校正后的超声造影图像序列，根据造影微泡在血液中表现出的高亮、闪烁，区别于其他背景组织的信号特点提取出灌注区域，实验结果表明，本方法计算量小，实现简单，应用于临床数据时，提取效果好。

Description

基于brox光流法的超声造影图像血管灌注区提取方法

技术领域

本发明涉及到医学图像处理技术，具体地说，是一种基于brox光流法的超声造影图像血管灌注区提取方法。

背景技术

超声造影成像是一种通过注入造影剂来动态、清晰地显示微细血管的新型成像技术，该技术可以获得组织丰富的供血及血流灌注信息。造影微泡具有较强的散射性并与周围血液形成高声阻抗差，注入血流后可使血液回声增强，达到良好的“血管显影”效果。从周围组织强背景中提取出血管灌注区域是对造影信号定量分析的基础,这对临床诊断组织病变,以及了解肿瘤的生长状况具有重要的临床意义。

目前对测定感兴趣区域(ROI)的时间—强度曲线(TICs)研究的比较多，TICs反映了造影剂注入后组织某点增强的回波信号(即图像中像素的灰度)随时间变化的曲线。基于TICs的特征参数进行参数成像，借助超声灌注成像来分析血管灌注区域以此来定量分析超声造影图像，然而由于超声图像噪声强、伪影多以及组织回波信号的多样性，难以确定合适的数学模型，大大增加了超声灌注参数成像的难度，因此该成像技术尚处于研究阶段。也有人提出了一种利用微泡破坏机制的多脉冲释放成像(multi-pulse release imaging)方法，对血管灌注区域进行成像。这种成像方法要求释放探头与成像探头共同聚焦于同一成像目标区，两个探头很难协调扫描，因此该方法很难成功实现；西安交通大学的宋延淳等人还提出一种基于高机械指数(MI)、单成像脉冲发射导致微泡破坏的射频回波解相关(DCR)检测方法，在仿体实验上取的较好的实验效果。

通过以上描述可以发现，现有技术中大多都是基于破坏性成像所产生的回波信号的基础上进行解相关运算，建模及分析过程复杂，仅仅是实现于离体灌注仿体模型中，难以应用于临床数据。

发明内容

针对现有技术的不足，本发明提出了一种直接对造影图像序列进行动态分析，通过造影微泡与周围组织成像的差别提取出血管灌注区域的方法。与普通图像分析不同，超声造影动态图像的分析方法通常以固定脏器和病变部位的造影信号的变化为基础，而采集图像的过程中ROI区域会受到呼吸及心跳的影响产生位移，因此超声造影图像的运动校正是保证动态分析造影图像序列结果正确的关键步骤。

为达到上述目的，本发明的具体技术方案如下：

一种基于brox光流法的超声造影图像血管灌注区提取方法，按照以下步骤进行：

步骤1：获取图像帧，形成超声造影图像序列；

步骤2：图像预处理，减少每帧超声造影图像的斑点噪声；

步骤3：基于brox光流法估计相邻两帧超声造影图像的运动场，并根据估计出的偏差量对相邻两帧超声造影图像进行校正；

步骤4：利用阈值分割法提取出血管灌注区域的点。

作为进一步描述，步骤1中的具体步骤为：

步骤1-1：截取3～8秒灌注峰值期的录像；

步骤1-2：按照5帧/秒的采样率获取图像帧，并进行连续编号；

步骤1-3：选取编号为奇数的图像，重新编号形成新的图像序列。

再进一步描述，步骤2中的图像预处理采用低通高斯滤波器，高斯滤波窗口为3*3。

更具体地，步骤4中的阈值分割法包括以下步骤：

步骤4-1：利用8邻域平均法对每帧图像的每个像素点进行平滑处理；

步骤4-2：按照计算相邻两帧图像间每个像素点灰度变化的平均值，并将所有像素点灰度变化的平均值构建成矩阵T；其中s为图像的帧序列，a为图像的总帧数，I(s)为第s帧图像对应像素点的灰度值；

步骤4-3：设定阈值t，并将矩阵T中高于阈值t的点标定为灰度变化剧烈点，即为造影微泡灌注的位置。

作为优选，所述阈值t按照以下方式设定：

步骤4-3-1：在医生指导下初步判断血管灌注区域，并选择三个样本区域作为观察对象；

步骤4-3-2：按照求取样本区域的像素点灰度变化的平均值，分别构建为m₁×n₁的矩阵T₁，m₂×n₂的矩阵T₂以及m₃×n₃的矩阵T₃，对应的m₁×n₁为第一样本区域的大小，m₂×n₂为第二样本区域的大小，m₃×n₃为第三样本区域的大小；

步骤4-3-3：按照

计算得到所述阈值t。

本发明的显著效果是：

本方法采用对噪声不敏感、鲁棒性强的brox光流法估计出相邻图像帧的运动位移并进行校正，基于校正后的超声造影图像序列，根据造影微泡在血液中表现出的高亮、闪烁，区别于其他背景组织的信号特点提取出灌注区域，实验结果表明，本方法计算量小，实现简单，应用于临床数据时，提取效果好。

附图说明

图1是本发明的方法步骤图；

图2是子宫肌瘤超声造影图像序列以及矩形框中的光流场，其中图2a为原始子宫肌瘤超声造影图像序列；图2b是第1帧与第2帧图像白色矩形框中像素的光流矢量图，箭头表示位移方向，线段的长短表示位移量的大小；

图3是样本区域的选择结果图；

图4是血管灌注区域提取结果示意图，其中图4a是用二值图像表示提取结果；图4b是在造影图像上用红色点标记出灌注区域，曲线内为ROI；

图5是运动校正后和灌注区域处理后图像序列的ROI的TICs。

具体实施方式

下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。

如图1所示，一种基于brox光流法的超声造影图像血管灌注区提取方法，按照以下步骤进行：

步骤1：获取图像帧，形成超声造影图像序列，具体为：

步骤1-1：截取3～8秒灌注峰值期的录像；

步骤1-2：按照5帧/秒的采样率获取图像帧，并进行连续编号；

步骤1-3：选取编号为奇数的图像，重新编号形成新的图像序列。

本例中的实验数据来源于重庆海扶医院的子宫肌瘤患者，实施过程中，先在高强度聚焦超声(HIFU)治疗前需要先静脉推注2.0ml(10.0mg)造影剂声诺维，以观察肌瘤边界及血流灌注情况。超声实时监控设备为意大利百盛公司生产的My-Lab70型B超，频率1.0～8.0MHz，超声探头置于聚焦超声换能器的中央，注射造影剂前将超声探头调至病人矢状位、子宫肌瘤最大切面处。有统计研究表明,经肘静脉注入造影剂至子宫肌瘤开始显影时间，平均约为15s，肌瘤开始灌注到达到峰值强度约为12s。本例中采集数据时，待注入造影剂10s后开始待录制造影过程，录制时间为53s，按照5帧/秒的采样率，可获得265帧图像。我们需要的是子宫肌瘤灌注达到峰值时期的图像，观看录像，发现录制至第30s的时候图像灌注达到峰值强度并持续超过10s，本例选取30s到34s共20帧图像，都处于灌注峰值期，编号为图像1～20，为了避免相邻两帧间造影剂的变化过于微小，取编号为奇数的图像，形成新的图像序列，重新编号为图像1～10。

步骤2：图像预处理，减少每帧超声造影图像的斑点噪声；

本例中，子宫肌瘤的超声造影图像序列来自于B超机的视频录像，要想分析一帧帧图像必须得先按照视频的录制帧率解码为连续的单帧图像。由于超声图像存在大量的斑点噪声，所以必须进行降噪处理。这里采用高斯滤波的方法进行降噪，将所有的超声造影图像序列通过一核大小为3*3的低通高斯滤波器，高斯核函数其中A为一常数，可以根据不同的设备调整，σ表征高斯滤波器的宽度，决定着平滑程度，本例取σ＝1.8，经过预处理减少了斑点噪声的影响以增强血管灌注区域提取算法的准确性与稳定性。

步骤3：基于brox光流法估计相邻两帧超声造影图像的运动场，并根据估计出的偏差量对相邻两帧超声造影图像进行校正；

由于子宫肌瘤属于盆腔肿瘤，人体的呼吸会影响超声图像运动。为了减少呼吸运动的影响，最常见的方法是控制呼吸，然而录制造影过程的时间往往需要20s以上，这对患者来说是很困难的，有些虚弱的患者甚至不能坚持几秒钟。还有一种方法是对超声造影图像序列进行运动校正，针对超声造影图像的运动校正的研究较少，很多配准算法无法取得满意的结果，主要原因是超声图像本身噪声多，伪影多，加上造影剂的不断摄入，使得相邻造影图像间出现明显的灰度差异，致使算法失效。

而本发明中采用brox光流法进行运动估计，光流法是分析运动图像的重要方法，光流表达了图像的变化，包含了图像运动的信息。光流场的计算最初由Horn和Schunk提出,在最近的20年得到了较大的发展，研究人员提出了很多不同的改进算法。本发明使用的是brox等人于2004年提出的一种新颖的变分光流法。传统的光流方程为:

I_xu+I_yv+I_t＝0 (1)

其中u和v分别表示此像素点在x和y方向的位移，I_x，I_y，I_t分别表示像素灰度沿x，y，t(时间)方向的梯度，表示的物理含义是图像上某一点灰度的变化率是亮度变化率与该点运动速度的乘积。传统光流法的计算方法主要是基于灰度连续的假设和光流场平滑的假设，而超声造影图像上存在斑点噪声，伪影，以及造影剂的地方会使得这些假设不成立。brox光流法在此基础上增加了梯度连续的假设，减小传统光流法对灰度值变化的敏感性，约束方程如下：

其中表示图像的空间梯度，和(1)式中一样,(u,v)是位移向量场。引入梯度连续的假设后，允许前后帧图像的灰度值发生小的变化的情况下，也能估计出位移。同时为了增加对大范围位移的估计，以及避免目标函数容易陷入局部最小值的缺陷，在算法中使用了多尺度技术。多尺度技术是采用一系列的向下采样滤波，降低图像的尺寸和分辨率，先将尺寸和分辨率最小的图像作为运动估计的初始对象，再一步步的提高图像的尺寸与分辨率，将前一层估计出来的结果作为下一层的初始值迭代至下一层估计中去，直至达到原始图片的大小。相对于传统的光流法无法估计位移较大情况而言，多尺度分解技术可以将超过几个像素的大位移分解为从粗糙到精细的迭代估计过程，得出准确的位移场。因为超声造影图像前后帧的位移能够达到5个像素以上，所以多尺度技术是必不可少的。

brox光流法的能量函数为：

E(u,v)＝E_Data+αE_Smooth (3)

展开如(4)，(5)式：

其中令x＝(x,y,t)^T,w＝(u,v,1)^T；表示时间—空间梯度，当只估计两帧图像时，可用空间梯度来代替。E_Data(u,v)是灰度连续假设与梯度连续假设的能量函数，γ为两种假设之间的权重，γ越大，算法对灰度守恒的依赖越小。E_smooth(u,v)为光流场平滑假设的能量函数，该假设是为了避免图像中出现图像梯度消失而造成算法估计出异常值或位移不连续的的情况，α为数据项(4)与平滑项(5)之间的权重。公式通过二次平方函数(ε为大于0的常数)加强了对孤立点的惩罚，提高了算法的鲁棒性。

本例中brox光流法中平滑项权重α设为30.0，梯度不变的权重γ设为80.0，针对预处理后的1～10帧图像，求出使得总能量函数(3)式最小值时的(u，v)即得图像运动位移场，最后根据(u，v)对相邻两帧图像每个像素进行校正。

运动估计结果如图2所示，图2b是图2a中白色矩形框中像素的光流场,为了显示方便，作了放大处理。白色矩形框选取在肌瘤的边界处，可见光流场的左上部分有明显的位移，与观看造影录像观察到的位移方位一致，右下部分光流的场的模很小，说明几乎没有位移。这跟实际图像的情况是一致的：左上部分是肿瘤部分，随着呼吸会有几个像素的位移，而右下角为贴近皮肤部分，位移很小。

为了验证校正运动算法的性能，本文采用计算互信息(MI)和互相关系数(CC)作为衡量指标，结果如表1所示：

表1两两相邻帧图像运动校正结果

互信息(MI)是两幅图像相关性的测度,互信息越大说明两幅图像越接近,互相关系数(CC)是用来表达两幅图像的线性相关性,该值越大表明图像校正的精确度越高。从表1可以看出校正后较于校正前，MI增大，CC增大，表明校正算法的有效性。

步骤4：利用阈值分割法提取出血管灌注区域的点。

超声造影成像是利用微泡造影剂经过外周静脉注射后进入全身血液循环，使得病灶处的回声和血流信号增强而提高超声诊断的分辨力、敏感性和特异性的技术。造影微泡的背向散射较红细胞强得多，而且微泡与组织不同的是它还是声波发射体，在入射声强的影响下，造影微泡会发生膨胀与收缩振动，当膨胀幅度超过收缩幅度时，导致波形的畸变，从而产生谐波效应；当超声功率达到某一临界点时，微泡局部受到的压力大于它本身可承受的压力,造影微泡会被破坏，产生强烈的背向散射，显示出该病灶区血管容量的信息，这些都是周围组织并不具备的特性。对子宫肌瘤内血供丰富区域的判断主要依赖于肌瘤相对于周围正常组织的回声情况。子宫肌瘤周围的组织有呈强回声的，有呈低回声的，而肌瘤内血供丰富的区域在灌注峰值的时候表现为强回声，单从孤立的一帧图像上无法区分出呈强回声的部分是正常组织还是肌瘤血供丰富的区域。当多帧图像连续播放时，可以发现与正常组织强回声区域在多帧图像中强度变化不大不同，造影微泡灌注的区域呈现出高增强并且强度会发生很大的变化，这是造影微泡在入射声强的影响下不断发生膨胀、压缩的甚至破灭的动态过程的表现。可以把超声成像束看成是有一定厚度的平面薄板，薄板上的造影气泡流进、流出、破灭、波动等变化在灌注峰值时达到一个动态平衡，在超声图像就表现为高亮、闪烁的似血液样的流动区域，明显区别于周围正常组织以及肿瘤内部乏血供的部分。因此方法根据图像序列相邻帧之间灌注区域与非灌注区域灰度值变化幅度的差异设计了一种提取血管灌注区域的算法，结合阈值分割的思想，具体包括以下步骤：

步骤4-1：利用8邻域平均法对每帧图像的每个像素点进行平滑处理；

步骤4-2：按照计算相邻两帧图像间每个像素点灰度变化的平均值，并将所有像素点灰度变化的平均值构建成矩阵T；其中s为图像的帧序列，a为图像的总帧数，本例中a＝10，I(s)为第s帧图像对应像素点的灰度值；

步骤4-3：设定阈值t，并将矩阵T中高于阈值t的点标定为灰度变化剧烈点，即为造影微泡灌注的位置。

阈值t的确定是分割出血管灌注区域的重要参数，t与具体图像密切相关，本发明提出利用采样平均算法来确定t，具体步骤为：

步骤4-3-1：在医生指导下初步判断血管灌注区域，并选择三个样本区域作为观察对象，所选区域如图3所示；

步骤4-3-2：按照求取样本区域的像素点灰度变化的平均值，分别构建为m₁×n₁的矩阵T₁，m₂×n₂的矩阵T₂以及m₃×n₃的矩阵T₃，对应的m₁×n₁为第一样本区域的大小，m₂×n₂为第二样本区域的大小，m₃×n₃为第三样本区域的大小；

步骤4-3-3：按照

计算得到所述阈值t。

本例中根据所选的三个样本区域计算得到t＝7.6996，提取结果如图4所示，从图4可以看出，提取出来的点几乎都处于肿瘤内部，超声远场以及肿瘤周边的点都排除在灌注区域的之外。这是因为肌瘤内部血管丰富，在灌注达到峰值的时候血管越丰富的区域像素值变化的越大，阈值t是在肿瘤内部血供丰富区域采样得到的，反映的是血管中造影微泡在超声图像中的变化，肌瘤外部以及内部乏血供的区域像素灰度变化较小，基于此可以将血管灌注区域从周围组织背景中提取出来。图4a所示的为用二值图像表示提取结果，白色区域表示血管灌注区域，将该区域在原超声造影图像上用红色点标记出来，如图4b所示。

分别从来自校正后的和对灌注区域处理后的造影图像序列的ROI(图4b曲线中的区域)中提取TICs，用于评估提取策略的性能。录制的53s造影录像中，第20s起造影微泡开始进入肌瘤。根据先前提取出的血管灌注区域的像素点坐标，将20s至53s采集到的图像中的这些相同位置点的灰度值用前10s图像序列对应位置点的像素值代替，也就是保留灌注区域的基础像素值，消除造影微泡对像素灰度值的影响。提取出的TICs如图5所示。

从两组图像序列的ROI提取出的TICs可以看出，仅仅是做过运动校正的图像序列的ROI的TICs符合肌瘤灌注的过程，造影剂未进入子宫肌瘤时信号强度先保持在一个较低的水平，造影剂开始进入肿瘤后，由于造影微泡的作用，信号强度开始增加，达到峰值时期后保持一段时间，随后随着造影剂的流出，信号强度慢慢减弱。而对灌注区域处理后的图像序列ROI的TICs则始终保持在较低的水平。从图5可以看出，前100帧图像的两组TICs是重合的，这是因为此时造影剂还未进入肌瘤，从100帧往后造影剂开始进入肌瘤内部，对灌注区域处理后的图像序列ROI的TICs信号强度并未发生剧烈变化，仅有轻微的增强，表明被替换过的点是有造影微泡存在的点，即本方法准确的识别出了血管灌注区域。

最后需要说明的是，运动估计与校正过程中所采用的brox光流法以及相应的校正方法均是图像处理技术中的成熟算法，具体步骤未在本说明书中做详细描述，作为本领域普通技术人员，如有难以理解之处，可以参考以下文献：

[1]Brox T,Bruhn A,Papenberg N,et al.High accuracy optical flowestimation based on a theory for warping[M]//Computer Vision-ECCV2004.Springer Berlin Heidelberg,2004:25-36.

[2]Mémin E,Pérez P.Hierarchical estimation and segmentation of densemotion fields[J].International Journal of Computer Vision,2002,46(2):129-155.

Claims (4)

1.一种基于brox光流法的超声造影图像血管灌注区提取方法，其特征在于按照以下步骤进行：

步骤1：获取图像帧，形成超声造影图像序列；

步骤2：图像预处理，减少每帧超声造影图像的斑点噪声；

步骤3：基于brox光流法估计相邻两帧超声造影图像的运动场，并根据估计出的偏差量对相邻两帧超声造影图像进行校正；

步骤4：利用阈值分割法提取出血管灌注区域的点；

步骤4中的阈值分割法包括以下步骤：

步骤4-1：利用8邻域平均法对每帧图像的每个像素点进行平滑处理；

步骤4-2：按照计算相邻两帧图像间每个像素点灰度变化的平均值，并将所有像素点灰度变化的平均值构建成矩阵T；其中s为图像的帧序列，a为图像的总帧数，I(s)为第s帧图像对应像素点的灰度值；

步骤4-3：设定阈值t，并将矩阵T中高于阈值t的点标定为灰度变化剧烈点，即为造影微泡灌注的位置。

2.根据权利要求1所述的基于brox光流法的超声造影图像血管灌注区提取方法，其特征在于，步骤1中的具体步骤为：

步骤1-1：截取3～8秒灌注峰值期的录像；

步骤1-2：按照5帧/秒的采样率获取图像帧，并进行连续编号；

步骤1-3：选取编号为奇数的图像，重新编号形成新的图像序列。

3.根据权利要求1所述的基于brox光流法的超声造影图像血管灌注区提取方法，其特征在于，步骤2中的图像预处理采用低通高斯滤波器，高斯滤波窗口为3*3。

4.根据权利要求1所述的基于brox光流法的超声造影图像血管灌注区提取方法，其特征在于：所述阈值t按照以下方式设定：

步骤4-3-1：在医生指导下初步判断血管灌注区域，并选择三个样本区域作为观察对象；

步骤4-3-2：按照求取样本区域的像素点灰度变化的平均值，分别构建为m₁×n₁的矩阵T₁，m₂×n₂的矩阵T₂以及m₃×n₃的矩阵T₃，对应的m₁×n₁为第一样本区域的大小，m₂×n₂为第二样本区域的大小，m₃×n₃为第三样本区域的大小；

步骤4-3-3：按照

$t=\frac{\left(\underset{1}{\overset{m_1}{\Σ}}\underset{1}{\overset{n_1}{\Σ}}{T}_1\right)\frac{1}{m_1\×n_1}+\left(\underset{1}{\overset{m_2}{\Σ}}\underset{1}{\overset{n_2}{\Σ}}{T}_2\right)\frac{1}{m_2\×n_2}+\left(\underset{1}{\overset{m_3}{\Σ}}\underset{1}{\overset{n_3}{\Σ}}{T}_3\right)\frac{1}{m_3\×n_3}}{3}$

计算得到所述阈值t。