CN102663709A - 一种x射线冠脉造影图像血管增强方法 - Google Patents

Abstract

本发明公开了一种X射线冠脉造影图像血管增强方法，包括以下步骤：对X射线冠脉造影图像I进行二维经验模态分解；分析各层经验模态分解的本征模态函数进行去噪，并利用能够代表冠脉造影图像的特征层来加权重构去除背景噪声的图像；对重构的图像进行Hessian矩阵分解，利用血管测度函数或血管特征拟合函数达到增强血管结构的目的。本发明利用经验模态分解X射线冠脉造影图像为一系列本征模态函数，利用噪声的分布规律达到去除背景噪声的目的，进而选择一些特定层本征模态函数进行加权构造冠状动脉图像，利用基于Hessian矩阵特征值的血管测度函数进一步增强图像中的血管结构。本发明能明显改善冠脉造影图像的视觉效果。

Description

一种X射线冠脉造影图像血管增强方法

技术领域

本发明属于医学图像处理、计算机医学辅助诊疗系统、数字信号处理、计算机人工智能及模式识别技术领域，特别涉及利用计算机对噪声条件下的X射线医学冠状动脉造影图像进行增强的方法。

背景技术

X射线冠状动脉造影(X射线冠脉造影图像：X-Ray Angiogram)检查是目前国际上公认的诊断冠心病的最佳手段，被称为冠心病诊断的“金标准”。由于成像机理原因，医学图像往往存在时间/空间分辨率与信噪比的矛盾。X射线冠脉造影成像时，X线经人体衰减后得到的图像包含了关于人体组织的各种相关信息，在X线透射、光的转换、视频信号的形成及A/D转换等过程中也引入了各种噪声，如：X线的散射噪声、影像增强器的量子噪声、摄像管及A/D转换等引起的系统噪声。而成像又受到各种实际因素的影响，如患者的舒适度、系统的要求等等，需要快速成像。快速成像会导致图像低的信噪比和对比度。在提高图像的时间分辨率的同时，图像中的噪声会大大降低图像的质量，使一些组织边界变得模糊，细微结构难以辨认，加大了对图像细节识别和分析的难度，影响医学诊断。

噪声的存在使获得的图像不清晰，尤其是掩盖和降低了造影图像中某些特征细节的可见度。可见度的损失对低对比度物体尤为重要，对以后的血管图像识别更为不利。因此需要对X射线冠脉造影图像进行增强，突出血管的线条结构以利于后期的进一步分析与处理。

图像增强采用各种数字处理技术改善图像的视觉效果，使图像结果更适合于观看，提高图像视觉质量及其可诊断性。对医学X射线冠脉造影图像进行增强，就是要压制(抑制)图像中的噪声同时突出冠状动脉血管结构，得到对医务人员来说视觉效果更好、更易于诊断的图像。

人们已经对图像去噪和增强研究了几十年，主要可分为两大类：传统的低通去噪方法(包括均值滤波、高斯平滑等)；非线性方法(中值滤波、小波变换、各向异性扩散方程、数学形态学等)。但这些方法在对X射线冠脉造影图像进行去噪时效果都不理想。X射线冠脉造影图像(简称XRA)具有如下特点：1)血管的形状复杂，而且容易发生变形；2)血管的密度和血管的直径变化很大，导致血管灰度的变化范围也很大；3)图像中的背景噪声比较高(由于不同组织对X射线的衰减系数不相同造成图像背景的灰度分布不均匀以及各种散射噪声、量化噪声等形成了大量背景噪声)。传统的低通去噪方法在进行X射线冠脉造影图像(XRA：X-Ray Angiogram)降噪会破坏图像中的重要组织(主要指血管)结构信息，使一些组织边界变得模糊，细微结构难以辨认，加大了对图像细节识别和分析的难度，影响医学诊断。基于小波变换去噪方法的效果依靠小波基的构造，而合适的小波基构造方法现在还没有一个成型的理论，本身就是一个难题。基于数学形态学的方法，去噪结果和结构元素的大小和形状密切相关。各向异性扩散方程在去噪过程中存在扩散尺度参数选择以及扩散终止条件确定问题，直接影响了去噪效果。这些方法最明显的缺点就是都需要冠脉造影图像相关的先验知识，根据不同的图像需要选择合理的参数才能处理，因此降噪能力受到很大限制。

经验模态分解(Empirical Mode Decomposition，EMD)是近年来引起人们关注的一种新的时频分析方法，适合于分析非线性、非平稳信号序列，具有很高的信噪比。它依据数据自身的时间尺度特征来进行信号分解，无须预先设定任何基函数。这一点与建立在先验性的谐波基函数和小波基函数上的傅里叶分解与小波分解方法具有本质性的差别。经验模态分解将复杂信号分解为有限个(层)本征模态函数(Intrinsic Mode Function，简称IMF)，所分解出来的各IMF分量包含了原信号的不同时间尺度的局部特征信号。因为基函数(IMF分量)是由数据本身所分解得到，同时由于分解是基于信号序列时间尺度的局部特性，所以与短时傅立叶变换、小波分解和各向异性扩散方程等方法相比，这种方法具有自适应性。本发明利用经验模态分解方法对X射线冠脉造影图像进行自适应分解，对得到的各IMF分量进行去噪，并选择能够代表血管结构的IMF高频分量进行造影图像重构即可达到去除噪声的目的，且在去噪过程中无须先验知识、能够根据造影图像特点进行自适应降噪。

二维经验模态分解可将一幅图像I分解为n层本征模态函数和一个残余的趋势分量，表示为：

$I=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i+R$

式中F_i为第i层本征模态函数，R为残余的趋势分量。

发明内容

为解决现有技术存在的上述问题，本发明要提出一种具有自适应降噪且具有血管结构增强的X射线冠脉造影图像增强方法。

为了实现上述目的，本发明的技术方案如下：一种X射线冠脉造影图像血管增强方法，包括以下步骤：

A、对X射线冠脉造影图像I进行二维经验模态分解

A1、利用输入的X射线冠脉造影图像I作为分解初始值，设定当前层为i＝1，第i层本征模态函数的估计为r_i＝I；

A2、令第i层本征模态函数的第j次估计量为h_j＝r_i，j＝1；

A3、求取h_j的局部极大值和极小值；

A4、对所获得的极大和极小值分别用三次样条插值得到上包络env_max和下包络env_min；

A5、计算上包络env_max和下包络env_min的均值，

A6、从第i层本征模态函数的第j次估计量h_j中减去m，h_j+1＝h_j-m；

A7、令第i层本征模态函数的估计次数j＝j+1；

A8、判断：如果||h_j‖＜ε或j＞＝4不成立，则转到步骤A3；如果||h_j‖＜ε或j＞＝4成立，则获得X射线冠脉造影图像I的第i层本征模态函数F_i＝h_j；

A9、得到第i+1层本征模态函数的初始估计

A10、令本征模态函数的层序数i＝i+1；

A11、如果r_i为单调函数，则R＝r_i为残余分量，分解过程结束；否则转到步骤A2；

B、分析各层经验模态分解的本征模态函数进行去噪，并利用能够代表冠脉造影图像的特征层来加权重构去除背景噪声的图像

B1、分析X射线冠脉造影图像I的各层本征模态函数F_i，计算强度直方图得到强度分布点数最大的强度值hist_i；

B2、计算各层本征模态函数F_i的强度方差σ_i；

B3、将方差σ_i或强度值hist_i的最大者作为第i层本征模态函数的阈值Th_i；

B4、对每一层本征模态函数利用阈值Th_i将小于阈值的点置0；

B5、设定进行重构的本征模态函数层数为M_ite＝3，迭代初始值ite＝1；

B6、设定过程变量temp为0；

B7、将去噪后的第ite层本征模态函数F_i取绝对值abs(F_ite)；

B8、过程变量

其中abs(·)为取绝对值操作；

B9、令ite＝ite+1；

B10、如果ite＜(M_ite+1)，则转到B8；否则重构过程结束，输出重构冠脉造影图像I_re＝temp；

C、对重构的图像进行Hessian矩阵分解，利用血管测度函数或血管特征拟合函数达到增强血管结构的目的

C1、构造多尺度高斯滤波器集

σ_n为高斯滤波器的方差，σ_n∈[σ₁，...，σ_N]；

C2、对多尺度高斯滤波器集

求二阶导数，用

的二阶导数与重构后的冠脉造影图像I_re进行卷积得到

和

C3、利用图像I_re的二阶导数计算Hessian矩阵的特征值；

C4、若Hessian矩阵的两个特征值|λ₁|＞|λ₂|，利用血管测度函数或血管特征拟合函数计算血管特征图

C5、得到血管特征图

C6、输出血管结构增强结果图

与现有技术相比，本发明具有以下有益效果：

1、本发明不依赖先验知识进行对冠脉造影图像自适应分解

经验模态分解EMD时频分析是最新发展起来的处理非线性非平稳信号的时频分析方法。它既吸取了小波变换多分辨的优点，又克服了小波变换中选择小波基的困难。经验模分解是一种完全由数据驱动的分解方法，不需要预设任何滤波器或小波函数。因此，可用来对非平稳信号进行滤波和去噪。且该法由信号本身的特征尺度出发，将信号分解为有限个具有不同特征尺度的本征模态函数(Intrinsic Mode Function，IMF)，利用这些本征模态函数可以分析冠脉造影图像信号能量在时间、空间等各种尺度上的分布规律。

2、本发明根据图像特点自适应产生随机噪声去除阈值

X射线冠脉造影图像中的噪声主要包含两部分：1)X射线散射、量化噪声等形成的随机背景噪声，主要分布在高频区域；2)不同组织(如肺等)对X射线的衰减系数不同造成在图像中区域灰度分布不均匀形成结构噪声，主要分布在低频区域。和小波变换类似，经验模态分解可以看作是一个滤波器组，而且其产生的各阶本征模态函数是从高频到低频逐渐展开的过程，即第一层本征模态函数包含图像中的高频信息，最后一层本征模态函数包含图像中的最低频率信息(或者说是趋势信息)。在X射线冠脉造影图像中冠状动脉血管是进行分析的目标，血管的边缘、结构等重要信息分布在高频区域。可利用阈值在各本征模态函数内去除随机噪声。对经验模态分解出的各层本征模态函数计算强度方差σ_i和强度直方图中强度分布点数最大的强度值hist_i，利用σ_i和hist_i构造去噪阈值，这一阈值只与图像中的含噪量有关(受信噪比影响)，这一去噪过程是一种自适应的过程。

3、本发明利用加权系数对重要的图像血管结构进行加重并去除结构噪声

在X射线冠脉造影图像中冠状动脉血管是进行分析的目标，血管的边缘、结构等重要信息分布在高频区域。冠脉造影图像中的结构噪声在低频区域。因此将低频区域本征模态函数丢弃并利用加权系数对高频区域的本征模态函数进行冠脉造影图像重构，不仅可以去除结构噪声而且可以加强血管结构的视觉效果。

4、本发明利用血管的解剖结构数学模型构造血管测度函数来加强血管结构

血管的解剖结构为管状，经X射线成像后血管区域灰度近似符合高斯函数分布规律，因此血管图像区域的导函数构成的Hessian矩阵的两个特征值λ₁和λ₂存在|λ₁|＞＞|λ₂|，利用特征值这种关系构造血管测度函数就可以将图像中的血管区域进一步加强，突出血管区域的视觉效果达到增强血管结构的目的，为进一步的血管图像识别与量化等处理打下良好基础。

5、本发明能够根据X射线冠脉造影图像的特点而自适应调整，具有很好抗噪性和血管增强效果。这种方法的原理主要是：1)利用经验模态分解X射线冠脉造影图像为一系列本征模态函数，而噪声在各层本征模态函数上的分布具有一定的规律，因此可以利用噪声的这种分布规律达到去除背景噪声的目的，进而用与图像的重要结构相关的本征模态函数对图像进行加权重构，因此可以选择一些特定层(本征模态函数)进行加权构造冠状动脉图像；2)对上一步去除噪声并重构的图像，针对理想血管形态为管状，利用基于Hessian矩阵特征值的血管测度函数进一步增强图像中的血管结构。本方法能明显改善冠脉造影图像的视觉效果

附图说明

本发明共有附图3张，其中：

图1是基于经验模态分解冠脉造影图像增强方法的流程图。

图2是利用本征模态函数重构去除背景噪声图像的流程图。

图3是基于Hessian矩阵的血管测度函数增强图像的流程图。

具体实施方式

下面参考附图，详细描述本发明提出的X射线冠脉造影图像冠脉造影图像增强方法。

如图1-3所示，本方法包括三个分离的阶段：(1)对X射线冠脉造影图像进行二维经验模态分解；(2)分析各层经验模态分解的本征模态函数(intrinsicmode functions：IMF)进行去噪，并利用能够代表冠脉造影图像的特征层(代表X射线冠脉造影图像血管结构的本征模态函数)来重构去除背景噪声的图像；(3)对重构的图像进行Hessian矩阵分解，利用血管测度函数(拟合函数)达到增强血管结构的目的，如图3所示。

A、第一阶段利用经验模态分解(the bidimensional empirical modedecomposition：BEMD)技术来分解X射线冠脉造影图像I(x，y)，得到物体的各层本征模态函数F_i(intrinsic mode functions：IMF)和残余分量R，包含以下步骤：

A1、利用输入的X射线冠脉造影图像I作为分解初始值，设定当前层为i＝1，第i层本征模态函数的估计为r_i＝I；

A2、令第i层本征模态函数的第j次估计量为h_j＝r_i，j＝1；

A3、求取h_j的局部极大值和极小值；

A4、对所获得的极大和极小值分别用三次样条插值得到上包络env_max和下包络env_min；

A5、计算上包络env_max和下包络env_min的均值，

A6、从第i层本征模态函数的第j次估计量h_j中减去m，h_j+1＝h_j-m；

A7、令第i层本征模态函数的估计次数j＝j+1；

A8、判断如果||h_j‖＜ε，ε＝1.5或j＞＝4不成立，则转到步骤A3；如果||h_j‖＜ε或j＞＝4成立，则获得X射线冠脉造影图像I的第i层本征模态函数F_i＝h_j；

A9、得到第i+1层本征模态函数的初始估计

A10、令本征模态函数的层序数i＝i+1；

A11、如果r_i为单调函数，则R＝r_i为残余分量，分解过程结束；否则转到步骤A2；

B、第二阶段对本征模态函数去噪与冠脉造影图像重构包含以下步骤：

B1、分析X射线冠脉造影图像I的各层本征模态函数F_i，计算强度直方图得到强度分布点数最大的强度值hist_i；

B2、计算各层本征模态函数F_i的强度方差σ_i；

${\mathrm{\σ}}_i=\sqrt{\frac{1}{p\×q}\underset{(x,y)\∈F_i}{\mathrm{\Σ}}{[F_i(x,y)-{\stackrel{\‾}{F}}_i(x,y)]}^2}$

${\stackrel{\‾}{F}}_i(x,y)=\frac{1}{p\×q}\underset{x=1}{\overset{p}{\mathrm{\Σ}}}\underset{y=1}{\overset{q}{\mathrm{\Σ}}}\left[F_i(x,y)\right]$

图像尺寸为p×q。

B3、将方差σ_i或强度值hist_i的最大者作为第i层本征模态函数的阈值Th_i；

${\mathrm{Th}}_i=\left\{\begin{array}{cc}{\mathrm{\&sigma;}}_i& \mathrm{if}{\mathrm{\&sigma;}}_i\&GreaterEqual;{\mathrm{hist}}_i\\ {\mathrm{hist}}_i& \mathrm{if}{\mathrm{\&sigma;}}_i<{\mathrm{hist}}_i\end{array}\right.$

B4、对每一层本征模态函数利用阈值Th_i将小于阈值的点置0；

$F_i=\left\{\begin{array}{cc}{F}_i& \mathrm{if}{F}_i\&GreaterEqual;{\mathrm{Th}}_i\\ 0& \mathrm{if}{F}_i<{\mathrm{Th}}_i\end{array}\right.$

B5、设定进行重构的本征模态函数层数为M_ite＝3，迭代初始值ite＝1；

B6、设定过程变量temp为0；

B7、将去噪后的第ite层本征模态函数F_i取绝对值abs(F_ite)；

B8、过程变量

其中abs(·)为取绝对值操作；

B9、令ite＝ite+1；

B10、如果ite＜(M_ite+1)，则转到B8；否则重构过程结束，输出重构冠脉造影图像I_re＝temp。

$I_{\mathrm{re}}=\underset{\mathrm{ite}=1}{\overset{3}{\mathrm{\&Sigma;}}}{5e}^{-\frac{{(\mathrm{ite}-1)}^2}{4}}\mathrm{abs}\left[F_{\mathrm{ite}}(x,y)\right]$

C、对重构的图像进行Hessian矩阵分解，利用血管测度函数或血管特征拟合函数达到增强血管结构的目的，包含以下步骤：

C1、构造多尺度高斯滤波器集σ_n为高斯滤波器的方差，σ_n∈[σ₁，...，σ_N]；

${\mathrm{\&sigma;}}_n\&Element;[1,\sqrt{2},\mathrm{2,2}\sqrt{2},\mathrm{4,4}\sqrt{2},\mathrm{8,8}\sqrt{2},16];$

$g_{{\mathrm{\&sigma;}}_n}(x,y)=\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_n^2}\exp [-\frac{(x^2+y^2)}{2{\mathrm{\&sigma;}}_n^2}]$

C2、对多尺度高斯滤波器集

求二阶导数，用

的二阶导数与重构后的冠脉造影图像I_re进行卷积得到和

$I_{\mathrm{xx}}^{\mathrm{re}}\left({\mathrm{\&sigma;}}_n\right)=I_{\mathrm{re}}\&CircleTimes;\frac{{\&PartialD;}^2{g}_{{\mathrm{\&sigma;}}_n}}{\&PartialD;x^2}$

$I_{\mathrm{xx}}^{\mathrm{re}}\left({\mathrm{\&sigma;}}_n\right)=I_{\mathrm{re}}\&CircleTimes;\frac{{\&PartialD;}^2{g}_{{\mathrm{\&sigma;}}_n}}{\&PartialD;x^2}$

$I_{\mathrm{xx}}^{\mathrm{re}}\left({\mathrm{\&sigma;}}_n\right)=I_{\mathrm{re}}\&CircleTimes;\frac{{\&PartialD;}^2{g}_{{\mathrm{\&sigma;}}_n}}{\&PartialD;x^2}$

式中代表卷积操作。

C3、利用图像I_re的二阶导数计算Hessian矩阵的特征值；

$H=\left[\begin{array}{cc}{I}_{\mathrm{xx}}^{\mathrm{re}}& I_{\mathrm{xy}}^{\mathrm{re}}\\ I_{\mathrm{xy}}^{\mathrm{re}}& I_{\mathrm{yy}}^{\mathrm{re}}\end{array}\right]$

C4、若Hessian矩阵的两个特征值|λ₁|＞|λ₂|，利用血管测度函数或血管特征拟合函数计算血管特征图

C5、得到血管特征图

C6、输出血管结构增强结果图

上面是本发明提出的方法的一种实施方式，但在某些步骤上，可以进行适当改变，以适应具体情况的需求。例如，在第一阶段的步骤A4利用三次样条进行插值，可以根据需要，进行适当的调整。例如，在第二阶段利用本征模态函数重构图像时，求和的层数可以适当调整。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种X射线冠脉造影图像血管增强方法，其特征在于：包括以下步骤：

A、对X射线冠脉造影图像I进行二维经验模态分解

A1、利用输入的X射线冠脉造影图像I作为分解初始值，设定当前层为i＝1，第i层本征模态函数的估计为r_i＝I；

A2、令第i层本征模态函数的第j次估计量为h_j＝r_i，j＝1；

A3、求取h_j的局部极大值和极小值；

A4、对所获得的极大和极小值分别用三次样条插值得到上包络env_max和下包络env_min；

A5、计算上包络env_max和下包络env_min的均值，

A6、从第i层本征模态函数的第j次估计量h_j中减去m，h_j+1＝h_j-m；

A7、令第i层本征模态函数的估计次数j＝j+1；

A8、判断：如果||h_j‖＜ε或j＞＝4不成立，则转到步骤A3；如果||h_j‖＜ε或j＞＝4成立，则获得X射线冠脉造影图像I的第i层本征模态函数F_i＝h_j；

A9、得到第i+1层本征模态函数的初始估计

A10、令本征模态函数的层序数i＝i+1；

A11、如果r_i为单调函数，则R＝r_i为残余分量，分解过程结束；否则转到步骤A2；

B、分析各层经验模态分解的本征模态函数进行去噪，并利用能够代表冠脉造影图像的特征层来加权重构去除背景噪声的图像

B1、分析X射线冠脉造影图像I的各层本征模态函数F_i，计算强度直方图得到强度分布点数最大的强度值hist_i；

B2、计算各层本征模态函数F_i的强度方差σ_i；

B3、将方差σ_i或强度值hist_i的最大者作为第i层本征模态函数的阈值Th_i；

B4、对每一层本征模态函数利用阈值Th_i将小于阈值的点置0；

B5、设定进行重构的本征模态函数层数为M_ite＝3，迭代初始值ite＝1；

B6、设定过程变量temp为0；

B7、将去噪后的第ite层本征模态函数F_i取绝对值abs(F_ite)；

B8、过程变量

其中abs(·)为取绝对值操作；

B9、令ite＝ite+1；

B10、如果ite＜(M_ite+1)，则转到B8；否则重构过程结束，输出重构冠脉造影图像I_re＝temp；

C、对重构的图像进行Hessian矩阵分解，利用血管测度函数或血管特征拟合函数达到增强血管结构的目的

C1、构造多尺度高斯滤波器集

σ_n为高斯滤波器的方差，σ_n∈[σ₁，...，σ_N]；

C2、对多尺度高斯滤波器集

求二阶导数，用的二阶导数与重构后的冠脉造影图像I_re进行卷积得到

和

C3、利用图像I_re的二阶导数计算Hessian矩阵的特征值；

C4、若Hessian矩阵的两个特征值|λ₁＞|λ₂|，利用血管测度函数或血管特征拟合函数计算血管特征图

C5、得到血管特征图

C6、输出血管结构增强结果图

Images