CN110084770B - 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法 - Google Patents

基于二维Littlewood-Paley经验小波变换的脑部图像融合方法 Download PDF

Info

Publication number
CN110084770B
CN110084770B CN201910160670.2A CN201910160670A CN110084770B CN 110084770 B CN110084770 B CN 110084770B CN 201910160670 A CN201910160670 A CN 201910160670A CN 110084770 B CN110084770 B CN 110084770B
Authority
CN
China
Prior art keywords
image
sub
images
brain
fused
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910160670.2A
Other languages
English (en)
Other versions
CN110084770A (zh
Inventor
金鑫
江倩
周维
储星
姚绍文
黄锦芳
黄珊珊
武丽雯
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Yunnan University YNU
Original Assignee
Yunnan University YNU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Yunnan University YNU filed Critical Yunnan University YNU
Priority to CN201910160670.2A priority Critical patent/CN110084770B/zh
Publication of CN110084770A publication Critical patent/CN110084770A/zh
Application granted granted Critical
Publication of CN110084770B publication Critical patent/CN110084770B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20064—Wavelet transform [DWT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明公开了一种基于二维Littlewood‑Paley经验小波变换的脑部图像融合方法，获取已经配准好的待融合的脑部CT图像和脑部MRI图像，分别对脑部CT图像和脑部MRI图像进行二维Littlewood‑Paley经验小波变换，得到残差分量子图像、固有模态函数子图像和各子图像对应的检测滤波器图像，分别对残差分量子图像、固有模态函数子图像和检测滤波器图像进行融合，根据融合后的残差分量子图像、融合后的固有模态函数子图像和融合后的检测滤波器图像进行二维Littlewood‑Paley经验小波逆变换，得到融合后的脑部图像。采用本发明可以提高脑部CT图像和脑部MRI图像的融合效果。

Description

基于二维Littlewood-Paley经验小波变换的脑部图像融合 方法
技术领域
本发明属于脑部图像处理技术领域，更为具体地讲，涉及一种基于二维Littlewood-Paley经验小波变换的脑部图像融合方法。
背景技术
医学成像技术可提供准确的人体组织和结构信息，给临床医学带来了极大的便利。但由于传感器成像机制的不同，不同模态的脑医学图像从不同角度表示了人脑的信息，具有一定互补性，且存在冗余信息。在脑医学图像中，以磁共振成像(Magnetic ResonanceImaging,MRI)和计算机断层扫描(Computed Tomography,CT)在医学诊断和治疗中最流行。MRI图像可以清晰显示含水率较高的组织，CT图像可以清晰显示高密度组织，将两种图像融合可以为医生提供更加可靠而丰富的人体生理与解剖信息。因此，图像融合技术在医疗领域的作用较为显著。
当前主要的医学图像融合方法大致可分为空间域和变换域。空间域融合方法可在图像空间域采用特定方法对多源图像直接进行融合处理，经典的方法有：基于神经网络、独立成分分析法、主成分分析法、模糊集的医学图像融合技术等。空间域融合方法通常通过产生与源图像相对应的决策图或权重图实现医学图像融合，此类方法的优点是简单易行，计算量小；然而，其在融合过程中常常会丢失大量图像细节特征，导致所得融合图像边缘和轮廓等细节特征损失较大。变换域融合方法在医学图像融合中较为流行，例如金字塔变换，离散小波变换、轮廓波变换、剪切波变换和tetrolet变换等。然而，大多数基于变换域的方法都有局限性，如，拉普拉斯金字塔换不能准确地描述图像的轮廓和对比度，离散小波变换、轮廓波变换和剪切波变换模型常常在融合后的图像中造成伪影和吉布斯效应，非下采样轮廓变换(NSCT)和非下采样剪切波变换(NSST)因产生相当数量与原始图像大小相同的子图像导致其融合计算量很大。
在经验模态分解(EMD)被提出后，也像其他传统变换方法一样被引入医学图像融合中，但其图像融合性能常常受到经验分解模型容量的影响。近年来提出了许多有前途的EMD模型，为提高基于EMD的医学图像融合的性能提供了很大的可能性。2013年和2014年，Gilles等人提出了基于EMD思想的提出了经验小波分解方法(empirical wavelettransform，EWT)，并迅速应用于青光眼图像的自动诊断、高光谱图像分类、故障诊断等。但对于不同的源图像，EWT分解所得IMF的个数不相等的，为基于EWT的医学图像融合带来了一定挑战，也是基于EWT的图像融合技术首先要解决的问题。
发明内容
本发明的目的在于克服现有技术的不足，提供一种基于二维Littlewood-Paley经验小波变换的脑部图像融合方法，提高脑部CT图像和脑部MRI图像的融合效果。
为实现上述发明目的，本发明基于二维Littlewood-Paley经验小波变换的脑部图像融合方法具体包括以下步骤：
S1：获取已经配准好的待融合的脑部CT图像和脑部MRI图像；
S2：分别对脑部CT图像和脑部MRI图像进行二维Littlewood-Paley经验小波变换，分别得到一组子图像，包括一个残差分量子图像和若干固有模态函数子图像，并将每个子图像对应的检测滤波器采用二维图像表示，记脑部CT图像的残差分量子图像为
Figure BDA0001984524050000021
固有模态函数子图像为Wi CT，i＝1,2,…,NCT，检测滤波器图像为
Figure BDA0001984524050000022
其中NCT表示脑部CT图像的固有模态函数数量；记脑部MRI图像的残差分量子图像为
Figure BDA0001984524050000023
固有模态函数子图像为
Figure BDA0001984524050000024
Figure BDA0001984524050000025
检测滤波器图像为
Figure BDA0001984524050000026
其中NMRI表示脑部CT图像的固有模态函数数量；
S3：分别提取残差子图像
Figure BDA0001984524050000027
和
Figure BDA0001984524050000028
的L2范数特征，基于L2范数特征进行两幅残差子图像的融合，得到融合后的残差子图像
Figure BDA0001984524050000029
S4：采用以下方法进行固有模态函数子图像融合：
如果NCT＞NMRI，则令待融合固有模态函数子图像数量R＝NMRI，融合后固有模态函数子图像总数T＝NCT，否则令待融合固有模态函数子图像数量R＝NCT，融合后固有模态函数子图像总数T＝NMRI；分别取脑部CT图像的前R个固有模态函数子图像和脑部MRI图像的前R个固有模态函数子图像，将对应固有模态序号的子图像作为一组，得到R组待融合固有模态函数子图像，根据以下公式进行加权融合，得到融合后的固有模态函数子图像
Figure BDA0001984524050000031
Figure BDA0001984524050000032
其中，r＝1,2,…,R，
Figure BDA0001984524050000033
表示融合后第r个固有模态函数子图像
Figure BDA0001984524050000034
中像素点(x,y)的像素值，
Figure BDA0001984524050000035
表示脑部CT图像的第r个固有模态函数子图像
Figure BDA0001984524050000036
中像素点(x,y)的像素值，
Figure BDA0001984524050000037
表示脑部MRI图像的第r个固有模态函数子图像
Figure BDA0001984524050000038
中像素点(x,y)的像素值，
Figure BDA0001984524050000039
和
Figure BDA00019845240500000310
分别表示预设的
Figure BDA00019845240500000311
和
Figure BDA00019845240500000312
的权值，且
Figure BDA00019845240500000313
如果NCT＞R，则直接将脑部CT图像的后NCT-R个固有模态函数子图像作为融合后的固有模态函数子图像
Figure BDA00019845240500000314
否则直接将脑部MRI图像的后NMRI-R个固有模态函数子图像作为融合后的固有模态函数子图像
Figure BDA00019845240500000315
从而得到T个融合后的固有模态函数子图像Wt F，t＝1,2,…,T；
S5：采用以下方法进行检测滤波器融合：
分别取脑部CT图像的前R+1个检测滤波器图像
Figure BDA00019845240500000316
和脑部MRI图像的前R+1个检测滤波器图像
Figure BDA00019845240500000317
将对应序号的检测滤波器图像作为一组，得到R+1组待融合检测滤波器图像，基于最大值选择法得到融合后的检测滤波器图像
Figure BDA00019845240500000318
Figure BDA00019845240500000319
其中，
Figure BDA00019845240500000320
分别表示融合后的检测滤波器图像
Figure BDA00019845240500000321
检测滤波器图像
Figure BDA00019845240500000322
检测滤波器图像
Figure BDA00019845240500000323
中像素点(x,y)的像素值；
如果NCT＞R，则直接将脑部CT图像的后NCT-R个检测滤波器图像作为融合后的检测滤波器图像
Figure BDA00019845240500000324
否则直接将脑部MRI图像的后NMRI-R个检测滤波器图像作为融合后的检测滤波器图像
Figure BDA00019845240500000325
从而得到T+1个融合后的检测滤波器图像
Figure BDA00019845240500000326
S6：根据融合后的残差分量子图像
Figure BDA00019845240500000327
融合后的固有模态函数子图像Wt F和融合后的检测滤波器图像
Figure BDA00019845240500000328
进行二维Littlewood-Paley经验小波逆变换，得到融合后的脑部图像。
本发明基于二维Littlewood-Paley经验小波变换的脑部图像融合方法，获取已经配准好的待融合的脑部CT图像和脑部MRI图像，分别对脑部CT图像和脑部MRI图像进行二维Littlewood-Paley经验小波变换，得到残差分量子图像、固有模态函数子图像和各子图像对应的检测滤波器图像，分别对残差分量子图像、固有模态函数子图像和检测滤波器图像进行融合，根据融合后的残差分量子图像、融合后的固有模态函数子图像和融合后的检测滤波器图像进行二维Littlewood-Paley经验小波逆变换，得到融合后的脑部图像。
本发明具有以下有益效果：
1)本发明将Littlewood-Paley EWT引入到医学图像融合中，同时定义了一种有效的计算操作解决了源图像具有不同数量固有模态函数子图像的问题；
2)本发明提出了一种新的二范数图像提取方法来提取残差分量的主要特性，从而实现残差子图像的有效融合；
3)本发明提出了一种基于显著性度量和匹配度量的方法来确定固有模态函数子图像的融合权重，从而实现固有模态函数子图像的有效融合；
4)本发明分别对残差子图像和固有模态函数子图像进行融合，可以同时反应Littlewood-Paley EWT域子图像系数振幅特征和区域特征，从而使最终的融合更加精确。经实验证明，与其他传统方法相比，本发明可将更多的脑部细节组织信息融合到最终图像中，从而有助于医学诊断和治疗。
附图说明
图1是本发明基于二维Littlewood-Paley经验小波变换的脑部图像融合方法的具体实施方式流程图；
图2是本实施例中脑部CT图像和脑部MRI图像二维Littlewood-Paley经验小波变换实例一的变换前后图像；
图3是本实施例中脑部CT图像和脑部MRI图像二维Littlewood-Paley经验小波变换实例二的变换前后图像；
图4是本发明中基于L2范数特征进行残差子图像融合的流程图；
图5是实例一中脑部CT图像残差子图像和脑部MRI图像残差子图像的基于L2范数的特征图；
图6是实例一中脑部CT图像残差子图像和脑部MRI图像残差子图像的融合结果图；
图7是实例二中脑部CT图像残差子图像和脑部MRI图像残差子图像的基于L2范数的特征图；
图8是实例二中脑部CT图像残差子图像和脑部MRI图像残差子图像的融合结果图；
图9是本实施例中权值
Figure BDA0001984524050000051
和
Figure BDA0001984524050000052
优选设置方法的流程图；
图10是实例一中脑部CT图像的前两幅固有模态函数子图像和脑部MRI图像的两幅固有模态函数子图像的融合结果图；
图11是实例二中脑部CT图像的两幅固有模态函数子图像和脑部MRI图像的两幅固有模态函数子图像的融合结果图；
图12是第一组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图；
图13是第二组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图；
图14是第三组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图；
图15是第四组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图；
图16是第五组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图；
图17是第六组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图。
具体实施方式
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
实施例
图1是本发明基于二维Littlewood-Paley经验小波变换的脑部图像融合方法的具体实施方式流程图。如图1所示，本发明基于二维Littlewood-Paley经验小波变换的脑部图像融合方法的具体步骤包括：
S101：获取已配准图像：
获取已经配准好的待融合的脑部CT图像和脑部MRI图像。
S102：二维Littlewood-Paley经验小波变换：
分别对脑部CT图像和脑部MRI图像进行二维Littlewood-Paley经验小波变换，分别得到一组子图像，包括一个残差分量子图像和若干固有模态函数子图像，并将每个子图像对应的检测滤波器采用二维图像表示，记脑部CT图像的残差分量子图像为
Figure BDA0001984524050000061
固有模态函数子图像为Wi CT，i＝1,2,...,NCT，检测滤波器图像为
Figure BDA0001984524050000062
其中NCT表示脑部CT图像的固有模态函数数量；记滤波器二维记脑部MRI图像的残差分量子图像为
Figure BDA0001984524050000063
固有模态函数子图像为
Figure BDA0001984524050000064
检测滤波器图像为
Figure BDA0001984524050000065
其中NMRI表示脑部CT图像的固有模态函数数量。
二维Littlewood-Paley经验小波变换(2D Littlewood-Paley empiricalwavelet transform，LPEWT)是一种常用的二维经验小波变换，其具体变换方法可参考文献“J.Gilles,G.Tran,S.Osher"2D Empirical transforms.Wavelets,Ridgelets andCurvelets Revisited",submitted in SIAM Journal on Imaging Sciences,2013.”
图2是本实施例中脑部CT图像和脑部MRI图像二维Littlewood-Paley经验小波变换实例一的变换前后图像。如图2所示，实例一中脑部CT图像转换得到的固有模态函数子图像有3幅，脑部MRI图像转换得到的固有模态函数子图像有2幅。
图3是本实施例中脑部CT图像和脑部MRI图像二维Littlewood-Paley经验小波变换实例二的变换前后图像。如图3所示，实例二中脑部CT图像和脑部MRI图像转换得到的固有模态函数子图像的数量均为2幅。
S103：基于L2范数特征进行残差子图像融合：
分别提取残差子图像
Figure BDA0001984524050000066
和
Figure BDA0001984524050000067
的L2范数特征，基于L2范数特征进行两幅残差子图像的融合，得到融合后的残差子图像
Figure BDA0001984524050000068
图4是本发明中基于L2范数特征进行残差子图像融合的流程图。如图4所示，本发明中基于L2范数特征进行残差子图像融合的具体步骤包括：
S401：提取行L2范数：
对于残差子图像
Figure BDA0001984524050000071
中的各个像素点(x,y)，ω∈{CT,MRI}，x＝1,2,...,W，y＝1,2,...,H，W×H表示图像大小，分别以该像素点(x,y)为中心提取出大小为p×q的图像块，p和q的大小根据需要确定，记该图像块左上顶点的坐标为(x0,y0)，根据以下公式计算得到像素点(x,y)的行L2范数
Figure BDA0001984524050000072
Figure BDA0001984524050000073
其中，imω()表示对应像素点在残差子图像
Figure BDA0001984524050000074
中像素值，当像素点(x0+m-1,y0+n-1)或像素点(x0+m-1,y0+n-2)不存在(即超出残差子图像
Figure BDA00019845240500000717
范围)时，令对应的imω(x0+m-1,y0+n-1)或imω(x0+m-1,y0+n-2)为0。
S402：提取列L2范数：
根据以下公式计算得到残差子图像
Figure BDA0001984524050000075
中各个像素点(x,y)的列L2范数
Figure BDA0001984524050000076
Figure BDA0001984524050000077
类似地，当像素点(x0+m-1,y0+n-1)或像素点(x0+m-2,y0+n-1)不存在(即超出残差子图像
Figure BDA0001984524050000078
范围)时，令对应的imω(x0+m-1,y0+n-1)或imω(x0+m-2,y0+n-1)为0。
S403：计算L2范数特征：
根据以下公式计算得到残差子图像
Figure BDA0001984524050000079
中各个像素点(x,y)的L2范数特征LNω(x,y)：
Figure BDA00019845240500000710
S404：残差子图像融合：
采用L2范数特征可以通过子图像像素或系数振幅表示出医学图像的结构信息，以便融合不同源图像的残差分量。基于L2范数特征的残差子图像融合的具体公式如下：
Figure BDA00019845240500000711
其中，
Figure BDA00019845240500000712
表示融合后残差子图像
Figure BDA00019845240500000713
中像素点(x,y)的像素值，
Figure BDA00019845240500000714
表示残差子图像
Figure BDA00019845240500000715
中像素点(x,y)的像素值，
Figure BDA00019845240500000716
表示残差子图像
Figure BDA0001984524050000081
中像素点(x,y)的像素值。
图5是实例一中脑部CT图像残差子图像和脑部MRI图像残差子图像的基于L2范数的特征图。图6是实例一中脑部CT图像残差子图像和脑部MRI图像残差子图像的融合结果图。
图7是实例二中脑部CT图像残差子图像和脑部MRI图像残差子图像的基于L2范数的特征图。图8是实例二中脑部CT图像残差子图像和脑部MRI图像残差子图像的融合结果图。
S104：固有模态函数子图像融合：
如果NCT＞NMRI，则令待融合固有模态函数子图像数量R＝NMRI，融合后固有模态函数子图像总数T＝NCT，否则令待融合固有模态函数子图像数量R＝NCT，融合后固有模态函数子图像总数T＝NMRI。分别取脑部CT图像的前R个固有模态函数子图像和脑部MRI图像的前R个固有模态函数子图像，将对应固有模态序号的子图像作为一组，得到R组待融合固有模态函数子图像，根据以下公式进行加权融合，得到融合后的固有模态函数子图像
Figure BDA0001984524050000082
Figure BDA0001984524050000083
其中，r＝1,2,...,R，
Figure BDA0001984524050000084
表示融合后第r个固有模态函数子图像
Figure BDA0001984524050000085
中像素点(x,y)的像素值，
Figure BDA0001984524050000086
表示脑部CT图像的第r个固有模态函数子图像
Figure BDA0001984524050000087
中像素点(x,y)的像素值，
Figure BDA0001984524050000088
表示脑部MRI图像的第r个固有模态函数子图像
Figure BDA0001984524050000089
中像素点(x,y)的像素值，
Figure BDA00019845240500000810
和
Figure BDA00019845240500000811
分别表示预设的
Figure BDA00019845240500000812
和
Figure BDA00019845240500000813
的权值，且
Figure BDA00019845240500000814
如果NCT＞R，则直接将脑部CT图像的后NCT-R个固有模态函数子图像作为融合后的固有模态函数子图像
Figure BDA00019845240500000815
否则直接将脑部MRI图像的后NMRI-R个固有模态函数子图像作为融合后的固有模态函数子图像
Figure BDA00019845240500000816
从而得到T个融合后的固有模态函数子图像Wt F，t＝1,2,…,T。
对于加权融合公式，权值
Figure BDA00019845240500000817
和
Figure BDA00019845240500000818
的取值对融合后的固有模态函数子图像具有较大影响，因此需要对权值
Figure BDA00019845240500000819
和
Figure BDA00019845240500000820
进行优选设置。图9是本实施例中权值
Figure BDA00019845240500000821
和
Figure BDA00019845240500000822
优选设置方法的流程图。如图9所示，本实施例中权值
Figure BDA00019845240500000823
和
Figure BDA00019845240500000824
优选设置方法的具体步骤包括：
S901：计算显著性度量：
显著性度量可以表示IMF系数的重要性，一般来说，重要系数往往比不重要系数更加显著。本发明中需要分别计算固有模态函数子图像
Figure BDA0001984524050000091
和固有模态函数子图像
Figure BDA0001984524050000092
的显著性。首先设置一个大小为p′×q′的全1矩阵P，p′和q′的大小根据需要确定，对于固有模态函数子图像
Figure BDA0001984524050000093
中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的显著性度量
Figure BDA0001984524050000094
Figure BDA0001984524050000095
式中，P(m′,n′)表示全1矩阵P中坐标为(m′,n′)的元素值，即P(m′,n′)＝1。
Figure BDA0001984524050000096
表示像素点(x-m′,y-n′)在固有模态函数子图像
Figure BDA0001984524050000097
中对应的像素值，当(x-m′,y-n′)不存在(即超出固有模态函数子图像
Figure BDA0001984524050000098
范围)时，对应
Figure BDA0001984524050000099
为0。
S902：计算匹配度量：
匹配度量不但可以表示固有模态函数系数的振幅信息，同时描述了两个固有模态函数子图像的区域特征。对于固有模态函数子图像
Figure BDA00019845240500000910
和固有模态函数子图像
Figure BDA00019845240500000911
中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的匹配度量ψr(x,y)：
Figure BDA00019845240500000912
式中，
Figure BDA00019845240500000913
Figure BDA00019845240500000914
为像素点(x-m′,y-n′)在固有模态函数子图像
Figure BDA00019845240500000915
中的像素值，
Figure BDA00019845240500000916
为像素点(x-m′,y-n′)在固有模态函数子图像
Figure BDA00019845240500000917
中的像素值；
Figure BDA00019845240500000918
为固有模态函数子图像
Figure BDA00019845240500000919
中像素点(x,y)的显著性度量；
Figure BDA00019845240500000920
为固有模态函数子图像
Figure BDA00019845240500000921
中像素点(x,y)的显著性度量；ε是一个预设常数，其值很小，用于避免分母为0。
S903：判断是否ψr(x,y)＞α，α表示预设的匹配度量阈值，本实施例中设置α＝0.75，如果是，进入步骤S904，否则进入步骤S905。
S904：设置权值：
令
Figure BDA00019845240500000922
根据
Figure BDA00019845240500000923
确定权值。
近年来，模糊集理论在图像融合中得到了广泛的应用，并取得了较好的效果。模糊集理论的隶属度函数(MF)量化了某个元素属于特定集合的隶属关系，MF的值区间为[0,1]；“0”表示完全不隶属，“1”表示完全隶属，“0”与“l”之间的值表示不同程度的隶属关系。由于固有模态函数系数分布服从高斯分布，因此本步骤中为了更加合理地确定权值，可以采用高斯隶属度函数权重，从而确定权值
Figure BDA0001984524050000101
和
Figure BDA0001984524050000102
其具体方法如下：
计算得到固有模态函数子图像
Figure BDA0001984524050000103
和固有模态函数子图像
Figure BDA0001984524050000104
中各像素点匹配度量ψr(x,y)的标准差σ和均值c，计算得到以下参数：
Figure BDA0001984524050000105
Figure BDA0001984524050000106
其中，e表示自然常数。
如果
Figure BDA0001984524050000107
则令
Figure BDA0001984524050000108
否则令
Figure BDA0001984524050000109
S905：判断是否
Figure BDA00019845240500001010
如果是，进入步骤S906，否则进入步骤S907。
S906：令
Figure BDA00019845240500001011
S907：令
Figure BDA00019845240500001012
图10是实例一中脑部CT图像的前两幅固有模态函数子图像和脑部MRI图像的两幅固有模态函数子图像的融合结果图。图11是实例二中脑部CT图像的两幅固有模态函数子图像和脑部MRI图像的两幅固有模态函数子图像的融合结果图。
S105：检测滤波器融合：
在二维Littlewood-Paley经验小波变换中，不同医学图像得到的检测滤波器是一个二进制图像，可以看作是一个决策图。检测滤波器有助于保持医学图像的细节特征，对医学图像重建具有重要意义。由于二维Littlewood-Paley经验小波变换中检测滤波器数量与医学图像对应的子图像(残差子图像和固有模态函数子图像)数量相同，因此本发明采用了一种简单而有效的方法来融合不同源图像的检测滤波器，即最大值选择法，其具体方法如下：
与固有模态函数子图像融合类似，分别取脑部CT图像的前R+1个检测滤波器图像
Figure BDA0001984524050000111
和脑部MRI图像的前R+1个检测滤波器图像
Figure BDA0001984524050000112
将对应序号的检测滤波器图像作为一组，得到R+1组待融合检测滤波器图像，基于最大值选择法得到融合后的检测滤波器图像
Figure BDA0001984524050000113
Figure BDA0001984524050000114
其中，
Figure BDA0001984524050000115
分别表示融合后的检测滤波器图像
Figure BDA0001984524050000116
检测滤波器图像
Figure BDA0001984524050000117
检测滤波器图像
Figure BDA0001984524050000118
中像素点(x,y)的像素值。
如果NCT＞R，则直接将脑部CT图像的后NCT-R个检测滤波器图像作为融合后的检测滤波器图像
Figure BDA0001984524050000119
否则直接将脑部MRI图像的后NMRI-R个检测滤波器图像作为融合后的检测滤波器图像
Figure BDA00019845240500001110
从而得到T+1个融合后的检测滤波器图像
Figure BDA00019845240500001111
S106：二维Littlewood-Paley经验小波逆变换：
根据融合后的残差分量子图像
Figure BDA00019845240500001112
融合后的固有模态函数子图像Wt F和融合后的检测滤波器图像
Figure BDA00019845240500001113
进行二维Littlewood-Paley经验小波逆变换，得到融合后的脑部图像。二维Littlewood-Paley经验小波逆变换具体变换方法同样可参考文献“J.Gilles,G.Tran,S.Osher"2D Empirical transforms.Wavelets,Ridgelets and CurveletsRevisited",submitted in SIAM Journal on Imaging Sciences,2013.”
为了更好地说明本发明的技术效果，采用一个具体实施例进行实验。本实施例中以http://www.med.harvard.edu/aanlib/home.html网站上的医学图像为例作为素材，对比本发明和目前流行的几种图像融合方法的融合效果进行对比，对比方法包括梯度金字塔(GAP)、离散小波变换(DWT)、滤波减法-抽取金字塔(FSDP)、引导图像滤波和图像统计(GFS)、稀疏表示(CSR)、自适应稀疏表示(ASR)、引导滤波融合(GFF)、多分辨率奇异值分解(MSVD)、双树复离散小波变换(DTDWT)、平稳小波变换(SWT)、曲波变换(CVT)、曲波变换与稀疏表示(DTCWT-SR)、曲波与稀疏表示(CVT-SR)。
图12是第一组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图。如图12所示，采用本发明生成的融合图像的亮度高于其他方法，其中头骨的尖锐边缘被完整保留，脑细胞的纹理也被有效地融合到最终图像中。因此，本发明生成的融合图像显然优于其他对比方法。
图13是第二组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图。如图13所示，采用本发明生成的融合图像清晰地描绘了颅骨和大脑的亮边，同时也保留了脑白质和脑灰质的纹理和边缘，比其他方法更为优越。
图14是第三组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图。如图14所示，与其他方法相比，采用本发明生成的融合图像与源图像的对比度最为接近，此外，本发明能够对脑部CT图像和脑部MRI图像中的人体组织信息进行有效描述。
图15是第四组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图。如图15所示，采用本发明生成的融合图像的对比度明显优于其他方法，可以为医生提供更多有价值的信息。本发明保留了脑灰质和脑白质的纹理和边缘，而且在融合图像中提供了眼部含水量信息。
图16是第五组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图。如图16所示，采用本发明以及GFF方法生成的融合图像的对比度、亮度和边缘细节都优于其他方法。然而，GFF方法生成的融合图像却丢失了一些重要的特征，例如脑部MRI图像含水量的信息。总体来看，本发明生成的融合图像还是明显优于其他的对比方法。
图17是第六组脑部CT图像和脑部MRI图像的原图像和采用本发明和对比方法得到的融合图像对比图。如图17所示，采用本发明生成的融合图像的对比度和亮度均优于其他方法，而且没有引入任何误差信息。而例如GFS方法生成的融合图像，其在脑灰质和脑白质区域有明显的伪影。因此本发明生成的融合图像将脑医学图像的细节特征和结构信息保存得很好，优于其他方法。
为了从客观角度评价图像融合质量，本发明采用一些常用的客观评价指标来检验本发明的性能，包括：基于相似性度量的边缘(Qabf)、互信息(MI)、标准差(STD)、熵(En)、空间频率(SF)、特征互信息(FMI)、平均梯度(AG)、平均值(MV)。表1为采用不同客观评价指标对本实施例所得到的融合图像的指标平均值。
Figure BDA0001984524050000131
表1
与其他传统方法相比，本发明的STD、SF、M值最高，更具竞争力。Qabf和MI是图像融合中最重要的评价指标，本发明所得融合图像的这两个指标明显高于其他方法，说明本发明可以保留更多的医学图像信息。虽然本发明的EN、FMI和AV值与其中一些对比方法较为接近，但仍略优于其他对比方法。综上可知，本实施例表明本发明生成的融合图像优于其他方法，是一种有效的脑部图像融合方法。
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于二维Littlewood-Paley经验小波变换的脑部图像融合方法，其特征在于，包括以下步骤：
S1：获取已经配准好的待融合的脑部CT图像和脑部MRI图像；
S2：分别对脑部CT图像和脑部MRI图像进行二维Littlewood-Paley经验小波变换，分别得到一组子图像，包括一个残差分量子图像和若干固有模态函数子图像，并将每个子图像对应的检测滤波器采用二维图像表示，记脑部CT图像的残差分量子图像为W0 CT，固有模态函数子图像为Wi CT，i＝1,2,…,NCT，检测滤波器图像为
Figure FDA0003959055870000011
其中NCT表示脑部CT图像的固有模态函数数量；记脑部MRI图像的残差分量子图像为
Figure FDA0003959055870000012
固有模态函数子图像为
Figure FDA0003959055870000013
Figure FDA0003959055870000014
检测滤波器图像为
Figure FDA0003959055870000015
其中NMRI表示脑部CT图像的固有模态函数数量；
S3：分别提取残差子图像
Figure FDA0003959055870000016
和
Figure FDA0003959055870000017
的L2范数特征，基于L2范数特征进行两幅残差子图像的融合，得到融合后的残差子图像
Figure FDA0003959055870000018
S4：采用以下方法进行固有模态函数子图像融合：
如果NCT＞NMRI，则令待融合固有模态函数子图像数量R＝NMRI，融合后固有模态函数子图像总数T＝NCT，否则令待融合固有模态函数子图像数量R＝NCT，融合后固有模态函数子图像总数T＝NMRI；分别取脑部CT图像的前R个固有模态函数子图像和脑部MRI图像的前R个固有模态函数子图像，将对应固有模态序号的子图像作为一组，得到R组待融合固有模态函数子图像，根据以下公式进行加权融合，得到融合后的固有模态函数子图像
Figure FDA0003959055870000019
Figure FDA00039590558700000110
其中，r＝1,2,…,R，
Figure FDA00039590558700000111
表示融合后第r个固有模态函数子图像
Figure FDA00039590558700000112
中像素点(x,y)的像素值，
Figure FDA00039590558700000113
表示脑部CT图像的第r个固有模态函数子图像Wr CT中像素点(x,y)的像素值，
Figure FDA00039590558700000114
表示脑部MRI图像的第r个固有模态函数子图像
Figure FDA00039590558700000115
中像素点(x,y)的像素值，
Figure FDA00039590558700000116
和
Figure FDA00039590558700000117
分别表示Wr CT(x,y)和
Figure FDA00039590558700000118
的权值，且
Figure FDA00039590558700000119
权值
Figure FDA00039590558700000120
和
Figure FDA00039590558700000121
按照以下方法设置：
S4.1：设置一个大小为p′×q′的全1矩阵P，p′和q′的大小根据需要确定，对于固有模态函数子图像Wr ω中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的显著性度量
Figure FDA0003959055870000021
Figure FDA0003959055870000022
式中，P(m′,n′)表示全1矩阵P中坐标为(m′,n′)的元素值，即P(m′,n′)＝1；
Figure FDA0003959055870000023
Figure FDA0003959055870000024
表示像素点(x-m′,y-n′)在固有模态函数子图像Wr ω中对应的像素值，当(x-m′,y-n′)不存在时，对应
Figure FDA0003959055870000025
为0；
S4.2：对于固有模态函数子图像Wr CT和固有模态函数子图像Wr MRI中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的匹配度量ψr(x,y)：
Figure FDA0003959055870000026
式中，
Figure FDA0003959055870000027
Figure FDA0003959055870000028
为像素点(x-m′,y-n′)在固有模态函数子图像Wr CT中的像素值，
Figure FDA0003959055870000029
为像素点(x-m′,y-n′)在固有模态函数子图像
Figure FDA00039590558700000210
中的像素值；
Figure FDA00039590558700000211
为固有模态函数子图像Wr CT中像素点(x,y)的显著性度量；
Figure FDA00039590558700000212
为固有模态函数子图像
Figure FDA00039590558700000213
中像素点(x,y)的显著性度量；ε是一个预设常数；
S4.3：判断是否ψr(x,y)＞α，α表示预设的匹配度量阈值，如果是，进入步骤S4.4，否则进入步骤S4.5；
S4.4：令
Figure FDA00039590558700000214
根据
Figure FDA00039590558700000215
确定权值；
S4.5：判断是否
Figure FDA00039590558700000216
如果是，进入步骤S4.6，否则进入步骤S4.7；
S4.6：令
Figure FDA00039590558700000217
S4.7：令
Figure FDA00039590558700000218
如果NCT＞R，则直接将脑部CT图像的后NCT-R个固有模态函数子图像作为融合后的固有模态函数子图像
Figure FDA00039590558700000219
否则直接将脑部MRI图像的后NMRI-R个固有模态函数子图像作为融合后的固有模态函数子图像
Figure FDA00039590558700000220
从而得到T个融合后的固有模态函数子图像Wt F，t＝1,2,…,T；
S5：采用以下方法进行检测滤波器融合：
分别取脑部CT图像的前R+1个检测滤波器图像
Figure FDA0003959055870000031
和脑部MRI图像的前R+1个检测滤波器图像
Figure FDA0003959055870000032
将对应序号的检测滤波器图像作为一组，得到R+1组待融合检测滤波器图像，基于最大值选择法得到融合后的检测滤波器图像
Figure FDA0003959055870000033
Figure FDA0003959055870000034
其中，
Figure FDA0003959055870000035
分别表示融合后的检测滤波器图像
Figure FDA0003959055870000036
检测滤波器图像
Figure FDA0003959055870000037
检测滤波器图像
Figure FDA0003959055870000038
中像素点(x,y)的像素值；
如果NCT＞R，则直接将脑部CT图像的后NCT-R个检测滤波器图像作为融合后的检测滤波器图像
Figure FDA0003959055870000039
否则直接将脑部MRI图像的后NMRI-R个检测滤波器图像作为融合后的检测滤波器图像
Figure FDA00039590558700000310
从而得到T个融合后的检测滤波器图像
Figure FDA00039590558700000311
S6：根据融合后的残差分量子图像
Figure FDA00039590558700000312
融合后的固有模态函数子图像Wt F和融合后的检测滤波器图像
Figure FDA00039590558700000313
进行二维Littlewood-Paley经验小波逆变换，得到融合后的脑部图像。
2.根据权利要求1所述的脑部图像融合方法，其特征在于，所述步骤S3中融合得到残差子图像
Figure FDA00039590558700000314
的具体步骤包括：
S3.1：对于残差子图像
Figure FDA00039590558700000315
中的各个像素点(x,y)，ω∈{CT,MRI}，分别以该像素点(x,y)为中心提取出大小为p×q的图像块，记该图像块左上顶点的坐标为(x0,y0)，根据以下公式计算得到像素点(x,y)的行L2范数
Figure FDA00039590558700000316
Figure FDA00039590558700000317
其中，imω()表示对应像素点在残差子图像
Figure FDA00039590558700000318
中像素值，当像素点(x0+m-1,y0+n-1)或像素点(x0+m-1,y0+n-2)不存在时，令对应的imω(x0+m-1,y0+n-1)或imω(x0+m-1,y0+n-2)为0；
S3.2：根据以下公式计算得到残差子图像
Figure FDA00039590558700000319
中各个像素点(x,y)的列L2范数
Figure FDA00039590558700000320
Figure FDA0003959055870000041
当像素点(x0+m-1,y0+n-1)或像素点(x0+m-2,y0+n-1)不存在，令对应的imω(x0+m-1,y0+n-1)或imω(x0+m-2,y0+n-1)为0；
S3.3：根据以下公式计算得到残差子图像W0 ω中各个像素点(x,y)的L2范数特征LNω(x,y)：
Figure FDA0003959055870000042
S3.4：根据以下公式进行残差子图像融合：
Figure FDA0003959055870000043
其中，
Figure FDA0003959055870000044
表示融合后残差子图像
Figure FDA0003959055870000045
中像素点(x,y)的像素值，
Figure FDA0003959055870000046
表示残差子图像
Figure FDA0003959055870000047
中像素点(x,y)的像素值，
Figure FDA0003959055870000048
表示残差子图像
Figure FDA0003959055870000049
中像素点(x,y)的像素值。
3.根据权利要求1所述的脑部图像融合方法，其特征在于，所述步骤S4.4中权值
Figure FDA00039590558700000410
和
Figure FDA00039590558700000411
按照以下方法设置：
计算得到固有模态函数子图像
Figure FDA00039590558700000412
和固有模态函数子图像
Figure FDA00039590558700000413
中各像素点匹配度量ψr(x,y)的标准差σ和均值c，计算得到以下参数：
Figure FDA00039590558700000414
Figure FDA00039590558700000415
如果
Figure FDA00039590558700000416
则令
Figure FDA00039590558700000417
否则令
Figure FDA00039590558700000418
CN201910160670.2A 2019-03-04 2019-03-04 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法 Active CN110084770B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910160670.2A CN110084770B (zh) 2019-03-04 2019-03-04 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910160670.2A CN110084770B (zh) 2019-03-04 2019-03-04 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法

Publications (2)

Publication Number Publication Date
CN110084770A CN110084770A (zh) 2019-08-02
CN110084770B true CN110084770B (zh) 2023-03-07

Family

ID=67413123

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910160670.2A Active CN110084770B (zh) 2019-03-04 2019-03-04 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法

Country Status (1)

Country Link
CN (1) CN110084770B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN118037560A (zh) * 2024-01-16 2024-05-14 北京长木谷医疗科技股份有限公司 基于同态滤波的多模态医学图像融合方法、装置及设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102129676A (zh) * 2010-01-19 2011-07-20 中国科学院空间科学与应用研究中心 一种基于二维经验模态分解的显微图像融合方法
CN104156918A (zh) * 2014-08-01 2014-11-19 西安电子科技大学 基于联合稀疏表示与残差融合的sar图像噪声抑制方法
CN107451595A (zh) * 2017-08-04 2017-12-08 河海大学 基于混合算法的红外图像显著性区域检测方法
CN109242888A (zh) * 2018-09-03 2019-01-18 中国科学院光电技术研究所 一种结合图像显著性和非下采样轮廓波变换的红外与可见光图像融合方法
CN109359607A (zh) * 2018-10-25 2019-02-19 辽宁工程技术大学 一种基于纹理的掌纹掌脉融合识别方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
TWI446897B (zh) * 2011-08-19 2014-08-01 Ind Tech Res Inst 超音波影像對齊裝置及其方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102129676A (zh) * 2010-01-19 2011-07-20 中国科学院空间科学与应用研究中心 一种基于二维经验模态分解的显微图像融合方法
CN104156918A (zh) * 2014-08-01 2014-11-19 西安电子科技大学 基于联合稀疏表示与残差融合的sar图像噪声抑制方法
CN107451595A (zh) * 2017-08-04 2017-12-08 河海大学 基于混合算法的红外图像显著性区域检测方法
CN109242888A (zh) * 2018-09-03 2019-01-18 中国科学院光电技术研究所 一种结合图像显著性和非下采样轮廓波变换的红外与可见光图像融合方法
CN109359607A (zh) * 2018-10-25 2019-02-19 辽宁工程技术大学 一种基于纹理的掌纹掌脉融合识别方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Multilinear commutators of Littlewood-Paley operator on herz-hardy spaces;Kaibin Zhao;《2011 International Conference on Multimedia Technology》;20110825;全文 *
图像融合及其性能评估若干问题研究;张小利;《中国博士学位论文全文数据库信息科技辑》;20160815;摘要，第27-114页 *
基于经验模态分解的图像融合研究;晁永国;《现代商贸工业》;20100215(第04期);第47-48页 *

Also Published As

Publication number Publication date
CN110084770A (zh) 2019-08-02

Similar Documents

Publication Publication Date Title
Xia et al. A novel improved deep convolutional neural network model for medical image fusion
CN106339998B (zh) 基于对比度金字塔变换的多聚焦图像融合方法
CN109242888B (zh) 一种结合图像显著性和非下采样轮廓波变换的红外与可见光图像融合方法
CN109389585B (zh) 一种基于全卷积神经网络的脑组织提取方法
Gao et al. A deep learning based approach to classification of CT brain images
CN108399611B (zh) 基于梯度正则化的多聚焦图像融合方法
Yang Multimodal medical image fusion through a new DWT based technique
Sapkal et al. Image fusion based on wavelet transform for medical application
Jeevakala Sharpening enhancement technique for MR images to enhance the segmentation
CN105894483A (zh) 一种基于多尺度图像分析和块一致性验证的多聚焦图像融合方法
CN116664462B (zh) 一种基于ms-dsc和i_cbam的红外和可见光图像融合方法
Chen et al. IOSUDA: an unsupervised domain adaptation with input and output space alignment for joint optic disc and cup segmentation
CN116452618A (zh) 一种三输入脊柱ct图像分割方法
CN114693682A (zh) 一种基于图像处理的脊椎特征识别方法
CN110084770B (zh) 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法
Zhao et al. A multi-module medical image fusion method based on non-subsampled shear wave transformation and convolutional neural network
Vani et al. Multi focus and multi modal image fusion using wavelet transform
CN117217997A (zh) 一种基于上下文感知边缘增强的遥感图像超分辨率方法
Teng et al. Wavelet-based texture fusion of CT/MRI images
CN115410032A (zh) 基于自监督学习的octa图像分类结构训练方法
Wan et al. Multi-focus color image fusion based on quaternion multi-scale singular value decomposition
Wang et al. Multi-modality anatomical and functional medical image fusion based on simplified-spatial frequency-pulse coupled neural networks and region energy-weighted average strategy in non-sub sampled contourlet transform domain
Abood Edges enhancement of medical color images using add images
Liang MRI image reconstruction based on artificial intelligence
Sulaiman et al. A Convolutional Neural Network Model for Image Enhancement of Extremely Dense Breast Tissue in Digital Breast Tomosynthesis Images

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant