基于混沌遗传算法和contourlet变换的医学图像融合

【摘要】  　　目的： 改善传统医学图象融合方法对细节信息的丢失. 方法： 利用contourlet的多尺度、方向性和各向异性等优点，提出了一种基于局域特性匹配度融合的改进算法. 首先对原图进行contourlet分解，对分解各子带构造局域特征的匹配度，设定匹配度的阈值，对匹配度在不同阈值范围内的系数进行不同的加权融合；最后对融合系数进行contourlet逆变换得到融合图像. 在阈值的选取上，引入混沌遗传算法求解阈值的全局最优解. 结果： 运用传统小波变换和本文提出对方法对两组医学图象进行融合处理，客观评价参数表明本文提出的算法效果更优. 结论： 该算法能够在保有原图信息的同时，有效的增强细节信息.

【关键词】  contourlet 医学图像 图像融合 遗传算法

　　【Abstract】AIM: To improve the traditional image fusion algorithm in order to avoid the loss of  the detailed information in the processe of image fusion. METHODS: Utilizing the contourlet?s advantages of multiscale,directionality and anisotropy,an advanced image fusion algorithm based on characteristic matching of region statistics in contourlet domain was proposed in this paper. First，the source images were transformed into contourlet domain. Then,we constructed the characteristic matching of each subband and set a threshold for it. The coefficients whose characteristic matching was in different ranges applied different fusion rules. Finally, the fused image was obtained by taking the inverse contourlet transformation. Chaos genetic algorithm was adopted to optimize fusion image by searching optimal solution of the threshold of characteristic matching. RESULTS: The comparison between the result of the traditional method and that of the proposed method showed that the later one was better. CONCLUSION: The proposed method can effectively preserve the information of the source images and enhance the details of the fused image.

　　【Keywords】 contourlet; medical image; image fusion; genetic algorithms

　　0   引言
   
　　传统的医学图像融合算法多采用基于小波变换，小波系数加权平均，或者阈值法融合，但二维可分小波仅能捕捉有限的方向信息，不能稀疏的表示含线或者面奇异的二维图像［1-3］. Do等［4］提出了一种稀疏的图像二维表示方法： contourlet变换. 该变换的“各向异性(anisotropy)”使得它能比小波变换更好的表现边缘特征，更加适合用来进行图像的处理［5］. 并且小波变换下的传统方法有很大的局限，比如加权平均法中加权因子和基于特征匹配度融合方法中的阈值大小都不易确定［6］. 本文提出了一种新的基于匹配度融合算法：在contourlet变换多方向的优势下，阈值的确定采用混沌遗传算法. 混沌遗传算法是近年提出的一种优化算法，能够有效的防止遗传算法的“早熟”问题，使得遗传算法能够以更快的速度收敛，得到阈值的全局最优解.

　　1    材料和方法

　　1.1   材料   选用两组融合图象，一组为实验标准图象，另一组为腰椎冠状位的MRI-T1和MRI-T2成像，编码格式为DICOM，来源于西安红十字会磁共振室.

　　1.2   方法

　　1.2.1   contourlet变换原理   contourlet变换也称塔型方向滤波器组(pyramidal directional filter bank,PDFB). contourlet 变换将多尺度分析和方向分析分拆进行：用拉普拉斯塔形滤波器结构(laplacian pyramid,LP)对图像进行多分辨率分解. 首先产生原始信号的一个低通采样逼近及原始图像与低通预测图像之间的一个差值图像，对得到的低通图像继续分解得到下一层的低通逼近和差值图像. 如此逐步滤波得到图像的多分辨率分解. 而后应用二维方向滤波器组(directional filter bank,DFB)对分解得到的每一级高频分量在任意尺度上再分解得到 个方向子带，将分布在同方向上的奇异点合成为一个系数. 图1A给出了离散Contourlet 变换的滤波器组结构. 原始图像经PDFB结构分解得到一个低通图像和分布于多尺度多方向上的高频分量. 图1B为Contourlet 频域分解图实际应用中方向数随着尺度增大而增多.   

　　A： 滤波器组结构图； B： 频率分解图.

　　图1   contourlet变换（略）

　　1.2.2   图像融合框架   本文的融合方法的改进之一就是将contourlet 引入图像融合，可以利用其优良特性更好地提取原始图像中的几何特征，为融合图像提供更多的信息.

　　1.2.2.1   contourlet分解   将精确配准的原图像进行Contourlet变换,得到相应的contourlet系数集合. 在融合中一般采用3层或4层分解. 设原图像为A和B，使用J级LP分解，每个尺度上的方向数分别为i(j)，则原图像的Contourlet分解过程为:
   
　　A?(aJA,bJA,bJ-1A,……,b2A,b1A)(1)
   
　　B?(aJB,bJB,bJ-1B,……,b2B,b1B)(2)
   
　　bJid2Ji,d2JI,……,d1i,   i=A,B(3)
   
　　其中aJi是低频子带; bJi为尺度j上的方向子带集合; dki为尺度j上第k个方向的高频子带,k=1,2,……2J.

　　1.2.2.2   图像融合   对于分解后的低频子带和所有高频子带，使用基于区域特性的融合规则进行判别和融合处理,得到各尺度上融合后的contourlet系数.

　　1.2.2.3   Contourlet重构   重构是分解的逆过程,对融合后的contourlet系数进行contourlet逆变换，得到重构的融合图像，该图像包含原有多幅图像中的信息. 设 为融合后的结果，此过程可以表示为
   
　　(aJF,bJF,bJ-1F,……,b2F,b1F)?F(4)

　　1.2.3   融合规则   不同的融合算法都有其相应的应用价值. 目前还没有一种融合算法适合所有的图像融合. contourlet对线条有高度敏感性，这对细节质量的提高很有好处. 由于contourlet变换系数与其邻域有极强的依赖性，在子带系数处理中，我们选用了基于区域特征的方法，所以对若干个方向的高频分量，本文采用如下方法进行融合：
   
　　（1） 选择一个矩形窗口L×K，如，contourlet分解各子带系数的均值. 子带中以点(n,m)为中心的区域均值为：
   
　　M-i(n,m)=1/L×K∑/n′∈L∑/m′∈Kfi(n+n′,m+m′)　(5)
   
　　i代表原图A或B，fi表示图A或图B的相应子带系数.
   
　　（2） 计算以点(n,m)为中心的区域方差δi2(n,m)
    
　　δ2i(n,m)=1/L×K∑/n′∈L∑/m′∈K(fi(n+n′,m+m′)-M-I(n,m))2(6)
   
　　（3） 计算以点(n,m)为中心图像A与图像B子带对应区域的协方差：
   
　　β2(n,m)=1/L×K∑/n′∈L∑/m′∈K｜(fA(n+n′,m+m′)-M-A(n,m))×fB((n+n′,m+m′)-M-B(n,m))｜　(7)
   
　　（4） 构造匹配度ρ(n,m)：
   
　　ρ(n,m)=2×β2(n,m)/σ2A(n,m)+σ2B(n,m)(8)

   
　　（5） 选定阈值α，其值大小一般在0.5～0.8之间. 最优值取决于待融合的图像. 本文用混沌遗传算法来确定阈值. 具体方法见后. 如果ρ(n,m)≥α，则融合结果中对应点值dKF(n,m)为（d2j为高频子带）：
  
　　dKF(n,m)=dKA(n,m)   if   σ2A(n,m)≥σ2B(n,m)

　　dKF(n,m)=dKB(n,m)   if   σ2A(n,m)≥σ2B(n,m)
 
　　(9.1)
           　　 
　　如果ρ(n,m)<α，则融合结果中对应点的值d2jF(n,m)为:
   
　　dKF(n,m)=dKA(n,m)×wmax(n,m)+dkB(n,m)×

　　wmin(n,m)   if   σ2A(n,m)≥σ2B(n,m)

　　dKF(n,m)=dKB(n,m)×wmax(n,m)+dkA(n,m)×

　　wmin(n,m)   if   σ2A(n,m)≥σ2B(n,m)(9.2)
   
　　其中     wmin(n,m)=1/2-1/2(1-ρ/1-α)

　　wmin(n,m)=1-wmin(n,m)(10)
   
　　对于低频分量，采用对应点取平均值的融合方法.

　　1.2.4   基于混沌遗传算法的阈值选取   下边讨论如何选取阈值. 阈值的选取是本文提出的算法的关键，也是本文对融合算法的改进之一. 对于不同的原图，能够取得最优融合效果的阈值也必然不同.

　　1.2.4.1   混沌理论   多目标函数的全局最优求解问题，由于存在目标之间的无法比较和冲突的问题，情况会复杂的多. 遗传算法在搜索过程中，为了保证算法能够高效快速的收敛到全局最优，必须保证种群的多样性. 在机上可以用混沌迭代方程来模拟混沌现象-Logistic迭代方程：
   
　　xl+1=uxl(1-xl)    l=1,2…… (11)

　　式中xl为混沌变量，u为控制参量,0

　　1.2.4.2   适应度函数   混沌遗传算法的核心所在是设立有效的适应度函数.  本文采用了交叉熵和平均梯度进行评价,这两个客观评价参数都是我们的目标函数.
   
　　设一幅图像的灰度分布为P={p(0),p(1),…,p(L-1)}，L 表示图像的总的灰度级数；图像大小为M×N.
   
　　①交叉熵：直接反映了两幅图像对应像素的差异,是对两幅图像X,Y所含信息的相对衡量,公式为：
   
　　CEX,Y=∑/L-1g=0 pXi(g)· log2pXi(g)/pYi(g)(12)
   
　　这里我们定义融合图像 与原图,的交叉熵为
   
　　CEN=(CE2A,F+CE2B,F)/2(13)
   
　　②平均梯度：反映图像质量的改进,及图像中微小细节反差和纹理变换特征. 其公式为
        
　　G=1/MN∑M  n=1∑N  m=1△xf(n,m)2+△yf(n,m)2 (14)
   
　　其中△xf(n,m)和△yf(n,m)分别为图像沿着x方向和y方向的差分.
   
　　由以上两个评价参数我们用妥协方法定义适应度函数：
   
　　Fitness(α)=［CEN2(Fα)+(Gmax-G(Fα))2］1/2 (15)
   
　　Fα表示阈值为α时，得到的融合图像. 注意这里适应度函数的取值越小，适应度越高.

　　1.2.4.3   实现步骤   ①编码：本文采用了实数编码. 实数编码对于函数优化问题最为有效，这一点已经得到了广泛的验证［8］. 实数编码的基因型空间的拓朴结构与其表现型空间中的拓朴结构一致，不存在Hamming悬崖的问题. ②初始群体的生成：由阈值的定义，α∈(0,1)；由先验知识，阈值在（0.5，0.8）时融合效果较佳，初始群体定义在（0.5，0.8），均匀、随机产生Ｎ＝16个初始个体；在约束条件中体现α∈(0,1). ③)选择、交叉、变异. ④混沌化：随机均匀的选取k个个体αi，将其映射到Logistic方程的定义域［0,1］，由于阈值本身的定义域就在［0,1］,所以映射与逆映射可以免除. 方程迭代500次，用得到的混沌变量替代之前的种群个体. 转到第三步继续执行，直到满足结束条件. 　
   
　　αl=4αl(1-αl)　(16)
   
　　k的选取原则是随着遗传循环次数的增加而减小的. 这样是为了让群体刚开始循环时充分利用混沌的遍历性来搜索，当群体向最优解靠近时，只需要在其周围小范围内搜索.

　　2   结果
   
　　本文把这种方法运用到了两组医学图像中. 图2A和B分别为MRI和CT的颅腔横断面，尺寸为；E和F分别为腰椎冠状位的MRI-T1和MRI-T2成像，尺寸为. 本文分别采用了小波变换法(采用Daubechies 双正交样条小波bior2)和本文提出的算法进行融合，分解的细节空间是相同的，都分解了四层. 我们知道阈值在（0.5，0.8）时融合效果较佳，所以小波变换法中的阈值这里选取的是0.65（图2）.  
   
　　这里选取了4个图像融合的性能评价参数，除了作为目标函数的交叉熵和平均梯度以外，本文还选取了熵和互信息作为客观评价准则. 熵(entropy)：图像的熵是衡量图像信息丰富程度的一个重要指标. 一幅图像的熵的大小反映了图像携带信息的多少，其定义为:
   
　　H=∑L-1  i=0pi·log2 pi(17)
   
　　互信息量MI（mutual information）也叫相关熵. 融合图像与原始图像的互信息（相关熵）越大越好，公式为：
    
　　MIX,Y=∑L  i=1∑L  j=0hX,Y(p,q)glog2hX,Y(i,,j)/hX(i)+hY(j)(18)
   
　　其中为hX,Y(i,,j)图像X,Y的归一化联合灰度直方图分布,hX(i)和hY(,j)分别为两幅图像的归一化边缘直方图分布. 这里我们定义融合图像F与源图A,B的互信息为:
   
　　MI=MI2A,F+MI2B,F/2 (19)

　　A：MRI; B： CT; C：小波融合结果; D：本文方法融合结果; E： MIR-T1 ; F：MRI-T2; G：小波融合结果; H：本文方法融合结果.

　　图2   融合图

　　式(19)表示融合图像F包含源图A，B的交互信息量的总和，这个值越大，表示融合图像从源图中获取的信息量越大，效果越好. 从融合结果的客观评价（表1），可以看出相对小波变换法，本文算法融合结果的细节更为清晰.

　　表1   客观评价（略）

　　3   讨论
   
　　本文提出了一种基于contourlet变换和遗传算法相结合的融合算法，并应用于医学图像的融合中. 这种算法有效的利用了contourlet变换多尺度、方向性和各向异性等优点，弥补了wavelet变换分解方向单一（只有三个方向）的不足，继而利用混沌遗传算法求取匹配度阈值的最优解，使得融合算法能够对不同输入原图像，自适应的给出最优阈值，从而得到这一组原图像的最优融合效果：融合图像既有效的保存了边缘能重要信息，又较好的体现了各原图的细节特征. 四种有效的客观评价参数证明了此方法的有效性和优越性.
   
　　但是本文提出的方法也有一定的不足，基于混沌遗传算法的阈值选取在处理一般大小的图像时，时间消耗是能够接受的，一般大约为20 min左右，就可以收敛到最优阈值，并得到最后的融合图像. 但是对于像素点很多的图像时，处理的速度会大幅下降，时间长度会达到若干小时. 在这个问题上，还有待进一步的研究和改进.

【】
  　　［1］ Zhu YM,Steven M. Cochoff,an object-oriented framework for medical image registration,fusion,and visualization［J］. Computer Meth Prog Biomed,2006,82: 258-267.

　　［2］ Nemec SF.CT-MR image data fusion for computer assisted navigated neurosurgery of temporal bone tumors［J］. Eur J Radiol,2006,10:10-16.

　　［3］ Petrovic V,Cootes T. Objectively adaptive image fusion［J］. Infor Fusion,2005,10:143-151.

　　［4］ Do MN,Vetterli M. Contourlets ［A］// Welland GV. Beyond wavelets［M］. New York: Academic Press,2001: 1-27.

　　［5］ Duncan DY,Do MN. Directional multiscale statistical modeling of images using the contourlet transform［J］. IEEE Trans Image Proc,2006,15(6):1610-1620.

　　［6］ 陶冰洁. 多原图像配准与融合方法研究［D］. 北京：院工程研究所,2003: 39-54.

　　［7］ Buccigrossi RW,Simoncelli E P,Image compression via joint statistical characterization in the wavelet domain［J］. IEEE Trans Image Proc,1999,8: 1688-1701.