2008-功能磁共振图像融合方法研究

东南大学
硕士学位论文
功能磁共振图像融合方法研究
姓名：王欢
申请学位级别：硕士
专业：生物医学工程
指导教师：鲍旭东
20080115
东南大学硕士学位论文
而如今，它作为一种常用的医学成像方式和检查手段，已经成为影像医学中不可或缺的一项重要内容。

ＭＲＩ是利用原子核在磁场内共振所产生的信号从而重建成像的一种技术。

磁共振是自旋的原子核在磁场中与电磁波相互作用的一种物理现象。

原予核由质子和中子组成，质子带有正电荷，自旋将产生一个小磁场，称为磁矩，而中子虽然是中性的，但由于内部电荷不均匀分布，自旋时也会产生磁矩（大约为质子磁矩的２／３）。

很多原子核中，质子和中子是成对存在的，磁矩能彼此抵消，对外磁矩很弱。

但对于一些原子核，质子和中子总数为奇数时。

就存在自旋磁矩，于是就存在磁共振现象。

氢原子核１Ｈ只有一个自旋的质子，结构最单纯，又能够提供最强的磁共振信号，目前磁共振成像主要是利用人体内的氢原子核．
磁共振成像利用投影重建图像。

成像过程中先把检查层面分成一
定数量的小体积（即体素），用接收器收集信息。

数字化后输入计算机
进行重建运算，获得每个体素的信号强度ｐＪ。

磁共振成像是一种安全、快速、准确的ｌ临床诊断方法，通过参数
调配可以针对具体病症、部位选择最佳成像方式。

它的成像质量，尤
其是针对人体软组织的成像效果远胜于其他技术。

图１．１是一张典型的解剖磁共振图像。

图ｌ－Ｉ解剖磁共振图像１．２．２功能磁共振图像（ｆＭＲｌ）
功能磁共振成像（ｆｕｎｃｔｉｏｎａｌｍａｇｎｅｔｉｃｒｅｓｏｎａｎｃｅｉｍａｇｉｎｇ，ｆＭＲＤ是近十余年来发展的一种成像技术，反映人脑在执行某项任务或受到某种刺激时皮质功能区的激活情况，以其无创性在脑的Ｉ临床及基础研究中得以广泛应用。

基于ＢＯＬＤ（ｂｌｏｏｄｏｘｙｇｅｎａｔｉｏｎｌｅｖｅｌｄｅｐｅｎｄｅｎ０技术的功能磁共振成像，由Ｏｇａｗａ等于１９９０年首先报道［６１。

其基本原理是当大脑执行一项任务或受到某种刺激时，某些脑区神经元的活动增强，其邻近静脉血管和毛细血管床的血流量和血流容增加，导致局部氧合血红蛋白含量增加，而这个增加额远大于局部耗氧量的增加，使得氧供应和氧消耗之间失衡，即该区域脱氧血红蛋白含量相对降低。

脱氧血红蛋白作为顺磁性物质，有明显的亿・缩短效应。

因此，在激活状态下该脑区的Ｔ２・的相对延长，ＭＲ信号强度增加，在脑功能成像时功能活动区的皮质表现为高信号，利用ＥＰＩ快速成像序列就可以将其检测出来１７１。

在ＢＯＬＤ法ｆＭＲＩ成像中，由于功能信号比较弱，所以需要在功能活跃和功能静止时反复多次扫描，然后将功能活跃时的扫描图像序列和功能静止时的扫描图像序列与功能激励曲线进行一定的相关分析，就可以得到功能活跃造成的信号变化，即功能信号。

至于功能信号参数的具体重建过程，将在下一章，功能磁共振参数图像重建中进行详细的介绍。

ＥＰＩ成像技术的时间分辨率较高，但问题在于空间分辨率较低。

因此，它需要与空间分辨率高的解剖磁共振图像结合起来，才具有临床诊断和研究意义。

２
——一．墨二！堕兰
ＥＰＩ圈像重建后鲍功能参数圉健
’ｌ：：三ｋ．：：：！！：：≥＞绪果显示圈像
解剖硅共振图像
闺１－２日ｘＡ豫Ｌ姻图像
１．３功能磁共振图像显示现状及存在的闯题
目前对于功能磁共振图像显示方面的研究还不是很多，大多数课题只将目光放在了功能区域的定位上面，也就是说・它们的结果图着重于表现功能区在什么位置，因此，通常的作法是将得到的功能图像直接叠加在解剖图像之上Ｅ蚺１”，结果如图１－３所示；
（丑＞英文字形实验结果ｌ司晰觉功能散发结果柳
圈ｌｏ目前常用的功能硅共振田像显示方法
这样的显示方法存在两个坷题。

第一，如图１－２中重建后的功能参数图像所示，每一个像素上都会有一个功能参数的系数。

功能区域与非功能区域的差别即在于功能区域中像素的功能参数值较大，而非功能区域中像素的功能参数值较小一且前常用的显示方法是确定一个闲值，大于这个闲值的像素即判定它属于功能区域，
被赋上色彩或显示为高亮，小于这个阈值的像素被认为属于非功稚区。

维持其原灰度值不变。

因此．
３
东南大学硕士学位论文
第二章功能磁共振参数图像重建
２．１实验图像
本次实验源ｆＭＲＩ数据来Ｉ！ｌｈｔｔｐ：／／ｗｗｗ．ｆｉｌ．ｉｏｎ．ｕｃｌ．∞．ｕｌ（，ｓｐＩｌ似ａｔａ，它是利用ＳｉｅｍｅｎｓＭＡＧＮＥＴＯＭ
Ｖｉｓｉｏｎ２Ｔ系统进行听觉功能激发实验获取的实验数据。

听觉实验包含静息和激发两种功能状态，静
息期不施加任何刺激，激发期以每分钟６０个双音节英文单词的速度施加听觉刺激。

这样的实验共重
复８次，每次实验中，对静息期和激发期分别进行６次采样，每次采样得到６４幅连续的ＢＯＬＤ／ＥＰＩ
横轴位断层图像，扫描矩阵为６４ｘ６４ｘ６４，体素大小为３ｍｍｘ３ｍｍｘ３ｍｍ，单次扫描采样重复时间
ＴＲ＝７ｓ，整个实验持续６分钟。

其相应的解剖图像为ｓＭ００２２３００２．ｉｍｇ，体素大小为ｌｍｍｘｌｍｍｘ
３ｍｍ０２１。

窗宽窗位。

可见此ＥＰＩ图像噪声较大，如图２－１ｂ所示。

噪声是影响医学图像质量的重要因素之一，影响有效信
息的提取，进而影响图像分析的精度与结果。

在进行功能参数图像重建之前，我们要对这些
ＢＯＬＤ／ＥＰｌ图像进行滤波去噪的预处理。

（ａ）ＥＰＩ断层图像（ｂ）断层图像中的噪声
图２．１ＥＰＩ断层图像
２．２数据预处理
去除噪声，提高信噪比是图像处理中的基本任务之一，通常采用低通滤波器滤除噪声。

滤波的
主要工作分为两步，一是确定卷积的核，二是利用确定的核对源图像进行卷积，得到结果图像。

高
斯滤波是目前应用最广泛的滤波方法之一。

本课题中也将采用高斯滤波对图像去噪。

２．２．１高斯卷积核参数选取
相距中心点距离为ｊ的单元，其高斯函数被定义为ｌ”１：
ｇ（Ｊ）；ｅ一』２“２／２ｘ／磊Ｔｍ２（２．１）
当进行三维的高斯卷积时，Ｊ２ｚｘ２＋ｙ２＋ｚ２．（２．１）式中，ｓ２为方差，ｓ为标准差。

ｓ的计算方法
６
第二章功能碰共振图像重建
如式（２．２）所示。

ｓ：！些
ｘ／８Ｉｎ（２）
（２．２）ｏ＇ⅢＨＭ是高斯核的半最大值所占宽度（ＦｕｎｗｉｄｔｈａｔＨａｌｆＭａｘｉｍｕｍ）Ｉ“１。

一维高斯函数如图２－２所示。

盯Ｍ值越大，意味着高斯核的宽度越宽，图像处理起来也就越慢，
去噪效果较好，但也使处理后的图像越模糊。

在肭Ｉｕ数据处理中，如何选择最优的平滑宽度仍存在争议。

根据标准的滤波理论，最优的滤波核宽度应与检测的激活区域的尺寸相匹配。

当图像信噪比较低时，可采用较宽的滤波核以较好地过滤噪声。

由于本实验中．ＥＰＩ体数据的分辨率为３ｍｍＸ３ｍｍＸ３ｒａｍ，可取其２～３个体素点的宽度为高斯核半最大值所占宽度，因此这里取ｏⅢⅢ＝８ｒａｍ
根据式（２．２），计算可得：ｓ＝ｉｚ丽８＝ｌ・１３５
２．２．２ＥＰＩ数据预处理结果／＼图２－２一维高斯函数
利用上述的方法对ＢＯＬＤ／ＥＰＩ图像进行处理，
处理结果如图２．３所示。

２—３ａ是未经高斯滤波前的源图像，可见图像
噪声较为严重，信噪比低。

２．３ｂ和２－３ｃ是经过高
斯卷积后的结果图像，其高斯函数中的标准差Ｓ
分别取ｌ（盯删＝７ｍｍ）和１．１３５（盯删＝Ｓｍｍ）。

由图可见，随Ｓ增大，噪声逐渐减小，但同时。

图
像也愈发模糊，图像细节信息丧失越多・
（ａ）滤波前的ＥＰＩ（ｂ）滤波后的ＥＰＩ（ｃ）滤波后的ＥＰＩ由此可见，在进行卷积滤波的时候应该选取图像图像睁１）图像（ｓ－Ｉ，１３５）
合适的盯日”ｎ。

值，即高斯核半最大值所占宽度。

图２－３８９１数据高斯滤波前后结果对比
７
东南大学硕士学位论文
表２．２改进前后方法得到的相关系数矩阵统计结果比较
０．４７Ｏ．∞ｏ．５５Ｏ．６Ｏ．６２修正前
相关系
ｌ∞６５３３１４ｌＯ数大于方法
ｎ的像修正后
素个数３０７２３５１５９１２７１１０方法
通过修正前后方法的统计学对比，可见改进后的方法令结果相关系数有了较大幅度的提升。

通过二次相关分析，得到了各像素被激活情况的参数映射数据，它是一个体数据，反映了各像素的被激活程度。

图２－６是其在横断面第３ｌ层上的截图。

（ａ）一次相关分折法所得重建结果（ｂ）改进后相关分析法所得重建结果
图２－６功能参数重建结果图
图２－６里高亮的像素表示了被激活程度高的区域。

也就是大脑在特定激励之下的功能区。

由图可见，重建得到的功能参数图像不包括解剖信息，因此，它需要与解剖图像联合使用，才能够知道究竟是大脑的什么位置在实验中被激活了。

１０
第三章功能磁共振闰像与解剖磁共振图像配准
图３＝２是功能参数体数据、ＥＰＩ数据和ＭＩｕ数据某三个层面上的截图。

ｆａ）功能参数甜－数捌截同
ｂ）ＥＰＩ融旧
（ｃ）ＭＲＩ截图
圈３－２配准源数据截面图
图３－３是配准前的ＥＰ！与Ｍｌｕ，及功能参数体数据和ＭＲＪ数据在某一截面上的横断面，冠状面、矢状面图。

１７
东南大学硕士学位论文
圈３－３配准前图像位置关系
经过几何变换、重采样、相似性测度和优化搜索的循环步骤后，图３－４ａ显示了ＥＰＩ和ＭＰｄ数据的配准结果，将配准矩阵应用于功能参数体数据，对其进行重采样后，功能参数图与ＭＲＩ对齐结果如图３－４ｂ所示。

图３－４配准后图像位置关系
目前在医学图像配准的问题上，还没有所谓的“金标准”，因此通常采用目标检测方法对结果进行评估。

由上图可见，原先在空闻位置上有很大差异的解剖磁共振图像与功能磁共振图像，在经由最大互信息配准处理后，达到了不错的配准效果。

最大互信息配准有以下优点：
第五章改进的区域能量融合算法投评价
其中，
ｒ‰寺三ｃ竿》
｛
Ｌ阡，二。

；１一野Ｚ。

－（占＝１，２’３）（５．５）这种方法相较于传统的功能磁共振图像与解割图像融合方法（直接叠加法），应用Ｊ＇泛而原理简单的加权平均融合方法，以及同样也是基于小波分析的小波域内最大系数法．是有明显优势的。

实验中，我们选取源解削磁共振体数据与功能参数体数据横断面上第２７层的图像（横断面上共５４层数据），利用直接彝加法、加权平均法闭、基于小波变换的最大系数法和基Ｔ－ｄ，渡变换的局部区域能量法分别进行融合实验，得到结果如图５－１所示。

（ａ）鱼描吾加，占站果ｆｂ｝加权平均珐结果
∞基十小政变挽的最大系数征结皋（ｄ）基十小敬业挟的局＿右ｌ；区域能量法站粜
脚５－！功能参教幽与ＭＲＩ融合结果幽
图５．１左侧的两张图像分别是解剖磁共振和功能参数体数据第２７层图像。

功能参数图中．像素的明瞎程度表示了大脑相应位置被激活程度的高低，也就是大脑响应功能信息的表现。

从此图也可以看出，重建得到的功能图像不包括解剖信息，因此，它需要与解剖图像联合使用。

才能够知道究竟是大脑的什么位簧在实验，ｌｔ被激活了。

５一ｌａ使用的是功能图像与解剖图像融合的传统方法，也就是直接叠加法。

由５－１ａ可见。

这种方ｊ去虽然能准确地反映出大脑中被撒活区域的位置，但由于解剐图像中原有的信息被覆盖，因而丧失了被激活区域详细的解剖信息，图像无法反映出此区域的大脑沟回情况，而沟回信息对丁大脑功能医学实验米说是重要的。

此外，５－ｌａ的方法需要一个阚值，用以判断功能参数图像中哪些区域被撤活了，而５．１ｂ、ｃ，ｄ的方法皆可使用功能参数图像的全部信息与解剖图像直接进行融合。

无需这样
３５
东南大学颈上学位论文
图５—２原局部区域能量法融合结果
豳５－３改进融合权重后局椰区域能量法融合结粜
由图５－２和５．３可见，与以往的显示或融合方法相比较，小波融合对于解剖图像中的沟回信息，以及功能图像中的强弱信息，都有了比较蟊好的表现。

实验结果可见，在颞上回等沟回处的激活响应明显高于在其余的灰质和白质部分的激活响应，且高亮（激活更明显）区域陷于沟回内部，这也很好地符合了听觉响应实验中的实际情况。

而改进前后的融合方法相比，改进后的方法对于细节显示更明显。

高亮的红色激活区域陷于沟回内部，且较连续，而改进前的方法所得的结果，高亮区域显示不是非常连续．因此从主观方面评价，基于小渡变换的局部区域能量法要优于以往的显示及融合方法，而对融合权重加以改进了的融合算法，又饔优于原算法。

东南大学硕士学位论文
（ａ）ＬＬ频带匹配度直方图（ｂ）ＬＨ频带匹配度直方图
（ｃ）ＨＬ频带匹配度直方圈（ｄ）ＨＨ频带匹配度直方图
图Ｓ４各频带匹配度直方图
本文选择Ｏｓｔｕ提出的类问方差法对匹配度集合膨；。

进行分析，自适应地选取最佳阈值。

由Ｏ吱ｕ提出的类间方差法【５ａ１是一种受到广泛关注的阈值选取方法。

该方法又称为大律闲值分割法，是在判决分析最小＝乘法原理的基础上推导得出的，算法较为简单，实时性好。

最大类间方差法的基本思路是：选取的最佳阈值Ｔ应当使得不同类间分离住最好。

首先对整幅图像的匹配度加以统计，计算出每一匹配度对于此组图像而言发生的概率，称为匹配度直方图，然后基于匹配度直方图得到备分割特性值的发生概率，并以阈值变量Ｔ将分割特征值分为两类。

再求出每一类的类内方差及类间方差，选取使类问方差最大，类内方差最小的Ｔ作为晟佳阌值。

即：
盯；＝口口２＋盯；，７（，）：嘻仃：Ｔ＝■曙ｍ钗叩（，）妊ｔｓＬ（５９）（５．１０）
（５．１１）
其中，口∥２为类内方差，矿口２为类问方差，仃ｒ２为总体方差，Ｌ为匹配度的范围，理想情况下为１。

方差是描述随机变量与其均值之间偏离程度的统计量。

假设匹配度统计数组为Ｍ（ｉ），即匹配度为ｉ的像素概率为Ｍｊ，而阈值ｔ将匹配度集合分为Ａ、Ｂ两类。

匹配度小于ｔ的为Ａ类，匹配度大于等于ｔ
４０
第五章改进的区域能量融合算法及评价
小波分解—＿Ｉ◆
ＬＬＨＬＬＨＨＨ
图５－５源图像小波分解示意图
现在所常用的基于局部区域能量的小波融合算法，是利用５．１式的方法计算窗口内的局域能量。

无论是ＬＬ、ＬＨ、ＨＬ还是ＨＨ子频带。

窗口权重都是相同的∞（ｍＩ，Ｖ／Ｉ）＝：
ｏ１
０１１４１Ｉ。

这相当于忽ｏ１ｏＩ
略了小波子频带的方向特性。

以ＬＨ子频带为例，该频带保持了图像竖直方向上的高频边缘信息，即在理想情况下，如果时域图像内有一条水平方向的边缘，那么，该边缘在小波域的ＬＨ频带内的相应位置也会表现出一条水平方向数值变化较平缓而均值相对周围其它区域的值来说较大的高频特征信息。

由于我们融合的目标是在融合结果图像上凸显两幅源图像的细节信息，也就是边缘信息，因此作为与窗Ｅ１中心点呈水平方向的像素理应比其它方位的像素占有更大的能量权重。

依此思路，将计算小波域内局部区域能量的窗口权矩阵调整如下：
ｒｏ１０１ｒｏ２０１
ＬＬ频一扩ｃ∥∥，‘都列蝴一矿ｃ耐∥，２弛４２：ｌ
ＬＨ频一∥ｔ∥，们。

轭。

４。

２０００ｌＨＨ频带一脚３（ｍｕ，／，／ｉ）＝；１０
１１ｏ４ｏＪ（＆１３）
１０１Ｊ
针对不同的频带，利用不同的窗口权矩阵根据５．１和５．２式计算局部区域能量和匹配度，利用这样改进后的方法进行图像融合，对融合结果的客观评价如表５４所示：。