CN102663710B - 一种x射线冠脉造影图像降噪方法 - Google Patents

Abstract

本发明公开了一种X射线冠脉造影图像降噪方法，包括以下步骤：对X射线冠脉造影图像I进行二维经验模态分解；分析经验模态分解得到的各层本征模态函数，进行去噪；利用能够代表冠脉造影图像血管结构的本征模态函数进行加权求和来重构去噪后的冠脉造影图像。本发明不依赖先验知识进行对冠脉造影图像自适应分解，根据图像特点自适应产生随机噪声去除阈值。本发明利用经验模态分解方法对X射线冠脉造影图像进行自适应分解，对得到的各IMF分量进行去噪，并选择能够代表血管结构的IMF高频分量进行造影图像重构即可达到去除噪声的目的，且在去噪过程中无须先验知识、能够根据造影图像特点进行自适应降噪。

Description

一种X射线冠脉造影图像降噪方法

技术领域

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

背景技术

冠状动脉造影检查(简称“冠脉造影”)是目前国际上公认的诊断冠心病的最佳手段，被称为冠心病诊断的“金标准”。冠脉造影由于成像过程要经过X射线管、光学成像镜头等许多设备，最后还要通过图像数字采集卡将视频信号转化为数字信号并存储在计算机中，这就不可避免的引入许多噪声，必然也给造影图像的处理增加了难度。特别是在X射线机透视状态下采集图像，X线的能量很低，产生图像的随机噪声较大，其图像背景噪声经常同血管造影的密度变化范围呈重叠状，图像的信噪比非常低。噪声的存在使获得的图像不清晰，尤其是掩盖和降低了造影图像中某些特征细节的可见度。可见度的损失对低对比度物体尤为重要，对以后的血管图像识别更为不利。因此需要尽量减少噪声的影响。

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

经验模态分解(Empirical Mode Decomposition,EMD)是近年来引起人们关注的一种新的时频分析方法，适合于分析非线性、非平稳信号序列，具有很高的信噪比。它依据数据自身的时间尺度特征来进行信号分解,无须预先设定任何基函数。这一点与建立在先验性的谐波基函数和小波基函数上的傅里叶分解与小波分解方法具有本质性的差别。经验模态分解将复杂信号分解为有限个(层)本征模态函数(Intrinsic Mode Function，简称IMF)，所分解出来的各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层即第i个本征模态函数的阈值Th_i，即Th_i＝max(σ_i,hist_i)，式中max(,)为取最大值操作；

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

C、利用能够代表冠脉造影图像血管结构的本征模态函数进行加权求和来重构去噪后的冠脉造影图像

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

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

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

C4、过程变量 $\mathrm{temp}=\mathrm{temp}+5e^{-\frac{{(\mathrm{ite}-1)}^2}{4}}\mathrm{abs}\left(F_{\mathrm{ite}}\right);$

C5、令ite＝ite+1；

C6、如果ite＜4，则转到C4；否则重构过程结束，输出重构冠脉造影图像I_en＝temp。

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

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

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

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

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

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

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

4、本发明的方法包括三个过程：一个图像分解过程、一个噪声估计与随机噪声消除过程以及一个加权重构过程；利用经验模态分解自适应的将X射线冠脉造影图像分解为本征模态函数，进而对本征模态函数进行噪声估计和去噪，选择三个高频分量(前三层本征模态函数)来重构去噪后的图像。

5、本发明利用经验模态分解方法对X射线冠脉造影图像进行自适应分解，对得到的各IMF分量进行去噪，并选择能够代表血管结构的IMF高频分量进行造影图像重构即可达到去除噪声的目的，且在去噪过程中无须先验知识、能够根据造影图像特点进行自适应降噪。

附图说明

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

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

图2是在本征模态函数中进行降噪的流程图。

具体实施方式

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

如图1-2所示，本方法包括三个分离的阶段：(1)利用二维经验模态分解方法对X射线冠脉造影图像进行分解，得到本征模态函数；(2)对各层本征模态函数确定阈值并去掉随机噪声；(3)利用去噪后的高频本征模态函数重构图像。具体步骤如下：

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

${\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层即第i个本征模态函数的阈值Th_i，即Th_i＝max(σ_i,hist_i)，式中max(,)为取最大值操作；

${\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.$

C、第三阶段利用去噪后的高频本征模态函数重构图像。利用去除随机噪声后的三个最高频本征模态函数来重构冠状动脉图像I_en：

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

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

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

C4、过程变量 $\mathrm{temp}=\mathrm{temp}+5e^{-\frac{{(\mathrm{ite}-1)}^2}{4}}\mathrm{abs}\left(F_{\mathrm{ite}}\right);$

C5、令ite＝ite+1；

C6、如果ite＜4，则转到C4；否则重构过程结束，输出重构冠脉造影图像I_en＝temp。输出的重构冠脉造影图像为：

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

上面是本发明提出的方法的一种实施方式，但在某些步骤上，可以进行适当改变，以适应具体情况的需求。例如，在第一阶段的步骤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层即第i个本征模态函数的阈值Th_i，即Th_i＝max(σ_i,hist_i)，式中max(,)为取最大值操作；

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

C、利用能够代表冠脉造影图像血管结构的本征模态函数进行加权求和来重构去噪后的冠脉造影图像

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

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

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

C4、过程变量 $\mathrm{temp}=\mathrm{temp}+5e^{-\frac{{(\mathrm{ite}-1)}^2}{4}}\mathrm{abs}\left(F_{\mathrm{ite}}\right);$

C5、令ite＝ite+1；

C6、如果ite＜4，则转到C4；否则重构过程结束，输出重构冠脉造影图像I_en＝temp。