CN110084770A - 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法 - Google Patents
基于二维Littlewood-Paley经验小波变换的脑部图像融合方法 Download PDFInfo
- Publication number
- CN110084770A CN110084770A CN201910160670.2A CN201910160670A CN110084770A CN 110084770 A CN110084770 A CN 110084770A CN 201910160670 A CN201910160670 A CN 201910160670A CN 110084770 A CN110084770 A CN 110084770A
- Authority
- CN
- China
- Prior art keywords
- image
- subgraph
- brain
- mode function
- intrinsic mode
- 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.)
- Granted
Links
- 210000004556 brain Anatomy 0.000 title claims abstract description 141
- 230000009466 transformation Effects 0.000 title claims abstract description 35
- 238000007500 overflow downdraw method Methods 0.000 title claims abstract description 17
- 238000001514 detection method Methods 0.000 claims abstract description 61
- 238000000034 method Methods 0.000 claims description 55
- 230000004927 fusion Effects 0.000 claims description 52
- 238000013421 nuclear magnetic resonance imaging Methods 0.000 claims description 11
- 238000005259 measurement Methods 0.000 claims description 9
- 239000011159 matrix material Substances 0.000 claims description 4
- 238000007689 inspection Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 8
- 238000002595 magnetic resonance imaging Methods 0.000 description 58
- 238000002156 mixing Methods 0.000 description 28
- 238000010586 diagram Methods 0.000 description 14
- 230000008859 change Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 210000001519 tissue Anatomy 0.000 description 5
- 238000003745 diagnosis Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 3
- 238000000354 decomposition reaction Methods 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 238000010008 shearing Methods 0.000 description 3
- 210000004885 white matter Anatomy 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000002591 computed tomography Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000011426 transformation method Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 208000010412 Glaucoma Diseases 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 210000000988 bone and bone Anatomy 0.000 description 1
- 210000004958 brain cell Anatomy 0.000 description 1
- 230000002860 competitive effect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000007499 fusion processing Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000000155 melt Substances 0.000 description 1
- 230000035479 physiological effects, processes and functions Effects 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 210000003625 skull Anatomy 0.000 description 1
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经验小波变换的脑部图像融合方法。
背景技术
医学成像技术可提供准确的人体组织和结构信息，给临床医学带来了极大的便利。但由于传感器成像机制的不同，不同模态的脑医学图像从不同角度表示了人脑的信息，具有一定互补性，且存在冗余信息。在脑医学图像中，以磁共振成像(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图像的残差分量子图像为固有模态函数子图像为Wi CT，i＝1,2,…,NCT，检测滤波器图像为其中NCT表示脑部CT图像的固有模态函数数量；记脑部MRI图像的残差分量子图像为固有模态函数子图像为 检测滤波器图像为其中NMRI表示脑部CT图像的固有模态函数数量；
S3：分别提取残差子图像和的L2范数特征，基于L2范数特征进行两幅残差子图像的融合，得到融合后的残差子图像
S4：采用以下方法进行固有模态函数子图像融合：
如果NCT＞NMRI，则令待融合固有模态函数子图像数量R＝NMRI，融合后固有模态函数子图像总数T＝NCT，否则令待融合固有模态函数子图像数量R＝NCT，融合后固有模态函数子图像总数T＝NMRI；分别取脑部CT图像的前R个固有模态函数子图像和脑部MRI图像的前R个固有模态函数子图像，将对应固有模态序号的子图像作为一组，得到R组待融合固有模态函数子图像，根据以下公式进行加权融合，得到融合后的固有模态函数子图像
其中，r＝1,2,…,R，表示融合后第r个固有模态函数子图像中像素点(x,y)的像素值，表示脑部CT图像的第r个固有模态函数子图像中像素点(x,y)的像素值，表示脑部MRI图像的第r个固有模态函数子图像中像素点(x,y)的像素值，和分别表示预设的和的权值，且
如果NCT＞R，则直接将脑部CT图像的后NCT-R个固有模态函数子图像作为融合后的固有模态函数子图像否则直接将脑部MRI图像的后NMRI-R个固有模态函数子图像作为融合后的固有模态函数子图像从而得到T个融合后的固有模态函数子图像Wt F，t＝1,2,…,T；
S5：采用以下方法进行检测滤波器融合：
分别取脑部CT图像的前R+1个检测滤波器图像和脑部MRI图像的前R+1个检测滤波器图像将对应序号的检测滤波器图像作为一组，得到R+1组待融合检测滤波器图像，基于最大值选择法得到融合后的检测滤波器图像
其中，分别表示融合后的检测滤波器图像检测滤波器图像检测滤波器图像中像素点(x,y)的像素值；
如果NCT＞R，则直接将脑部CT图像的后NCT-R个检测滤波器图像作为融合后的检测滤波器图像否则直接将脑部MRI图像的后NMRI-R个检测滤波器图像作为融合后的检测滤波器图像从而得到T+1个融合后的检测滤波器图像
S6：根据融合后的残差分量子图像融合后的固有模态函数子图像Wt F和融合后的检测滤波器图像进行二维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是本实施例中权值和优选设置方法的流程图；
图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图像的残差分量子图像为固有模态函数子图像为Wi CT，i＝1,2,...,NCT，检测滤波器图像为其中NCT表示脑部CT图像的固有模态函数数量；记滤波器二维记脑部MRI图像的残差分量子图像为固有模态函数子图像为检测滤波器图像为其中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范数特征进行残差子图像融合：
分别提取残差子图像和的L2范数特征，基于L2范数特征进行两幅残差子图像的融合，得到融合后的残差子图像图4是本发明中基于L2范数特征进行残差子图像融合的流程图。如图4所示，本发明中基于L2范数特征进行残差子图像融合的具体步骤包括：
S401：提取行L2范数：
对于残差子图像中的各个像素点(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范数
其中，imω()表示对应像素点在残差子图像中像素值，当像素点(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。
S402：提取列L2范数：
根据以下公式计算得到残差子图像中各个像素点(x,y)的列L2范数
类似地，当像素点(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。
S403：计算L2范数特征：
根据以下公式计算得到残差子图像中各个像素点(x,y)的L2范数特征LNω(x,y)：
S404：残差子图像融合：
采用L2范数特征可以通过子图像像素或系数振幅表示出医学图像的结构信息，以便融合不同源图像的残差分量。基于L2范数特征的残差子图像融合的具体公式如下：
其中，表示融合后残差子图像中像素点(x,y)的像素值，表示残差子图像中像素点(x,y)的像素值，表示残差子图像中像素点(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组待融合固有模态函数子图像，根据以下公式进行加权融合，得到融合后的固有模态函数子图像
其中，r＝1,2,...,R，表示融合后第r个固有模态函数子图像中像素点(x,y)的像素值，表示脑部CT图像的第r个固有模态函数子图像中像素点(x,y)的像素值，表示脑部MRI图像的第r个固有模态函数子图像中像素点(x,y)的像素值，和分别表示预设的和的权值，且
如果NCT＞R，则直接将脑部CT图像的后NCT-R个固有模态函数子图像作为融合后的固有模态函数子图像否则直接将脑部MRI图像的后NMRI-R个固有模态函数子图像作为融合后的固有模态函数子图像从而得到T个融合后的固有模态函数子图像Wt F，t＝1,2,…,T。
对于加权融合公式，权值和的取值对融合后的固有模态函数子图像具有较大影响，因此需要对权值和进行优选设置。图9是本实施例中权值和优选设置方法的流程图。如图9所示，本实施例中权值和优选设置方法的具体步骤包括：
S901：计算显著性度量：
显著性度量可以表示IMF系数的重要性，一般来说，重要系数往往比不重要系数更加显著。本发明中需要分别计算固有模态函数子图像和固有模态函数子图像的显著性。首先设置一个大小为p′×q′的全1矩阵P，p′和q′的大小根据需要确定，对于固有模态函数子图像中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的显著性度量
式中，P(m′,n′)表示全1矩阵P中坐标为(m′,n′)的元素值，即P(m′,n′)＝1。表示像素点(x-m′,y-n′)在固有模态函数子图像中对应的像素值，当(x-m′,y-n′)不存在(即超出固有模态函数子图像范围)时，对应为0。
S902：计算匹配度量：
匹配度量不但可以表示固有模态函数系数的振幅信息，同时描述了两个固有模态函数子图像的区域特征。对于固有模态函数子图像和固有模态函数子图像中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的匹配度量ψr(x,y)：
式中， 为像素点(x-m′,y-n′)在固有模态函数子图像中的像素值，为像素点(x-m′,y-n′)在固有模态函数子图像中的像素值；为固有模态函数子图像中像素点(x,y)的显著性度量；为固有模态函数子图像中像素点(x,y)的显著性度量；ε是一个预设常数，其值很小，用于避免分母为0。
S903：判断是否ψr(x,y)＞α，α表示预设的匹配度量阈值，本实施例中设置α＝0.75，如果是，进入步骤S904，否则进入步骤S905。
S904：设置权值：
令根据确定权值。
近年来，模糊集理论在图像融合中得到了广泛的应用，并取得了较好的效果。模糊集理论的隶属度函数(MF)量化了某个元素属于特定集合的隶属关系，MF的值区间为[0,1]；“0”表示完全不隶属，“1”表示完全隶属，“0”与“l”之间的值表示不同程度的隶属关系。由于固有模态函数系数分布服从高斯分布，因此本步骤中为了更加合理地确定权值，可以采用高斯隶属度函数权重，从而确定权值和其具体方法如下：
计算得到固有模态函数子图像和固有模态函数子图像中各像素点匹配度量ψr(x,y)的标准差σ和均值c，计算得到以下参数：
其中，e表示自然常数。
如果则令否则令
S905：判断是否如果是，进入步骤S906，否则进入步骤S907。
S906：令
S907：令
图10是实例一中脑部CT图像的前两幅固有模态函数子图像和脑部MRI图像的两幅固有模态函数子图像的融合结果图。图11是实例二中脑部CT图像的两幅固有模态函数子图像和脑部MRI图像的两幅固有模态函数子图像的融合结果图。
S105：检测滤波器融合：
在二维Littlewood-Paley经验小波变换中，不同医学图像得到的检测滤波器是一个二进制图像，可以看作是一个决策图。检测滤波器有助于保持医学图像的细节特征，对医学图像重建具有重要意义。由于二维Littlewood-Paley经验小波变换中检测滤波器数量与医学图像对应的子图像(残差子图像和固有模态函数子图像)数量相同，因此本发明采用了一种简单而有效的方法来融合不同源图像的检测滤波器，即最大值选择法，其具体方法如下：
与固有模态函数子图像融合类似，分别取脑部CT图像的前R+1个检测滤波器图像和脑部MRI图像的前R+1个检测滤波器图像将对应序号的检测滤波器图像作为一组，得到R+1组待融合检测滤波器图像，基于最大值选择法得到融合后的检测滤波器图像
其中，分别表示融合后的检测滤波器图像检测滤波器图像检测滤波器图像中像素点(x,y)的像素值。
如果NCT＞R，则直接将脑部CT图像的后NCT-R个检测滤波器图像作为融合后的检测滤波器图像否则直接将脑部MRI图像的后NMRI-R个检测滤波器图像作为融合后的检测滤波器图像从而得到T+1个融合后的检测滤波器图像
S106：二维Littlewood-Paley经验小波逆变换：
根据融合后的残差分量子图像融合后的固有模态函数子图像Wt F和融合后的检测滤波器图像进行二维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为采用不同客观评价指标对本实施例所得到的融合图像的指标平均值。
表1
与其他传统方法相比，本发明的STD、SF、M值最高，更具竞争力。Qabf和MI是图像融合中最重要的评价指标，本发明所得融合图像的这两个指标明显高于其他方法，说明本发明可以保留更多的医学图像信息。虽然本发明的EN、FMI和AV值与其中一些对比方法较为接近，但仍略优于其他对比方法。综上可知，本实施例表明本发明生成的融合图像优于其他方法，是一种有效的脑部图像融合方法。
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
Claims (4)
1.一种基于二维Littlewood-Paley经验小波变换的脑部图像融合方法，其特征在于，包括以下步骤：
S1：获取已经配准好的待融合的脑部CT图像和脑部MRI图像；
S2：分别对脑部CT图像和脑部MRI图像进行二维Littlewood-Paley经验小波变换，分别得到一组子图像，包括一个残差分量子图像和若干固有模态函数子图像，并将每个子图像对应的检测滤波器采用二维图像表示，记脑部CT图像的残差分量子图像为固有模态函数子图像为Wi CT，i＝1,2,…,NCT，检测滤波器图像为i′＝0,1,…,NCT，其中NCT表示脑部CT图像的固有模态函数数量；记脑部MRI图像的残差分量子图像为固有模态函数子图像为j＝1,2,…,NMRI，检测滤波器图像为j′＝0,1,…,NMRI，其中NMRI表示脑部CT图像的固有模态函数数量；
S3：分别提取残差子图像和的L2范数特征，基于L2范数特征进行两幅残差子图像的融合，得到融合后的残差子图像
S4：采用以下方法进行固有模态函数子图像融合：
如果NCT＞NMRI，则令待融合固有模态函数子图像数量R＝NMRI，融合后固有模态函数子图像总数T＝NCT，否则令待融合固有模态函数子图像数量R＝NCT，融合后固有模态函数子图像总数T＝NMRI；分别取脑部CT图像的前R个固有模态函数子图像和脑部MRI图像的前R个固有模态函数子图像，将对应固有模态序号的子图像作为一组，得到R组待融合固有模态函数子图像，根据以下公式进行加权融合，得到融合后的固有模态函数子图像
其中，r＝1,2,…,R，表示融合后第r个固有模态函数子图像中像素点(x,y)的像素值，表示脑部CT图像的第r个固有模态函数子图像中像素点(x,y)的像素值，表示脑部MRI图像的第r个固有模态函数子图像中像素点(x,y)的像素值，和分别表示和的权值，且
如果NCT＞R，则直接将脑部CT图像的后NCT-R个固有模态函数子图像作为融合后的固有模态函数子图像r′＝R+1,R+2,…,T，否则直接将脑部MRI图像的后NMRI-R个固有模态函数子图像作为融合后的固有模态函数子图像从而得到T个融合后的固有模态函数子图像Wt F，t＝1,2,…,T；
S5：采用以下方法进行检测滤波器融合：
分别取脑部CT图像的前R+1个检测滤波器图像和脑部MRI图像的前R+1个检测滤波器图像g＝0,1,…,R，将对应序号的检测滤波器图像作为一组，得到R+1组待融合检测滤波器图像，基于最大值选择法得到融合后的检测滤波器图像
其中，分别表示融合后的检测滤波器图像检测滤波器图像检测滤波器图像中像素点(x,y)的像素值；
如果NCT＞R，则直接将脑部CT图像的后NCT-R个检测滤波器图像作为融合后的检测滤波器图像g′＝R+1,R+2,…,T，否则直接将脑部MRI图像的后NMRI-R个检测滤波器图像作为融合后的检测滤波器图像从而得到T个融合后的检测滤波器图像t′＝0,1,…,T；
S6：根据融合后的残差分量子图像融合后的固有模态函数子图像Wt F和融合后的检测滤波器图像进行二维Littlewood-Paley经验小波逆变换，得到融合后的脑部图像。
2.根据权利要求1所述的脑部图像融合方法，其特征在于，所述步骤S3中融合得到残差子图像的具体步骤包括：
S3.1：对于残差子图像中的各个像素点(x,y)，ω∈{CT,MRI}，分别以该像素点(x,y)为中心提取出大小为p×q的图像块，记该图像块左上顶点的坐标为(x0,y0)，根据以下公式计算得到像素点(x,y)的行L2范数
其中，imω()表示对应像素点在残差子图像中像素值，当像素点(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：根据以下公式计算得到残差子图像中各个像素点(x,y)的列L2范数
当像素点(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：根据以下公式计算得到残差子图像中各个像素点(x,y)的L2范数特征LNω(x,y)：
S3.4：根据以下公式进行残差子图像融合：
其中，表示融合后残差子图像中像素点(x,y)的像素值，表示残差子图像中像素点(x,y)的像素值，表示残差子图像中像素点(x,y)的像素值。
3.根据权利要求1所述的脑部图像融合方法，其特征在于，所述步骤S4中权值和按照以下方法设置：
S4.1：设置一个大小为p′×q′的全1矩阵P，p′和q′的大小根据需要确定，对于固有模态函数子图像中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的显著性度量
式中，P(m′,n′)表示全1矩阵P中坐标为(m′,n′)的元素值，即P(m′,n′)＝1。 表示像素点(x-m′,y-n′)在固有模态函数子图像中对应的像素值，当(x-m′,y-n′)不存在(即超出固有模态函数子图像范围)时，对应为0；
S4.2：对于固有模态函数子图像和固有模态函数子图像中的各个像素点(x,y)，根据以下公式计算得到像素点(x,y)的匹配度量ψr(x,y)：
式中， 为像素点(x-m′,y-n′)在固有模态函数子图像中的像素值，为像素点(x-m′,y-n′)在固有模态函数子图像中的像素值；为固有模态函数子图像中像素点(x,y)的显著性度量；为固有模态函数子图像中像素点(x,y)的显著性度量；ε是一个预设常数；
S4.3：判断是否ψr(x,y)＞α，α表示预设的匹配度量阈值，如果是，进入步骤S4.4，否则进入步骤S4.5；
S4.4：令根据确定权值；
S4.5：判断是否如果是，进入步骤S4.6，否则进入步骤S4.7；
S4.6：令
S4.7：令
4.根据权利要求3所述的脑部图像融合方法，其特征在于，所述步骤S4.4中权值和按照以下方法设置：
计算得到固有模态函数子图像和固有模态函数子图像中各像素点匹配度量ψr(x,y)的标准差σ和均值c，计算得到以下参数：
如果则令否则令
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 true CN110084770A (zh) | 2019-08-02 |
| CN110084770B 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) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118037560A (zh) * | 2024-01-16 | 2024-05-14 | 北京长木谷医疗科技股份有限公司 | 基于同态滤波的多模态医学图像融合方法、装置及设备 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102129676A (zh) * | 2010-01-19 | 2011-07-20 | 中国科学院空间科学与应用研究中心 | 一种基于二维经验模态分解的显微图像融合方法 |
| US20130044929A1 (en) * | 2011-08-19 | 2013-02-21 | Industrial Technology Research Institute | Ultrasound image registration apparatus and method thereof |
| 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 | 辽宁工程技术大学 | 一种基于纹理的掌纹掌脉融合识别方法 |
-
2019
- 2019-03-04 CN CN201910160670.2A patent/CN110084770B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102129676A (zh) * | 2010-01-19 | 2011-07-20 | 中国科学院空间科学与应用研究中心 | 一种基于二维经验模态分解的显微图像融合方法 |
| US20130044929A1 (en) * | 2011-08-19 | 2013-02-21 | Industrial Technology Research Institute | Ultrasound image registration apparatus and method thereof |
| 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)
| Title |
|---|
| KAIBIN ZHAO: "Multilinear commutators of Littlewood-Paley operator on herz-hardy spaces", 《2011 INTERNATIONAL CONFERENCE ON MULTIMEDIA TECHNOLOGY》 * |
| 张小利: "图像融合及其性能评估若干问题研究", 《中国博士学位论文全文数据库信息科技辑》 * |
| 晁永国: "基于经验模态分解的图像融合研究", 《现代商贸工业》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118037560A (zh) * | 2024-01-16 | 2024-05-14 | 北京长木谷医疗科技股份有限公司 | 基于同态滤波的多模态医学图像融合方法、装置及设备 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110084770B (zh) | 2023-03-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Xia et al. | A novel improved deep convolutional neural network model for medical image fusion | |
| Luo et al. | Image fusion with contextual statistical similarity and nonsubsampled shearlet transform | |
| Ma et al. | Infrared and visible image fusion technology and application: A review | |
| Bin et al. | Efficient image fusion with approximate sparse representation | |
| Vijayakumar et al. | Segmentation and grading of brain tumors on apparent diffusion coefficient images using self-organizing maps | |
| Tirupal et al. | Multimodal medical image fusion based on yager’s intuitionistic fuzzy sets | |
| DE102004057026A1 (de) | Erfassung zeitlicher Bildveränderungen und Darstellungsverfahren und -vorrichtung | |
| DE102006041309A1 (de) | System und Verfahren zur 3D-CAD unter Verwendung von Projektionsbildern | |
| CN103793711A (zh) | 一种基于脑部核磁共振图像的多维度纹理提取方法 | |
| Bhutto et al. | CT and MRI medical image fusion using noise-removal and contrast enhancement scheme with convolutional neural network | |
| Pei et al. | Two-scale multimodal medical image fusion based on guided filtering and sparse representation | |
| CN110570387A (zh) | 一种基于特征级Copula模型相似性的图像融合方法 | |
| Sivakumar et al. | Implementation of differential evolution algorithm to perform image fusion for identifying brain tumor | |
| Saleh et al. | A brief analysis of multimodal medical image fusion techniques | |
| Kavita et al. | Study of image fusion optimization techniques for medical applications | |
| CN105809650A (zh) | 一种基于双向迭代优化的图像融合方法 | |
| Ramana | Alzheimer disease detection and classification on magnetic resonance imaging (MRI) brain images using improved expectation maximization (IEM) and convolutional neural network (CNN) | |
| Das et al. | Multimodal classification on PET/CT image fusion for lung cancer: a comprehensive survey | |
| Zhang et al. | Image fusion using online convolutional sparse coding | |
| CN110084770A (zh) | 基于二维Littlewood-Paley经验小波变换的脑部图像融合方法 | |
| Venkatesan et al. | Integrated fusion framework using hybrid domain and deep neural network for multimodal medical images | |
| Vani et al. | Multi focus and multi modal image fusion using wavelet transform | |
| CN116664462A (zh) | 一种基于ms-dsc和i_cbam的红外和可见光图像融合方法 | |
| Aishwarya et al. | An image fusion framework using morphology and sparse representation | |
| CN108022230A (zh) | 肾脏多期ct图像融合系统 |
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 |