基于SPM的功能磁共振成像图像头动校正配准方法的分析与改进

陈 文 CHEN Wen 郁 芸 YU Yun 刘宏毅 LIU Hongyi 胡新华 HU Xinhua 杨 坤 YANG Kun 赵 宇 ZHAO Yu

基于SPM的功能磁共振成像图像头动校正配准方法的分析与改进

陈 文 CHEN Wen 郁 芸 YU Yun 刘宏毅 LIU Hongyi 胡新华 HU Xinhua 杨 坤 YANG Kun 赵 宇 ZHAO Yu

作者单位 南京医科大学 江苏 南京 210029

功能磁共振成像作为一种新兴的神经影像学成像技术，其重要性受到神经、认知和心理学等领域的广泛关注。本文针对功能磁共振成像数据预处理过程中的头动校正环节，较为全面地阐述了国际知名脑功能成像软件包SPM的配准原理，并在此基础上提出了一种改进的配准方法，通过由粗到细的处理有效地提高了配准的准确度。此外，将差值平方和与归一化互信息的比值作为配准测度的思想弥补了单用一种测度的局限性，增加了配准过程的鲁棒性和可靠性。研究结果显示该方法具有一定的实际意义和应用价值。

磁共振成像；图像配准；图像处理，计算机辅助；SPM；头动校正

功能磁共振成像（functional magnetic resonance imaging，fMRI）是一种新的无侵入、无创伤的脑功能成像技术，为研究脑的认知活动、生理功能及病理状态提供了有利的工具，已成为脑科学和生命科学研究的重要工具。

在fMRI技术中，影像图像的后期处理具有重要意义，图像的配准（registration）是其中的一个关键环节。在扫描采样过程中，受试者不由自主的轻微头动仍然不可避免，包括物理头部运动（轻微的左右摆动或点头等）和生理头部运动（心跳、呼吸等），导致伪影出现［1］。为了提高脑功能图像统计分析的准确性，需对图像进行配准处理，以确保不同时间序列图像空间位置的完全对应。

目前用于fMRI图像处理的软件包有AFNI、SPM、FSL、MRIcro等，其中SPM是一个专为脑功能成像数据分析而设计的通用软件包，是国际上应用最普遍、最权威的脑功能成像处理分析软件。本文基于SPM分析了fMRI图像头动校正的配准方法，并对其进行优化和改进。

1 SPM图像头动校正配准的算法原理与结果分析

1.1 概述 为了适应图像数据分析的需要，fMRI数据必须经过多种图像预处理过程，包括头动校正、图像融合、图像分割、标准化、空间平滑等［1-3］。SPM中的头动校正功能可以确定时间序列图像与其参考图像配准的参数，并在此基础上将一个实验序列中的每一帧图像均与这个序列的参考图像进行对齐，以矫正头动［4］。本文讨论的配准对应于SPM的Realign（estimate）过程。

1.2 配准方法 作为一种理想的单体（subject）单模态（modality）配准，头部运动的校正常采用刚体运动模型。SPM配准算法的基本思想是通过迭代计算平移、旋转参数，使参考图像（通常为时间序列的第1幅）与后续序列图像的不匹配度最小化，实现所有序列的配准［1］。

1.2.1 图像配准原理 以2幅图像为例，头动校正是在相似性测度下寻找一个空间变换关系，使得经过该变换后2幅图像的相似性达到最大。例如，对参考图像g和资源图像f进行配准，则需要寻找一个变换Mr，使得：

其中，S是相似性测度，Mr为所求的最佳空间变换［5］。

因为需保证各向同性，需将体素坐标转换为笛卡尔坐标。笛卡尔坐标系下的参考图像ge和资源图像fe为：

其中，Mg、Mf分别代表将参考图像g和资源图像f从体素坐标变换到笛卡尔坐标的转换矩阵。

由公式（1）及公式（2）可得，变换过程为：

1.2.2 刚体变换（rigid transformation） 对于同一受试者，由于未产生结构间相对位置的改变，头部运动即可分解为旋转和平移两部分的刚体变换。定义刚体变换向量=（q1，q2，q3，q4，q5，q6），其中q1，q2，q3依次代表头部沿x，y，z 坐标轴的平移量大小；q4，q5，q6依次代表头部绕x，y，z坐标轴的旋转角度大小，则：

1.2.3 Gauss-Newton迭代算法 建立刚体变换模型之后，头动校正的配准即可化归为刚体变换向量= （q1，q2，q3，q4，q5，q6）的确定，在SPM中求解这6个独立未知参数使用Gauss-Newton迭代算法来实现。

考虑将参考图像g的坐标xi经过刚体变换逆变换为资源图像的坐标要使得两者的相似性达到最大，SPM采用差值平方和作为测度，应使其最小化，定义目标函数为：

其中，xi表示第i个体素的坐标，为变换矩阵，f（）表示资源图像的灰度值，g（）表示参考图像的灰度值，u为2幅图像的灰度平衡系数。

当差值平方和函数趋于平缓或达到指定的迭代次数Nth时迭代结束，此时的即为最佳的刚体变换参数。

1.3 实验结果及分析 实验所用为100张已做完时间校准的img图像（a*.img），配准共用时27.2 s，所得结果为Graphics窗口中的2幅图像，清晰地显示了受试者在检查时的平动和转动情况（图1）。

图1 配准结果显示的平动（A）和转动（B）情况

经分析发现，SPM的配准过程仍存在一些不足，有待进一步优化：①SPM迭代可能出现过早收敛至局部最优解的情况，具有一定的局限性［8］。伍亚军等［9］提出了“安全区域”的概念，即配准目标函数最值的附近存在一个特殊区域，在此区域内函数较为光滑，局部极值要缓和的多甚至消失，如果在迭代之前进行一次粗配准，使迭代的初始值进入安全区域，即可避开局部极值的影响。②SPM采用的配准测度是差值平方和最小化，然而这仅是一种相似性测度。目前基于像素灰度值的配准测度有很多，最常用的是图像的差值平方和与互信息［10］。适当调整配准测度，将互信息最大化纳入以弥补单用一种测度的局限性，将使配准的准度得到一定的提高。

2 SPM图像头动校正配准的算法改进与结果分析

2.1 配准方法的改进 从配准的精度和速度出发，本文在SPM的基础上做了以下改进：为了避免收敛至局部极值，采用力矩主轴法在迭代之前进行一次粗配准；为了减少粗配准增加的时间，先提取图像轮廓，省去不必要的计算负担，相对加快了配准速度；为了减弱测度不全面性对配准精确度的影响，在原先配准测度的基础上加入互信息测度，增加配准过程的鲁棒性和可靠性。优化后的算法流程见图2。

图2 优化后的算法流程

2.1.1 轮廓提取 轮廓提取主要涉及形态学在边缘提取、区域填充和移除对象方面的应用。首先对Canny算子提取的边缘进行连接运算，然后进行区域填充、细化、腐蚀，用细化后图像减去腐蚀过的图像，移除干扰轮廓，即可得到去除内部边缘线的外轮廓图［11］。

2.1.2 主轴矩法 力矩主轴法借用经典力学中刚体质量分布的概念，计算2幅图像的质心位置和主轴方向，进一步得到质心和主轴对齐的平移、旋转参数，从而限制位置偏差在一个较小范围内，达到粗略配准的目的。

首先，在笛卡尔坐标体系下，通过经典力学方法计算质心：

其中，x，y，z代表像素点的坐标，p（x，y，z）代表像素点的灰度值（视为像素的质量）。

其次，根据“主轴是惯量矩阵的特征向量”，计算惯量矩阵：

其中，Ixx=μ020+μ002，Iyy=μ200+μ002，Izz=μ200+μ020，Ixy=Iyx=μ110，Ixz=Izx=μ101，Iyz=Izy=μ011。μpqr=∑x，y，z（x－xc）p（y－yc）q（z－zc）r，表示p+q+r阶中心矩。

由于惯量矩阵I的标准化特征向量E等于旋转矩阵R。

其中，α，β，γ分别为图像绕x，y，z轴旋转的角度。

最后求得：

由公式（8）、公式（11）可分别算出2幅图像的xc，yc，zc，α，β，γ，即质心和主轴方向得到确定。经过作差，可粗略地得到平移和旋转的配准参数［11-13］。

2.1.3 SPM迭代 完成粗配准后，将得到的配准参数作为SPM迭代的初始值，以达到避免误差函数过早收敛的目的，同时可一定程度地减少迭代次数和缩短配准时间。

2.1.4 寻优 经过SPM迭代，在差值平方和最小化的原则下获取了较为准确的配准参数，但这只是相似性测度的一种表现形式。为了提高准确度，弥补单用一种配准测度的局限性，考虑将互信息（mutual information）作为测度的一部分。

互信息具有精度高、适应性强、自动化程度高的优点，在配准多模图像和具有不同成像参数的同模图像领域均有较好的效果［10，14］。在信息统计学中，互信息通常用于描述两个系统之间的统计相关性，把2幅图像的灰度值看作随机变量A、B，其概率密度分布函数分别为pA（a）和pB（b），两者的联合概率密度分布函数为pAB（a，b），则随机变量A和B的互信息可表示为：

其中，H（A），H（B）和H（A，B）为随机变量A和B的个体熵和联合熵，其定义为：

既往研究表明，当两组基于共同解剖结构的图像的空间位置达到一致时，由于图像之间存在的相关性最大，灰度联合概率密度分布p（a，b）最集中，联合熵H（A，B）最小，互信息量I（A，B）最大。互信息测度对于2幅图像的重叠覆盖程度较为敏感，Studholme等提出了归一化互信息（normalized mutual information）的概念［5］，并将其作为一种更为准确的配准测度，其定义为：

综合两种测度，本文选择差值平方和与归一化互信息的比值作为相关性测度：

当S值取得最小值时，则近似地认为相关程度达到最大。

该步骤在SPM所得配准参数的基础上运用交替方向法进行进一步寻优。

2.2 实验结果分析 本文针对前面所述的100幅图像进行头动校正处理。根据实验结果，100幅图像中76%的图像数据提示该方法较改进前具有较好的改善，其中17%显著改善，其余24%的数据表示算法无改进。表1中为几组典型的数据。

由表1可见，经过配准后第9、26、43、65、87幅图像的SSD减小、NMI增大、SSD/NMI减小，其配准的准确度得到较为全面的提高；第55幅和第73幅图像尽管SSD有所增加，但NMI增大、SSD/NMI减小，即互信息水平得到提高，且互信息对配准效果的影响较大；第35幅图像NMI有一定的增大，但意义不如前两者明显。

表1 改进前和改进后的相关典型数据

以第43幅图像的第15层为例，改进前后的图像见图3。

图3 改进前后配准结果比较

本文算法中，轮廓提取有效地省去了不必要的计算负担，相对加快了粗配准的速度；基于主轴矩法的粗配准可限制位置偏差在一个较小范围内，避免迭代时收敛至局部极值，同时减少了迭代次数，一定程度地提高了运算速度；互信息测度的加入及其交替方向法判优，弥补了单用一种测度的局限性，因而增加了配准过程的鲁棒性和可靠性，但由于计算量较大，又延长了配准时间。

因此，此算法对SPM配准算法的改进有一定的意义，对fMRI图像预处理中的配准操作具有一定的实用价值。然而，在配准时间上，此算法仍存在缺陷，有待进一步优化。

3 结束语

本文较为全面地阐述了SPM头动校正配准的算法原理，客观地分析了SPM处理的优、缺点，并在此基础上提出了一种改进的配准方法，加入了轮廓提取、主轴矩法和基于SSD/NMI新测度的交替方向法判优的步骤，由粗到精地进行处理，避免了局部极值问题，有效地提高了配准过程的准确度。后续研究将针对算法运算时间较长的问题进一步深入分析，探讨并提出优于SSD/NMI的配准测度。

［1］ 张志强，王世杰，卢光明.功能磁共振数据处理分析的原理及应用.中国医学影像技术，2004，20（10）：1632-1635.

［2］ 赵健仪，尹建忠，纪东旭，等.嗅觉脑功能MRI对阿尔茨海默病的早期诊断.中国医学影像学杂志，2013，21（5）：321-325.

［3］ 曾统军，赖新生，黄泳，等.针刺外关穴和非穴fMRI脑功能成像的比较研究.中国医学影像学杂志，2010，18（6）：532-537.

［4］ 刘永，张莉，马竟锋，等.基于fMRI的立体视觉数据分析辅助诊断系统.中国医疗器械信息，2015，（7）：66-70.

［5］ 王丽，孙丰荣，王奕琨，等.基于互信息的颅脑MR影像序列的三维配准.计算机工程与应用，2011，47（31）：160-163.

［6］ 吕晓琪，李娜，张宝华，等.基于体素相似性的三维多模脑图像配准研究.中国医学影像学杂志，2013，21（2）：146-151.

［7］ 李京娜，王国宏，孙少燕，等.基于修正的结构相似度为测度的三维脑图像配准.中国医学影像学杂志，2013，21（8）：618-623.

［8］ 王超，杜剑平，闫镔，等.SPM软件包头动校正原理的分析与评价.中国组织工程研究与临床康复，2010，14（52）：9772-9776.

［9］ 伍亚军，周正东，戴耀东.基于特征点互信息预配准的医学图像配准技术.南京航空航天大学学报，2007，39（5）：622-627.

［10］ 孙少燕，李斌，唐一源.fMRI时间序列图像的配准测度及其相关优化方法.应用基础与工程科学学报，2004，12（3）：241-249.

［11］ 钱慧，胡志忠.基于粗配准和混合互信息的医学图像配准算法.生物医学工程研究，2007，3（3）：266-270.

［12］ 吴锋，钱宗才，杭洽时，等.基于轮廓的力矩主轴法在医学图像配准中的应用.第四军医大学学报，2001，22（6）：567-569.

［13］ 胡久乡，丁国治.主轴法在三维图像匹配中的应用.华中科技大学学报（自然科学版），2002，30（1）：38-40.

［14］ 杜俊俐，郭清宇.互信息医学图像配准的实时性研究.计算机工程与应用，2011，47（28）：162-165.

（本文编辑 张春辉）

Analysis and Improvement of Functional MRI Image Registration Method Based on SPM

As a new neuroimaging method，functional magnetic resonance （fMRI） with its importance of data processing has been widely recognized by neurology and cognitive psychology.Focusing on the realignment section in fMRI image preprocessing，this paper comprehensively describes the registration principle of SPM，an internationally-known software package which is specially designed for cerebral function imaging.An improved registration method is presented which effectively increases the accuracy.In addition，choosing the ratio of SSD and NMI as the registration measure can compensate for the limitation of using single measurement，which improves the robustness and reliability of the registration process.Experimental results prove the feasibility of this method.

Magnetic resonance imaging； Image registration； Image processing，computer-assisted； SPM； Realignment

10.3969/j.issn.1005-5185.2015.10.017

郁 芸

Nanjing Medical University，Nanjing 210029，China

Address Correspondence to：YU Yun E-mail：yuyun@njmu.edu.cn

江苏省大学生实践创新训练计划项目（201410312057X）；南京市医学科技发展项目（YKK12125）；南京市医学科技发展项目（201303009）。

R445.2

2015-03-17

2015-04-20

中国医学影像学杂志2015年 第23卷10期：780-784

Chinese Journal of Medical Imaging 2015 Volume 23（10）：780-784