CN101799918B - 基于脊波变换的医学数字减影图像融合方法 - Google Patents

Abstract

本发明公开了一种基于脊波变换的医学数字减影图像融合方法，包括下列步骤：(1)对2个图像分别进行脊波变换，获得脊波变换矩阵；(2)进行融合处理；(3)对融合后的脊波变换矩阵进行脊波逆变换，获得的重构图像即为融合后的图像；其特征在于：步骤(1)中，设定初始判断阈值和步长，据此进行脊波变换，利用反变换重构融合图像，计算融合图像的信息熵，用动态模糊方法根据步长改变判断阈值重复上述操作，用对应最大信息熵的判断阈值作为最终判断阈值，获得步骤(1)所述的脊波变换矩阵。本发明能有效提高了融合图像的信息熵，降低了均方根误差率，且性能优于其他传统融合方法。算法的时间复杂度低，得到的结果较为精确，极大地丰富了医学图像的纹理细节。

Description

基于脊波变换的医学数字减影图像融合方法

技术领域

本发明属于医学图像处理技术中的图像融合领域，涉及一种医学数字减影(DSA)图像的自适应融合方法，具体涉及动态模糊技术与Ridgelet变换理论。

背景技术

当今世界，心脑血管疾病正严重威胁着人类的健康和生命安全。但血管类软组织的诊断却相当复杂，极大地依赖于计算机与X线血管造影技术等高科技手段。

数字减影血管造影(DSA)是20世纪80年代兴起的一种医疗图像学新技术，是计算机与常规X线血管造影相结合的一种新的检查方法，它集数字电子学、计算机技术与血管造影于一体的新技术，它具有实时成像、分辨率高、减影效果及绘制血管路径图等特点，可清晰地显示血管形态及分支，是研究血管影像特征的有效方法。DSA图像是通过将造影剂注入到心脑血管等感兴趣部位，将相同部位造影剂注入前后的两幅X射线成像图像相减，去除非感兴趣器官对X射线吸收形成的背景图像，以得到的清晰的血管造影图像。

由于脑部血管形态结构复杂，粗细不均，个体差异又客观存在，导致造影剂在血管里的流动速度不一致，从而产生了差异极大的DSA序列图像。单幅序列图像只能代表局部血管影像状况。医生观察整个序列图像费时费力且容易产生漏诊。因此将每一幅DSA序列图像融合到一幅完整的图像中，使之包含血管整体结构，又不乏细节表述，对医疗辅助诊断具有重要意义。

传统的图像融合技术可分为像素级的图像融合、特征级的图像融合和决策级的图像融合。其中，像素级和特征级的图像融合的研究和使用最广泛。像素是图像的基本元素，像素间灰度值的差异显现出图像中所包含的结构信息，以像素为基础的方法由于直接对图像进行逐点处理，所以，用到的数学原理易于理解，算法实现也比较简单，不过实现的效果相对较差。而以图像特征为基础的方法，由于要对图像进行特征提取、目标分割等处理，算法复杂，但是实现的效果基本能满足诊断的要求。

由于心脑血管瘤的特殊性和复杂性，这些传统的融合技术在医学图像的融合过程中都不能发挥很好的作用。首先，由于心脑脑血管图像大多是灰度图像、比较模糊，将所有序列图像按同固定阈值进行直接融合，其对比度分辨率相对就低。其次，由于人体血管组织比较复杂，结构很不规则，融合时滤波器的选择至关重要。

因而，需要一种新的方法对DSA图像进行融合，以获得较好的融合图像。利用获得的融合图像可以对图像中的血管实现检测，以构建DSA图像计算机辅助诊断系统，对心脑血管的诊断起到辅助作用。

近年来，美国斯坦福大学E.J.Candes和D.L.Donoho建立了一种非常适合于表示方向差异性的多尺度方法—Ridgelet变换。由于Ridgelet变换本质是通过对小波基函数添加一个表征方向参数得到的，所以它不但和小波一样具有局部时频分析能力，而且还具有很强的方向选择和识别能力，可以非常有效的表示信号中具有方向性的奇异性特征，如图像的线性轮廓等，这是传统小波变换方法无法做到的。大量实验表明，Ridgelet变换在直线特征的表示和提取中非常有效.[参见文献：E.J.Candes，Ridgelets.Theory and Applications[D]，Departmentof Statistics，StanfordUniversity，1998.]

由于Ridgelet变换与小波变换类似，但是Ridgelet函数引入了表示直线的参数，而小波函数采用了表示点的参数，因此小波变换可以刻画点(零维)的奇异性，但是无法刻画图像中线、面(一维或更高维)的奇异性，这一性质直接影响小波变换在表示图像边缘等几何机构方面的有效性；而Ridgelet函数的横截面是一条类似小波的曲线，Ridgelet沿着脊线是一条直线。正是这样的几何结构使得Ridgelet变换可以有效的处理图像中直线状和超平面状的奇异性。

李国新在“基于Ridgelet变换的图像融合”(通信技术2009年第6期P.144-146)一文中公开了一种基于Ridgelet变换的图像融合方法，用于照片类图像的融合。Ridgelet变换所产生的高频及低频系数是进行图像滤波处理的关键，但高频和低频的划分标准是一个模糊概念，传统的阈值处理已不能很好的处理包含复杂纹理信息的图像。因此，当将上述方法用于医学数字减影(DSA)图像融合时，难以获得较好的、均衡的融合效果。

发明内容

本发明目的是提供一种基于Ridgelet变换的DSA图像自适应融合方法，用以解决现有DSA图像无法在保持血管纹理细节的前提下实现自适应融合的缺点。

为达到上述目的，本发明首先针对DSA图像的灰度，纹理等特性进行了研究，发现DSA图像大多是灰度图像、比较模糊，将所有序列图像在不分析其所携带信息特点的情况下，按同一阈值直接融合，其对比度分辨率相对就低。其次，由于人体血管组织比较复杂，结构很不规则，简单的叠加融合会带来更多的噪声，不利于医生观察主要部位。

基于上述脑血管的形态特点，本发明考虑采用将Ridgelet变换与动态模糊理论相结合的方法进行图像融合。

本发明采用的技术方案是：一种基于脊波变换的医学数字减影图像融合方法，包括下列步骤：

(1)对待处理的序列医学数字减影图像中的2个图像分别进行脊波变换，获得各图像的脊波变换矩阵；

(2)进行融合处理，获得融合后的脊波变换矩阵；

(3)对融合后的脊波变换矩阵进行脊波逆变换，获得的重构图像即为融合后的图像；

其中，所述步骤(1)中，首先设定高低频系数的初始判断阈值和步长，根据判断阈值确定高低频系数进行脊波变换，并利用反变换重构融合图像，计算融合图像的信息熵，用动态模糊方法根据步长改变判断阈值重复上述操作，用对应最大信息熵的判断阈值作为最终判断阈值，据此获得步骤(1)所述的脊波变换矩阵；

所述信息熵计算公式为：

$E=\underset{k=0}{\overset{L-1}{\mathrm{\Σ}}}{P}_k\log \left(P_k\right)$

式中，E指信息熵，P_k是指灰度值等于K的概率，L指图像包含的灰度值量。

优选的技术方案，所述初始判断阈值为0.5，所述步长为0.1。

一种基于脊波变换的医学数字减影图像融合方法，对不少于3幅图像的序列进行图像融合，每次取序列中的2幅图像，采用前述方法进行融合，用融合后的图像替代该2幅图像，构成新的图像序列，重复上述过程，直至完成所有图像的融合。

图像融合的目标是将多幅图像合成一幅单一图像，该融合DSA图像所包含的细节信息比任何一单一源DSA图象的信息都要丰富。本方法针对多幅图像进行融合时，可以将第1张与第2张进行融合，将得到的融合图像再与第3张融合，不断迭代直至最后一张DSA图像。此过程中保持每次融合2张图像。

Ridgelet由于同时具有时、频局域性，适于表示瞬变信号，能够有效地处理高维直线或超平面奇异性问题，因而在图像处理中得到了广泛的应用。

Ridgelet处理后得到的高低频系数分别代表了图像中不同的细节信息，然而处理过程中高低频系数的划分又是一个动态模糊问题，将动态模糊逻辑理论应用到图像融合领域非常必要。因此本发明构建了一个动态模糊成员函数来对高低频系数进行优化配置，从而提升融合图像的各项质量指标。

由于上述技术方案的运用，本发明与现有技术相比具有下列优点：

1.本发明给出了一种全新的DSA图像自适应、智能融合方法，提高了融合图像的信息熵，降低了均方根误差率，具有很好的鲁棒性。

2.本发明使用图像的信息熵作为图像所携带信息量的判断标准，不但符合图像领域的科学原理，而且容易实现，为产品化奠定基础。

3.使用动态模糊技术，在处理图像等非结构化数据的时候，比传统的(静态)模糊技术具有更强的计算优势，它能判断数据变化的趋势，减少不必要的计算量，提高了算法效率。

4.实验证明，本发明的基于动态模糊Ridgelet变换的DSA图像融合方法是可行的，通过与多种其他融合方法的比较，验证了该方法性能优于平均权重法、拉普拉斯金字塔法、离散小波法等传统融合方法。提高了融合图像的信息熵，降低了均方根误差率，极大得丰富了医学图像的纹理细节，为医疗辅助诊断提供了一种新的可靠的图像融合方法。

5.将本发明应用于心脑血管的辅助诊断系统中，算法的时间复杂度低，得到的结果较为精确，可以很好地辅助医生诊断心脑血管疾病。

附图说明

图1是实施例一中基于脊波变换的DSA图像自适应融合方法总体流程图；

图2是实施例一中针对具体的2幅DSA序列图像进行自适应融合的方法示意图。

具体实施方式

下面结合附图及实施例对本发明作进一步描述：

实施例一：图1为基于Ridgelet变换的DSA图像自适应融合方法的总体流程图，数据文件(图片文件)是符合DICOM标准的脑部DSA血管图片。

步骤一：读入DSA序列图片，并且对其进行Ridgelet变换。具体步骤如下(参见图2)：

步骤1：读入DSA序列图片中的2张图片，图片大小统一由系统设置为1024像素(长)×1024像素(宽)，标记为图片1和图片2。这样即可避免由于图像过大所带来的图像分给问题。

步骤2：对图片1、图片2分别进行数字Ridgelet变换。

即通过Radon变换，一幅1024×1024的图像的像素点变为1024×2048的阵列，再对1024×2048阵列进行一维小波变换就得到了2048×2048阵列Ridgelet变换的结果，矩阵中每个单元的值都是0-1之间的一个实数。

步骤3：Ridgelet变换后分别得到图片1和图片2所对应的高频系数和低频系数(矩阵)。这里的高频和低频是一个模糊的概念，并且随着阈值的选择动态变化着。一般通过设定一个固定阈值来划分高频系数和低频系数。由于高低频系数都分布于[0，1]这个区间内，因此可以通过建立动态模糊逻辑(DFL)优化模型来优化分布高频系数和低频系数。

具体方法如下：

步骤二：建立动态模糊逻辑(DFL)优化模型实现图像的自适应融合：

首先要定义介绍动态模糊逻辑的基本概念。首先要明确动态模糊数据的几个公约.

公约1.把既包括“动态的”又包括“模糊的”的数据特性称作动态模糊性。

公约2.把具有动态模糊性的数据称作动态模糊数据。

公约3.把动态模糊数据构成的集合称作动态模糊数据集U.

建立动态模糊逻辑(DFL)优化模型具体步骤：

步骤1：在论域U上定义映射：

$\stackrel{\→}{A}:\stackrel{\→}{a}\→\left[\mathrm{0,1}\right],\stackrel{\→}{a}\→\stackrel{\→}{A}\left(\stackrel{\→}{a}\right)$

或者

$\stackrel{\←}{A}:\stackrel{\←}{a}\→\left[\mathrm{0,1}\right],\stackrel{\←}{a}\→\stackrel{\←}{A}\left(\stackrel{\←}{a}\right).$

记为 $(\stackrel{\→}{A},\stackrel{\←}{A})=\stackrel{\→}{A}\mathrm{or}\stackrel{\←}{A},$ 则称

为U上的动态模糊数据集(dynamic fuzzy data sets)，简称DFDS，称为隶属函数对

的隶属度。

任何一个动态模糊数据a∈[0，1]，都可动态模糊化为：

且 $\max (\stackrel{\←}{a},\stackrel{\→}{a})\stackrel{\mathrm{\Δ}}{=}\stackrel{\→}{a},$ $\min (\stackrel{\←}{a},\stackrel{\→}{a})\stackrel{\mathrm{\Δ}}{=}\stackrel{\←}{a}.$

步骤2：对DSA序列图片进行基于动态模糊的Ridgelet变换。

具体步骤如下：

步骤2-1：本发明定义 $(\stackrel{\←}{U},\stackrel{\→}{U})=[\stackrel{\←}{0},\stackrel{\→}{0}],[\stackrel{\←}{1},\stackrel{\→}{1}]=\left[\mathrm{0,1}\right]\×[\←,\→],$ 动态模糊数据集为和

分别表明频率“增高”和频率“降低”.其中“←”表明增高方向的动态变化趋势；“←”表明降低方向的动态变化趋势。

步骤2-2：假设动态模糊数据成员α∈[0，1]，则α属于动态模糊数据集，定义如下：

且 $\max (\stackrel{\←}{a},\stackrel{\→}{a})\stackrel{\mathrm{\Δ}}{=}\stackrel{\→}{a},$ $\min (\stackrel{\←}{a},\stackrel{\→}{a})\stackrel{\mathrm{\Δ}}{=}\stackrel{\←}{a}$

步骤2-3：计算两个动态模糊数据子集(分别对应两幅被融合图片高低频信息)可以理解为两个成员函数之间的计算。图片1对应使用成员函数A，图片2对应使用成员函数B。本发明构建的新的成员函数如下：

$\stackrel{\←}{A}\left(\stackrel{\←}{u}\right)=\left\{\begin{array}{cc}0& \mathrm{if}0\≤\stackrel{\←}{u}\≤\stackrel{\←}{0.5}\\ {(1+{\left(\frac{\stackrel{\←}{u}-0.5}{0.05}\right)}^{-2})}^{-1}& \mathrm{if}\stackrel{\←}{0.5}\≤\stackrel{\←}{u}\≤\stackrel{\←}{1}\end{array}\right.$

$\stackrel{\→}{A}\left(\stackrel{\→}{u}\right)=\left\{\begin{array}{cc}0& \mathrm{if}0\≤\stackrel{\→}{u}\≤\stackrel{\→}{0.5}\\ {(1+{\left(\frac{\stackrel{\→}{u}-0.5}{0.05}\right)}^{-2})}^{-1}& \mathrm{if}\stackrel{\→}{0.5}\≤\stackrel{\→}{u}\≤\stackrel{\→}{1}\end{array}\right.$

$\stackrel{\←}{B}\left(\stackrel{\←}{u}\right)=\left\{\begin{array}{cc}0& \mathrm{if}0\≤\stackrel{\←}{u}\≤\stackrel{\←}{0.5}\\ {(1+{\left(\frac{\stackrel{\←}{u}-0.5}{0.05}\right)}^{-2})}^{-1}& \mathrm{if}\stackrel{\←}{0.5}\≤\stackrel{\←}{u}\≤\stackrel{\←}{1}\end{array}\right.$

$\stackrel{\→}{B}\left(\stackrel{\→}{u}\right)=\left\{\begin{array}{cc}0& \mathrm{if}0\≤\stackrel{\→}{u}\≤\stackrel{\→}{0.5}\\ {(1+{\left(\frac{\stackrel{\→}{u}-0.5}{0.05}\right)}^{-2})}^{-1}& \mathrm{if}\stackrel{\→}{0.5}\≤\stackrel{\→}{u}\≤\stackrel{\→}{1}\end{array}\right.---\left(6\right)$

本实施例设置对高低频系数的判断阈值的步长为0.1，即阈值最多可以取0.1，0.2，0.3，……0.9这9个值。动态模糊技术可以计算阈值的变化趋势，减少计算量，提升图像融合效率。设置阈值初值取0.5，公式(6)列出了阈值取0.5时成员函数的构建方式。

步骤三：用Ridgelet反变换重构DSA图像，计算图像信息熵，不断重复步骤二，直到判断得到最大熵。具体如下：

步骤1：将图片1和图片2通过Ridgelet变换和动态阈值分析后生成的高频系数一一比较取最大值保留；同样将图片1和图片2通过Ridgelet变换和动态阈值分析后生成的低频系数一一比较取最小值保留。

步骤2：将优化所得的高、低频系数通过下式进行ridgelet反变换，重建图像，即得初步的融合图像：

$f\left(x\right)={\∫}_0^{2\mathrm{\π}}{\∫}_{-\∞}^{\∞}{\∫}_0^{\∞}{R}_f(a,b,\mathrm{\θ}){\mathrm{\ψ}}_{a,b,\mathrm{\θ}}\left(x\right)f\left(x\right)\frac{\mathrm{da}}{a^3}\mathrm{db}\frac{\mathrm{d\θ}}{4\mathrm{\π}}---\left(7\right)$

步骤3：通过信息熵计算公式计算图像信息熵，不断重复以上，直到判断得到最大熵。信息熵计算公式如下：

$E=\underset{k=0}{\overset{L-1}{\mathrm{\Σ}}}{P}_k\log \left(P_k\right)---\left(8\right)$

其中，E指信息熵，P_k是指灰度值等于K的概率，L指图像包含的灰度值量。

例如：当阈值＝0.5时，计算Ridgelet反变换后生成的融合图像的信息熵E1。

根据成员函数，求当阈值＝0.4时，计算Ridgelet反变换后生成的融合图像的信息熵E2。

若E1＜E2，下一个阈值取0.3，得信息熵与E2比较。不断递进计算得最大信息熵。

若E1＞E2，根据成员函数，求阈值变化趋势的反方向，即求当阈值＝0.6时，计算Ridgelet反变换后生成的融合图像的信息熵E3。

若E1＜E3，下一个阈值取0.7，得信息熵与E3比较。不断递进计算得最大信息熵。

若E1＞E3，则取信息熵E1为最大信息熵。

步骤四：输出包含最大信息熵的DSA图像，即得到最终包含丰富细节的融合图像。

表1表示取不同阈值时融合图像所对应的信息熵。经计算，阈值取0.6时所对应的图片(c)的信息熵最大，所以输出图片(c)为最优融合图片。

表1各阈值所对应信息熵

注：“←”表明增高方向的动态变化趋势；

“→”表明降低方向的动态变化趋势。

说明：本实施例所用DSA实例图片都遵守DICOM标准。典型的DSA医学序列图片一般由20-22张构成。

通过构建动态模糊成员函数，我们可以充分利用Ridgelet高低频系数的动态模糊特性。对系数分布进行优化配置，提升融合性能。Ridgelet是一种能够表示方向选择能力的多尺度变换方法，促进着计算和应用谐波分析的领域的发展。与传统小波相比Ridgelet有着很好的方向选择它能有效地处理二维空间中具有直线奇性的信号。动态模糊逻辑理论在处理此类非结构化数据时具有极大的优越性，本发明将二者有机结合成基于Ridgelet变换的DSA图像自适应融合方法，其融合性能有了很大的提高，且具有很好的鲁棒性。

Claims (3)

1.一种基于脊波变换的医学数字减影图像融合方法，包括下列步骤：

(1)对待处理的序列医学数字减影图像中的2个图像分别进行脊波变换，获得各图像的脊波变换矩阵；

(2)进行融合处理，获得融合后的脊波变换矩阵；

(3)对融合后的脊波变换矩阵进行脊波逆变换，获得的重构图像即为融合后的图像；

其特征在于：

所述步骤(1)中，首先设定高低频系数的初始判断阈值和步长，根据判断阈值确定高低频系数进行脊波变换，并利用反变换重构融合图像，计算融合图像的信息熵，用动态模糊方法根据步长改变判断阈值重复上述操作，用对应最大信息熵的判断阈值作为最终判断阈值，据此获得步骤(1)所述的脊波变换矩阵；

所述信息熵计算公式为：

$E=\underset{k=0}{\overset{L-1}{\mathrm{\Σ}}}{P}_k\log \left(P_k\right)$

式中，E指信息熵，P_k是指灰度值等于K的概率，L指图像包含的灰度值量。

2.根据权利要求1所述的基于脊波变换的医学数字减影图像融合方法，其特征在于：所述初始判断阈值为0.5，所述步长为0.1。

3.一种基于脊波变换的医学数字减影图像融合方法，其特征在于：对不少于3幅图像的序列进行图像融合，每次取序列中的2幅图像，采用权利要求1的方法进行融合，用融合后的图像替代该2幅图像，构成新的图像序列，重复上述过程，直至完成所有图像的融合。

Images