本申请一般涉及磁共振图像处理技术领域，尤其涉及自适应小动物脑功能磁共振成像数据分析方法。
背景技术：
小动物实验在神经科学研究中具有不可替代的作用，例如新药物研发，重大脑疾病发病机理研究等，都需要利用适当的动物模型进行试验。在体小动物脑功能磁共振成像(functionalmagneticresonanceimaging，fmri)技术使研究人员能够动态观测小动物脑功能活动，为小动物脑功能活动及脑功能网络研究提供了强有力的手段。目前常用基于手工勾画感兴趣区(regionofinterest，roi)的分析方法，即由研究人员逐一且逐层的勾画出各个roi，工作量繁重，效率低，且对操作者的背景知识要求较高，受人为主观因素的影响较大，且不能满足高通量的小动物fmri脑功能成像数据的分析。因此，小动物图像的数据分析是限制其发展和应用的主要瓶颈。
为了实现客观且高效的小动物fmri全脑功能网络分析，首先将所有被试动物的fmri脑功能像标准化到统一空间，即通常所说的空间标准化，其中标准脑模板为其提供了统一的参考标准。与人脑功能成像不同，小动物脑远远小于人脑，因此，为了提高其成像的空间分辨率及降低信噪比，小动物脑功能成像的时间通常要远远长于人脑成像。因此，为了提高实验效率，降低实验成本，研究人员通常针对特定的研究任务设定不同的成像视野(fieldofview,fov)，例如关注恐惧记忆情感环路功能的研究人员就仅对大鼠的大脑部分进行成像；关注抑郁症疾病动物模型脑功能网络异常情况的研究人员就会对小动物的全脑进行成像，包括嗅球，大脑和小脑等。
目前在人脑功能成像研究中使用最为广泛的单一fov的标准脑模板就难以满足小动物脑功能成像中不同研究的需要。
技术实现要素：
鉴于现有技术中的上述缺陷或不足，期望提供一种适应性广、不受人为主观因素影响的自适应小动物脑功能磁共振成像数据分析方法。
本申请提供一种自适应小动物脑功能磁共振成像数据分析方法，包括以下步骤：
s100：多次获取每只被试动物的全脑扫描图像，并将每只所述被试动物的每次所述全脑扫描图像进行时间片校正，以消除数据采集周期内各扫描层的时间偏差；
s200：以每只所述被试动物成像时间序列中的第一个全脑扫描图像为头动校正参考图像，分别以每只所述被试动物的头动校正参考图像为参考，对该只所述被试动物的剩余各全脑扫描图像分别进行刚体变换以获得头动校正图像，对每只被试动物的所述头动校正参考图像及所述头动校正图像进行平均，获得每只被试动物的平均图像；
s300：对所述平均图像迭代的进行图像配准及图像平均，获得用户数据驱动的自适应脑模板；
s400：以所述自适应脑模板为参考标准，将所述每只被试动物的平均图像与所述参考标准进行图像配准，并将相同的配准参数应用于所述头动校正图像，以完成每只所述被试动物的空间标准化；
s500：使用多通道高斯滤波的方式提取出所述自适应脑模板的脑内像素作为初始的脑内像素掩膜图像；
s600：以标准脑图谱所提供的脑轮廓信息为先验知识，通过图像形态学算法及统一化分割的方法对初始的脑内像素掩膜图像进行边缘修正，获得最终的与自适应脑模板空间统一的颅内像素掩膜图像；
s700：以所述颅内像素掩膜图像，去除自适应脑模板中的脑外像素，获得仅包含颅内信息的标准自适应脑模板；
s800：将标准自适应脑模板与标准脑图谱进行配准，根据标准脑图谱中的感兴趣区的信息，获取所有被试动物的所有感兴趣区的时间序列；
s900：通过皮尔逊相关算法，计算被试动物的所有感兴趣区两两之间的相关系数，获得所有被试动物的全脑网络。
进一步地，上述步骤s300具体为：
s301：从所有所述平均图像中随机选取一副图像，对其进行空间直角三坐标系旋转，使该幅图像的冠状面与标准脑图谱的切片方位一致，将进行空间直角三坐标系旋转后的该幅图像作为初始参考图像；
s302：以初始参考图像与平均图像的图像强度间的残差小于5％为目标函数，对所述平均图像进行空间变换，以完成平均图像与初始标准图像的空间对齐，所述目标函数表示为：
其中，imgt表示初始参考图像和imgi表示平均图像，下标ijk表示该像素点在图像矩阵中的坐标位置，n表示图像中包含的所有的像素个数；
s303：将空间对齐后的所有平均图像进行空间平均，获得参考平均图像；
s304：判断参考平均图像与初始参考图像之间的图像强度间残值是否小于5％，若是则将参考平均图像作为由用户数据驱动的自适应脑模板，否则将参考平均图像作为初始参考图像，重新执行步骤s302-s304。
进一步地，所述对所述平均图像进行空间变换，具体为：首先通过12参数仿射变换对所述初始参考图像与所述平均图像进行一次对齐操作，然后通过非线性形变对已对齐操作后的所述初始参考图像与所述平均图像进行再次对齐。
进一步地，上述步骤s500，具体为：
s501：建立所有元素值为零的三维积分图像，所述三维积分图像的矩阵大小与所述自适应脑模板的矩阵大小一致；
s502：通过二维高斯滤波器迭代对所述自适应脑模板进行逐层的高斯滤波，所述二维高斯滤波器的均值为零，方差均匀递增；每次迭代后将三维积分图像中与自适应脑模板高斯滤波后剩余像素点对应的像素点位置的图像强度值增加1，直至达到预定迭代次数后，将图像强度值增加后的所述三维积分图像作为所述掩膜图像。
进一步地，上述步骤s600，具体为：
s601：采用图像形态学算法去除所述掩膜图像中的独立区域；
s602：按照掩膜图像的分辨率对标准脑图谱进行三维降采集，获取与掩膜图像矩阵、空间分辨率均相当的脑图谱图像；
s603：根据所述脑图谱图像所提供的图像强度及位置信息，建立所述脑图谱图像的概率密度图，以此概率密度图为先验知识，通过统一化分割的方法将所述掩膜图像进行迭代的配准及分割，获得所述掩膜图像各层的图像强度阈值；
s604：将所述掩膜图像中图像强度值大于所述图像强度阈值的像素点赋值为1，否则赋值为0，获得二值化的颅内像素掩膜图像。
本申请提供的自适应小动物脑功能磁共振成像数据分析方法，采用基于用户数据本身建立脑模板后再完成所有被试的空间对齐。该方法自适应性强，精度高，几乎适用于所有不同成像参数的小动物fmri脑功能成像的空间标准化。本方案根据小动物脑功能成像特点，利用多通道的高斯滤波及形态学先验信息，建立了全自动小动物脑功能成像颅骨剥离算法，能够完成对自适应脑模板的颅骨剥离进而实现其与标准脑图谱的空间对齐，较人工勾画法具有工作强度低，效率高，且不受人为主观因素的影响，其剥离结果具有客观性强，可重复性高的优点。该技术方案一定程度上为小动物fmri脑功能数据的批量快速分析，减少数据分析的主观性及对全脑网络的深度挖掘和分析提供了强有力的工具。
具体实施方式
下面结合实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释相关发明，而非对该发明的限定。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合实施例来详细说明本申请。
本发明实施例提供的自适应小动物脑功能磁共振成像数据分析方法，包括以下步骤：
s100：多次获取每只被试动物的全脑扫描图像，并将每只所述被试动物的每次所述全脑扫描图像进行时间片校正，以消除数据采集周期内各扫描层的时间偏差；
采用磁共振成像仪分别对每只被试的小动物顺次进行多次扫描，获得每只被试小动物的多个全脑扫描图像，对每次扫描结果均进行时间片校正，以消除数据采集周期内各扫描层的时间偏差，使个扫描层采用同一时间参考点。
s200：以每只所述被试动物成像时间序列中的第一个全脑扫描图像为头动校正参考图像，分别以每只所述被试动物的头动校正参考图像为参考，对该只所述被试动物的剩余各全脑扫描图像分别进行刚体变换以获得头动校正图像，对每只被试动物的所述头动校正参考图像及所述头动校正图像进行平均，获得每只被试动物的平均图像；
在进行头动校正时，以头动校正参考图像为基准，对每只被试小动物的其余图像在空间三坐标系内进行平动、旋转，以消除扫描过程中由呼吸、心跳等因素引起的微小的头动。头动校正后，对每只被试动物的所有图像的各像素强度进行相加在求平均值，得到每只被试动物的平均图像。例如，每只被试小动物对应40张全脑扫描图像，则需要将每只小动物40张全脑扫描图像对应像素处的像素强度相加后再除以40，即得到每只被试动物的平均图像。
s300：对所述平均图像迭代的进行图像配准及图像平均，获得用户数据驱动的自适应脑模板；
s400：以所述自适应脑模板为参考标准，将每只被试动物的平均图像与所述参考标准进行图像配准，并将相同的配准参数应用于所述头动校正图像，以完成每只所述被试动物的空间标准化。
s500：使用多通道高斯滤波的方式提取出所述自适应脑模板的脑内像素作为初始的脑内像素掩膜图像；
s600：以标准脑图谱所提供的脑轮廓信息为先验知识，通过图像形态学算法及统一化分割的方法对初始的脑内像素掩膜图像进行边缘修正，获得最终的与自适应脑模板空间统一的颅内像素掩膜图像；
s700：以所述颅内像素掩膜图像，去除自适应脑模板中的脑外像素，获得仅包含颅内信息的标准自适应脑模板；
s800：将标准自适应脑模板与标准脑图谱进行配准，根据标准脑图谱中的感兴趣区的信息，获取所有被试动物的所有感兴趣区的时间序列；
s900：通过皮尔逊相关算法，计算被试动物的所有感兴趣区两两之间的相关系数，获得所有被试动物的全脑网络。
计算获得的相关系数，可以保存为excel文件，以方便用户进行后续的统计学分析。
上述方案的自适应小动物脑功能磁共振成像数据分析方法，采用基于用户数据本身建立脑模板后再完成所有被试的空间对齐。该方法自适应性强，精度高，几乎适用于所有不同成像参数的小动物fmri脑功能成像的空间标准化。此外，由于是完全由用户数据驱动生成标准脑模板，对数据扫描参数和成像fov没有特定要求，因此是适应非常广，几乎可以使用与所有不同成像参数的小动物fmri脑功能成像的空间标准化。本方案根据小动物脑功能成像特点，利用多通道的高斯滤波及形态学先验信息，建立了全自动小动物脑功能成像颅骨剥离算法，能够完成对自适应脑模板的颅骨剥离进而实现其与标准脑图谱的空间对齐，较人工勾画法具有工作强度低，效率高，且不受人为主观因素的影响，其剥离结果具有客观性强，可重复性高的优点。一般计算一只大鼠全脑网络的运行时间可以控制在10分钟之内，大大提高了数据分析的效率。该技术方案一定程度上为小动物fmri脑功能数据的批量快速分析，减少数据分析的主观性及对全脑网络的深度挖掘和分析提供了强有力的工具。
此外，上述方案数据处理便捷，全脑网络的相关系数可以直接保存为excel文件，供研究人员后续使用。
进一步地，上述步骤s300具体为：
s301：从所有所述平均图像中随机选取一副图像，对其进行空间直角三坐标系旋转，使该幅图像的冠状面与标准脑图谱的切片方位一致，将进行空间直角三坐标系旋转后的该幅图像作为初始参考图像；
例如，该处所指的冠状面可以是从嗅球向小脑方向所作的切片。
s302：以初始参考图像与平均图像的图像强度间的残差小于5％为目标函数，对所述平均图像进行空间变换，以完成平均图像与初始标准图像的空间对齐，所述目标函数表示为：
其中，imgt表示初始参考图像和imgi表示平均图像，下标ijk表示该像素点在图像矩阵中的坐标位置，n表示图像中包含的所有的像素个数；
s303：将空间对齐后的所有平均图像进行空间平均，获得参考平均图像；
s304：判断参考平均图像与初始参考图像之间的图像强度间残值是否小于5％，若是则将参考平均图像作为由用户数据驱动的初始自适应脑模板，否则将参考平均图像作为初始参考图像，重新执行步骤s302-s304。
进一步地，所述对所述平均图像进行空间变换，具体为：首先通过12参数仿射变换对所述初始参考图像与所述平均图像进行一次对齐操作，该次对齐实现初始参考图像与平均图像的整体对齐，然后通过非线性形变对已对齐操作后的所述初始参考图像与所述平均图像进行再次对齐，该次对齐实现初始参考图像与平均图像的细节对齐。
进一步地，上述步骤s500，具体为：
s501：建立所有元素值为零的三维积分图像，所述三维积分图像的矩阵大小与所述自适应脑模板的矩阵大小一致；
s502：通过二维高斯滤波器迭代对所述自适应脑模板进行逐层的高斯滤波，所述二维高斯滤波器的均值为零，方差均匀递增；每次迭代后将三维积分图像中与自适应脑模板高斯滤波后剩余像素点对应的像素点位置的图像强度值增加1，直至达到预定迭代次数后，将图像强度值增加后的所述三维积分图像作为所述掩膜图像。
例如但不限于，所建立的二维高斯滤波器方差均匀递增，即|σ1|2＜|σ2|2＜…＜|σk|2，且|σ2|2-σ1|2＝|σ3|2-σ2|2＝…＝|σk|2-σk-1|2，迭代次数为5，通过二维高斯滤波器迭代对初始自适应脑模板进行逐层的高斯滤波，每次迭代后将三维积分图像中与初始自适应脑模板高斯滤波后剩余像素点对应的像素点位置的图像强度值增加1，用关系是表示为：
imagec{voxel{x,y,z}}＝imagec{voxel{x,y,z}}+1；
其中，voxel{x,y,z}为自适应脑模板中没有被滤掉的像素点，imagec为三维积分图像。
当然，根据需要也可以是其他数量的迭代次数。
进一步地，上述步骤s600，具体为：
s601：采用图像形态学算法去除所述掩膜图像中的独立区域；
这里的独立区域一般是高亮的小区域，其一般为颅骨或颅外像素区。
s602：按照掩膜图像的分辨率对标准脑图谱进行三维降采集，获取与掩膜图像矩阵、空间分辨率均相当的脑图谱图像；
s603：根据所述脑图谱图像所提供的图像强度及位置信息，建立所述脑图谱图像的概率密度图，以此概率密度图为先验知识，通过统一化分割的方法将所述掩膜图像进行迭代的配准及分割，获得所述掩膜图像各层的图像强度阈值；
s604：将所述掩膜图像中图像强度值大于所述图像强度阈值的像素点赋值为1，否则赋值为0，获得二值化的颅内像素掩膜图像。
为进一步描述上述方法，下面以用户孤独症模型大鼠为例进行说明。
以每组数据包括40只实验大鼠，例如但不限于在brukeravance7.0-t磁共振成像仪上对实验大鼠的大脑和嗅球区域进行高分辨的脑功能成像，成像过程采用epi扫描序列和表面线圈。其中，第一只实验大鼠的脑成像次数多于其他大鼠，且其余大鼠的脑成像次数相等。例如但不限于，第一只实验大鼠的脑成像次数比其余大鼠多十次。在进行成像前，将实验大鼠放入麻醉诱导盒内(异氟烷浓度控制在2％-3％之间)，待大鼠彻底麻醉后，立即放在磁共振动物床上采取俯卧姿势摆上扫描床固定，采取俯卧姿势放置，同时利用进行呼吸麻醉方式(异氟烷浓度控制在1％-1.5％之间)维持浅麻状态，进行扫描。扫描参数为：图像矩阵大小为64×64×30，图像分辨率为0.44×0.44×0.6mm3，扫描层厚为0.6mm，采用无间隔隔层扫描，tr＝2000ms，时间序列长度为320volumes，图像存储格式为dicom。
在扫描获得过图像后，为使基于matlab软件开发的分析程序可以对采集到的全脑扫描图像进行分析处理，首先通过bruker2analyzeparavisionconverter(bru2anz)将文件2dseq(dicom数据格式)批量dicom数据转存为可分析的nifti格式，其中2dseq为包含时间域的4d数据格式。
然后对格式转化后的各全脑扫描图像进行时间片校正，时间片校正可以采用以下数据采集参数：tr＝2s；slices＝30；图像采集顺序：隔层采集；参考层：slice002。
时间片校正后，剔除第一只实验大鼠前十次脑成像的全脑扫描图像，以消除磁共振成像仪前期磁场不均匀对数据分析结果带来的不利影响。
然后，以每只被试大鼠成像时间序列中的第一个全脑扫描图像为头动校正参考图像，分别以每只被试大鼠的头动校正参考图像为参考，对该只被试大鼠的剩余各全脑扫描图像分别进行刚体变换以获得头动校正图像，对每只被试大鼠的头动校正参考图像及头动校正图像进行平均，获得每只被试动物的平均图像。
对40只被试大鼠的平均图像进行迭代，生成标准脑模板，迭代方法参见上述表述，这里不再赘述。任选一平均图像进行空间直角三坐标系旋转，旋转后作为初始参考图像，迭代使用图像培训准和图像平均的方法获得适用于该样本的由用户数据驱动的初始自适应脑模板。
然后，通过高斯滤波及形态学卷积的方法，去除用户数据驱动的初始自适应脑模板的脑外像素，获取仅包含脑内像素的掩膜图像，并将掩膜图像与标准脑图谱(paxinosandwatson，大鼠脑立体坐标图谱)进行空间对齐，并获得空间变换参数。
通过上述空间变换参数，将标准自适应脑模板统一至标准脑图谱，并通过高斯平滑处理，进一步去除随机噪声，再通过线性漂移和带通滤波，仅保留0.01hz—0.08hz的图像信息，根据标准脑图谱中的感兴趣区的信息，获取所有被试大鼠的所有感兴趣区的时间序列。
最后，通过皮尔逊相关算法，计算所有被试大鼠两两roi之间的相关性，仅保留显著性p>0.95的相关系数，也即获得所有被试大鼠的全脑网络，并保存为excel文件，供用户进行后续的统计分析。
以上描述仅为本申请的较佳实施例以及对所运用技术原理的说明。本领域技术人员应当理解，本申请中所涉及的发明范围，并不限于上述技术特征的特定组合而成的技术方案，同时也应涵盖在不脱离所述发明构思的情况下，由上述技术特征或其等同特征进行任意组合而形成的其它技术方案。例如上述特征与本申请中公开的(但不限于)具有类似功能的技术特征进行互相替换而形成的技术方案。