CN108198158A - 一种基于瑞利噪声分布的超声与磁共振影像融合方法 - Google Patents

一种基于瑞利噪声分布的超声与磁共振影像融合方法 Download PDF

Info

Publication number
CN108198158A
CN108198158A CN201711454377.4A CN201711454377A CN108198158A CN 108198158 A CN108198158 A CN 108198158A CN 201711454377 A CN201711454377 A CN 201711454377A CN 108198158 A CN108198158 A CN 108198158A
Authority
CN
China
Prior art keywords
matrix
mri
eog
magnetic resonance
negative
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
Application number
CN201711454377.4A
Other languages
English (en)
Other versions
CN108198158B (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.)
Jiangsu Normal University
Original Assignee
Jiangsu Normal University
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 Jiangsu Normal University filed Critical Jiangsu Normal University
Priority to CN201711454377.4A priority Critical patent/CN108198158B/zh
Publication of CN108198158A publication Critical patent/CN108198158A/zh
Application granted granted Critical
Publication of CN108198158B publication Critical patent/CN108198158B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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 by the use of more than one image, e.g. averaging, 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/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/10—Image acquisition modality
    • G06T2207/10132—Ultrasound image
    • 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

Abstract

本发明公开了一种基于瑞利噪声分布及非负矩阵分解的超声与磁共振影像融合方法，本方法具体步骤如下：步骤1：基于自适应加权非负矩阵分解A＝WUV，将超声图像和磁共振图像分别转化为列向量US和MRI，构建非负矩阵A，并令A＝[US,MRI]，构建非负的权重矩阵W＝[WUS,WMRI]，W的每一个元素值都被初始化为0.5，U为按照行优先转换后的列向量；步骤2：计算超声和磁共振图像对应的梯度能量矩阵并记为EOGUS和EOGMRI，对应两个梯度能量矩阵元素的较大值构成梯度能量矩阵EOGUS～MRI；步骤3：对U和V进行随机初始化；步骤4：对U和V进行迭代求解；步骤5：更新非负的权重矩阵W；步骤6：如果迭代满足结束条件，那么转到步骤7；如果迭代不满足结束条件，转到步骤4；步骤7：将向量U转换为矩阵，即为融合结果。

Description

一种基于瑞利噪声分布的超声与磁共振影像融合方法
技术领域
本发明涉及医学影像图像处理技术领域，具体是一种基于瑞利噪声分布及非负矩阵分解的超声与磁共振影像融合方法。
背景技术
超声成像由于其实时性好、成本低等优点，在临床上有着广泛的应用。和磁共振和CT成像磁共振比较，超声成像分辨率不高。磁共振是临床上另一种被广泛使用的成像模式，对软组织成像效果好，并有着较高的空间分辨率。除此之外，磁共振还能够提供临床上需要的功能信息。但是磁共振成像实时性不好。目前，没有任何一种成像模式能够满足临床上得所有成像需求。各种成像模式在临床上有着各自的有点和缺点，在很多情况下，超声和磁共振在临床上是互相补充的，为了能够更好的辅助临床医生进行诊断，有必要对超声和磁共振进行融合。
由于超声影像斑点噪声突出、对比度低，直接对超声影像与其他模态影像(包括磁共振)进行融合将导致融合后影像噪声突出、清晰度降低。目前融合算法的研究集中于以下4个方面:基于数学/统计学的融合算法，包括加权平均法、NMF方法、PCA方法、贝叶斯估计融合和EM融合方法；基于多分辨率分析的图像融合算法，包括多分辨金字塔法、小波变换法；基于彩色空间的图像融合方法，包括HSI变换法、视觉神经生理学法；智能图像融合算法，主要包括模糊理论方法、神经网络法等。上述融合算法本质都是构建和优化数学模型，对两幅或多幅参考图像的信息进行综合，使融合后的图像具有更高的清晰度和信息量。上述融合算法通过降低图像清晰区域的图像质量为代价提升图像模糊区域的图像质量，以得到目标均清晰的融合结果，导致融合图像中一些重要边缘信息的丢失。
发明内容
为了克服上述现有技术的缺点，本发明提供一种基于瑞利噪声分布及非负矩阵分解的超声与磁共振影像融合方法，提高融合图像的对比度。
本发明是以如下技术方案实现的：一种基于瑞利噪声分布的超声与磁共振影像融合方法，其具体步骤如下：
步骤1：基于自适应加权非负矩阵分解A＝WUV，首先根据行优先，将超声图像和磁共振图像分别转化为列向量US和MRI，构建非负矩阵A，并令A＝[US,MRI]，构建非负的权重矩阵W＝[WUS,WMRI]，并且非负的权重矩阵W的每一个元素值都被初始化为0.5；W和A均为k×2矩阵，这里k＝n×m；U、V分别为非负的k×r矩阵和非负的r×2矩阵，其中r＝1，即U为按照行优先转换后的列向量；
步骤2：计算超声图像和磁共振图像对应的梯度能量矩阵并记为EOGUS和EOGMRI，并对上述两个梯度能量矩阵对应的元素进行比较，对应元素的较大值构成梯度能量矩阵EOGUS～MRI；
步骤3：对U和V进行随机初始化；
步骤4：对U和V进行迭代求解；
步骤5：更新非负的权重矩阵W；
步骤6：如果迭代满足结束条件||A-WUV||2＜0.01，那么转到步骤7；如果迭代不满足结束条件，转到步骤4；
步骤7：基于行优先，将向量U转换为矩阵，即为融合结果，融合过程结束；
步骤5其更新非负的权重矩阵W具体方法如下：
步骤5.1：迭代一次后，将向量U按照行优先转化为矩阵U，并计算U对应的梯度能量矩阵，记为EOGU；
步骤5.2：将WUS和WMRI按照行优先转化为矩阵；
步骤5.3对于超声和磁共振影像的每一个元素，如果EOGU(i,j)≧EOGUS～MRI(i,j)，进入步骤5.4，如果EOGU(i,j)<EOGUS～MRI(i,j)，则对应的W自适应更新算法如下：
SUM(i,j)＝WUS(i,j)+WMRI(i,j)；
步骤5.4：将WUS和WMRI按照行优先转换分别为列向量，并令W＝[WUS,WMRI]。
本发明的有益效果：将自适应非负矩阵算法用于本发明中，获得的融合图像具有更强的对比度。
附图说明
图1是本发明融合流程图；
图2是基于不同融合方法的融合结果比较图。
具体实施方式
一种基于瑞利噪声分布的超声与磁共振影像融合方法，具体步骤如下：
步骤1：基于自适应加权非负矩阵分解A＝WUV，首先根据行优先，将超声图像和磁共振图像分别转化为列向量US和MRI，构建非负矩阵A，并令A＝[US,MRI]，构建非负的权重矩阵W＝[WUS,WMRI]，并且非负的权重矩阵W的每一个元素值都被初始化为0.5；W和A均为k×2矩阵，这里k＝n×m；U、V分别为非负的k×r矩阵和非负的r×2矩阵，其中r＝1，即U为按照行优先转换后的列向量；
步骤2：计算超声图像和磁共振图像对应的梯度能量矩阵并记为EOGUS和EOGMRI，并对上述两个梯度能量矩阵对应的元素进行比较，对应元素的较大值构成梯度能量矩阵EOGUS～MRI；
步骤3：对U和V进行随机初始化；
步骤4：对U和V进行迭代求解；
步骤5：更新非负的权重矩阵W；
步骤5.1：迭代一次后，将向量U按照行优先转化为矩阵U，并计算U对应的梯度能量矩阵，记为EOGU；
步骤5.2：将WUS和WMRI按照行优先转化为矩阵；
步骤5.3对于超声和磁共振影像的每一个元素，如果EOGU(i,j)≧EOGUS～MRI(i,j)，进入步骤5.4，如果EOGU(i,j)<EOGUS～MRI(i,j)，则对应的W自适应更新算法如下：
SUM(i,j)＝WUS(i,j)+WMRI(i,j)；
步骤5.4：将WUS和WMRI按照行优先转换分别为列向量，并令W＝[WUS,WMRI]。
步骤6：如果迭代满足结束条件||A-WUV||2＜0.01，那么转到步骤7；如果迭代不满足结束条件，转到步骤4；
步骤7：基于行优先，将向量U转换为矩阵，即为融合结果，融合过程结束；
本发明原理如下：
设A为非负的待分解的n×m矩阵，U、V分别为非负的n×r矩阵和非负的r×m矩阵，非负矩阵分解的计算模型为：A＝UV。
非负矩阵分解即为求解一个非负的n×r矩阵U和非负的r×m矩阵V，在此基础上构建代价函数L(A|U,V)。L(A|U,V)作为能够度量A和UV逼近程度的代价函数，当非负矩阵分解算法迭代到一定次数，代价函数L(A|U,V)最小，迭代结束。本发明拟基于超声影像的噪声分布，构建一种新型的代价函数，基于上述代价函数，提出了自适应加权非负矩阵分解算法并进行超声与磁共振的融合。下面给出新型代价函数及自适应加权非负矩阵分解算法(A＝WUV)的构建过程，其中W为权重矩阵。
超声影像的噪声分布符合瑞利分布，瑞利分布定义如下：
在(1)式中x≥0，a为常数且a＞0.为了简化后续计算，a被设置为常数1.对(1)式进行求偏导：
求得
当x等于时，px(x)获得最大值。由于超声的斑点噪声符合噪声分布，本发明提出了一种新的架构如下：
在式(2)中，A为一个n×m非负矩阵。A被近似的分解为两个非负的矩阵U(n×r)和V(r×m)。
设
那么(3)式对应的最大似然解，即为对如下损失函数进行求最小极值：
为了对(4)式进行化简，将预先加到矩阵A的所有元素上，(4)式可以化简为：
通过最小化(6)式构建的加权代价函数，我们能够推导出自适应加权非负矩阵分解算法。
在(6)式中W＝{Wij}＞0为非负的权重矩阵，XοY代表矩阵X和矩阵Y的Hadamard积(即矩阵X和矩阵Y点乘)，这里我们定义(6)式为加权瑞利测度。
在本发明中我们仅给出V的迭代公式的推导过程，因为我们可以用同样的方法获得U的迭代式。A,V和W中任意一列向量分别表示为a,v和w，我们引入部分代价函数F(v)：
在(7)式中：
DW＝diag(w), (8)
其中DW为基于w构建的对角方阵。对(9)式和(10)式分别进行一阶和二阶求偏导，分别得到式(11)、(12)、(13)和(14)：
式(13)和(14)中的Ui代表U的第i列。对(7)式进行迭代，可以得到v的迭代式(15)
同理能够得到u的迭代式：
(16)式中，AT,UT和WT分别为A,U和W对应的转置矩阵,aT,uT和wT分别为AT,UT和WT对应的列向量，其中的一阶偏导和二阶偏导定义见式(17)、(18)、(19)和(20)：
在式(17)、(18)、(19)和(20)中，(VT)i代表VT的第i列,DwT代表基于wT构建的对角方阵。
在图像融合过程中，随着迭代的进行，融合图像像素灰度也是动态变化的，图像各区域的灰度分布也会发生变化，因此图像各区域的权重也应该发生相应的变化。采用固定权重矩阵不能体现这种动态变化过程。针对以上问题，本文提出一种动态迭代策略，随着迭代的进行，根据最新迭代结果对加权系数进行动态调整。对于迭代结果的特定指标改善已经达到目标阈值的情况，则停止对对应加权系数的调整；反之，加权系数朝有利于特定指标改善的方向进行调整。WNMF的加权系数体现了对应元素的重要程度，通过对加权系数进行持续动态的调整，即可实现融合图像特定指标的改善。
本发明采用平均梯度(average gradient，AG)、边缘强度(edge intensity，EI)、信息熵(entropy，EN)和空间频率(space frequency，SF)作为融合图像的评价标准。AG也称为清晰度，反映了图像的细节反差与纹理变化特征，融合图像的AG值越大，表明融合图像的清晰度越高；EN表示图像所包含信息量的多少，EN值越大，表明融合图像包含的信息量越大；EI反应了边缘与邻近像素的对比强度，EI值小意味着强度变化不大而难以分辩，值越大表示边缘越清晰；SF反应了图像空间域的活跃程度，目标与背景对比明显，图像中存在剧烈的突变等情况下，SF值较大，反之较小。
为了对比本发明提出的融合方法的有效性，另外分别基于DWT(discrete wavelettransform)的融合算法、基于PCA的融合算法、基于SIDWT(shift invariance discretewavelet transform)的融合算法对两只裸鼠的超声和磁共振影像进行了融合实验，实验结果如图2所示，直观观察，本发明提出方法获得的融合图像具有更强的对比度。为了定量比较，基于四种方法分别进行了融合指标计算，计算结果见表1和表2。分析表1和表2，使用本发明提出的方法对应的AG、EI、EN和SF值更大，说明本发明或的的融合结果在定量指标方面表现更好。
表1.基于第一只裸鼠的融合结果指标比较.
表2.基于第二只裸鼠的融合结果指标比较.
图2
图2(a)～(d)分别为基于第一只裸鼠使用本发明方法、DWT、PCA以及SIDWT进行融合的结果；
图2(e)～(h)分别为基于第二只裸鼠使用本发明方法、DWT、PCA以及SIDWT进行融合的结果。

Claims (2)

1.一种基于瑞利噪声分布的超声与磁共振影像融合方法，其具体步骤如下：
步骤1：基于自适应加权非负矩阵分解A＝WUV，首先根据行优先，将超声图像和磁共振图像分别转化为列向量US和MRI，构建非负矩阵A，并令A＝[US,MRI]，构建非负的权重矩阵W＝[WUS,WMRI]，并且非负的权重矩阵W的每一个元素值都被初始化为0.5；W和A均为k×2矩阵，这里k＝n×m；U、V分别为非负的k×r矩阵和非负的r×2矩阵，其中r＝1，即U为按照行优先转换后的列向量；
步骤2：计算超声图像和磁共振图像对应的梯度能量矩阵并记为EOGUS和EOGMRI，并对上述两个梯度能量矩阵对应的元素进行比较，对应元素的较大值构成梯度能量矩阵EOGUS～MRI；
步骤3：对U和V进行随机初始化；
步骤4：对U和V进行迭代求解；
步骤5：更新非负的权重矩阵W；
步骤6：如果迭代满足结束条件||A-WUV||2＜0.01，那么转到步骤7；如果迭代不满足结束条件，转到步骤4。
步骤7：基于行优先，将向量U转换为矩阵，即为融合结果，融合过程结束；
2.根据权利要求1所述的一种基于瑞利噪声分布的超声与磁共振影像融合方法，其更新非负的权重矩阵W具体方法如下：
步骤5.1：迭代一次后，将向量U按照行优先转化为矩阵U，并计算U对应的梯度能量矩阵，记为EOGU；
步骤5.2：将WUS和WMRI按照行优先转化为矩阵；
步骤5.3对于超声和磁共振影像的每一个元素，如果EOGU(i,j)≧EOGUS～MRI(i,j)，进入步骤5.4，如果EOGU(i,j)<EOGUS～MRI(i,j)，则对应的W自适应更新算法如下：
SUM(i,j)＝WUS(i,j)+WMRI(i,j)；
步骤5.4：将WUS和WMRI按照行优先转换分别为列向量，并令W＝[WUS,WMRI]。
CN201711454377.4A 2017-12-28 2017-12-28 一种基于瑞利噪声分布的超声与磁共振影像融合方法 Active CN108198158B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711454377.4A CN108198158B (zh) 2017-12-28 2017-12-28 一种基于瑞利噪声分布的超声与磁共振影像融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711454377.4A CN108198158B (zh) 2017-12-28 2017-12-28 一种基于瑞利噪声分布的超声与磁共振影像融合方法

Publications (2)

Publication Number Publication Date
CN108198158A true CN108198158A (zh) 2018-06-22
CN108198158B CN108198158B (zh) 2021-06-01

Family

ID=62584955

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711454377.4A Active CN108198158B (zh) 2017-12-28 2017-12-28 一种基于瑞利噪声分布的超声与磁共振影像融合方法

Country Status (1)

Country Link
CN (1) CN108198158B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102499711A (zh) * 2011-09-28 2012-06-20 无锡祥生医学影像有限责任公司 三维或四维超声图像自动优化调节方法
CN102542547A (zh) * 2011-12-29 2012-07-04 北京航空航天大学 一种基于光谱约束的高光谱图像融合方法
EP2759817A1 (en) * 2013-01-25 2014-07-30 General Electric Company Ultrasonic holography imaging system and method
CN104091355A (zh) * 2014-06-06 2014-10-08 华南理工大学 一种采样角度受限下的ct图像重建方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102499711A (zh) * 2011-09-28 2012-06-20 无锡祥生医学影像有限责任公司 三维或四维超声图像自动优化调节方法
CN102542547A (zh) * 2011-12-29 2012-07-04 北京航空航天大学 一种基于光谱约束的高光谱图像融合方法
EP2759817A1 (en) * 2013-01-25 2014-07-30 General Electric Company Ultrasonic holography imaging system and method
CN104091355A (zh) * 2014-06-06 2014-10-08 华南理工大学 一种采样角度受限下的ct图像重建方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
侯漠 等: "基于双模态造影剂的超声影像与磁共振影像的配准", 《中国科学》 *
侯漠 等: "基于磁性微泡对比剂的超声和磁共振影像融合算法", 《中国科学》 *

Also Published As

Publication number Publication date
CN108198158B (zh) 2021-06-01

Similar Documents

Publication Publication Date Title
Yin et al. A novel infrared and visible image fusion algorithm based on shift-invariant dual-tree complex shearlet transform and sparse representation
Chai et al. Image fusion scheme using a novel dual-channel PCNN in lifting stationary wavelet domain
Yin et al. Simultaneous image fusion and super-resolution using sparse representation
CN101882304B (zh) 一种sar图像自适应去噪和特征增强方法
Chai et al. Multifocus image fusion scheme based on features of multiscale products and PCNN in lifting stationary wavelet domain
CN101216939B (zh) 基于量子行为粒子群算法的多分辨率医学图像配准方法
CN101847257B (zh) 基于非局部均值与多级定向图像的图像降噪方法
CN111563915B (zh) 一种融合运动信息检测和Radon变换的KCF目标跟踪方法
CN105631820A (zh) 基于小波变换和三边滤波器的医学超声图像去噪方法
NL2025814B1 (en) Precise positioning system for three-dimensional position of tumor by multi-image fusion
CN104021536B (zh) 一种自适应的sar图像和多光谱图像融合方法
CN103236046A (zh) 基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法
CN102567973A (zh) 基于改进的形状自适应窗口的图像去噪方法
CN110322404A (zh) 一种图像增强方法及系统
CN114792287A (zh) 一种基于多图像融合的医学超声图像超分辨率重建方法
Liu et al. An effective wavelet-based scheme for multi-focus image fusion
CN103077507B (zh) 基于Beta算法的多尺度SAR图像降噪方法
CN103198456B (zh) 基于方向波域隐马尔可夫树模型的遥感图像融合方法
CN102096913A (zh) 压缩感知框架下的多策略图像融合方法
CN103839225A (zh) 基于模糊密度权的支持向量场景图像去噪算法
CN104408695B (zh) 直方图模糊加权调整与均衡的红外图像增强方法
CN108198158A (zh) 一种基于瑞利噪声分布的超声与磁共振影像融合方法
CN106530277A (zh) 一种基于小波方向相关系数的图像融合方法
Bhattacharya et al. Brain image segmentation technique using Gabor filter parameter
Chaitanya et al. PET and MRI medical image fusion using STDCT and STSVD

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