基于泊松方程和互信息分解的多模态医学图像融合方法

1.本发明涉及一种医学图像融合方法，特别涉及一种基于泊松方程和互信息分解的多模态医学图像融合方法，该方法主要将功能代谢图和解剖结构图融合起来，提供互补的医学信息，属于医学图像处理领域。


背景技术：

2.现代医学影像技术为医学研究提供了多种模态的医学图像。根据各模态图像所提供的信息内涵，可以将这些图像分为解剖结构图(如ct、mri)和功能代谢图(如spect、pet)。由于成像机制的差异，不同模态的图像往往呈现不同的特色。例如，ct图像中骨骼、植入物等致密结构成像清晰，但是软组织的分辨率较低；而mri图像则可以高分辨率地反映血管等软组织或器官的解剖结构；pet和spect图像可以提供功能代谢信息，但是缺乏必要的解剖结构。为了综合利用不同模态所提供的医学信息，通常将多个模态的医学图像融合在一起进行分析。例如，将解剖结构图和功能代谢图融合起来，利用两者信息的互补性，提高病灶定位与评估的准确性。
3.目前，针对一般图像的融合问题，已经提出了很多方法。按照融合时数据的表征层次，这些方法可以分为像素级融合、特征级融合和决策级融合。在多模态医学图像融合领域，应用得比较多的是像素级融合，因为像素级融合直接对原始图像进行处理，可以有效保留源图像的细节，准确性较高。按照融合时数据的表现形态，像素级融合又分为空间域融合和变换域融合。空间域融合直接在图像的像素值空间进行融合。变换域融合利用人眼视觉系统多分辨率的特点，通常先对待融合图像进行多尺度变换，然后在变换域进行融合。与空域融合相比，变换域融合得到的图像质量更高，细节更丰富；但是变换域融合在保留不同模态细节特征方面，存在特征边界模糊的现象。
4.在图像编辑领域，为了实现图像的无缝克隆(seamless cloning)，文献1(patrick p
é
rez,michel gangnet,and andrew blake.poisson image edting.acm transaction on graphics,vol.22,no.3,2003.)提出了基于泊松方程的图像融合方法。该方法为了把感兴趣的源图像区域融入目标图像，首先根据源图像和目标图像的梯度建立一个引导场(guidance field)，然后在引导场的约束下构建泊松方程，通过求解泊松方程得到融合后的图像。泊松图像融合效果取决于引导场的构造。为同时保留源图像和目标图像的特征，通常选用模值更大的梯度或者是梯度的线性组合来构造引导场。理论上，这种构造方法可以直接用于不同模态的医学图像融合。然而，跟一般图像融合不同，不同模态的医学图像并不是完全独立的，对其中一种模态的观察往往会获取一些有关另一种模态的信息。传统的泊松图像融合并没有考虑不同模态之间的相关性，存在丢失某些重要信息的可能。
5.不同模态之间的相关性，可以借助互信息(mutual information，mi)来刻画。理论上，互信息是指一个随机变量中包含的关于另一个随机变量的信息量。对于图像融合而言，可以将两个图像的像素取值看作两个不同的随机变量，那么两个图像的互信息就是两个图像之间共享的信息量。对互信息进行分解，可以进一步得到与单一像素值关联的特定信息
(specific information)，如单一像素值的意外性(surprise，i1)、预见性(predictability，i2)和卷入度(entanglement，i3)等。文献2(r.bramon,i.boada,et al.mutlimodal data fusion based on mutual information.ieee transactions on visualization and computer graphics,2012,18(9)：1574
‑
1587.)给出了分解互信息mi而得到i1、i2和i3的具体方法。给定随机变量x和y，以及各自的观测值x和y，从互信息i(x；y)中分解出的i2(x；y)刻画了观测值x的预测能力，x的i2值越高，y的不确定性就越小；i3(x；y)则在i2的基础上进一步展示了观测值x卷入观测值y的程度，i3值越大，y中与x有关的y的信息量越大。
6.另外，在医学图像处理领域，功能代谢图一般是包含rgb三个分量的彩色图像，而解剖结构图则是只包含强度信息的灰度图像。为了避免颜色失真，在融合功能态的彩色图像与解剖态的灰度图像时，通常将功能态的彩色图像变换到ihs(intensity
‑
hue
‑
saturation)颜色空间，提取其强度分量与解剖态灰度图像进行融合；然后再对融合结果进行ihs逆变换，得到保留功能态彩色信息以及解剖态结构信息的融合图像。
7.鉴于上述情况，本发明针对功能代谢图与解剖结构图的融合问题，提出一种基于泊松方程和互信息分解的多模态医学图像融合方法。在给定功能代谢图和解剖结构图的情况下，该方法首先对功能代谢图进行ihs变换得到其强度分量x；然后计算x和解剖结构图y的互信息i(x；y)，以及互信息的分解结果i2(x；y)和i3(x；y)；再根据互信息分解结果构造泊松方程，求解泊松方程得到融合了功能态和解剖态信息的强度分量f；最后对强度分量f进行ihs逆变换，得到包含功能态彩色信息以及解剖态结构信息的融合图像。该方法利用功能态图像与解剖态图像的互信息分解结果构造泊松方程，通过求解泊松方程而得到融合结果。这不仅使融合结果包含了尽可能多的信息，还提高了融合结果中不同模态特征的清晰度，特别是解剖态结构特征的细节清晰度。


技术实现要素：

8.本发明的目的是提供一种基于泊松方程和互信息分解的多模态医学图像融合方法，该方法是将功能代谢图与解剖结构图融合起来，提供尽可能多的不同模态的信息，同时提高融合结果中不同模态特征的细节清晰度，特别是解剖态结构特征的细节清晰度。
9.本发明的目的是通过以下技术方案实现的。
10.基于泊松方程和互信息分解的多模态医学图像融合方法，包括以下步骤：
11.步骤1：输入两个已配准的功能代谢图和解剖结构图，分别令其为图像a和图像b。
12.步骤2：对处于rgb模式的彩色功能代谢图像a，令其rgb分量为(r
a
,g
a
,b
a
)，然后按照下面的公式计算a的强度分量i
a
以及相关参数v1和v2：
[0013][0014]
步骤3：令图像a的强度分量i
a
为x，令处于灰度图模式的图像b为y，计算x和y的互信息i(x；y)的分解结果i2(x；y)和i3(x；y)，具体计算过程如下：
[0015]
步骤3.1：令图像x中像素取值的集合为s
x
，x(i,j)表示x中第(i,j)个像素的取值，图像y中像素取值的集合为s
y
，y(i,j)表示y中第(i,j)个像素的取值；然后按照下面的公式，统计图像x和y的像素取值概率p)x)和p(y)，以及图像x和y之间像素取值的条件转移概率p(x|y)和p(y|x)：
[0016][0017][0018][0019][0020]
其中n(x)为图像x中像素取值为x的像素个数，n(y)为图像y中像素取值为少的像素个数；n
x
为图像x的像素总数，n
y
为图像y的像素总数；n(x，y)是将图像x的像素x(i，j)和y的像素y(i，j)一一配对后，像素对中取值为(x，y)的像素对数量。
[0021]
步骤3.2：针对每一个x∈s
x
和y∈s
y
，计算x和y的互信息分解结果i2(x；y)和i3(x；y)，计算公式如下：
[0022][0023][0024][0025]
步骤4：根据x和y互信息分解结果，构造泊松方程，求解泊松方程，得到x和y的融合结果f，具体计算过程如下：
[0026]
步骤4.1：按照下面的公式构造x和y的特定信息映射图m
i2
和m
i3
：
[0027]
m
i2
(i，j)＝i2(x(i，j)；y)，
[0028]
m
i3
(i，j)＝i3(x(i，j)；y)。
[0029]
步骤4.2：创建与x大小相同的掩码m
c
，在m
c
中将与x的待融合区域对应的位置设置为1，其他位置都设为0。
[0030]
步骤4.3：令图像x和y的融合结果为f，按照下面的公式对f进行初始化：
[0031][0032]
步骤4.4：计算图像x和y的梯度场g
x
和g
y
，然后按照下面的公式构造融合结果f的梯度场g
f
：
[0033][0034]
步骤4.5：根据f的梯度场g
f
，计算该梯度场的散度场，令其结果为δf，即
[0035][0036]
其中和分别表示g
f
(i，j)的x分量和y分量。
[0037]
步骤4.6：令m
c
中取值为1的元素编号分别为(m
k
，n
k
)，k＝1，2，...k，根据前面的计算结果δf，构造如下以f(m
k
，n
k
)为未知数的泊松方程组：
[0038][0039]
其中f(m
k
±
1，n
k
±
1)是f(m
k
，n
k
)相邻位置的像素值，它们的值可能是已知的，也可能是待求解的；如果它们所对应的掩码m
c
(m
k
±
1，n
k
±
1)＝1则是待求解的像素值，否则其值为步骤4.3的初始化结果y(m
k
±
1，n
k
±
1)。
[0040]
步骤4.7：求解步骤4.6构造的泊松方程组，得到融合结果f。
[0041]
步骤5：将融合结果f与步骤2得到的参数v1和v2结合起来，按照下面的公式计算得到新图像c的rgb分量(r
c
，g
c
，b
c
)；新图像c就是功能代谢图像a和解剖结构图像b的最终融合结果：
[0042][0043]
有益效果
[0044]
本发明所述的基于泊松方程和互信息分解的多模态医学图像融合方法，主要针对功能代谢图与解剖结构图的融合，在实现技术和融合效果方面具有以下几个方面的优势和特点：
[0045]
(1)本发明方法不针对源图像的像素值直接进行融合，而是基于源图像的梯度场构造泊松方程，通过求解泊松方程得到融合结果。这种融合方式可以更好地保留细节特征，提高特征边界的清晰度。
[0046]
(2)本发明方法考虑了不同模态之间的相关性，利用不同模态的互信息分解结果指导融合过程，可以使融合结果包含尽可能多的信息。
[0047]
(3)本发明方法基于his颜色模型，将灰度解剖结构图与彩色功能代谢图的强度分量进行融合，避免了融合结果中功能代谢内容颜色失真的问题。
[0048]
(4)与变换域的多尺度融合方法相比，本发明方法的计算量更小，效率更高，因为本发明方法不涉及复杂的多尺度变换，也不需要多层级融合处理。
[0049]
(5)本发明方法的融合效果更好，融合图像不仅保留了更多解剖态和结构态信息，
而且提高了模态特征(特别是解剖态细节特征)的清晰度。
附图说明
[0050]
图1基于泊松方程和互信息分解的多模态医学图像融合方法的流程图；
[0051]
图2待融合的pet图像a；
[0052]
图3待融合的mri t2图像b；
[0053]
图4按照本发明方法融合图像a和b得到的结果图像c；
[0054]
图5按照传统泊松融合方法融合图像a和b得到的结果；
[0055]
图6按照多尺度变换方法融合图像a和b得到的结果；
[0056]
图7待融合图像以及不同融合方法得到的结果图像的局部放大对比。
[0057]
具体实施方法
[0058]
下面结合附图和实施例阐述本发明的具体实施方式。
实施例
[0059]
如图1所示，本发明所述的基于泊松方程和互信息分解的多模态医学图像融合方法，包括如下步骤：
[0060]
步骤1：输入两个已配准的功能代谢图和解剖结构图，分别令其为图像a和图像b。
[0061]
本实施例从harvard medical school公开的the wholebrain atlas数据集中(http：//www.med.harvard.edu/aanlib/home.htm)，选取了一个阿尔兹海默症患者的脑部pet图像和mri t2图像作为输入，分别令其为a和b。这两个图像的大小都是256
×
256，其中rgb模式的彩色pet图像提供功能代谢信息，如图2所示；灰度模式的mri t2图像提供高分辨率的解剖结构信息，如图3所示。
[0062]
步骤2：对处于rgb模式的彩色功能代谢图像a，令其rgb分量为(r
a
,g
a
,b
a
)，然后按照下面的公式计算图像a的强度分量i
a
以及相关参数v1和v2：
[0063][0064]
上面从(r
a
,g
a
,b
a
)到(i
a
,v1,v2)的变换，实际上是从rgb颜色空间变换到ihs颜色空间的中间结果，可以在此基础上根据参数v1和v2进一步计算得到与分量i
a
对应的h
a
分量和s
a
分量。本发明以及本实施例中不涉及h
a
分量和s
a
分量，因此没有进一步计算h
a
分量和s
a
分量。之所以采用ihs颜色模型提取图像的强度分量是因为ihs颜色模型更符合人眼视觉特性。
[0065]
步骤3：令图像a的强度分量i
a
为x，令处于灰度图模式的图像b为y，计算x和y的互信息i(x；y)的分解结果i2(x；y)和i3(x；y)，具体计算过程如下：
[0066]
步骤3.1：令图像x中像素取值的集合为s
x
，x(i，j)表示x中第(i，j)个像素的取值，图像y中像素取值的集合为s
y
，y(i，j)表示y中第(i，j)个像素的取值；然后按照下面的公式，统计图像x和y的像素取值概率p(x)和p(y)，以及图像x和y之间像素取值的条件转移概
率p(x|y)和p(y|x)：
[0067][0068][0069][0070][0071]
其中n(x)为图像x中像素取值为x的像素个数，n(y)为图像y中像素取值为y的像素个数；n
x
为图像x的像素总数，n
y
为图像y的像素总数；n(x，y)是将图像x的像素x(i，j)和y的像素y(i，j)一一配对后，像素对中取值为(x，y)的像素对数量。
[0072]
本实施例中n
x
＝n
y
＝256
×
256，且图像x和y中强度值的取值范围是0到255，即0≤x，y≤255。
[0073]
步骤3.2：针对每一个x∈s
x
和y∈s
y
，计算x和y的互信息分解结果i2(x；y)和i3(x；y)，计算公式如下：
[0074][0075][0076][0077]
其中i2(x；y)实际上是x和y的互信息i(x；y)的一部分，因为
[0078][0079]
因此，i2(x；y)和i3(x；y)是分别采取不同的分解策略，将x和y的互信息分解到不同像素取值上的结果。i2(x；y)刻画了x的预测能力，i2(x；y)值越高，y的不确定性就越小；i3(x；y)则在i2(x；y)的基础上展示了x和y卷入度，i3(x；y)值越大，y中与x有关的y的信息量越
大。
[0080]
由于不同位置的像素可能取相同的像素值，所以本发明先针对每一个x∈s
x
和y∈s
y
，计算出i2(x；y)，i2(y；x)和i3(x；y)，以免在后续步骤中重复计算。
[0081]
步骤4：根据x和y互信息分解结果，构造泊松方程，求解泊松方程，得到x和y的融合结果f，具体计算过程如下：
[0082]
步骤4.1：按照下面的公式构造x和y的特定信息映射图m
i2
和m
i3
：
[0083]
m
i2
(i，j)＝i2(x(i，j)；y)，
[0084]
m
i3
(i，j)＝i3(x(i，j)；y)。
[0085]
在本实施例中，由于输入图像的大小是256
×
256，因此i和j的取值范围0≤i，j≤255，从而特定信息映射图m
i2
和m
i3
的大小也是256
×
256。
[0086]
步骤4.2：创建与x大小相同的掩码m
c
，在m
c
中将与x的待融合区域对应的位置设置为1，其他位置都设为0。
[0087]
在本实施例中，根据pet图像的特点，采用阈值法确定x的待融合区域，将x中强度值大于阈值的非边界区域设置为待融合区域。令阈值为t，那么掩码m
c
的设置结果是
[0088][0089]
在本实施例中，阈值t＝0.28
×
(x
max
‑
x
min
)+x
min
，其中x
max
和x
min
分别是x的最大和最小强度值。
[0090]
步骤4.3：令图像x和y的融合结果为f，按照下面的公式对f进行初始化
[0091][0092]
步骤4.4：计算图像x和y的梯度场g
x
和g
y
，然后按照下面的公式构造融合结果f的梯度场g
f
：
[0093][0094]
步骤4.5：根据f的梯度场g
f
，计算该梯度场的散度场，令其结果为δf，即
[0095][0096]
其中和分别表示g
f
(i，j)的x分量和y分量。
[0097]
步骤4.6：令m
c
中取值为1的元素编号分别为(m
k
，n
k
)，k＝1，2，...k，根据前面的计算结果δf，构造如下以f(m
k
，n
k
)为未知数的泊松方程组：
[0098][0099]
其中f(m
k
±
1，n
k
±
1)是f(m
k
，n
k
)相邻位置的像素值，它们的值可能是已知的，也可能是待求解的；如果它们所对应的掩码m
c
(m
k
±
1，n
k
±
1)＝1则是待求解的像素值，否则其值
为步骤4.3的初始化结果y(m
k
±
1，n
k
±
1)。
[0100]
步骤4.7：求解步骤4.6构造的泊松方程组，得到融合结果f。
[0101]
步骤5：将融合结果f与步骤2得到的参数v1和v2结合起来，按照下面的公式计算得到新图像c的rgb分量(r
c
，g
c
，b
c
)；新图像c就是功能代谢图像a和解剖结构图像b的最终融合结果：
[0102][0103]
本实施例在步骤5得到的最终融合结果图像c如图4所示。与图2和图3所示的融合前的图像(图2为功能态的pet图，图3为解剖态的mri t2图)对照，可以看出图4包含图2用彩色表示的功能代谢信息，同时也呈现了图3高分辨率的脑沟回结构信息，彼此的位置关系也一目了然，实现了功能态与解剖态信息的互补融合。
[0104]
为了说明本发明方法的优势与特色，分别采用传统的泊松图像融合方法(参见patrick p
é
rez,michel gangnet,and andrew blake.poisson image edting.acm transaction on graphics,vol.22,no.3,2003.)，以及基于非下采样剪切波变换的融合方法(参见m.yin,x.liu,et al.medical image fusion with parameter
‑
adaptive pulse coupled neural network in nonsubsampled shearlet transform domain.ieee transactions on instrumentation and measurement,vol.68,n0.1,2018.)，对本实施例中待融合的图像a和b进行了融合处理，得到的融合结果如图5和图6所示。
[0105]
图5虽然也是按照构造泊松方程求解泊松方程的方式得到的融合结果，但是图5对功能代谢图pet的色彩保留程度不如图4。与图2对照，图5的颜色失真比较大。
[0106]
在图6所采取的基于非下采样剪切波变换的融合方法中，为避免颜色失真，先对功能代谢图进行yuv变换得到y分量，再对y分量和解剖结构图进行非下采样剪切波变换，然后分别用能量加权和自适应脉冲耦合神经网络融合低频和高频系数，最后对融合系数进行非下采样剪切波逆变换和yuv逆变换得到彩色融合图像。就融合过程而言，图6对应的融合方法需要经过多次数据变换和处理，比本发明提出的方法更复杂。就融合效果而言，图4和图6对功能代谢信息的表现都比较好，色调和饱和度都非常接近融合前的图像。但是，图6在解剖结构特征的表现上存在一些不足。对比图3和图6可以看到，解剖结构的一些细节特征在图6中模糊了或者是消失了。而图4则比较好地保留很多微小的细节特征，解剖结构的边界更清晰，不同模态之间对比度更明显。从图7的局部放大对比结果中，可以更清楚地观察到这一点。
[0107]
此外，采用标准互信息(normalized mutual information，nmi)评价指标(参见hossny m,nahavandi s,creighton d.comments on information measure for performance of image fusion.electronics letters,2008,44(18)：1066
‑
1067)对图4、图5和图6的融合效果进行评价，得到的nmi值依次是0.7782、0.7146和0.6485。从这个评价指标看，图4包含的信息最多，融合效果最好。
[0108]
上述步骤及实例说明了本发明所述的基于泊松方程和互信息分解的多模态医学
图像融合方法的全部过程。
[0109]
应该理解的是，本实施方式只是本发明实施的具体实例，不应该是本发明保护范围的限制。在不脱离本发明的精神与范围的情况下，对上述内容进行等效的修改或者变更均应包含在本发明所要求保护的范围之内。