专利名称：一种心肌声学造影图像定量分析的方法
技术领域：
本发明涉及一种心肌声学造影图像定量分析的方法，尤其涉及一种基于光流场计算散斑 跟踪技术的心肌声学造影图像定量分析方法。
(二)
背景技术：
在心肌声学造影(Myocardial Contrast Echocardiography, MCE)图像定量分析的临 床应用中，人们经由MCE时间-强度曲线得到某心肌节段的冠脉微循环血流动力学参数，如 血容量、血流速度、血流量等供临床医师定量评价心肌灌注。MCE图像定量分析的一般原理 是测量某心肌节段的所选择感兴趣(Region of Interest, R0I)区域的取样点dB值。 这里，R01区域的中心称作取样点，ROI区域像素值的平均称作取样点dB值。国内外学者 以往的研究都证明超声切面内的微泡浓度随触发时间延长呈指数上升，其变化规律可
用Y二^(卜e-w) (^>0，">0)数学模型描述。模型中的参数^反映了心肌微血管的血
容量，参数/ 反映了心肌微血管的血流速度，i/ 则反映了心肌微血管的血流量。由所
选择ROI区域的取样点dB值测量数据，采用非线性回归分析方法，求解以上数学模型， 便可获得反映某心肌节段冠脉微循环各项临床指标的MCE时间-强度曲线。上述MCE图像 定量分析的基本前提是某心肌节段所选择的ROI区域，在存储的MCE图像序列不同帧 上，是完全相同的心肌区域。但在临床实践中，诸如超声探头位置或方位的细微变化、 患者体动、心脏跳动、呼吸及其他脏器蠕动之类的因素是不可能完全避免的，致使MCE 图像序列存在较严重的帧间失准问题，并最终导致所选择的R0I区域在MCE图像序列不 同帧上的位置发生"漂移"。若不考虑上述MCE图像序列的帧间失准问题，换言之，若 不跟踪所选择R0I区域在MCE图像序列不同帧上的位置变化，MCE图像的定量分析便会失 去准确性与客观性，分析结果的临床有效性也难以保证。
据査，目前已投入临床应用的MCE定量分析软件，大致分为两类。 一类如孙丰荣、张 明强等在《生物医学工程学杂志》(心肌造影超声心动图定量分析软件系统的设计与实 现，2005; —种新的基于MCE的心肌微循环定量分析系统的研制，2008)等学术刊物中 报道的MCE图像定量分析软件系统。该类软件在其MCE图像定量分析中，某心肌节段所 选择的R0I区域，在MCE图像序列不同帧上的位置是静止不变的。S卩，该类软件没有考 虑MCE图像序列的帧间失准所导致的R0I区域在MCE图像序列不同帧上的位置变化。显 然，应用该类软件对MCE图像进行定量分析，其准确性和临床有效性有待商榷。另外的 一类MCE定量分析软件则采用手动或交互式处理方法以应对所选择R0I区域在MCE图像序 列不同帧上的位置变化，比如GE公司的Vivid7型超声设备，该设备在超声序列图像播放 过程中，允许观察者对所选择R0I区域位置"漂移"比较明显的帧手动改变所选择ROI 区域的位置，这种人工调整所选择ROI区域位置的方法操作繁琐、耗时，MCE图像定量 分析的准确性很大程度上取决于临床医师的经验和知识。另据査，目前为止，也有学者 关注MCE图像序列的帧间失准所导致的R0I区域在MCE图像序列不同帧上的位置变化问 题，如M. L- Lope博士 (2003)在文章/卿r。w'/7g7fcy。carc/jW C。"tra" ^b力ocaro7ograp/ y 7/Hages y egistraWo/7 i/7crea57'/7g i/asges 5Y肌7arj'ty中，综合基于分害U与基于像素属 性的图像配准技术提出的一种半自动的MCE超声序列图像的弹性配准算法，该算法的配准精度和有效性取决于分割阶段人工选择的模板，无法实现配准过程的自动化，且计算效率 和鲁棒性不理想。
发明内容
针对背景技术中所述的(1)或者不考虑MCE图像序列的帧间失准所导致的ROI区域在
MCE图像序列不同帧上的位置变化问题，致使分析结果的准确性和临床有效性难以保证，(2) 或者考虑该问题但依赖于人工调整，致使分析过程缺乏自动化、分析结果缺乏客观性，(3) 或者考虑该问题并采用序列图像配准的解决问题思路，但方法大都准确性差、鲁棒性低，本 发明提出一种基于光流场计算散斑跟踪技术的MCE图像定量分析方法，该方法通过计算MCE 图像序列的光流场，进而获得所选择ROI区域的位移，通过位移准确定位所选择ROI区域在 MCE图像序列每一帧中的位置，继而由所选择ROI区域的取样点dB值测量数据，采用非线 性回归分析方法，求解相关数学模型并绘制MCE时间-强度曲线。 本发明的技术方案如下
一种基于光流场计算散斑跟踪技术的MCE定量分析方法,步骤如下
51) 确定光流计算中的能量方程
根据光流约束方程、梯度约束方程与平滑约束条件得到总的能量方程，其公式如下
五(",v)二五o她+ o^s謂^
v) =+『)—/(x)|2 + y| v/(x +『)—v/(x)|2)a
v) = J。(|V3"|2 +1 V3v|2) 其中Z:二(x,;Mf ,『:=(",v,lf ， y是权重系数，V为梯度算子，V33/为时空
梯度算子，a —O亦为权重系数。五O,v)为总的能量方程，五SM。。,""，力为平滑约束项，
^z^a(w,"为光流约束方程与梯度约束方程结合项；
52) 为MCE图像序列的各帧构造其金字塔表示
对MCE图像序列中的每一帧图像都构造一个金字塔表示，金字塔的层数可由操作者设 定；MCE图像序列中的一帧图像作为该帧图像金字塔表示中的最低层，图像层数设定为零， 然后对最低层图像进行平滑下采样，其结果作为金字塔表示中的倒数第二层，图像层数设定 为一，以此类推，分解到金字塔表示预先设定的层数为止， 一帧MCE图像的金字塔表示便构 造完毕，在金子塔表示中较高的层是它相邻下面层的平滑后下采样形式；
53) 采用由粗到精的搜索策略计算图像序列的光流场
计算两幅图像之间的光流场时，按照图像的金字塔表示，采取由高层到低层的形式进 行，即从两幅图像的金字塔表示中的最高层(设最高层层数为M)开始计算，设定与M层相 对应的光流初始值为零，计算两幅图像金字塔表示中相对应的M层图像之间的光流增量值
Ap^，将其与该层的光流初始值相加，结果为Ap^，作为金字塔表示中M-l层图像之间光
流场的初始值；再计算两幅图像金字塔表示中相对应的M-1层图像之间的光流增量值.
Arfi，将其与该层的光流初始值相加，结果为(么^^一+ AT")作为M-2层图像之间光流
5场的初始值，以此类推，直到计算到最低层为止，即可估计出原始图像的光流；
54) 选择某心肌节段的ROI区域
由用户选定初始帧，并在其上手动选取某个特定心肌节段的ROI区域作为散斑跟踪 的目标；
55) 确定该心肌节段所选择R01区域在MCE图像序列每一帧中的位移
由S3)得到的图像序列的光流场，计算所选择ROI区域在MCE图像序列每一帧中的位 移，从而实现对所选择ROI区域的跟踪； -
56) 确定该心肌节段所选择R01区域在MCE图像序列每一帧中的位置，并测量其取样点 dB值
根据所选择ROI区域在初始帧中的位置及其在MCE图像序列每一帧中的位移，得到所选 择ROI区域在MCE图像序列每一帧中的准确位置，进而计算其取样点dB值；
57) 确定MCE定量分析数学模型中的未知参数
采用非线性最小二乘回归方法计算数学模型中的未知参数，其中数学模型表示如下
其中参数A反映出心肌血容量，参数/ 反映出心肌血流平均速度，A./ 综合反映出心肌血
流量； 、
58) 绘制MCE时间-强度曲线
依据上述S7)得出的结果，绘制MCE时间-强度曲线。
上述步骤S4)中所述散斑是超声图像的固有属性，是一种结构性信号，反应人体组织器 官的局部解剖结构特性。散斑模式具有唯一性和稳定性，据此可以实现对人体组织器官(或 其感兴趣区域)运动的跟踪。
上述MCE是英文Myocardial Contrast Echocardiography的縮写，即心肌声学造影。ROI 是英文Region of Interest的縮写，即感兴趣区域。
一种使用本发明方法的MCE定量分析系统，主要包括有影像源(超声诊断仪)、视频传 输装置(光端机)、视频采集装置(图像采集卡)、工作站PC机。 本发明MCE定量分析系统的工作过程如下
101. 工作站PC机采集并以DIC0M格式存储心肌显影稳定后的MCE图像序列；
102. 打开己存储的MCE图像序列文件，可自动连续播放或者手动逐帧顺序播放该序列图
像；
103. 确定光流计算中的能量方程
根据光流约束方程、梯度约束方程与平滑约束条件得到总的能量方程，其公式如下 五(",v)二五zto" + a五w。础
五一("，v) = Hl"i +『)一 ""l2+H v/(z +『)—w(i)「)血
其中Z^(x,;;力、『:=(",v,l)T， Z是权重系数，V为梯度算子，v3:= (&, ,&f为时空梯度算子，《 —0亦为权重系数。五("，力为总的能量方程，E^。。^0,v)为平滑约束项，
fz^a(",力为光流约束方程与梯度约束方程结合项；
104. 为MCE图像序列的各帧构造其金字塔表示
对MCE图像序列中的每一帧图像都构造一个金字塔表示，金字塔的层数可由操作者设 定。MCE图像序列中的一帧图像作为该帧图像金字塔表示中的最低层，图像层数设定为零， 然后对最低层图像进行平滑下采样，其结果作为金字塔表示中的倒数第二层，图像层数设定 为一，依此类推，分解到金字塔表示预先设定的层数为止， 一帧MCE图像的金字塔表示便构 造完毕，在金子塔表示中较高的层是它相邻下面层的平滑后下采样形式；
105. 采用由粗到精的搜索策略计算图像序列的光流场
计算两幅图像之间的光流场时，按照图像的金字塔表示，釆取由高层到低层的形式进 行，即从两幅图像的金字塔表示中的最高层(设最高层层数为M)开始计算，设定与M层相 对应的光流初始值为零，计算两幅图像金字塔表示中相对应的M层图像之间的光流增量值
△f^，将其与该层的光流初始值相加，结果为Af^，作为金字塔表示中M-l层图像之间光
流场的初始值。再计算两幅图像金字塔表示中相对应的M-1层图像之间的光流增量值 Arf、将其与该层的光流初始值相加，结果为4 A{/M)作为M-2层图像之间光流
场的初始值，以此类推，直到计算到最低层为止，即可估计出原始图像的光流；
106. 选择某心肌节段的ROI区域
在MCE图像序列播放暂停状态下，由用户选定初始帧，并在其上手动选取某个特定心肌 节段的ROI区域作为散斑跟踪目标，ROI区域可以是矩形、椭圆或者手绘任意形状的区域；
107. 确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位移 由得到的图像序列的光流场，计算所选择ROI区域在MCE图像序列每一帧中的位移，
从而实现对所选择ROI区域的跟踪；
108. 确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位置 根据操作者所选择ROI区域在初始帧上的位置及其在MCE图像序列每一帧中的位移，获
得所选择ROI区域在MCE图像序列每一帧中的准确位置；
109. 测量该心肌节段所选择ROI区域的取样点dB值
通过获取的心肌节段所选择ROI区域在MCE图像序列每一帧中的位置测量该心肌节段 所选择ROI区域的取样点dB值，将ROI区域的中心视为取样点，ROI区域像素值的平均作 为取样点dB值；
110. 确定MCE定量分析数学模型中的未知参数
确定所选择ROI区域在MCE图像序列每一帧中的位置后，采用非线性最小二乘回归方法 计算数学模型中的未知参数，其中数学模型表示如下
7 =邓-￡—,(J>0，/ >0)
其中参数A反映出心肌血容量，参数yff反映出心肌血流平均速度，A^综合反映出心肌血
流量；
7111.绘制MCE时间-强度曲线
依据上述110)得出的结果，绘制MCE时间-强度曲线。
本发明方法基于光流场计算散斑跟踪技术完成MCE图像的定量分析，它首先计算得到 MCE图像序列的光流场，选择要跟踪的ROI区域之后，便根据MCE图像序列的光流场确定所 选择ROI区域在MCE图像序列每一帧中的位置，继而由所选择ROI区域的取样点dB值测量 数据，采用非线性回归分析方法，求解相关数学模型并绘制MCE时间-强度曲线，这实现了 MCE图像定量分析的全自动化，分析结果的准确，分析方法临床有效性和客观性髙。
(四)
图1是本发明计算方法的流程图。
图2是步骤S2中的图像序列金字塔表示。其中imageO为原始图像，image 1为image0 的平滑下采样形式，vZ值(li, 1，2)为金字塔中每一层图像相对应的光流值。
图3是步骤S3中的由粗到精搜索策略计算图像序列的光流场的示意图。其中V2 = 0表
示光流的初始值，其值为0，光流计算之后可以得到对应此层图像的光流增量AV2，将这一
层的初始值与光流增量值相加，得到的结果作为下一层的初始值，以此类推。
具体实施方式
下面结合附图和实施例对本发明做进一步说明，但不限于此。 实施例
一种基于光流场计算散斑跟踪技术的MCE定量分析方法,步骤如下
51) 确定光流计算中的能量方程
根据光流约束方程、梯度约束方程与平滑约束条件得到总的能量方程，其公式如下
v) =+『)-/(义)|2+H v/(x+『)—V/(义)|2)血
￡^"(",力=_f"(|v3"|2 + |v3v|2)
其中Z:-(x,少力。^是权重系数，V为梯度算子，V3^(3x,3^，^)T为时空 梯度算子，a》0亦为权重系数。五("，v)为总的能量方程，五&。。说(",V)为平滑约束项，
五"。to(w,力为光流约束方程与梯度约束方程结合项；
52) 为MCE图像序列的各帧构造其金字塔表示
对MCE图像序列中的每一帧图像都构造一个金字塔表示，金字塔的层数可由操作者设 定；MCE图像序列中的一帧图像作为该帧图像金字塔表示中的最低层，图像层数设定为零， 然后对最低层图像进行平滑下采样，其结果作为金字塔表示中的倒数第二层，图像层数设定 为一，以此类推，分解到金字塔表示预先设定的层数为止， 一帧MCE图像的金字塔表示便构 造完毕，在金子塔表示中较高的层是它相邻下面层的平滑后下采样形式；
53) 采用由粗到精的搜索策略计算图像序列的光流场计算两幅图像之间的光流场时，按照图像的金字塔表示，采取由高层到低层的形式进 行，即从两幅图像的金字塔表示中的最高层(设最高层层数为M)开始计算，设定与M层相 对应的光流初始值为零，计算两幅图像金字塔表示中相对应的M层图像之间的光流增量值
△FM，将其与该层的光流初始值相加，结果为AFM，作为金字塔表示中M-l层图像之间光
流场的初始值；再计算两幅图像金字塔表示中相对应的M-1层图像之间的光流增量值
AFM_'，将其与该层的光流初始值相加，结果为(Af^—^ 作为M-2层图像之间光流
场的初始值，以此类推，直到计算到最低层为止，即可估计出原始图像的光流；
54) 选择某心肌节段的ROI区域
由用户选定初始帧，并在其上手动选取某个特定心肌节段的ROI区域作为散斑跟踪 的目标；
55) 确定该心肌节段所选择R01区域在MCE图像序列每一帧中的位移
由S3)得到的图像序列的光流场，计算所选择ROI区域在MCE图像序列每一帧中的位 移，从而实现对所选择ROI区域的跟踪；
56) 确定该心肌节段所选择R01区域在MCE图像序列每一帧中的位置，并测量其取样点 dB值
根据所选择ROI区域在初始帧中的位置及其在MCE图像序列每一帧中的位移，得到所选 择ROI区域在心肌图像序列每一帧中的准确位置,进而计算其取样点dB值；
57) 确定MCE定量分析数学模型中的未知参数
采用非线性最小二乘回归方法计算数学模型中的未知参数，其中数学模型表示如下 F = j(l-e一,(h0》0)
其中参数A反映出心肌血容量，参数/ 反映出心肌血流平均速度，A.A综合反映出心肌血
流量参数；
58) 绘制MCE时间-强度曲线
依据上述S7)得出的结果，绘制MCE时间-强度曲线。
权利要求
1、一种基于光流场计算散斑跟踪技术的MCE定量分析方法，步骤如下S1)确定光流计算中的能量方程根据光流约束方程、梯度约束方程与平滑约束条件得到总的能量方程，其公式如下E(u，v)＝EData+αESmooth<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msub> <mi>E</mi> <mi>Data</mi></msub><mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo></mrow><mo>=</mo><msub> <mo>∫</mo> <mi>Ω</mi></msub><mrow> <mo>(</mo> <mo>|</mo> <mi>I</mi> <mrow><mo>(</mo><mi>X</mi><mo>+</mo><mi>W</mi><mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow><mo>(</mo><mi>X</mi><mo>)</mo> </mrow> <msup><mo>|</mo><mn>2</mn> </msup> <mo>+</mo> <mi>γ</mi> <msup><mrow> <mo>|</mo> <mo>▿</mo> <mi>I</mi> <mrow><mo>(</mo><mi>X</mi><mo>+</mo><mi>W</mi><mo>)</mo> </mrow> <mo>-</mo> <mo>▿</mo> <mi>I</mi> <mrow><mo>(</mo><mi>X</mi><mo>)</mo> </mrow> <mo>|</mo></mrow><mn>2</mn> </msup> <mo>)</mo></mrow><mi>dx</mi> </mrow>]]></math></maths><maths id="math0002" num="0002" ><math><![CDATA[ <mrow><msub> <mi>E</mi> <mi>Smooth</mi></msub><mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo></mrow><mo>=</mo><msub> <mo>∫</mo> <mi>Ω</mi></msub><mrow> <mo>(</mo> <msup><mrow> <mo>|</mo> <msub><mo>▿</mo><mn>3</mn> </msub> <mi>u</mi> <mo>|</mo></mrow><mn>2</mn> </msup> <mo>+</mo> <msup><mrow> <mo>|</mo> <msub><mo>▿</mo><mn>3</mn> </msub> <mi>v</mi> <mo>|</mo></mrow><mn>2</mn> </msup> <mo>)</mo></mrow> </mrow>]]></math></maths>其中X＝(x，y，t)T，W＝(u，v，1)T，γ是权重系数， id="icf0003" file="A2009100207450002C3.tif" wi="2" he="3" top= "81" left = "113" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>为梯度算子，<maths id="math0003" num="0003" ><math><![CDATA[ <mrow><msub> <mo>▿</mo> <mn>3</mn></msub><mo>:</mo><mo>=</mo><msup> <mrow><mo>(</mo><msub> <mo>∂</mo> <mi>x</mi></msub><mo>,</mo><msub> <mo>∂</mo> <mi>y</mi></msub><mo>,</mo><msub> <mo>∂</mo> <mi>t</mi></msub><mo>)</mo> </mrow> <mi>T</mi></msup> </mrow>]]></math> id="icf0004" file="A2009100207450002C4.tif" wi="31" he="6" top= "80" left = "143" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>为时空梯度算子，α＞0亦为权重系数。E(u，v)为总的能量方程，ESmooth(u，v)为平滑约束项，EData(u，v)为光流约束方程与梯度约束方程结合项；S2)为MCE图像序列的各帧构造其金字塔表示对MCE图像序列中的每一帧图像都构造一个金字塔表示，金字塔的层数可由操作者设定；MCE图像序列中的一帧图像作为该帧图像金字塔表示中的最低层，图像层数设定为零，然后对最低层图像进行平滑下采样，其结果作为金字塔表示中的倒数第二层，图像层数设定为一，以此类推，分解到金字塔表示预先设定的层数为止，一帧MCE图像的金字塔表示便构造完毕，在金子塔表示中较高的层是它相邻下面层的平滑后下采样形式；S3)采用由粗到精的搜索策略计算图像序列的光流场计算两幅图像之间的光流场时，按照图像的金字塔表示，采取由高层到低层的形式进行，即从两幅图像的金字塔表示中的最高层开始计算，设定与M层相对应的光流初始值为零，计算两幅图像金字塔表示中相对应的M层图像之间的光流增量值ΔVM，将其与该层的光流初始值相加，结果为ΔVM，作为金字塔表示中M-1层图像之间光流场的初始值；再计算两幅图像金字塔表示中相对应的M-1层图像之间的光流增量值ΔVM-1，将其与该层的光流初始值相加，结果为(ΔVM-1+ΔVM)作为M-2层图像之间光流场的初始值，以此类推，直到计算到最低层为止，即可估计出原始图像的光流；S4)选择某心肌节段的ROI区域由用户选定初始帧，并在其上手动选取某个特定心肌节段的ROI区域作为散斑跟踪的目标；S5)确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位移由S3)得到的图像序列的光流场，计算所选择ROI区域在MCE图像序列每一帧中的位移，从而实现对所选择ROI区域的跟踪；S6)确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位置，并测量其取样点dB值根据所选择ROI区域在初始帧中的位置及其在MCE图像序列每一帧中的位移，得到所选择ROI区域在MCE图像序列每一帧中的准确位置，进而计算其取样点dB值；S7)确定MCE定量分析数学模型中的未知参数采用非线性最小二乘回归方法计算数学模型中的未知参数，其中数学模型表示如下Y＝A(1-e-βt)(A＞0，β＞0)其中参数A反映出心肌血容量，参数β反映出心肌血流平均速度，A·β综合反映出心肌血流量；S8)绘制MCE时间-强度曲线依据上述S7)得出的结果，绘制MCE时间-强度曲线。
全文摘要
一种心肌声学造影(MCE)图像定量分析的方法，属医学图像分析处理技术领域，本发明方法步骤为确定光流计算中的能量方程；为MCE图像序列的各帧构造其金字塔表示；采用由粗到精的搜索策略计算图像序列的光流场；选择某心肌节段的ROI区域；确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位移；确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位置，并测量其取样点dB值；确定MCE定量分析数学模型中的未知参数；绘制MCE时间-强度曲线。本发明实现了MCE图像定量分析的全自动化，分析结果准确，分析方法的临床有效性和客观性高。
文档编号A61B8/06GK101536919SQ20091002074
公开日2009年9月23日 申请日期2009年4月27日 优先权日2009年4月27日
发明者孙丰荣, 张明强, 王晓婧, 贾晓波 申请人:山东大学