基于分阶段多原子正交匹配跟踪的多模图像融合方法.pdf

本发明提供了一种基于分阶段多原子正交匹配跟踪的多模图像融合方法，该方法将待融合的源图像划分为相同大小的图像块，并将每个图像块转换成为一维列向量，从而构成源图像的矩阵表达形式，再对每个图像块进行稀疏分解，然后对不同源图像的图像块稀疏分解系数采用一定的融合规则进行融合，得到融合图像块系数矩阵，最后对融合图像块系数矩阵进行重构，生成最终的融合图像，在对图像块信号进行稀疏分解时，在每个阶段中采用多原子选择策略，再通过自组织搜索裁剪算法去除冗余原子，降低支撑集的相干性，从而保证有更有效地获取原始图像块的稀疏表达，采用本发明提出的图像融合方法，能够极大地提高了算法执行速度并改善融合图像的质量。

1.  一种基于分阶段多原子正交匹配跟踪的多模图像融合 方法，其特征在于，该方法采用包括自组织搜索裁剪亚步骤(2-5) 的稀疏分解过程，所述自组织搜索裁剪亚步骤(2-5)包括以下 子步骤：
子步骤(2-5-1)，初始化子支撑集为空集；
子步骤(2-5-2)，在候选集中搜索与当前残差相干性最大 的原子，将搜索到的上述原子并入子支撑集，并将之从候选集 中剔除，更新候选集；
子步骤(2-5-3)，在子步骤(2-5-2)更新的候选集中继续 搜索出与子步骤(2-5-2)中搜索到的原子相关性大于第一阈值 的原子；
子步骤(2-5-4)，所述子步骤(2-5-3)中搜索到的原子构 成冗余原子集合，将所述冗余原子集合从候选集中剔除；
子步骤(2-5-5)，判断：候选集是否满足返回子步骤(2-5-2) 的迭代条件，如果满足，返回子步骤(2-5-2)，如果不满足， 则迭代终止。

2.  根据权利要求1所述的多模图像融合方法，其特征在于，
子步骤(2-5-3)中，所述第一阈值取值为0.1；
子步骤(2-5-5)中，判断是否返回子步骤(2-5-2)的迭代 条件为：当候选集不为空时，返回子步骤(2-5-2)。

3.  根据权利要求1或2所述的多模图像融合方法，其特征在 于，图像块像素矩阵都采用包括自组织搜索裁剪亚步骤的分阶 段多原子正交匹配跟踪算法进行稀疏分解，所述包括自组织搜 索裁剪亚步骤的分阶段多原子正交匹配跟踪算法包括以下亚步 骤：
亚步骤(2-1)，初始化残差信号为待分解信号，初始化信 号的支撑集为空集，初始化当前阶段数为零；
亚步骤(2-2)，利用过完备字典的转置与当前残差信号做 内积计算得到内积系数向量；
亚步骤(2-3)，定义第二阈值为常数与空原子系数绝对值 的最大值之积；
亚步骤(2-4)，在图像块像素矩阵中，选择原子系数绝对 值大于第二阈值的原子，构成候选集；
亚步骤(2-5)，对候选集进行自组织搜索裁剪，生成子支 撑集；
亚步骤(2-6)，将亚步骤(2-5)迭代过程产生的子支撑集 并入到信号的支撑集当中；
亚步骤(2-7)，根据当前信号的支撑集构成的子字典，利 用最小二乘法计算当前稀疏向量；
亚步骤(2-8)，更新当前残差信号，当前残差信号为待分 解信号减去过完备字典与当前稀疏向量之内积所得的差；
亚步骤(2-9)，判断当前残差信号或当前阶段数是否满足 迭代过程返回亚步骤(2-2)的迭代条件，如果满足，返回亚步 骤(2-2)，如果不满足，迭代终止。

4.  根据权利要求1～3之一所述的多模图像融合方法，其特 征在于，
亚步骤(2-3)中，常数取值为2.38；
亚步骤(2-4)中，所述原子系数为当前残差信号与过完备 字典中原子的内积；
亚步骤(2-9)中，判断是否返回亚步骤(2-2)的迭代条 件为：当前残差信号的范数大于逼近误差，或者当前阶段数小 于最大阶段数时，返回亚步骤(2-2)，其中，所述逼近误差是 指每次迭代的误差残余，所述最大阶段数是指最大迭代次数。

5.  根据权利要求1～4之一所述的多模图像融合方法，其特 征在于，该方法包括以下步骤：
步骤1，输入两幅待融合的原始图像A和B，将两幅原始图 像分别划分成图像块，每块图像块的大小一致、任何相邻的两 块图像块都具有重叠部分，并将每块图像块的像素表示成列向 量，得到每块图像块的像素矩阵，即得到图像块的像素矩阵表 示；
步骤2，对上述每块图像块像素矩阵进行分离，得到直流分 量和残余图像块像素矩阵，对得到的残余图像块像素矩阵进行 分阶段多原子正交匹配跟踪自组织搜索裁剪稀疏分解，获得残 余图像块的稀疏表示系数矩阵；
步骤3，根据能量最大原则，对步骤2获得的残余图像块的 稀疏表示系数矩阵进行融合，得到该位置的残余图像块的融合 图像块系数矩阵；对于步骤2分离出的直流分量，采取加权平均 方法进行融合，得到融合后的直流分量；
步骤4，对步骤3融合后的系数矩阵进行重构得到残余图像 块的融合图像块像素矩阵，对所述残余图像块的融合图像块像 素矩阵与所述融合后的直流分量进行合并，得到融合后的图像 块像素矩阵；
步骤5，将步骤4获得的融合图像块像素矩阵还原，得到最 终的融合图像；
其中，待融合的两幅原始图像是针对同一部位的，不同成 像设备提供的多模医学图像；
所述分离出直流分量包括计算各个图像块的平均像素值， 然后从各个图像块的像素值中减去该平均像素值，即直流分量 是各个图像块的平均像素值；
所述能量最大原则是将两幅原始图像的对应位置图像块的 稀疏分解系数绝对值最大者作为融合后图像块对应位置的稀疏 分解系数；
所述重构是重新构成各个图像块的空间像素表示；
所述还原是将图像块的一维表示还原为二维表示。

6.  根据权利要求1～5之一所述的多模图像融合方法，其特 征在于，步骤1中，
所述两幅原始图像为CT图像与PET图像，或者CT图像与 MRI图像，或者CT图像与SPET图像，或者MR图像与SPEC图像；
所述划分为滑动处理，所述滑动处理为在长为W、宽为H 的图像中，按一定规律移动一个长为宽为的晶格窗 口，将该窗口对应的像素块转换成块像素矩阵表示中的一列， 然后窗口向右或向下移动一步，直到完成对整幅图像的处理， 优选每次平移两个像素点，即滑动步长为2。

7.  根据权利要求1～6之一所述的多模图像融合方法，其特 征在于，步骤3中，图像块I_A的稀疏表示系数矩阵用X_A表示， 图像块I_B的稀疏表示系数矩阵用X_B表示，融合后图像块对应位 置的稀疏表示系数矩阵为：
X F = X A , | X A | ≥ | X B | X B , else ]]>
其中，X_F表示融合后图像块对应位置的稀疏表示系数矩阵；
即，当|X_A|≥|X_B|时，X_F＝X_A；
否则，X_F＝X_B。

8.  根据权利要求1～7之一所述的多模图像融合方法，其特 征在于，步骤4中，所述重构的方法为：由融合后的稀疏系数值 恢复出图像块的空间像素表示，对于相邻图像块间的重叠部分， 计算其同一位置的像素值平均值作为最终像素矩阵的计算结果。

9.  根据权利要求1～8之一所述的多模图像融合方法，其特 征在于，步骤5中，所述还原包括对于图像块像素矩阵滑动处理 的重叠部分进行取平均的操作。

10.  根据权利要求1～9之一所述的多模图像融合方法，其特 征在于，该方法中所用到的字典为过完备字典，优选为由小波 域字典和DCT冗余字典混合得到的过完备字典。

基于分阶段多原子正交匹配跟踪的多模图像融合方法
技术领域
本发明属于图像融合技术领域，涉及盲源信号分离，匹配 跟踪，信号分阶段稀疏分解表示，具体涉及医学图像分析、恢 复以及多模图像融合等领域。
背景技术
在医学影像领域，各种先进的医疗成像设备的不断涌现， 为临床诊断提供了更直观的依据，成为现代医疗手段不可或缺 的重要部分，但是不同的成像设备所呈现的不同的多模医学图 像，所提供的某些人体组织或者结构的图像，往往都有各自优 点的缺陷。如何利用各种单模医学图像的优势，获得更全面更 准确的病灶图像信息，对于医生的诊断结果的准确性显得尤为 重要，基于此提出的医学图像融合技术，正是通过获取多模医 学图像的互补信息进行融合，使得融合后的医学图像更清晰， 更完整，更可靠。
医学图像分为解剖图像和功能图像两个部分，常见的解剖 图像有X射线投射成像、计算机断层(CT)、核磁共振成像(MRI) 以及内窥镜等获取的序列图像等；功能图像主要包括正电子发 射计算机断层扫描(PET)、单光子发射计算机断层扫描(SPECT) 等，其中，CT图像具有良好的集合特性和较高的空间密度分辨 率，利于观察骨骼结构信息，但是对于人体组织等表现出较低 的对比度；MRI图像对于人体组织的分辨率较高，可以清晰地 展示出入血管、器官等结构，但是对于骨骼结构的观察敏感度 低。因此，要解决诸如此类的问题就需要将这些不同的图像信 息全部综合在一幅图像上，就是利用医学图像融合技术。
医学图像融合的原理是将不同医学成像设备采集到的关于 同一器官或者组织的图像进行处理，捕捉感兴趣的互补信息， 利用某种融合规则合成包含更多细节信息的融合图像，因此， 采用何种的融合算法、融合规则成为了解决这个问题的关键。
医学图像融合技术作为信息融合技术的一个重要分支，是 90年代才发展起来的医学图像后处理技术之一，将医疗成像手 段与数字图像处理技术结合，成为当代医学图像处理领域的前 沿课题。根据医学图像成像方式的不同，可分为同类方式融合 和交互方式融合，同类方式融合(即单模融合)是指同种成像 方式的图像融合，例如多幅MR图像融合，SPECT图像融合等； 交互方式融合(也称多模融合)是指不同成像方式获得的图像 进行融合，例如CT和PET图像融合，MR和SPECT图像融合等。
在信号处理领域众多应用中，人们希望找到一种稀疏的信 号表示方法，采用稀疏逼近替代原始信号，使得原始信号的能 量集中到少部分系数上，从而获得有效的信号表达方式。
传统的信号表示理论通常是基于正交变换，如离散余弦变 换，离散小波变换，Curvelet变换，Ridgelet变换等等，而这些 变换往往只能捕捉到信号当中某种特定的成分，如离散余弦变 换适合表达周期性的纹理信号，小波变换适合表达边缘或奇异 特征，而Curvelet变换则适合表达曲线轮廓特征，但是，自然界 的许多信号是各种自然现象的混合体，这些混合起来的信号在 单一的正交变换下，其独特的特征往往不能获得正确的表达， 因此，人们希望寻找到一种更为有效的信号分解方式，它能够 有效的提取出信号中所包含的各种特征。
在这种背景下，基于超完备字典的信号稀疏表示方法就诞 生了，1993年Mallat和Zhang首次提出了应用超完备字典对信号 进行稀疏分解的思想，并给出了求解该问题的一种算法，即匹 配跟踪(Matching Pursuit,MP)算法，该问题可以简单描述为： 对于给定信号，从一个超完备字典集合中寻找少数原子构成的 支撑集，在这个支撑子集生成的子空间内对该信号进行分解。
随着近十几年研究，信号分解大致可分为5类，贪婪跟踪算 法(Greedy Pursuit)，凸优化(Convex Optimization)，贝叶斯 框架(Bayesian Framework)，非凸优化(Nonconvex Optimization) 和穷尽搜索(Brute Force)。
在凸优化解决方案当中，基跟踪(Basis Pursuit,BP)算法 是其中最突出的一个，它将稀疏问题的求解转换称为范数约束 条件下的凸优化问题，然后采用线性规划的算法加以解决。然 而，该算法的计算复杂度很高，阻碍了其在大规模系统问题求 解上的应用。
最近，基于贪婪跟踪思想的这一类算法受到广泛的关注， 包括正交匹配跟踪OMP(Orthogonal Matching Pursuit)，硬门 限IHT(Iterative Hard Thresh-holding)，子空间跟踪SP(Subspace  pursuit)，压缩采样匹配跟踪CoSaOMP(Compressive Sampling  Matching Pursuit)和正交匹配跟踪算法OMP(Regularized OMP)。
在贪婪跟踪算法家族当中，正交匹配跟踪(OMP)算法应用 是最为广泛，OMP算法的计算复杂度与原始信号的稀疏度有关， 如果原始信号并不那么稀疏，那么该算法的计算性能就会下降。
目前OMP算法每次只允许一个原子进入OMP算法，的计算 复杂程度高，不适于大规模多模医学图像融合，因此，亟待开 发一种有助于减少迭代次数，提高算法执行效率的图像融合方 法。
发明内容
为了解决上述问题，本发明人进行了锐意研究，结果发现： 将图像分解为固定大小、部分重叠的图像块，并将每块图像块 像素表示成列向量，得到每块图像块像素矩阵，再对每块图像 块像素矩阵进行分离，得到直流分量和残余图像块像素矩阵， 以残余图像块像素矩阵作为待分解信号，以小波变换和离散余 弦变换混合基作为冗余字典集合，在进行正交匹配跟踪信号分 解前，设置第二阈值，在残余图像块的像素矩阵中，选择原子 系数绝对值大于第二阈值的原子构成候选集，在候选集中搜索 与当前残差相干性最大的原子，并将搜索到的原子并入子支撑 集中，再将之从候选集中剔除，在候选集中继续搜索出与上述 原子相关性大于第一阈值的原子，将搜索到的原子从候选集中 剔除，即，对候选集进行自组织搜索裁剪，迭代上述过程，至 候选集为空，获得子支撑集，将子支撑集并入信号的支撑集中， 根据当前信号的支撑集构成的子字典，利用最小二乘法计算当 前稀疏向量，进而更新残差信号，当前残差信号的范数大于逼 近误差，或者当前阶段数小于最大阶段数时，继续上述迭代过 程，否则停止，再根据最终得到的信号的支撑集对两幅图像进 行融合、重构和还原，从而完成了本发明。
本发明的目的在于提供以下方面：
第一方面，本发明提供一种基于分阶段多原子正交匹配跟 踪的多模图像融合方法，其特征在于，该方法采用包括自组织 搜索裁剪亚步骤(2-5)的稀疏分解过程，所述自组织搜索裁剪 亚步骤(2-5)包括以下子步骤：
子步骤(2-5-1)，初始化子支撑集为空集；
子步骤(2-5-2)，在候选集中搜索与当前残差相干性最大 的原子，将搜索到的上述原子并入子支撑集，并将之从候选集 中剔除，更新候选集；
子步骤(2-5-3)，在子步骤(2-5-2)更新的候选集中继续 搜索出与子步骤(2-5-2)中搜索到的原子相关性大于第一阈值 的原子；
子步骤(2-5-4)，所述子步骤(2-5-3)中搜索到的原子构 成冗余原子集合，将所述冗余原子集合从候选集中剔除；
子步骤(2-5-5)，判断：候选集是否满足返回子步骤(2-5-2) 的迭代条件，如果满足，返回子步骤(2-5-2)，如果不满足， 则迭代终止。
第二方面，本发明还提供上述多模图像融合方法，其特征 在于，
子步骤(2-5-3)中，所述第一阈值取值为0.1；
子步骤(2-5-5)中，判断是否返回子步骤(2-5-2)的迭代 条件为：当候选集不为空时，返回子步骤(2-5-2)。
第三方面，本发明还提供上述多模图像融合方法，其特征 在于，图像块像素矩阵都采用包括自组织搜索裁剪步骤的分阶 段多原子正交匹配跟踪算法进行稀疏分解，所述包括自组织搜 索裁剪步骤的分阶段多原子正交匹配跟踪算法包括以下亚步骤：
亚步骤(2-1)，初始化残差信号为待分解信号，初始化信 号的支撑集为空集，初始化当前阶段数为零；
亚步骤(2-2)，利用过完备字典的转置与当前残差信号做 内积计算得到内积系数向量；
亚步骤(2-3)，定义第二阈值为常数与空原子系数绝对值 的最大值之积，
亚步骤(2-4)，在图像块像素矩阵中，选择原子系数绝对 值大于第二阈值的原子，构成候选集；
亚步骤(2-5)，对候选集进行自组织搜索裁剪，生成子支 撑集；
亚步骤(2-6)，将亚步骤(2-5)迭代过程产生的子支撑集 并入到信号的支撑集当中；
亚步骤(2-7)，根据当前信号的支撑集构成的子字典，利 用最小二乘法计算当前稀疏向量x_s；
亚步骤(2-8)，更新当前残差信号，当前残差信号为待分 解信号减去过完备字典与当前稀疏向量之内积所得的差；
亚步骤(2-9)，判断当前残差信号或当前阶段数是否满足 迭代过程返回亚步骤(2-2)的迭代条件，如果满足，返回亚步 骤(2-2)，如果不满足，迭代终止。
第四方面，本发明还提供上述多模图像融合方法，其特征 在于，
亚步骤(2-3)中，常数取值为2.38；
亚步骤(2-4)中，所述原子系数为当前残差信号与过完备 字典中原子的内积；
亚步骤(2-9)中，判断是否返回亚步骤(2-2)的迭代条 件为：当前残差信号的范数大于逼近误差，或者当前阶段数小 于最大阶段数，返回亚步骤(2-2)，其中，所述逼近误差是指 每次迭代的误差残余，所述最大阶段数是指最大迭代次数。
第五方面，本发明还提供上述多模图像融合方法，其特征 在于，该方法包括以下步骤：
步骤1，输入两幅待融合的原始图像A和B，将两幅原始图 像分别划分成图像块，每块图像块的大小一致、任何相邻的两 块图像块都具有重叠部分，并将每块图像块的像素表示成列向 量，得到每块图像块的像素矩阵，即得到图像块的像素矩阵表 示；
步骤2，将上述每块图像块像素矩阵进行分离，得到直流分 量和残余图像块像素矩阵，对得到的残余图像块像素矩阵进行 分阶段多原子正交匹配跟踪自组织搜索裁剪稀疏分解，获得残 余图像块的稀疏表示系数矩阵；
步骤3，根据能量最大原则，对步骤2获得的残余图像块的 稀疏表示系数矩阵进行融合，得到该位置的残余图像块的融合 图像块系数矩阵；对于步骤2分离出的直流分量，采取加权平均 方法进行融合，得到融合后的直流分量；
步骤4，对步骤3融合后的系数矩阵进行重构得到残余图像 块的融合图像块像素矩阵，对所述残余图像块的融合图像块像 素矩阵与所述融合后的直流分量进行合并，得到融合后的图像 块像素矩阵；
步骤5，将步骤4获得的融合图像块像素矩阵还原，得到最 终的融合图像；
其中，待融合的两幅原始图像是针对同一部位的，不同成 像设备提供的多模医学图像；
所述分离出直流分量包括计算各个图像块的平均像素值， 然后从各个图像块的像素值中减去该平均像素值，即直流分量 是各个图像块的平均像素值；
所述能量最大原则是将两幅原始图像的对应位置图像块的 稀疏分解系数绝对值最大者作为融合后图像块对应位置的稀疏 分解系数；
所述重构是重新构成各个图像块的空间像素表示；
所述还原是将图像块的一维表示还原为二维表示。
第六方面，本发明还提供上述多模图像融合方法，其特征 在于，步骤1中，
所述两幅原始图像为CT图像与PET图像，或者CT图像与 MRI图像，或者CT图像与SPET图像，或者MR图像与SPEC图像；
所述划分为滑动处理，所述滑动处理为在长为W、宽为H 的图像中，按一定规律移动一个长为宽为的晶格窗 口，将该窗口对应的像素块转换成块像素矩阵表示中的一列， 然后窗口向右或向下移动一步，直到完成对整幅图像的处理， 优选每次平移两个像素点，即滑动步长为2。
第七方面，本发明还提供上述多模图像融合方法，其特征 在于，步骤3中，图像块I_A的稀疏表示系数矩阵用X_A表示，图 像块I_B的稀疏表示系数矩阵用X_B表示，融合后图像块对应位置 的稀疏表示系数矩阵为：
X F = X A , | X A | ≥ | X B | X B , else ]]>
其中，X_F表示融合后图像块对应位置的稀疏表示系数矩阵；
即，当|X_A|≥|X_B|时，X_F＝X_A；
否则，X_F＝X_B。
第八方面，本发明还提供上述多模图像融合方法，其特征 在于，步骤4中，所述重构的方法为：由融合后的稀疏系数值恢 复出图像块的空间像素表示，对于相邻图像块间的重叠部分， 计算其同一位置的像素值平均值作为最终像素矩阵的计算结果。
第九方面，本发明还提供上述多模图像融合方法，其特征 在于，步骤5中，所述还原包括对于图像块像素矩阵滑动处理的 重叠部分进行取平均的操作。
第十方面，本发明还提供上述多模图像融合方法，其特征 在于，该方法中所用到的字典为过完备字典，优选为由小波域 字典和DCT冗余字典混合得到的过完备字典。
根据本发明提供的基于分阶段多原子正交匹配跟踪的多模 图像融合方法，具有以下有益效果：
(1)本发明采用基于分阶段多原子正交匹配跟踪进行图像 融合的算法，该算法在原始OMP算法的基础上，分阶段进行多 原子正交匹配跟踪过程，不仅极大的提高了OMP算法执行速度， 而且获得了更有效的原始信号稀疏表示的结果；
(2)针对多模医学图像融合的特点，提出针对性的融合策 略，获得了更准确、更可靠、更丰富的融合图像。
附图说明
图1示出根据本发明一种优选实施方式的图像融合方法流 程图；
图2示出根据本发明一种优选实施方式的对图像块像素 矩阵的分阶段多原子正交匹配跟踪稀疏分解的流程图；
图3示出根据本发明一种优选实施式的自组织搜索裁剪 的流程图；
图4a示出MR图像；
图4b示出SPEC图像；
图5示出实施例1融合得到的图像；
图6示出对比例1融合得到的图像；
图7示出对比例2融合得到的图像；
图8示出混合字典的原子函数示意图。
具体实施方式
下面通过对本发明进行详细说明，本发明的特点和优点将 随着这些说明而变得更为清楚、明确。
根据本发明提供的一种基于分阶段多原子正交匹配跟踪的 多模图像融合方法，其特征在于，该方法流程如图3所示，具体 来说包括以下步骤：
步骤1，输入两幅待融合的原始图像A和B，将两幅原始图 像分别划分成图像块，每块图像块的大小一致、任何相邻的两 块图像块都具有重叠部分，并将每块图像块的像素表示成列向 量，得到每块图像块的像素矩阵，即得到图像块的像素矩阵表 示。
在本发明的一个优选实施方式中，所选取的原始图像为针 对同一部位的，不同成像设备提供的多模医学图像，如CT图像 与PET图像，或者CT图像与MRI图像，或者CT图像与SPET图像， 或者MR图像与SPEC图像等。
在进一步的优选实施方式中，输入的两幅原始图像大小相 等，划分方式优选为滑动处理，滑动处理得到的图像块的大小 与冗余字典的原子数目相互匹配，即，当冗余字典大小为 Φ∈R^N×M时，图像块的大小为其中，N是指图像块按 行扫描生成的列向量的维数，M表示字典中列向量的个数，即 字典中原子的个数。
本发明采用的滑动处理为滑动窗口技术，即从原始图像像 素左上角遍历到右下角，将每一块图像块对应像素按照字典顺 序重新排列成矩阵的一列矢量，得到图像块像素矩阵V_A和V_B， V_A和V_B的大小都可以通过如下公式表示:
V A = V B = N × ( W - N + 1 ) ( H - N + 1 ) S 2 ]]>
其中，S表示滑动步长，本发明中优选设定滑动步长为2； W表示待融合的源图像的长，H表示待融合的源图像的宽。
本发明中滑动处理是指：在长为W、宽为H的图像中，按 一定规律移动一个长为宽为的晶格窗口，将该窗口 对应的像素块转换成块像素矩阵表示中的一列，然后窗口向右 或向下移动一步，直到完成对整幅图像的处理；本发明中所述 的滑动步长是指：上述移动一步所滑过的像素点个数，在本发 明中优选每次平移两个像素点，即滑动步长为2。
作为一个优选的实施例，步骤1为，输入两幅大小相等的原 始图像A和B，将所述两幅原始图像划分成固定大小、相互重叠 的图像块并将每块图像块的像素表示成列向量，得 到图像块的像素矩阵表示，可以设定N等于64，即每个原始图 像都划分成若干组8×8像素大小的图像块，可以设定M等于128， 即，冗余字典大小为128；
步骤2，对上述每块图像块像素矩阵进行分离，得到直流分 量和残余图像块像素矩阵，对得到的残余图像块像素矩阵进行 分阶段多原子正交匹配跟踪自组织搜索裁剪稀疏分解，获得残 余图像块的稀疏表示的系数矩阵。
本发明中所述分离出直流分量包括计算各个图像块的平均 像素值，然后从各个图像块的像素值中减去该平均像素值，即 直流分量是指各个图像块的平均像素值。
本发明中所述分阶段多原子正交匹配跟踪算法中，自组织 搜索裁剪稀疏分解是指，对于每个残余图像块像素矩阵对应的 列向量都采用包括自组织搜索裁剪亚步骤的分阶段多原子正交 匹配跟踪算法进行稀疏分解，所述包括自组织搜索裁剪亚步骤 的分阶段多原子正交匹配跟踪算法的流程如图2所示，具体来说， 包括以下亚步骤：
亚步骤(2-1)，初始化残差信号为待分解信号，初始化信 号的支撑集为空集，初始化当前阶段数为零；
亚步骤(2-2)，利用过完备字典的转置与当前残差信号做 内积计算得到内积系数向量；
亚步骤(2-3)，定义第二阈值为常数与空原子系数绝对值 的最大值之积，此处常数取值为2.38，本发明中，利用K-means 聚类算法将内积系数向量的各个原子特征系数分为两类，一类 为空原子系数集合另一类为非空原子系数集合 其中，
表示信号的方差，
μ_s表示非空原子系数绝对值的平均值，
亚步骤(2-4)，在残余图像块像素矩阵中，选择原子系数 绝对值大于第二阈值的原子，构成候选集，其中，所述原子系 数为当前残差信号与过完备字典中原子的内积；
亚步骤(2-5)，对候选集进行自组织搜索裁剪，生成子支 撑集；
本发明人发现，在亚步骤(2-4)中形成的候选集中产生大 量冗余数据，不仅降低运算速度，而且融合图像还原后，图像 失真；而对候选集进行自组织搜索裁剪后，不仅运算速度明显 提高，而且融合图像还原后几乎不失真。
本发明所述自组织搜索裁剪的流程如图1所示，具体来说， 包括以下子步骤：
子步骤(2-5-1)，初始化子支撑集为空集；
子步骤(2-5-2)，在候选集中搜索与当前残差相干性最大 的原子，将搜索到的上述原子并入子支撑集，并将之从候选集 中剔除，更新候选集；
子步骤(2-5-3)，在子步骤(2-5-2)更新的候选集中继续 搜索出与子步骤(2-5-2)中搜索到的原子相关性大于第一阈值 的原子，其中，所述第一阈值取值与算法所采用的冗余字典的 相干性大小有关，本发明设定为0.1；
子步骤(2-5-4)，所述子步骤(2-5-3)中搜索到的原子构 成冗余原子集合，将所述冗余原子集合从候选集中剔除；
子步骤(2-5-5)，判断：候选集是否满足返回子步骤(2-5-2) 的迭代条件，如果满足，返回子步骤(2-5-2)，如果不满足， 则迭代终止；
在一个优选的实施方式中，步骤(2-5-5)中判断是否返回 子步骤(2-5-2)的迭代条件为：当候选集不为空时，返回子 步骤(2-5-2)。
亚步骤(2-6)，将亚步骤(2-5)迭代过程产生的子支撑集 并入到信号的支撑集当中；
亚步骤(2-7)，根据当前信号的支撑集构成的子字典，利 用最小二乘法计算当前稀疏向量x_s；
亚步骤(2-8)，更新当前残差信号，当前残差信号为待分 解信号减去过完备字典与当前稀疏向量之内积所得的差；
亚步骤(2-9)，判断当前残差信号或当前阶段数是否满足 迭代过程返回亚步骤(2-2)的迭代条件，如果满足，返回亚步 骤(2-2)，如果不满足，迭代终止。
在本发明的一个优选实施方式中，亚步骤(2-9)中判断是 否返回亚步骤(2-2)的条件为：当前残差信号的范数大于逼近 误差，或者当前阶段数小于最大阶段数时，返回亚步骤(2-2)， 其中，所述逼近误差是指每次迭代的误差残余，所述最大阶段 数是指最大迭代次数，即，当‖r_s‖>ε或者当前阶段数s<maxS 时，返回亚步骤(2-2)，否则停止，其中，逼近误差为一个预 定值，最大阶段数也为一个预定值，在一个优选的实施方式中， 预定逼近误差为10，预定最大阶段数为5。
如果至少满足上述两种情况中的一种，则判断为可以继续 迭代，返回亚步骤(2-2)，否则迭代终止。
在进一步优选的实施方式中，一个完整的分阶段多原子正 交匹配跟踪算法的迭代过程如下：
定义y为待分解信号，且y∈R^N，定义Φ(N×M)表示为过完 备字典，定义ε表示为逼近误差，定义maxS为最大阶段数，定 义x为稀疏向量，且x∈R^M：
亚步骤(2-1)，初始化残差信号r₀＝y，初始化信号的支撑 集初始化当前阶段数s＝0；
亚步骤(2-2)，利用公式c_s＝Φ^Tr_s计算得到内积系数向量， 其中，c_s表示内积系数向量，Φ^T表示过完备字典的转置，r_s表 示当前残差信号；
亚步骤(2-3)，利用K-means聚类算法将内积系数向量c_s的 各个原子特征系数分为两类，一类为空原子系数集另一类为非空原子系数集定义第二阈值 T＝C·argmax|Z_s|，此处C＝2.38，其中，
表示信号的方差，
μ_s表示非空原子系数绝对值的平均值，
亚步骤(2-4)，在残余图像块像素矩阵中，选择原子系数 绝对值大于第二阈值T的原子，构成候选集 K &OverBar; s = ]]> { k : | C s ( k ) | > T } ; ]]>
亚步骤(2-5)，对候选集进行自组织搜索裁剪，构成子 支撑集K_s，定义自组织搜索裁剪的运算符号为prune_atomset， 其方法如子步骤(2-5-1)～子步骤(2-5-5)所述；
亚步骤(2-6)，将亚步骤(2-5)迭代过程产生的子支撑集 K_s并入到信号的支撑集I_s当中，即，I_s＝I_s-1∪K_s；
亚步骤(2-7)，根据当前信号的支撑集构成的子字典，利 用最小二乘法计算当前稀疏向量其中，
x_s表示当前稀疏向量，
表示由索引属于I_s的原子函数构成的子字典；
亚步骤(2-8)，更新残差信号，公式为r_s＝y-Φx_s；
亚步骤(2-9)，判断‖r_s‖>ε，或者s<maxS，返回亚步骤 (2-2)继续执行，否则停止迭代。
本发明提供的多原子正交匹配跟踪自组织搜索裁剪稀疏分 解与本领域中常用的正交匹配跟踪算法的区别主要集中在于在 算法的每个阶段，通过设置合适的第二阈值，使多个原子同时 选入候选集，相比于传统每次只允许一个原子进入的正交匹配 跟踪算法，该算法允许原子间存在一定的冗余，对此，本发明 通过定义一种自组织搜索裁剪算法prune_atomset，剔除候选集 中相干性较大的原子，获得精简后的子支撑集，再将子支撑集 并入信号的支撑集中，从而保证更有效地获取原始信号的逼近， 并且极大的提高了算法的执行速度。
在一个优选实施方案中，亚步骤(2-5)中自组织搜索裁剪 算法prune_atomset一个完整实施过程如下：
定义网络N＝<V,W,E>，其中，
表示的是网络中的各个结点，即 所有位于候选集当中的原子初始；
w_i＝<r_s,v_i>，表示候选集中原 子与残差信号的相似性度量；
表示的是连接任意两个原子结点之间 边的集合，
其中边的权重e_ij＝<v_i,v_j>，表示两个原子函数之间的相干 程度，n表示结点的最大数目，v_i表示第i个结点，表示初始的 候选集。
定义子支撑集为K_s，定义初始的候选集为
子步骤(2-5-1)，初始化
子步骤(2-5-2)，根据在候选集中搜 索与当前残差相干性最大的原子，令K_s＝K_s∪{v_m}，其中，v_m表示结点m，
子步骤(2-5-3)，根据在候选集中 搜索与子步骤(2-5-2)中搜索到的原子相关性大于第一阈值的 原子；
其中，
H表示搜索到的与结点m相干性大于第一阈值的原子构成 的集合；
theta为度量原子之间相干性的第一阈值，在本子步骤中取 值为0.1，
v_kk表示结点k，；
子步骤(2-5-4)，所述子步骤(2-5-3)中搜索到的原子构 成冗余原子集合H，并将所述冗余原子集合从候选集中剔除， 即令 K &OverBar; s = K &OverBar; s \ H ; ]]>
子步骤(2-5-5)，判断：候选集是否满足返回子步骤(2-5-2) 的迭代条件，如果满足，返回子步骤(2-5-2)，如果不满足， 则迭代终止，
具体为，当候选集不为空时，返回子步骤(2-5-2)。
步骤3，根据能量最大原则，对步骤2获得的残余图像块的 稀疏分解的系数矩阵进行融合，得到该位置的残余图像块的融 合图像块系数矩阵；对于步骤2分离出的直流分量，采取加权平 均方法进行融合，得到融合后的直流分量。
本发明中所述能量最大原则是将两幅原始图像的对应位置 图像块的稀疏分解系数绝对值最大者作为融合后图像块对应位 置的稀疏分解系数，说明在此处绝对值最大的系数所对应的像 素对融合图像的贡献最大。
本发明中所述的对应位置是指两幅待融合原始图像上相对 位置相同的区域，进一步地理解为在上述两幅原始图像上相同 区域处划分出的两个图像块，在两幅原始图像上互称为对应位 置的两个位置处描所显示的图像内容是相同的。
在一个优选的实施方式中，图像块I_A的稀疏表示系数矩阵 用X_A，图像块I_B的稀疏表示系数矩阵用X_B表示，则融合后图像 块对应位置的稀疏表示系数矩阵为：
X F = X A , | X A | &GreaterEqual; | X B | X B , else ]]>
其中，X_F表示融合后图像块对应位置的稀疏表示系数矩阵；
即，当|X_A|≥|X_B|时，X_F＝X_A；
否则，X_F＝X_B。
在进一步优选的实施方式中，对于直流分量部分，所述加 权平均方法如下：
定义图像块I_A的直流分量加权系数为ω_A，图像块I_B的直流分 量加权系数为ω_B，I_A的直流分量为dc_A，I_B的直流分量为dc_B，对 两幅图像I_A和I_B的直流分量加权系数进行融合，则图像块I_A的直 流分量加权系数为：
ω A = | | I A | | 2 2 | | I A | | 2 2 + | | I B | | 2 2 ]]>
图像块I_B的直流分量加权系数为：
ω B = | | I B | | 2 2 | | I A | | 2 2 + | | I B | | 2 2 ; ]]>
则融合后的直流分量部分为：
dc_F＝ω_Adc_A+ω_Bdc_B，其中，
dc_F表示融合后的直流分量。
步骤4，对步骤3融合后的系数矩阵进行重构得到残余图像 块的融合图像块像素矩阵，对所述残余图像块的融合图像块像 素矩阵与所述融合后的直流分量dc_F进行合并，得到融合后的图 像块像素矩阵。
本发明中所述重构是指对融合图像块系数矩阵进行稀疏分 解的逆运算；即，重新构成各个图像块的空间像素表示，重构 方法为：
由融合后的稀疏系数值恢复出图像块的空间像素表示，对 于相邻图像块间的重叠部分，计算其同一位置的像素值平均值 作为最终像素矩阵的计算结果。
在本发明中图像块在像素矩阵中表现为某列向量，融合图 像块系数矩阵的每一列就是分解后的稀疏系数值。
在进一步优选的实施方式中，将融合系数矩阵经过稀疏分 解的逆运算，获得融合图像块像素矩阵，其运算公式如下：
V_F＝Φ·X_F
其中，V_F表示融合图像块像素矩阵，Φ表示稀疏分解所采 用的过完备字典。
步骤5，将步骤4获得的融合图像块像素矩阵还原，得到最 终的融合图像，其中，还原包括对于图像块像素矩阵滑动处理 的重叠部分进行取平均的操作，还原得到的融合图像块像素矩 阵包括：还原融合图像块重叠部分的像素矩阵和直流分量部分 的像素矩阵，所述还原是将图像块的一维表示还原为二维表示； 其中，对于重叠部分的像素矩阵，通过取其像素矩阵平均值的 方式获得融合图像的重叠分量像素矩阵；对于直流分量部分的 像素矩阵，通过对图像块像素矩阵的直流分量进行加权操作的 方式获得融合图像的直流分量像素矩阵。
在一个优选的实施方式中，发明中所用到的字典为过完备 字典，该字典是混合字典，所述混合字典是由小波域字典和DCT 冗余字典混合得到的；其中，字典可以用Φ(N×M)表示，在本 发明中，混合字典如图8中所示，是由小波域字典(Bior6.8) 和DCT冗余字典的混合得到的。
实施例
将图4a与图4b进行融合，得到多模图像，其中，图4a中示 出的图像为MR图像，图4b中示出的图像为SPEC图像，利用相 同的计算机对上述两幅图像用不同的方法进行图像融合。
实施例1
利用本发明提供的多原子正交匹配跟踪、自组织搜索裁剪 进行综合控制，用时9.3秒得到如图5所示的融合图像，该图像 质量指标如下表1所示：
表1融合图像质量指标

对比例
对比例1
利用多原子正交匹配跟踪进行控制，与实施例1的区别在于 图像融合过程中不进行自组织搜索裁剪，用时8.63秒得到如图6 所示的融合图像，该图像质量指标如下表2所示：
表2融合图像质量指标

对比例2
本对比例利用普通正交匹配跟踪进行控制，与实施例1的区 别在于，利用单个原子进行迭代计算进行图像融合，用时25.3 秒得到如图7所示的融合图像，该图像质量指标如下表3所示：
表3融合图像质量指标

其中，QABF是度量从两幅源图像中转移到融合图像中信 息量的大小，LABF是度量融合过程中丢失的信息，其值大小反 应丢失信息的多少，NABF是度量融合过程产生的人工失真效 应，其值越小表示失真效应越小。
QABF为评价融合图像质量的重要指标，在本领域中， QABF每提升0.01个单位都需要克服许多技术难题，而本发明实 施例1相比于对比例1提升了0.1个单位，为显著进步，取得了意 想不到的效果。
通过以上实施例1与对比例1，其在QABF指标上取得了惊 人的效果，从图5及图6也可以明显看出，实施例1融合得到的图 像与原始图像相符度大，图像质量好；而对比例1得到的图像在 图像中部左侧和图像中部右侧均存在白色斑点，这是在图像融 合过程中图像失真造成的，图像质量差。
实施例1和对比例2对比可知，实施例1与对比例2得到的图 像质量相差不大，但实施例1仅用9.3秒，而对比例2却用了25.3 秒，实施例1节约时间近200％。
以上结合具体实施方式和范例性实例对本发明进行了详细 说明，不过这些说明并不能理解为对本发明的限制。本领域技 术人员理解，在不偏离本发明精神和范围的情况下，可以对本 发明技术方案及其实施方式进行多种等价替换、修饰或改进， 这些均落入本发明的范围内。本发明的保护范围以所附权利要 求为准。