admin 管理员组

文章数量: 887021


2024年2月25日发(作者:urlgot官网)

转载者cell注: 本文件转存自蓝宝社区/, 并经作者同意首转发于心心水滴论坛.

本教程适合05年afni的版本, 软件的新增功能请登陆afni官方网下载查询,或者登陆心心水滴论坛浏览其相关介绍.

登陆蓝宝或心心水滴论坛,共同交流.

AFNI使用指南(连载)--精心整理的AFNI教程

特别奉献精心整理的AFNI使用指南,以飨读者。

目 录

0 前言 ………………………………………………………………… 1

1 如何获得并安装AFNI …………………………………………… 1

2 基本概念 …………………………………………………………… 1

2.1 数据集 …………………………………………………………… 1

2.2 数据集的存储 ………………………………………………… 2

2.3 数据集块和子数据集块 ……………………………………… 2

2.4 集合或路径 ……………………………………………………… 2

2.5 文件类型 ……………………………………………………… 2

3 AFNI的交互界面 ………………………………………………… 2

3.1 启动AFNI ……………………………………………………… 2

3.2 AFNI界面 ……………………………………………………… 3

3.3 在Batch Mode下使用AFNI命令 …………………………… 4

3.4 AFNI的Plugins ………………………………………………… 4

4 数据分析基本步骤 ……………………………………… 插页1~13

5 数据处理基本步骤 ………………………………………………… 5

5.1 原始功能像的重建 …………………………………………… 5

5.2 转换成AFNI格式 …………………………………………… 5

5.3 功能像层面时间校正和运动校正 …………………………… 7

5.4 功能像的时间域滤波 …………………………………………… 12

5.5.a 功能像的空间平滑 …………………………………………… 12

5.5.b 补充内容 ……………………………………………………… 16

包括去除头皮外伪影、去除异常时间点、时间序列时间点数据的标准化

(difference/mean)、多个run的数据连接及去线性漂移。

5.6.a 单个被试通用线性模型分析及结果显示 ……………………… 22

包括1D 参考波创建、t-检验及相关分析(及Cluster Analysis)、3dDecon-

volve反卷积分析。

5.6.b 激活图的clustering, montage, and render ……………………… 31

5.7 空间标准化(切换至Talairach坐标系) ……………………… 34

5.8 统计图的空间平滑 …………………………………………… 39

5.9 被试统计图的组分析 …………………………………………… 39

6 ROI(Region of Interesting)绘制及分析 …………………………… 45

6.1 手动绘制ROI及ROI分析 ……………………………………… 45

6.2 从激活图创建ROI …………………………………………… 51

附录1 DICOM图像标准 …………………………………………… 52

附录2 rnxpc源程序 ………………………………………………… 52

附录3 应用于3dDeconvlve分析的脚本源程序 ………………… 53

附录4 AFNI Programs & Plugins ……………………………………… 56

0 前言(Preface)

MCW AFNI是“Medical College of Wisconsin Analysis of Functional NeuroImage”的缩写形式。由美国Wisconxin医学院生物物理研究所开发研制,主要开发者为Cox博士。

AFNI是一个交互式的脑功能成像数据分析软件。它可以将低分辨率的脑功能成像的实验结果叠加在具有较高分辨率的脑结构成像上进行三维显示。

AFNI程序分为两种,一种是利用AFNI界面本身直接运行的程序(GUI mode),另一种是脱离AFNI界面执行的辅助程序(Batch mode)。另外还提供了可供扩展功能的Plugins。

1 如何获得并安装AFNI(How to download and install AFNI)

AFNI可运行在多个操作系统下,推荐使用Linux系统。可在/afni下载源程序linux_ ,并解压至/usr/local/bin目录下,可以直接使用。如提示无执行权限,可使用chmod 777 *命令修改文件权限。

(Linux的使用技巧:输入命令后,长文件名可以输入部分字母,然后按Tab键自动补上;对于较长的命令行,可在行尾加上空格和””回车再续行)

2 基本概念(Fundamental AFNI concepts)

数据集(datasets) 和集合(session) 是AFNI 中的两个基本概念,下面作简要的介绍。

2.1 数据集(datasets)

AFNI 中的基本数据单位是数据集(datasets),它是指由一个或者多个图像的3D 数据块(bricks)以及与之相关的附加信息所组成的数据集。 数据集有两种基本的类型:解剖数据集和功能数据集。AFNI 在进行数据集处理时在任何时候都是以解剖数据集作为背景,再将功能数据集叠加到解剖数据集上。

对于功能数据集,它有5 种类型:fim、fith、fico、fitt和fift。 其中各个类型的具体解释如下:

(注:f为function功能,i为intensity强度,th表示threshold阈值,co表示correlation相关,tt表示t-test,ft表示F-test)

fim 表示功能强度,每个体素用一个值表示。

fith 表示功能强度及阈值,每个体素用2 个值表示,第一个值表示“强度”,第二个值表示阈值,用来表示哪一个点是激活点。

fico 表示功能强度及相关性,每个体素用2 个值表示,第一个值表示“强度”,第二个值表示相关系数( -1.0 到1.0 之间) ,当给出一个显著性值p 时,确定激活点。

fitt 表示功能强度及t-test,每个体素用2 个值表示,第一个值表示“强度”,第二个值是t 检验值,当给出一个显著性值时,确定激活点。

fift 表示功能强度及F-test ,每个体素用2 个值来表示,第一个值表示强度,第二个值是F检验值,当给出一个显著性值时,确定激活点。

对应于每个体素的值,可以有4 种数据类型:byte、short、float和complex。

Byte型为8位有符号整数,表示范围为0~255,最多有3位有效数字;

Short型为16位有符号整数,表示范围为-32768~32767,最多有5位有效数字;

Float型为32位实数,最多有7位有效小数。

在AFNI 中,数据集存储为两种类型的文件:头文件(header) 和块文件(brick) 。

所有的数据集都有规定的命名方式,其基本形式是:prefix + view. NAME。

prefix 由用户指定;“+ view”表示数据集的显示形式(坐标系),由AFNI自动生成,通常有3 种类型即:+orig, +acpc,+tlrc,分别表示原始数据、ACPC和Talairach坐标系。NAME 可为BRIK或HEAD。

块文件存储了所有的3D 原始数据(Brick)及由AFNI程序计算衍生的统计参数数据。

头文件包含了所有的辅助信息,它提供了解释块文件的所有信息,以ASCII 的形式存储。通常包含以下信息:

(1) 每个体素x,y,z方向的大小(mm);

(2) 数据集的轴向:例如,X-axis为R-L,

Y-axis为AP,Z-axis为I-S,则为水平层面;

(3) 数据集在扫描坐标的定位;

(4) 3D数据集每个sub-brick之间的时间间隔

3D+time数据集为fMRI的基本数据集;

(5) 统计结果衍生的功能数据集,头文件包括与统计方法有关的参数,如t-test、F-test的自由度。

2.2 数据集的存储

然而,并不是每一个数据集都必须包含一个BRIK文件。当需要时图像的显示可以由真正含有.BRIK文件的父数据集变换而来。这种功能称为“Warp-on-Demand”,如空间标准化时生成的Talairach图像。

在程序设计时,图像数据集数组可以有两种方式实现存储。一种是使用malloc函数分配内存空间,另一种是使用UNIX

的mmap 函数。mmap 函数直接将.BRIK文件映射到内存地址空间。这种映射是以只读方式实现。

2.3 数据集块和子数据集块

一个数据集可以含有一个或多个3D 数据集子块(sub-bricks)。例如,一个3D+time数据集本质上是由包含3D 数据集子块构成的数组,每一个时间点的数据可以是一个数据集子块(sub-brick);又如,一个bucket数据集也是由多个3D 数据集子块组成。

2.4 集合或路径(session)

包含一系列数据集的路径称为一个session,AFNI 只读取形式为*. HEAD 和*. BRIK 的文件。也就是说,只能在某一个实际路径中运行AFNI。

所有存放在同一session的数据集,如果显示格式相同,则认为x,y,z坐标是配准的。所以可将一数据集(常为功能数据集)重叠在另一数据集上(常为解剖数据集),即使其轴向和体素大小不一致。

通常,在同一session目录中,是从一个被试一次扫描session中获得及其衍生的数据,通常包括:

(1) 解剖参考数据集(SPGR 或MP-RAGE);

(2) 10~20个3D+time EPI function runs;

(3) 从3D+time数据集计算获得的统计数据集,用来显示激活;

(4) 从orig转换至tlrc的数据集,用于被试间的比较。

2.5 文件类型

AFNI除能识别BRIK和HEAD文件外,还可以识别下列文件格式:

ANALYZE (.hdr/.img file pairs): SPM, FSL所使用的格式;

MINC (.mnc): MNItools所使用的格式;

CTF (.mri, .svl): MEG analysis volumes;

NIfTI-1 (.nii): 一个新的由AFNI, SPM, FSL, and BrainVoyager达成协议的标准格式;

ASCII text (.1D): 按列排放的数字,如刺激处理参数。

3 AFNI的交互界面

3.1 启动AFNI

将当前目录切换至数据集所在目录,直接键入AFNI启动界面。AFNI从当前目录读取数据集。

AFNI –dir1 dir2 则AFNI可从列出的目录中读取数据;

AFNI –R 则AFNI可从当前目录及其所有子目录读取数据集。(-Recursively,递归地)

3.2 AFNI界面

(1) 主界面(下图)

通常,坐标系都是按RAI的次序(DICOM标准),x = Right (negative) to Left (positive),y = Anterior (negative) to

Posterior (positive),z = Inferior (negative) to Superior (positive)。

通常,通过Switch Underlay

选择高分辩率的解剖数据集,而

通过Switch Overlay选择低分辩

率的功能数据集(并且根据需要

插值至解剖像分辩率和翻转至

解剖像轴像)。

(2) Define OverLay窗口(右图)

(3) Define DataMode窗口(下图)

Plug-ins 功能(具体功能见后面的数据分析步骤)

(4) 图像和图表窗口

图像(Image)窗口

Disp用于改变显示方式,Mont (Montage)用于对Slice进行剪切(顾名思义,同电影剪切一样),如同时显示2行4例,再设置一定的间距(spacing),可以美观地显示结果。

图表(Graph)窗口

其中,Opt用于控制显示选项,FIM用来提供对功能数据集的交互的统计计算(详见后述)。

(5) 其它窗口(如to3d, Draw Dataset, Render等),分别见后面数据分析步骤。

3.3 在Batch模式下使用AFNI命令

3.4 AFNI的Plug-ins(见上述)

4. 数据分析基本原则

(请参考PDF文件)

5数据处理基本步骤

(1) 原始功能像的重建(reconstruction of raw functional images) (P-files)

(2) 转换成AFNI格式(conversion to AFNI format)

(3) 功能像层面时间校正和运动校正(slice timing correction and motion correction of functional images)

(4) 功能像的时间域滤波(temporal filtering of functional images)

(5) 功能像的空间平滑(spatial blurring of functional images) <可选>

去除头皮外伪影、时间序列时间点数据的标准化(difference/mean)、多个run的数据连接及去线性漂移可在该步骤后进行。

(6) 单个被试通用线性模型分析(single-subject analysis with GLM)及结果显示

可以对得到的图像进行提取clusters、渲染(render)及剪辑(montage),以达到最美观的效果。

(7) 变换至Talairach坐标系(/空间归一化或标准化)(Talairach transforming /Spatial normalization)

(8) 统计图的空间平滑 (spatial blurring of statistical map) <可选>

(9) 各被试统计图的组分析 (Group analysis of individual statistical maps)

在整个分析过程中,我们应该利用AFNI优秀的可视化功能对数据进行visual check。

5.1 原始功能像的重建

Scanner保存的功能数据通常以K空间形式。也就是说,数据实际上是图像数据的Fourier变换。可使用epirecon将K空间的图像进行Fourier变换生成图像。Epirecon也可以进行对不同类型的图像变形做校正。另外除了直接调用epirecon外,还可以调用sip程序,通过sip再转而调用epirecon。Sip程序同时建立一些标准的目标和日志文件,这些在以后的分析和调试中都十分用用。

其它一些与图像重建及格式转换有关的软件包括:MRIcro, ezdicom, dicom2, etc.

有关DICOM图像标准,参见附录1。

5.2 转换成AFNI格式(with program to3d)

记录实验日志很重要,这对数据分析很有帮助。

(1) 复制图像文件、排序

将光盘中的数据复制到硬盘并分类。

(笔者注:一般可根据文件大小确定图像类别,如通常最大的文件为功能像,因为它同时包含多个层面,并且文件数通与时间点数一致;其次中等大小的文件为3D像,因为其分辨率较高,通常为128个;最小的文件为解剖像,通常为20个。将三种类别的图像分别分类放到3个目录中,分别起直观简洁的目录名,如xx-epi, xx-3d, xx-ana,xx为被试名称或ID。)

3D像(spoiled grass, SPGR):128 items;

解剖像(FSE ? ):20 items

(注:以扫描所得的层面数而定)

功能像(EPI): 如果是多个run就可以命名为run1, run2, … …,文件数通常与TRs数目一致。

举例,在记录fMRI研究中,实验设计为:

提示实验开始(2s)--系统饱和>[提示(2s)-->呈现记忆材料(40s)-->提示(2s)-->回忆(30s)-->提示(2s)-->基线(20s)]-->结束(2~6s)

“[……]”中步骤重复6次。

TR为2s,则共有299个时间点,则有299个图像文件。

由于描扫方式不同(如隔行扫描),所以生成的图像文件需要进行重新排序,可以用dicom2文件查看图像并生成包括头信息的同名txt文件,再用rnxpc程序进行排序。(rnxpc源程序参见附录2)

(2) 使用t3d创建解剖数据集(SPGR的3D像和FSE解剖像)

(笔者注:各厂商生产的Scanner,原始文件名称不同,如GE的I.*文件 和Simense 的*.IMA文件)

进入已经排序的3D或解剖像目录(cd xx-3d或cd xx-ana)后,直接键入to3d *,出现to3d程序界面:

该例中,to3d已从图像文件的头信息获得必需的信息,所以不需再进行参数调整,只需在窗口的底部填写session

directory(session目录)和prefix(前缀)。(笔者注:可点击View Images按钮查看图像)

对于无头信息的裸图像(naked images),则需进行参数设置,这时记录翔实的实验日志非常关键。执行to3d命令后出现界面:

上图黄色矩形框中的警告文字表示发现图像中有负值,需点击Byte Swap [2]按钮进行转换(未转换前的图像如下图示)。

出现这种现象是因为Intel/Linux和Sun/SGI系统使用的数字内码的不同。

另外,在窗口上部需填写一些信息,左面为数据集的轴向信息,中间为体素(像素)和视野的大小,右面为每一轴向第一层面的位移。轴向信息可由查看图像窗口获知(横跨屏幕的为x轴,上下方向为y轴, z轴可通过滑动图像窗口下部的slider确定);体素和体素大小及层面位移需根据实验日志填写。

对于3D(SPGR)像,选择irregular(即x,y,z各不同);对于功能像(EPI)和解剖(SE)像,选择square(即x=y,需填写z方向体素大小)。

(3) 使用to3d创建功能数据集(EPI)

功能数据集是EPI 3D+time数据集。需要指定一些参数才能运行,命令行如下:

to3d –time:zt nz nt TR tpattern (e.g. to3d –time:zt 20 299 2s altplus)

*笔者注:层面通常先以space (z)排序,再以time (t)排序;-time:zt就表示如此,20为全脑的层面数,299为时间序列的帧数(时间点数或TRs);根据实验具体参数略有不同。2s为重复时间(TR),如果填写0,则为从头信息读取。altplus为数据获取方式,可分为正向间隔获取(altplus/alt+z)、逆向间隔获取(altminus/alt-z)、正向顺序获取(seqplus/seq+z)、逆向顺序获取(seqminus/seq-z)和同时获取(simult)等。其它的包括zero(3D扫描)和@filename时间模式包含在文件中)。其中alt为alternatively, seq为sequentially。(对于排列方式未知的图像,可用aiv命令查看)

另X轴方向应设为右到左(R--L);Y轴设为从前到后(A—P);Z轴设为从下到上(I—S)(RAI坐标系)。AFNI并不能识别左右侧,因此对新机器或左右侧信息对研究具有重要意义的数据,在采集时做标记,如在左侧贴上几颗鱼肝油丸。

设置好各参数后,选择Save Dataset将出现下面的窗口:

窗口中指出奇异值(outliers, 即与同一时间序列中其它数据点显然不同的值);通常早期的outlier可能是长轴磁场饱和至稳态值之前的瞬时效果导致;而以后的则可能由头动、scanner故障导致。

5.3 功能像层面时间校正和运动校正

(1) 时间校正

对于一个给定的脑体积,每个层面图像获取时间稍有不同。层面的精确时间依赖于成像的设置。通常,层面以顺序获取,或以隔行获取(偶数行先获取,然后再奇数行)的方式进行。尤其是在后者,相邻层面的时间位移将导致相邻层面的信号变化显著,尤其是在事件相关实验设计中。虽然,这样的差异不会影响个别的体素分析,但任何涉及层面间功能时间序列平均、插值的处理过程将会受层面时间差异的影响。

这些处理过程包括3维空间平滑、cluster平均、运动校正以及空间标准化(e.g. to Talairach space)。当对功能数据集进行任一上述处理时,就需要对层面时间差异进行校正。在AFNI中,可以使用命令3dTshift实现,但一个更简单地方法就是在进行运动校正时使用-tshift选项(见下面的例子)。时间校正在运动校正前进行,因为运动校正涉及邻近层面的插值,如果邻近层面存在时间差异会使插值不准确。

Usage: 3dTshift [options] dataset

将输入数据集的体素时间序列进行移位,从而使各层面对齐到相同时间原点(temporal origin)。默认地,使用头文件中的层面间移位信息(slice-wise shifting information)(由to3d程序输入的”tpattern”信息)。

 插值(interpolate) 方法: 去除趋势(detrend) 恢复趋势(retrend)

输入数据集可以接一个子数据块选择器(参考3dcalc –help)

输出数据集时间序列将重新插值至新的时间格(temporal grid),这也许不是分析数据的最好方法,但非常便利。

注意:

* 请注意混叠(aliasing)现象:超过1/(2*TR)以上的频率不能被正确的插值。对于绝大多数3D fMRI数据,这意味心率和呼吸频率不能在该程序中被正确的处理。

* 对于高速fMRI 成像的开始的图像通常与较后的图像具有不同的质量,因为纵向磁化至稳态值之间的瞬时效果(如前所述的outlier)。这些图像应该不要包括在插值范围内。举例,如果你希望排除刚开始的4个图像,那么输入的数据集应按’prefix+orig[4..$]’形式指定。或者,可以使用’-ignore ii’选项。

* 最好在3dvolreg之前使用3dTshift。

Options:

-verbose = 在程序运行的时候显示一些信息

-TR ddd = 用'ddd' 作为TR值, 而不使用包含在数据集头信息内的值。可以附后缀 's' 或 'ms'

分别表示秒和毫秒

-tzero zzz = 配准每个层面的时间们移为 'zzz'; 'zzz' 值必须介于最小和最大位移之间

注: 默认配准时间为'tpattern' 值均值(来自数据集头信息或-tpattern选项)

-slice nnn = 配准每个层面的位移为 'nnn' 层面的时间位移

注:选项-tzero或-slice只能用1个

-prefix ppp = 用 'ppp' 作为输出文件的前缀。默认是'tshift'

-ignore ii = 忽略最前面'ii' 个时间点(points)(默认ii=0)。最初的ii个值在输出文件中不发生改变

(不管-rlt选项);也不会用于detrending 和time shifting

-rlt or –rlt+ = 在位移前,去除每个时间序列的均值和线性趋势。默认在位移后再加回这些值。

-rlt表示在输出中去除这些值;而-rlt+表示在输出中只加回均值。

-Fourier|linear|cubic|-quintic|heptic表示采用傅立叶或一、三、五、七次拉格郎日多插项插值方法

-tpattern ttt = 使用'ttt' 作为时间模式, 而不采用输入的数据文件头信息中包含的时间模式。

tpatter的定义可以为altplus(=alt+z), altminus(=alt-z), seqplus(seq+z), seqminus(seq-z),

@filename

举例:如果nz = 5, TR = 1000, 那么inter-slice time为dt = TR/nz = 200。在这个例子中,层面位移为下列数字:

S L I C E N U M B E R

Tpattern 0 1 2 3 4 Comment

————————————————————————————

altplus 0 600 200 800 400 Alternating in the +z direction

alt+z2 400 0 600 200 800 Alternating, but starting at #1

altminus 400 800 200 600 0 Alternating in the -z direction

alt-z2 800 200 600 0 400 Alternating, starting at #nz-2

seqplus 0 200 400 600 800 Sequential in the -z direction

seqplus 800 600 400 200 0 Sequential in the -z direction

——————————————————————————————

如果使用@filename作为tpattern, 那么nz个以ASCII码存储的数字从文件中读取,并作为每个层面的时间位移(offsets)。

输入数据集可以是:

'r1+orig[3..5]' {sub-brick selector}

'r1+orig<100.200>' {sub-range selector}

'r1+orig[3..5]<100..200>' {both selectors}

'3dcalc( -a r1+orig -b r2+orig -expr 0.5*(a+b) )' {calculation}

(2) 运动校正

将不同方式和不同时间获取的图像进行对齐,从而利于体素-体素(voxel-by-voxel)比较。这样功能时间序列将更少地受被试运动的影响;如果图像正确的校正的话,可以比较不同session的结果。

绝大多数图像校正都使用pairwise alignment方法。即给出一个基准图像(base image)J(x)和一个要进行校正的图像(target image)I(x),寻找一种几何变形(geom. transformation)T[x],使得I(T[x])≈J(x)。

T[x]依赖于一些参数:目标是寻找一些参数使得变换后的I与J最拟合。为校正整个时间序列,每个3D Volume In(X)都通过自己的变换Tn[X]进行与J(x)配准,n=0,1,…;所以结果也是一个时间序列In(Tn[X]),用户必须选择基准图像J(x)。

绝大多数图像校正都需进行3种算法选择:如何测量I(T[x])和J(x)之间的误差E?如何调正T[x]的参数使得E最小?如何对I(T[x])进行插值至J(x)的网格(grid)?

在fMRI图像处理过程中,刚体模型配准问题可以分为两步:(1)刚体模型:首先通过迭代法估计出描述空间坐标之间变换参数的最佳值,然后用这些参数对需要配准的图像进行空间变换和必要的内插处理。(2) 图像重取样(re-sample):图像重取样是决定变换到新的空间坐标之后每个体元的值。经过变换之后的体元位置大多数情况下不是正好一个体元位置,所以需要通过插值法重新取样。方法包括取最相邻体元的值(0阶重取样)或多点线性插值(一阶重取样)。几种插值方法可以使用,默认的是Fourier,最准确但速度最慢。其它的有1, 3, 5, 7次拉格郎日多项式插值(linear, cubic, quintic

and heptic)。

目前的AFNI程序通过灰度(强度)值进行图像的配准。

E=平方差的加权和= Sx w(x) • {I(T[x]) - J(x)}2

EPI, 但SPGR和EPI之间不行。SPGR, EPI所以只对类似的图像可以进行配准,如SPGR

□ 3D RegistrationPlugins3dvolreg或Define Datamode

用来对3D volume(sub-brick)进行对齐。T[x]有6个参数:R-L, A-P, I-S轴的位移及沿I-S, R-L, A-P轴的旋转(分别对应于Roll, Pitch, Yaw)。常用于session内(intra-session)和session间(inter-session)的对齐。对于发生在单个TR(2-3s)内的运动则无法校正。

Usage: 3dvolreg [optinos] dataset

e.g. 3dvolreg -base 4 -heptic -clipit -zpad 4 -prefix fred1_epi_vr -dfile fred1_vr_dfile fred1_epi+orig

-base 4 Þ 选择输入数据集(fred1_epi+orig)的子数据块(sub-brick) #4作为基准图像J(x)

也可以写成-base fred1_epi+orig[4]。可以使用不同的图像作为基准图像,但大多数情况下,

但最好使用最靠近解剖像扫描的时间点,因为这时描扫参数相似。

-heptic Þ 选择7次拉格郎日多项式插值方法(见前述)

-clipit Þ 将负的体素值设为0(注:负值是由于高次插值方法导致的伪影)

-zpad 4 Þ 在进行shift/rotation前,将每个耙图像(target image, 即I(x))垫0四层(Zero padding),最后再去掉。对于Fourier插值法,Zero padding非常需要。同样多项值插值方法也适合,因为如果有较大的旋转,如果不垫0的话将会有数据丢失

-prefix fred1_epi_vr Þ 指定输出文件的前缀

-dfile fred1_vr_dfile Þ 将估计的运动参数输出至指定的1D 文件(以后可使用1dplot绘图显示)

为查看是否有较大的平移和旋转,可以查看运动参数文本文件:

1dplot -volreg -dx 5 -xlabel Time ‘fred1_vr_dfile[1..6]’

( [1..6] 指出运动参数文件中6 列包括平移和旋转的估计值)

可以看出,在160s附近有最大的运动,因为被试在此时刻头动了一下。

3dvolreg可以相当好地处理小的运动(motion),但是较大的运动 (> 1mm) 不能被正确地校正。可以用AFNI查看运动校正后的数据集看是否是这样。时间序列会显示什么样(i.e. 是否有运动有关的数据不连续点)? 沿着时间轴是否还有明显的运动?

3DPlugins同样的操作,可以通过点击Define DateMode Registration完成,界面如下:

Datasets用于选择输入和输出数据集;Parameters用于选择基准图像和Resampling插值方法;Outputs用于指出1D输出文件名称。填写完毕,点击Run+Close即可。

□ 2DPlugins2dImReg或Define Datamode Registration

2dImReg用于对齐2D层面图像。T[x]有3个参数(x-, y-轴的位移和z-轴的旋转)。 对于矢状面的EPI扫描,如果被试点头的运动快于TR,则无法被3dvolreg校正,所以只能用2dImReg校正。如果有可能的话,在3dvolreg之后再运行2dImReg去除点头运动非常有意义。

2dImReg使用比较简单:

2dImReg -input fred2_epi+orig -basefile fred1_epi+orig -base 4 -prefix fred2_epi_2Dreg

2dImReg的功能同样可由Plugins来完成,点击2D Registration, 出现界面如下:

同3D Registration,分别选择输入和输出数据集有基准图像,点击Run+Close即可。

□ 同被试不同session之间的校正(对同一被试进行持续多日的研究)

如果进行纵向研究或对长时期神经行为(如Learning)的研究,则需进行inter-session registration。Session之间的变换通过计算每个Session高分辩率的解剖像来实现。因为to3d可以定义1个session内EPI和SP-GR之间关系;而3dvolreg可以计算不同session间的关系;所以可以将EPI数据集从session 2转换至sessoin 1的轴向。

Inter-session校正存在的问题:

a. 被试的头放置位置不同(方向和定位都不同);所以xyz坐标和解剖位置都不对应。(如下图,Day 1和Day 2的层面位于同一体素的并非同一组织,需进行旋转,然后进行平移。

b. 扫描覆盖的解部结构不同。

c. EPI和SPGR之间的几何关系在不同session之间不同。

d. 层面厚度也可能不同(不过可以尽量避免此点)。

如下图,进行旋转时还需要注意一个问题,即旋转的原点应该相同(以SPGR 中心),否则需要先进行中心平移。

上述问题的解决方案:

a. 在旋转后再加上适当的平移(on top of熟练掌握, 在...之上, 另外, 紧接着);

允许不同天进行的实验(E1-E2)之间进行xyz方向的平移,以及EPI和SPGR之间(E1-S1和E2-S2)的中心平移(center

shifts)

b. 将EPI数据集垫加(pad)额外的层面(zeros slices),这样进行对齐的数据集可以包含所有session的所有数据

c 对数据集进行Zero padding可在to3d(创建数据集)时进行,也以后使用3dZeropad进行

d 3dvolreg和3drotate可以进行zero pad 从而使输出符合网格父数据集(“grid parent” dataset)的大小和位置。

进行intra-subject S2-to-S1变换的步骤

a. 计算S2-to-S1的变换关系

3dvolreg –twopass –clipit –zpad 4 –base S1+orig –prefix S2reg S2+orig

b. 旋转/位移参数保存在S2reg+文件中

c. 如果以前没有做(e.g. in to3d),对E1数据集进行Zero pad

3dZeropad –z 4 –prefix E1pad E1+orig

d. 对E1数据集进行intra-session校正

3dvolreg –clipit –base ‘E1pad+orig[4]’ –prefix E1reg E1pad+orig

e. 对E2数据集进行intra-session校正,同时进行大的旋转/位移校正至Session的坐标系(需要的信息存在S2reg+)

3dvolreg –clipit –base ‘E2+orig[4]’ –rotparent S2reg+orig –gridparent E1reg+orig –prefix E2reg E2+orig

注:-rotparent 告知inter-session变换关系从哪获得

-gridparent 定义输出的新数据集的网格位置及大小

输出数据集将取决于E1reg+orig的需要进行平移和Zero padded

上面讨论的步骤没有考虑不同层面厚度(EPI and/or SPGR)时的问题。最好的解决方案是在扫描时尽量对同一类型图像使用相同的层面厚度;当然,如果没有做到这一点,可以使用3dZregrid对数据集进行线性插值,从而具有相同的厚度。

上面的步骤也没有考虑存在头信息内(from to3d, e.g. ‘alt+z’)的层面依赖的时间位移。在插值至新的旋转后的网格时,体素的值再也不能说来自某一特定的时间位移,因为来自不同层面的数据集将集合在一起。在做空间插值之前,首先对数据集进行时间位移至共同的时间原点非常有意义。时间位移可以使用3dTshift或3dvolreg –tsihft完成。

5.4 功能像的时间域滤波

数据可以在线性模型拟合前使用 3dFourier进行时间域的过滤,将会对时间序列进行相当严格的频率域过滤。

Usage: 3dFourier [options] dataset (Afni 3d+time 数据集)

-prefix 新的输出3d+time数据集名称[默认为fourier]

-lowpass f 低通滤波的界值(f Hz)

-highpass f 高通滤波的界值(f Hz)

-ignore n 忽略最开始的n个图像[默认为1]

-retrend 任何均值和线性趋势在滤波前去除,在滤波后再恢复

注: lowpass和highpass可以一起用,构建成band-pass。

例:3dFourier -prefix run_1_vr_lp -lowpass 0.05 -retrend run_1_vr+orig

低通(low pass)滤波界值(cutoff value)是基于输入刺激模型的频率谱(i.e. the 1D files),如下图(图像经过反转色彩)。

时间滤波在线性模型拟合前平滑数据,而提供更好的模型拟合。

在选择频率界值时应该保守一点,并在时间滤波后用AFNI查看时间序列。

5.5.a 功能像的空间平滑

在一些时间点,你可能想对你的数据进行空间平滑(smooth)或模糊(blur)。因为两个原因。一是许多平滑操作及空间平滑趋向于平均图像中的较高空间频率噪声,因此,具有较大范围的激活区域仍然保存,但一些较小范围的激活被去除。二是集合被试间的数据之前,需要进行空间模糊(blurring),因为脑解剖结构即使在标准化至标准脑空间(Talairach)后也存在差异,因此它们的功能激活区有可能并不精确重叠。

在数据处理过程中主要有三个步骤可以进行空间平滑。数据在重建(如使用hamming或fermi filter选项)时可被平滑。这样做具有优势:在它们进行进一步处理前消除不需要的噪声或伪影。然而,在重建过程中的空间平滑是层面内的(in-plane)。也就是说,它为每个脑层面独立应用平滑,因此不考虑层面间的不需要的较高空间频率噪声,也没有对被试间层面方向的功能激活配准产生帮助。

第2个可以进行空间平滑的时间是恰在第1次各被试内数据统计分析前。这是(必须)使用原始数据时间序列的最后阶

段。在各被试内统计分析后,进一步的处理是应用统计汇总数据,而对统计汇总数据进行空间平滑具有较低的有效性。

第3个可能进行空间平滑的时间点是由各被试内分析提供的统计图,尽管存在有效性问题(将在后面简单讨论)。

空间平滑使用3dmerge程序:

Usage: 3dmerge [options] datasets ...

□ 为每个输入数据集指定编辑选项:

-1thtoin = 拷贝阈值数据覆盖(over)强度数据(这样,是否即对阈值数据进行平滑 ?)

仅在数据集附有阈值统计值时有效。所有后继操作仅对这些取代的数据操作。

-2thtoin = 同-1thtoin,但处理时不改变阈值(from shorts to floats)。该选项仅为早期的AFNI

‘3d*’程序包提供。

-1noneg = 将负值改为0

-1abs = 取强度值的绝对值

-1clip val = 将在(-val val)范围内的强度值改为0

-2clip v1 v2 = 将在(v1,v2)范围内的强度值改为0

-1uclip val = 这些选项与上述的相似,但不在数据集后附上任何解剖的缩放因子(scaling fac-

or -2uclip v1 v2 tors)。这些选项仅在特殊的环境下使用。(‘u’意为‘unscaled’,3dinfo程序可以用

来查看缩放因子)

注:这些clip选项只能使用1个。

-1thresh thr = 使用阈值数据删减(censor)强度数据 censor the intensities

注:仅对'fith', 'fico' 或 'fitt' 数据集有效。阈值'thr'为浮点型,'fith' 和 'fico'数据

集范围为0.0 < thr < 1.0,'fitt'数据集范围为0.0 < thr < 32.7

-1blur_sigma bmm = Gaussian blur with sigma(总体标准差) = bmm (in mm)

-1blur_rms bmm = Gaussian blur with rms(均方根) deviation = bmm

-1blur_fwhm bmm = Gaussian blur with FWHM(半高全宽窗) = bmm

-t1blur_* bmm,同上述选项,只不过是对阈值数据进行blur操作。

-1zvol x1 x2 y1 y2 z1 z2 = 通过x1 <= x <= x2, y1 <= y <= y2, z1 <= z <= z2定义3D volume,包含

在该3D volume内的体素(entries) 改为0

注:数据集x,y,z的范围可以通过3dinfo程序查看(单位为mm)。

CLUSTERING

rmm为连接距离(connection distance)

Clustering的步骤:先发现一些非0体素,与指定像素距离小于rmm(中心-中心距离)的像素被包括在cluster内,然后再向外扩展。

假设voxels grid size为L mm:

L < rmm < L Þ 将共享一个面的体素连接

L < rmm < L Þ 将共享一个边的体素连接

L < rmm < 2L Þ 将共享一个角的体素连接

较大的rmm 值将跳过一些0的值把非0值连成一片。

你可以忽略真实的体素大小(也许体素并非立方体),使用-dxyz=1选项,这样可将体素大小伪装成size =1。

-dxyz=1 = 在clustering中,空间的clusters通过连通性(connectivity)来确定。使用记录在数

据集头信息内的体素的维度。-dxyz=1选项强制3个体素维度grid size为1mm。

下面各选项中,'rmm'用来确定所连接体素最大距离,'vmul'是cluster含voxel的

最小数目。(下面各选项互相排斥,即只能使用1个)

-1clust rmm vmul = 形成clusters, 以rmm 连接距离,生成的cluster至少包含 vmul个体素

-1clust_mean rmm vmul = 同-1clust,但cluster内体素强度全部替换为cluster的平均强度

[同样,-1clust_后接max, amax, smax, size, order分别表示将cluster内体素强度全部替换为clust内的最大强度、绝对值最大强度、最大有符号强度、cluster的size、cluster的size索引值 (e.g. 最大的cluster索引值为1,其次为2,……)]

* 如果rmm 为0,这意味着使用6个最近的相邻体素形成非0体素的cluster

* 如果vmul为0,那么所有size的cluster都被接受(可能结果没有用)

* 如果vmul为负数,那么abs(vmul) 是cluster保持的体素最小数目

下面的命令导致3D clusters的腐蚀和扩张(erosion and dilation)。这些命令假定1种-1clust命令已被使用。其目的是避免形成具有2(或更多)个主体部分通过细长部分('necks')连接的cluster。Erosion可以切离neck;Dilation可以最小化对主体部分的erosion。

注:对cluster(-1clust命令)内的值的操作发生在下面两个命令执行之后

-1erode pv 对于每个体素,除了rmm半径范围内pv%的体素非0,否则设为0

-1dilate 如果rmm范围内还保持有1个非0体素,恢复被前面的步骤去除的体素

下面这些filter选项互相排斥:

-1filter_mean rmm = 设置每个体素值为rmm半径内体素的平均强度值

[同样,-1filter_*后接nzmean, max, amax, smax, aver分别表示将每个体素值设为rmm半径内 非0体素的平均强度值、最大值、绝对值最大值、最大有符号值、平均强度值(与mean相同,但使用了新的代码,所以运算更快)]

对阈值的filter选项(-t1filter_* rmm)也互相排斥:各选项作用与上述filter相似。

-1mult factor = 将强度值与指定的factor相乘

-1zscore = 如果sub-brick被标记为一种已知分布的统计值,它将被转换为等价的N(0,1)

偏差(或 'z score');如果sub-brick没有被标记,则什么都不错。

上述 '-1' 选项以上面给出的先后顺序执行,而不管在命令行中输入的顺序。

注:上面3个'-1blur'选项只是提供不同的途径指定半径用于模糊(blurring)功能,其关系如下:

sigma = 0.57735027 * rms = 0.42466090 * fwhm

必不可少的卷积 (convolution) 运算是通过快速傅立叶变换(FFTs)完成的,该步骤是到目前为止所有编辑选项中最慢的操作。

其它选项:

-datum type = 强迫输出集存储为指定的类型,如byte, short或float

注:Byte类型不能表示负值,如果使用该类型,数据集内任何负值将设为0

-keepthr = 当使用3dmerge编辑附带有阈值的功能数据集,通常结果数据集为'fim'(in-

tensity only)类型,该选项告诉3dmerge拷贝阈值(threshold data)至输出数据集

(尽管阈值没有进行编辑过)

注:如果3dmerge用来连接2个或更多的数据集,该选项将被忽略

注:-datum选项对阈值的存储没有效果。(可使用-thdatum type选项)

-doall = 对所有sub-bricks一律地应用editing和merging选项

注:使用-doall选项时,所有输入数据集必须有相同数目的sub-bricks

注:针对阈值的选项(如-1thresh, -keepthr, -tgfisher, etc)不兼容-doall选项

(对-1dindex和1-tindex也不兼容)

注:各个sub-bricks的所有标签及统计参数均从第1个数据集拷贝。正由于此,

用户必须检查使用这些选项是否合适。

注:Sub-bricks辅助数据可通过3drefit程序修改

-1dindex j = 使用sub-brick #j 作为数据源,使用sub-brick#k作为阈值数据源。这样,你可以

-1tindex k = 操作输入数据集任何给定的sub-brick来产生一个输出1 brick数据集。如果愿意,

多个1 brick数据集的集合可以使用'3dbucket'整合成multi-brick bucket数据集,

注:如果这些选项没有使用,j=0 和 k=1 为默认值

下面这些选项允许你指定遮罩(mask)数据集限制'filter' 选项在遮罩内非0体素操作

-1fmask mset = 读取数据集 'mset' (可能包括sub-brick指定器)并用非0体素用为filter选项的遮

罩。Filtering将不对遮罩外的数据进行运算。如果输出体素rmm半径范围内不含

有任何被遮罩的体素,输出体素将设为0

注:仅-1filter_*以及-t1filter_*选项受-1fmask影响。

* 在线性平均filters(_mean, _nzmean, and _expr)中,不包括在mask中的体素不被

使用或计算在内。这会导致意想不到的结果。如果mask设计用来排除脑外的体

素,将会出现一些问题。

* 因此,如果-1fmask有用来去除filtering的结果中非脑组织数据,在其后应该

再使用3dcalc进行masking操作

3dmerge -1filter_aver 12 -1fmask mask+orig -prefix x input+orig

3dcalc -a x -b mask+orig -prefix y -expr 'a*step(b)'

rm -f x+orig.*

下面的选项允许你为3D 线性filtering指定任意的权重函数:

-1filter_expr rmm expr 定义一个线性filter

Filter的权重是在rmm相邻范围内每个体素的位移按表达式的估计值计算出的比例。表达式中只能使用下述符号:

r = radius from center

x = dataset x-axis offset from center

y = dataset y-axis offset from center

z = dataset z-axis offset from center

i = x-axis index offset from center

j = y-axis index offset from center

k = z-axis index offset from center

Example: -1filter_expr 12.0 'exp(-r*r/36.067)'

这进行一个以12mm为半径的Gaussian filter。在该例中,FWHM为10mm。(通常,指数的分母为:0.36067 * FWHM

* FWHM。这是结合-1fmask和高斯模糊的唯一途径)

Another example:-1filter_expr 20.0 'exp(-(x*x+16*y*y+z*z)/36.067)'

这是一个非球形的Gaussian filter。

-1filter_winsor rmm nw 选项对数据进行'Winsor' filter:

□ 形成输出数据集的MERGING OPTIONS:

即整合结果的不同途径,下面的'-g' 选项相互排斥:

-gmean = 通过平均强度值(包括0) combine数据集(默认)

(同样,-g后接nzmean, max, amax, smax, count, order, fisher分别表示通过使用非0值的均值、最大强度值、最大

绝对值、最大有符号值、统计'hits'的数目*、输入顺序**、tanh(mean(arctanh))***来整合数据集)

* 不等于0(!=0)为hits

** 如果一个体素在dataset #1中为非0,则采用#1中的值;如果在dataset#1中为0而在dataset#2中为非0,则采用#2中的值;依此类推

*** 计算每个输入的arctanh(反双曲正切函数),取平均值,再计算该平均值的tanh(双曲正切)。如果输入数据类型为'short',则输入值乘以0.0001,然后输出再乘以10000。该选项用于merging相关系数的bricks

-nscale = 如果输出数据集为short类型,不对其进行缩放(scale)

应用于阈值数据(Threshold data)的merging选项:

即整合阈值数据的不同途径。如果不使用这些选项,阈值将不进行merge,输出数据集中也将不附有阈值。注意下面'-tg' 选项也互相排斥,但与上面的merging强度数据的'-g'选项互相独立。

-tgfisher = 该选项只对'fico' 或 'fith'输入数据集有效——功能强度加上相关或阈值(在后者,

阈值看成相关系数)。作用同上述'-gfisher'。

对COMBINED RESULTS进行后处理的可选途径

-ghits count = 删除为0的一些像素,统计hits(!=0 is a 'hit')数目

-gclust rmm vmul = 以rmm为连接半径,以vmul为最小voxel数目,去除不包含在cluster内的像素

注:'-g' 和 '-tg' 选项作用于整组输入数据集

指定输出数据集名称的选项:

-session dirname = 输出到指定的目录(default=./)

-prefix pname = 用'pname' 作为输出数据集前缀 (default=mrg)

注意事项:

** 如果仅一个数据集读入, '-g*'选项将无效,输出数据集仅为'-1*'选项作用于输入数据集(i.e., editiing options)

** 你可以使用3dbucket程序将3dmerge的输出数据集同其它sub-bricks整合

** Complex-valued 类型的数据集不能被3dmerge

** 该程序如无'-doall'选项不能处理时间依赖(time-dependent)数据集

** 输入数据集通过它们的.HEAD 文件指定,但它们的.BRIK文件也必须存在。

例选择进行空间平滑的第2种途径,我们可使用3dmerge命令对功能数据集进行空间平滑:

3dmerge -1blur_fwhm 5 -doall -session . -prefix run_1_vr_lp_5 run_1_vr_lp+orig

这将对数据集应用FWHM为5的高斯平滑。

5.5.b 补充内容

包括去除头皮外伪影、去除异常时间点、时间序列时间点数据的标准化(difference/mean)、多个run的数据连接及去线性漂移可在该步骤后进行

(1) 去除头皮伪影

为帮助理解3dAutomask程序,先介绍一下3dClipLevel程序:

Usage: 3dClipLevel [options] dataset

估计1个值,当用此值对解剖数据集进行clip时可将背景区域全部高为0。

方法:寻找所有正值>=clip value的值,取其中间值,重设clip vlaue为该中值*0.5,再重复这两步,直至clip value不改变

选项:

-mfrac ff = 用ff值代替上述算法中的0.5

-verb = 使用该选项选择冗余(verbose)输出,除在stdout上显示clip值,同时也显示stderr

注:该程序仅对byte-和short-型数据集操作,如果输入数据集中有任何体素为负值,将会出现警告。

如果输入数据集包含1个以上sub-brick,所有sub-bricks被用来建立histogram

Usage: 3dAutomask [options] dataset

输入数据集是一个EPI 3D+time数据集(不同于3dintracranial),输出数据集是一个仅含脑组织的遮罩数据集(brain-only mask dataset)

方法:

+ 使用3dClipLevel 算法寻找clipping level.

+ 在erosion/dilation步骤后仅保持最大的阈上体素的连接部

+ 保存结果为功能数据集'fim' 型

选项:

-prefix ppp = 输入数据集以'ppp'为前缀[default='automask']

-q = 不输出程序执行过程中的信息(i.e., be quiet).

-eclip = 在创建mask后,去除低于clip阈值的外部体素

-dilate nd = 向外扩大mask 'nd' 次

-SI hh = 在创建mask后,寻找最上面的体素,然后将低于该体素超过'hh' 毫米的任何体素设

为0。对于人类大脑hh=130比较合适。

例:3dAutomask –prefix r1mask r1blur+orig

3dAutomask –prefix r2mask r2blur+orig

对于有几个run的实验,需将多个run求得的mask混合在一起,可使用3dcalc程序(具体用法见下):

3dcalc –prefix combmask –a r1mask –b r2mask –expr 'step(a+b)'

(2) 去除异常时间点

在创建功能数据集时,to3d便对时间序列进行outliers的检查。

我们在进行单个被试数据统计分析前,可利用3dToutcount 和 3dTqual程序对时间序列进行检测,可剔除异常的时间点。(如在censor.1D中将相应的时间点设为0,即有效点设为1,无效点设为0)

Usage: 3dToutcount [options] dataset

计算3D+times数据集中'outliers'的数目,并将各个时间点写输出到标准输出(stdout)

选项:

-mask mset = 仅包括mask数据集之内的体素

-qthr q = 使用'q' 替换alpha计算中的0.001

-autoclip = Clip掉'小(small) ' 的体素(as in 3dClipLevel);

-automask = 与-mask选项互相排斥

-range = 输出每个时间点outlier的median+3.5*MAD(中值绝对偏差)

-save ppp = 输出新的数据集,保存outlier Q至每个体素,其中Q通过体素值v 计算而来:

Q = -log10(qg(abs((v-median)/(sqrt(PI/2)*MAD)))) 或Q=0 (如果v接近中值,即

不是outlier)。也就是说,10**(-Q)是值v服从v’s iid normal假设p值的粗略估

计。新的数据集以'ppp'为前缀。 (注:MAD和qg详见程序3dcalc)

-polort nn = 在outlier估计前,先对每个体素时间序列进行'nn' 次多项式(polynominals of

order 'nn')去势。默认nn=0,即只去除中值(常量)。

OUTLIERS are defined as follows:

* 计算出每个时间序列的trend和MAD

(MAD 即median absolute deviation, 中值绝对误差,等于时间序列的绝对值减去趋势

* 对于每个时间序列,与trend相差甚远('far away')的时间点被称为outliers, 而'far'定义为:

alpha * sqrt(PI/2) * MAD

其中alpha = qginv(0.001/N) (inverse of reversed Gaussian CDF),N = length of time series

* 一般都会出现几个outliers,但如果出现比例较大的outliers,应该认真地检查数据集。

因为结果输出至stdout,因此你可能想重定向这些结果生成文件或至另一个程序,可以如下例:

3dToutcount -automask v1+orig | 1dplot -stdin

Usage: 3dTqual [options] dataset

计算3D+time数据集每个sub-brick的`quality index'。输出为由每个sub-brick的index组成的1D时间序列。结果输出至stderr。

注:小的index值为'好',意味着该sub-brick与norm无非常显著的差异。该程序的目的就是为筛选fMRI时间序列零星的异常图像提供粗糙的方法,例如由于被试头部较大运动或scanner故障导致的伪影。

不要完全照搬该程序的结果。只是提供GUIDE帮助研究者发现数据问题。

具有更高的index值的sub-bricks应该检查问题。如何定义'much higher'表示多少由研究者决定。研究者应该用图形来检测index的值。作为guide, 程序将输出中值质量指数(median quality index)和范围median-3.5*MAD ~

median+3.5MAD(MAD为中绝绝对偏差)。超过该范围的values值得怀疑,如果quality index服从正态分布,超过该范围的值发生的概率为1%。

选项

-spearman = Quality index值为1-Spearman (rank) correlation coefficient (每个sub-brick与

median sub-brick的Spearman相关系数) [default method.]

-quadrant =与-spearman相似,但用1-(1/4)*Spearman 相关系数

-autoclip = 去除(clip off) median sub-brick中低强度的区域,这样,相关系数只在高强度的

-automask 体素(大概脑组织范围内)之间计算。强度水平的决定与3dClipLevel相同。这就

防止脑组织以外的接近0的体素使相关系数的计算发生偏斜。

-clip val = 去除(clip off) median sub-brick中低于'val' 的体素。

-range = 与每个quality index一起输出-3.5*MAD 和median+3.5*MAD值,这样可以进

行plot绘图。

注:* 这些值在任何情况下输出至stderr。仅在用1dplot时有用。较低的值median-3.5*MAD 不允许小于0。

举例:3dTqual -range -automask fred+orig | 1dplot -one -stdin

(3) 3dcalc程序

该程序对3D数据集进行voxel-by-voxel运算(limited to inter-voxel computation)

该程序假设voxel-by-voxel运算在占有相同空间并有相同轴向的数据集进行。

Usage: 3dcalc -a dsetA [-] -expr EXPRESSION [options]

举例:

a. 基于voxel-by-voxel,将几个数据集进行平均

3dcalc -a fred+tlrc -b ethel+tlrc -c lucy+tlrc -expr '(a+b+c)/3' -prefix subjects_mean

b. 建立一个简单的mask,其中只包含sub-brick中大于某一指定值的体素

3dcalc -a 'func+orig[0]' -expr 'ispositive(a-3.14159)' -prefix mask

c 在准备组分析前,对被试的时间序列进行标准化

下面的例子中,将时间序列中(epi_r1+orig)的每个强度值除以平均(mean+orig)得到改变值百分比。

3dcalc -a epi_run1+orig -b mean+orig -expr '100 * a/b ' -prefix percent_chg

d. 从统计数据集建立一个混合的mask,其中3个刺激呈现激活

注:'step' 和 'ispositive' 是等价的,可交替使用

3dcalc -a 'func+orig[12]' -b 'func+orig[15]' -c 'func+orig[18]'

-expr 'step(a-4.2)*step(b-2.9)*step(c-3.1)' -prefix compound_mask

e. 创建3维球体ROI(region-of-interest) mask。ROI内的值

将被标记为1,ROI mask外的值将被标记为0。这样,就

可以对ROI内的体素进行统计分析:

3dcalc -a anat+tlrc -expr 'step((9-(x-20)*(x-20)-(y-30)*(y-30)-(z-70)*(z-70))' -prefix

ball

f. 数据集 'a' 通常被用作输出数据集的模板(template)。下面的例子中,假设d1+orig和d2+orig由小的整数组成:

a) 当a除以b时,结果应进行缩放,这样结果值2.4不会被四舍五入成2。为避免这种小数部分被删减现象,使用-fscale选项强迫进行缩放

3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -fscale

b) 最好的方法是结果采用'float'型,可设置输出数据集为float型:

3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -datum float

c) 也许有时只想取商数的整数部分,即9/4=2, 而不是2.24,可使用-nscale选项强迫结果不进行缩放

3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -nscale

3dcalc的参数(必须包含在命令行内):

-a dnam = 读取数据集'dname' ,并在表达式expression (-expr)中称其体素值为'a'。3dcalc可以最

多使用24个数据集名(-a, -b, -c, ... -z)。

** 如果在表达式中使用某字母,但在数据集选项中没有指定,则其值设为0

** 如果字母后接一数字,那么数字被用来选择数据集中某一sub-brick(indexes start at 0)。

-expr = 表达式以单引号(quotes)括信。将表达式应用于输入数据集。

3dcalc的选项:

-verbose = 让程序在执行过程中显示冗余程序

-datum type = 强迫数据集以指定的类型(byte, short, or float)存储。[默认为第1个输入数据集类型]

-fscale = 强迫输出数据集进行缩放至最大整数范围。仅在输出数据集为byte或short时有效。

(默认即进行缩放) [默认在计算的值看上去需要缩放时进行缩放,如所有均小于1.0或

有某一值超过整型的上限]

-gscale = 同于'-fscale',介同时强迫每个输出的sub-brick有相同的缩放比例(scaling factor),例

如,这也许对3D+time数据集有用

-nscale = 不对byte或short型输出数据集进行缩放。当对mask数据(仅包含0或1)进行操作时

非常有用。

-prefix pname = 使用 'pname' 作为输出数据集的前缀[default='calc']

-session dir = 使用'dir' 作为输出数据集目录[default='./',即当前目录]

-dt tstep = 使用 'tstep'作为人为生成的3D+time数据集的TR

-TR tstep = 如果没有指定的话,默认为1s

-taxis N = 如果仅3D数据集输入,通常只有一个3D数据集被计算。使用该选项,可以强迫建

*OR* 立长度为N的时间轴。Step可选。例如,你可能想使用在你的表达式中预先定义好的

-taxis N:tstep 时间变量't' 和/或 'k' ,否则每个结果sub-brick将相同。

例如:'-taxis 121:0.1' 将产生121个时间点,并以TR 0.1为间隙。

注:也可以使用-dt选项指定TR。

注:可以通过使用'1D:n@val,n@val'指定 1D输入数据。如:-dt 0.1 -w '1D:121@0'

注:对于'-dt' 和 '-taxis'选项,'tstep' 都是以s为单位。可以后缀'ms'指定以ms为单位。

-rgbfac A B C = 对一起RGB输入数据集,3个通道(channels r,g,b) 分开以使用3dcalc运算:value = A*r + B*g

+ C*b。默认值为:A=0.299 B=0.587 C=0.114,给出了灰度强度值。例如,挑选绿通道,使用'-rgbfac 0 1 0'。注意,RGB数据集中每个通道的值为byte型,范围为0..255,因此,'-rgbfac 0.001173 0.002302 0.000447'将重新缩放到范围0..1.0 (i.e., 0.001173=0.299/255, etc.)

数据集类型:

最常用的AFNI 数据集类型为 'byte', 'short' 和 'float'。

数据集可以在每个sub-brick后附1个向量。其主要用途是允许short型数据表示非整数值,但却只有float数据的一半大小。

例如,考虑一个short数据集后附1个标量0.0001,这样就可以表示-32.768~+32.767,精度为0.001。

3D+time数据集:

新版本的3dcalc可以操作3D+time 数据集。每个输入数据集可以是下列情况的一种:

a) 是普通的3D数据集(no time);b) 3D+time数据集并指定sub-brick索引(如 '-b3');c) 3D+time数据集不使用sub-brick索引(如 '-b')。也可使用multi-brick 'bucket' 数据集。

输入数据集的名称:

可以以下列形式'prefix+view', 'prefix+', 或 'prefix+'指定。并可以附上1个sub-brick选择列表。例如:

fred+orig[5..8] or [5-8] ==> use #5, #6, #7, and #8

fred+orig[5..13(2)] or [5-13(2)] ==> use #5, #7, #9, #11, and #13

fred+orig[0..$(3)] 可以用字符'$'表示数据集中最后一个sub-brick

同时,也可以在输入数据集后使用语法(syntax)限定数值的范围,例如:

fred+orig[5..7]<100..200>

注:字符 '$ ( ) [ ] < >' 都是Linux中shell专用字符,可以将整个输入数据集和选择列表一起括在前向单引号(')内避免冲突。例如:'fred+orig[5..7,9]'。(也可以使用双引号)

1D 时间序列:

也可以输入1个'*.1D' 时间序列来代替数据集。这种情况下,每个空间的体素值在n时间点将相同,并且为时间序列文件中的第n个值。至少必须输入1个真实数据集。如果所有输入数据集为3D(single sub-brick)或是来自multi-brick的一个sub-brick,输出数据集将是一个人为生成的3D+time数据集。

例如,假设'a3D+orig'为3D数据集

3dcalc -a a3D+orig -b b.1D -expr "a*b"

坐标和预定义的值:

如果你在指定数据集时没有使用'-x', '-y', or '-z' ,体素的空间坐标将会载入这些变量。例如:表达式

'a*step(x*x+y*y+z*z-100)' 将会把所有与距离半径10mm范围内的体素清0。

同样,'-t' 如果没有被数据集或*.1D输入使用,将会作为变量保存体素的时间坐标。请指定该变量的单位,可以是ms, s 或Hertz。另外,数据集的各层面之间可能具有不同的时间位移,这就允许对't'进行运算。使用3dinfo程序可以找出数据集的结构。

同样,如果没有被使用,'-i', '-j'和 '-k'可作为存储体素空间索引坐标的变量;'-l' 作为时间索引坐标变量。

否则,没有指定的字母初值设为0。

使用下标:

3dcalc通常的运算严格控制在per-voxel的基础上。没有不同空间或时间值的'cross-talk'。使用下标特性允许相对于基准体素,指向不同位置的变量。例如:

-a fred+orig -b 'a[1,0,0,0]' -c 'a[0,0,0,2]'

这表示'a' 代替fred+orig数据集中的1个值,'b' 代表x-方向下1坐标体素, 'c' 代表相同位置time serials 下1个体素。

使用该性质,首先必须定义基准数据集(e.g., 'a')。然后基准数据集后接不同的下标(整数),表示x-, y-, z-和-t方向的位移。

为方便使用,下面的简写可代替下标表示坐标平移:

[1,0,0,0] == +i [-1, 0, 0, 0] == -i [0,1,0,0] == +j [ 0,-1, 0,

0] == -j

[0,0,1,0] == +k [ 0, 0,-1, 0] == -k [0,0,0,1] == +l [ 0, 0, 0,-1]

== -l

当位移到达边界并引用数据集范围外的体素时,可以(假设x range 0..99,+2i):

STOP => 在数据集边缘停止平移(返回边缘值 voxel 99)

WRAP => 绕至数据集相反的边缘(返回相反方向体素值 voxel 1)

ZERO => 返回0

可以通过开关选项指定: '-dsSTOP', '-dsWRAP', 或 '-dsZERO',默认为STOP。

问题:

* 不能对Complex-valued类型数据集进行运算

* 使用Differential subscripts降低程序执行速度

表达式:

运算符 + - * / **及圆括号。指向数据集的变量及选择列列,以及未指定的预定义的变量。

内建的函数包括:

sin, cos, tan, asin, acos, atan, atan2, sinh, cosh, tanh, asinh, acosh, atanh, exp, log, log10, abs, int, sqrt,

max, min, J0, J1, Y0, Y1, erf, erfc, qginv, qg, rect, step, astep, bool, and, or, mofn, sind, cosd, tand, median,

lmode, hmode, mad, gran, uran, iran, eran, lran, orstat, mean, stdev, sem, Pleg

比较有用的函数说明如下:

* qg(x) 标准正态分布的翻转累积分布函数(reversed cdf)

* qginv(x) qg函数的逆函数

* mean(a,b,c,...) 计算各输入参数的均值

* stdev(a,b,c,...) 计算各输入参数的标准差

* sem(a,b,c,...) 计算各输入参数的标准误[sem(n arguments) = stdev(same)/sqrt(n)]

下面的这些函数设计用来实验逻辑功能:

step(x) = {1 if x>0 , 0 if x<=0},

rect(x) = {1 if abs(x)<=0.5, 0 if abs(x)>0.5},

bool(x) = {1 if x != 0.0 , 0 if x == 0.0},

equals(x,y = 1-bool(x-y) = { 1 if x == y , 0 if x != y },

ispositive(x) = { 1 if x > 0; 0 if x <= 0 },

and(a,b,...,c) = {1 if all arguments are nonzero, 0 if any are zero}

or(a,b,...,c) = {1 if any arguments are nonzero, 0 if all are zero}

mofn(m,a,...,c) = {1 if at least 'm' arguments are nonzero, 0 otherwise}

argmax(a,b,...) = index of largest argument; = 0 if all args are 0

还有一些用来转变统计值的函数(略)。

(4)3dTstat程序及数据的标准化

Usage: 3dTstat [options] dataset

对3D+time数据集计算1个或多个 voxel-wise统计量,并存储在1个bucket数据集中

Options:

-mean = 计算输入数据集的均值[默认]

-slope = 计算输入数据集的平均坡(mean slope) vs. time

-stdev = 计算输入数据集的标准差

注:该值在平均坡被去除后计算

-cvar = 计算输入数据集的变异系数[cv=stdev/fabs(mean)]

** 注:可在上述2选项后加NOD关闭detrending操作。例如-stdevNOD或-cvarNOD

-MAD = 计算输入数据集的中值绝对偏差(median absolute deviation, MAD) [undetrended]

MAD = median(|voxel-median(voxel)|)

-DW = 计算输入数据集的杜宾-沃森统计量(Durbin-Watson Statisti)[detrended]

**根据统计学原理,对时间序列数据而言,如果自相关存在,就意味着一种有显著影响的因素

——时间序列没有在回归模型的考虑之中,从而使误差平方和不是最小值,这样就不能进行有效的判断。因此,应对

时间序列数据进行自相关检验。最常见的自相关检验是杜宾—沃森检验,简称DW检验,DW检验统计量的计算公式为:

et表示某一项按时间序列排列的因变量的残差项,et-1表示前一项因变量的残差项。

查Durbin-Watson表可得到两个DW界值(上、下限):Du和Dl,如果Dl < DW < 4-Du ,认为无自相关。

还有选项-median, -min, -max, -absmax, -argmin, -argmax, -argabsmax分别表示计算输入数据集的中值、最小值、最大值、最大绝对值、最小值的索引值、最大值的索引值、最大绝对值的索引值。全部都不需要进行detrend操作[undetrended]。

-prefix p = 使用字符串 'p' 作为输出数据集的前缀[DEFAULT = 'stat']

datum d = 使用指定的类型d('byte', 'short', 或'float')作为输出数据集存储类型[DEFAULT=float]

-autocorr n = 计算自相关函数并返回最前n个系数值

-autoreg n = 计算自回归系数并返回最前n个系数值

输出数据集为1个bucket 数据集。

□ 数据的标准化(vs. time)

3dTstat –prefix r1mean r1blur+orig

3dTstat –prefix r2mean r2blur+orig

3dcalc –prefix r1norm –a r1blur+orig. –b r1mean+orig. –c combmask+orig –expr '(a/b*100)*c'

3dcalc –prefix r2norm –a r2blur+orig. –b r2mean+orig. –c combmask+orig –expr '(a/b*100)*c'

(5) 3dTcat 程序和多个run连接(去线性漂移)

将输入的数据集sub-bricks连接成一个大的3D+time数据集。

Usage: 3dTcat options

选项:

-prefix pname = 使用 'pname'作为输出数据集的前缀(或-output pname, 默认为'tcat')

-session dir = 指定'dir' 输出数据集的路径(default='./'=current working directory)

-glueto fname = 添加sub-bricks至'fname' 数据集的末尾(与-prefix可相互替换)

-dry = 执行 'dry run'; 也就是说只输出显示将要做什么。

-verb = 在程序执行过程中显示一些冗余信息。(-dry 意味着 -verb).

-rlt = 分别地为从每个输入数据集载入的每个体素的时间序列去除线性超势(remove

linear trends, RLT)。也就是说,从每个数据集的数据被分别地detrended。

注:1) -rlt 为每个体素去除'a+b*t'最小二乘法拟合。这意味着mean和trend都被

去除。这样就使得用AFNI内部的FIM无法计算强度的%改变。

2) 为使得每个数据时间序列的mean在执行后都加回,使用选项'-rlt+'。这样,只

有斜坡(线性超势)'b*t'被去除。

3) 加回所有数据集时间序列的总体均值(overall mean),使用选项'-rlt++'。这样,

各个输入数据集的'a+b*t'被分别去除,在最后,再加上所有输入数据集的均值。

4) -rlt 可对shorts或floats数据集处理,但不能处理complex或byte类型数据集。

注意事项:

* TR以及其它时间轴信息从第1个输入数据集获得。如果没有输入数据集包含这样的信息,TR则设为1.0s。可以使用3drefit程序改变。

* Sub-bricks以指定的顺序输出,并且不一定是原始数据集的顺序。例如使用fred+orig[0..$(2),1..$(2)] 会导致fred+orig内的sub-bricks以间隔的顺序输出至新的数据集。如果使用了-rlt选项,从输入数据集选择的sub-bricks将重新排序至输出数据集,并且对该顺序进行detrended。例:

3dTcat -prefix normlink r1stded+orig r2stded+orig

注:concat.1D,用于有多个run相接的3dDeconvolve处理,这样会使结果更准确。

假设一个run包含200个时间点,那么该1D文件内容为:

0

200

5.6.a 单个被试通用线性模型分析及结果显示

(1) 用waver生成1D 文件

在进行线性分析前,我们需要生成1个或多个刺激1D文件,该文本文件包含了我们假想的每1类型实验刺激(事件)的假设激活曲线。每个1D文件包含数字,对应于在每个可测量的时间点因给定的实验情况fMRI信号改变的设想值(如每1行代表1个时间点,每1个文件代表1个实验情况)。最简单的例子是为on-off block设计实验设计的boxcar 1D文件。这里与on-阶段对应的每个时间点为1,而在off-阶段的时间点为0。1个更合适的模型可以通过将简单的boxcar刺激函数与假想的血流动力学反应函数进行卷积计算,程序Waver可以完成。

对于每个体素,我们都有一个BOLD信号时间序列,并且我们希望根据实验处理(事件)来解释这一BOLD信号时间序列。

功能脑成像中所有体素都伴随着时间变化。每个体素随着时间呈现一定的模式(强度变化),与特定体素特定时间的血氧水平一致。这些模式可以用来跟与任务有关的预测的血流动力学反应比较。Afni 1-D files是包含预测血流动力学反应的文件,用来和真实的血流动力学反应(实验获得的体素的时间序列)进行相关分析。然后用它们之间的相关性鉴定出在任务(或事件)中被激活的位置/体素。

Waver是用来生成这些1-D files的程序。基本上,如果你有80个TRs(举例),你想生成一个含有一列80个数字文件。Waver通常在输出文件末生成额外的TRs。所以你需要检查和编辑生成的文件(如vi, emacs等文本编辑器)以确信每个TR只有一个值。

用法: waver [options] > output_filename

创建ideal waveform时间序列文件。通常输出至标准输出(stdout),可使用重定向符(>)输出至文件。

Options: (# 表示数字; [xx] 指默认值)

-WAV = 将波形设成 Cox special [default] (cf. AFNI FAQ list for formulas)

-GAM = 将波形设成 t^b * exp(-t/c) (cf. Mark Cohen)

-EXPR "expression" = 将波形设置成指定的表达式,其必须依赖于变量 't'

e.g.: -EXPR "step(t-2)*step(12-t)*(t-2)*(12-t)"

下面这些选项为-WAV 波形设置参数:

-delaytime # = 将延迟时间设置成 # 秒 [2]

-risetime # = 将上升时间设置成 # 秒 [4]

-falltime # = 将下降时间设置成 # 秒 [6]

-undershoot # = 将下冲设置为 # times the peak [0.2] (必须是非负数)

-restoretime # = 将恢复时间设置为下冲后 # 秒 [2]

下面这些选项为 –GAM 波形设置参数:

-gamb # = 设置参数 'b' 为 # [8.6]

-gamc # = 设置参数 'c' 为 # [0.547]

这些选项为所有波形设置参数

-peak # = 设置峰值(peak value)为 # [100]

-dt # = 设置输出和输入的时间步长 [0.1]

默认是输出由上述参数定义的波形。但如果通常下面其中一个选项指定输入文件,则由该文件定义的时间序列将和上述参数定义的ideal waveform进行卷积计算。

-xyout = 以2列输出数据, 第1列是时间,第2列是波形。如果不使用-xyout选项, waver

产生1列,波形。2列输出在绘图时很有用。

-input infile = 从 *.1D 格式化的 'infile' 文件中读取时间序列; 并与波形卷积生成输出。

注: 你可以使用下标向量(sub-vector selector)选择infiles文件中特定的列,如-input 'fred.1D[3]'

-inline DATA = 从命令行 DATA中读取时间序列; 与波形卷积生成输出。DATA 以数字count

@value的形式,如-inline 20@0.0 5@1.0 30@0.0 1.0 20@0.0 2.0,表示时间序列由:

20个0, 然后5个1, 然后30个0, 然后1个1, 然后20个0,最后1个2组成。

-tstim DATA = 从命令行DATA中读取离散的刺激时间,并以这些时间点的冲激函数(delta-

functions)与波形进行卷积。在这种输入格式中,时间并不需要间隔’-dt’。如:

-dt 2.0 -tstim 5.6 9.3 13.7 16.4 指定TR为2秒,并刺激4次(5.6 s, 等)并不与TR

的整数倍对应。DATA值不能为负。如果DATA存储在文件中,你可以将其读

命令行,如:-tstim `cat filename` ,这里使用后向单引号(`)操作符。

如果你有’xmgr’ 绘形程序,你可以以下列方式预览输出结果:

>waver -dt 0.25 -xyout -inline 16@1 40@0 16@1 40@0 | xmgr -source stdin

使用 AFNI 包程序 1dplot, 可以同样实现:

>waver -GAM -tstim 0 7.7 | 1dplot -stdin

如果需要设计方波, 可使用 'sqwave' 程序。

编辑有时需要涉及不考虑特定的时间点,可以在相应的位置插入数字99999。

假设实验有3种情况(conditions): rest, control and hum。想像你有100个TRs,每个情况有一定数目的TRs。假设TRs

1-10发生在rest condition, 11-20发生在control condition, 21-30发生在hum condition,并且该模式重复2次,最后接一个rest condition.下面是用Waver程序生成该实验的ideal waveform:

>waver -WAV -delaytime 0.0 -risetime 2.0 -dt 1.0 -inline 10@0 10@1 10@1 10@0 10@1 10@1 10@0 10@1 10@1

10@0 > ContHum_Rest.1D

>waver -WAV -delaytime 0.0 -risetime 2.0 -dt 1.0 -inline 10@0 10@99999 10@1 10@0 10@99999 10@1 10@0

10@9999 10@1 10@0 > Hum_Rest.1D

waver -GAM -gamb 8.6 -gamc 0.547 -peak 100 -dt 3 -input visBox.1D > visGamm.1D

你可以用1dplot查看所生成的1D文件。

1dplot visBox.1D

1dplot visGamm.1D

结果如右图所示。

注: 通过卷积生成的1D文件通常比boxcar文件多含几个时间点,这是因为卷积过程导致。可用vi或emacs等文本编辑程序修改(Open office功能更强大)。

FIM对于生成的1D文件,也可以通过Graphick 选择该1D文件来查看。Ideal

一个更复杂的例子在事件相关实验设计(event-related design)中可能含有两类刺激(如2种conditions)。我们应该首先创建两个1D文件,每个文件与对应刺激相应的时间点为1。每个1D文件然后再和设想的HRF进行卷积创建相应的功能激活1D文件。

在建立1D刺激文件时有许多变化。你可以与任何不同的假想HRF进行卷积,包括从被试先前的任务中测量的真实的HRF(如简单的运动或视觉任务)。你也可以为每个实验状态创建不止1个刺激1D文件,这样你可以将不同的HRF加权组合与数据进行拟合。这种类型的回归模型与事件相关分析更有关。

(2) t-检验和相关分析

□ 3D数据集的t-检验(Student’s t-test)s

Usage 1: 3dttest [options] -set1 datasets ... -set2 datasets ...

比较2个数据集的均值(voxel-by-voxel)

Usage 2: 3dttest [options] -base1 bval -set2 datasets ...

比较数据集的均值与常量

输出:创建1个数据据,为voxel-by-voxel计算所得的set2的均值减去set1的均值(或常量bval)的差异。输出数据集类型为intensity+Ttset('fitt')型。每个体素的t-统计量可作为AFNI中(OverLay)的阈值。

t-testing选项:

-set1 datasets ... = 指定数据集集合作为set1。对set1的均值和set2的均值进行2-sample t-test。

注:-set1 和 -base1互相排斥

-base1 bval = 'bval' 是一常量数值,将对set2和bval进行1-sample t-test.

-set2 datasets ... = 指定数据集集合作为set2。每个set1(if used)和set2必须至少有2个datasets

-paired = 指定使用配对t-检验(aired-sample t-test)来比较set1和set2。如果使用该选项,

set1和set2必须具有相同的基数。

注:如果set1和set2数据集功能值具有配对关系,可使用配对t-检验。如果

它们之间没有配对关系,统计力(statistical power)将比非配对t-test降低。

-unpooled = 指定单独对set1和set2进行方差估计(not pooled together)。只有在未使用

-paired选项时有意义。

注:如果使用该选项,每个体素的自由度将是一个变量,而不是一个常数

-dof_prefix ddd = 如果同时使用'-unpooled' 选项,将创建包含每个体素自由度(DOF)的前缀为

'ddd'的数据集。可以使用-dof_prefix数据集将t-value转变成z-score。

命令行如下:

3dcalc -a 'pname+orig[1]' -b ddd+orig -datum float -prefix ddd_zz -expr 'fitt_t2z(a,b)'

3drefit -substatpar 0 fizt ddd_zz+orig

目前,AFNI可以直接处理其DOF参数在各体素间变化的数据集。转变至

z-score(with no parameters)是克服(get around)这种困难的一个途径。

-workmem mega = 'mega' 指定多少MB的RAM用来统计运算。默认为12,如果设置大点可回

快程序运行速度。(见下面的注意事项)

输入编辑选项(input editing options)同3dmerge程序。

输出选项也同于大多数程序。包括-session dirname, -prefix pname, -datum type三个选项。

注意事项:

** 为了节省内存占用资源,3dttest将对其输入数据集进行多次变化(passes)。为提高效率,可先用3dmerge程序生成临时的数据集,再用3dttest对这些临时数据集进行计算;尤其使用'blurring'选项时。另外,对数据集进行editing 时需要读入内存,所以将会增加内存的占用。可以通过-workmem选项增加指定内存。

** 输入数据集通过它们的.HEAD文件指定,但BRIK文件也必须存在。

** 该程序不能处理时间依赖(time-dependent)或复杂数据类型的数据。默认情况下(使用-datum option),如果第1个输入数据集为byte-或short-型,输出数据集的功能值为short型,否则为float型。

□ 相关分析

在该阶段你将应用AFNI检测哪个区域的激活与生成的1D文件定义的血流动力学反应波形相关。

选择Switch Anatomy(UnderLay),下拉窗口会出现,选择合适的数据集(预处理后的EPI数据集)。单击后选择,再点击主界面Image按钮,会出现层面图像。

点击"Graph"按钮将会显示相应体素的信号强度的时间序列曲线 (下图,黄色框中的体线是你当前指向的焦点)

在Graph窗口右下角,选择FIM => Pick Ideal =>选择合适的 *.1D file (Note: 当你选择 ideal后,它会在graph窗口中显示独立的一条线,并且在x轴每点都有对应的值)。可以对Ideal进行编辑,如进行平滑(smooth)、平移(shift)、忽略初始时间点(ignore)等操作。

FIM => Compute FIM =>fico (Graph窗口顶部将出现绿色的程序条)

Define Function(OverLay)按钮将会打开另外的GUI窗口,一个相关系数的Bar(条带)和一个Color Bar,选择4个色块,从橙色至红色。可以通过滑杆设置相关性系数阈值,将看到激活图像。

(* 该分析使用的AFNI版本较老)

重命名该分析结果: Define Data Mode=>Plugins=>Dataset Rename (e.g., normlink050912a)。

□ Cluster Analysis:

a) Define Function(OverLay)=>选择p值 (位于相关系数阈值滑杆下)。记录 p值 和r值。

b) 通常选择p值可简单地改变哪个功能激活可见。如果是第2次clustering,重新检测Switch Fucntion(OverLay)以确定使用了正确的数据集。

c) Define Data Mode =>Plugins=>3D Cluster

数据集(Dataset) : 前面分析中得到的数据集 (e.g., normlink050912a).

半径(Radius, rmm): 3.8

最小体积(Min vol): 150

阈值(Threshold): 从相关系数阈值滑杆中得到的r值(double check it).

输出(Output): Condition_ideal_norm_pvalue_minvol (e.g.,normlink050912a_3dclus)

Run and Close(运行并关闭)

然后重新查看激活图,选择合适的解剖像后,再选择"switch function(OverLay)" ,并选择3dclust 输出的数据集(e.g.

normlink050912a _3dclus)将会看到如下图的效果:

(3) 反卷积分析

Simple Regression Model VS. Multiple Regressor

Simple Regression Mode : 假设固定的HRF波形:e.g., h(t ) = t 8.6 exp(-t/0.547) [MS Cohen, 1997]

与刺激时间(timing)进行卷积计算(e.g., AFNI program Waver),获得理想的反应函数参考波r(t)。

假设基线的形状:e.g. a + b*t 表示常量加上线性趋势。

对于某个体素,将Z(t)进行拟合:

Z(t ) ≈ a + b*t +β∙r(t )

a, b为每个体素计算的未知参数。a,b 为“nuisance” 参数。

β是数据中r(t ) 的幅度 = “how much” BOLD change

Multiple Stimuli = Multiple Regressor :

通常在实验中有不止一种刺激或激活,e.g., 想看看“face activation”, “house activation” or “what” vs. “where” activity的大小,需要对每类刺激与特定的反应函数r1(t ), r2(t ), r3(t ), ….进行拟合。

每个rj(t ) 基于第j类的刺激时间,针对第j类刺激。

计算βj的幅度 = 时间序列Z(t)中rj(t )的量。

对比βs来看哪个体素在这些不同的刺激条件下具有不同的激活水平,如对β1–β2 = 0的检验。

简单回归步骤概要:选择HRF模型h(t) 为每个刺激类型建立模型反应rn(t) [使用h(t) 和刺激时间(stimulus[fixed-model regression] 选择基线模型时间序列[constant(a) + linear(b*t) + quadratic (c*t2)+timing)]

对每个voxel时间序将模型和基线时间序列整合成R maxtrix的各列movement(volreg.1D)]

^

列,求解Z=R∙β得β GroupIndividual Analysis Analysis

3dDeconvolve程序用法:

3dDeconvolve计算测量的3D+time数据集与一个给定的输入刺激时间序旬(*.1D)的反卷积。指定多个输入刺激时间序列,该程序也执行多元线性回归(multiple linear regression)。输出为AFNI 'bucket'型数据集,包括对线性回归系数(β)的最小二乘法估计、对系数的t-statistics、对某个输入刺激显著性的partial F-statistics以及对整个回归显著性的F-statistics。额外的输出为3D+time数据集,包括估计的系统脉冲响应函数。

Usage: 3dDeconvolve

**** 输入数据集和控制选项

-input fname fname = 输入3D+time 数据集名称(这里可以输入不止1个文件名称,这些数据集

将按时间顺序连接,如此这样做,则不需要并忽略'-concat'选项)

[-input1D dname] dname = single.1D 时间序列文件名称

[-nodata [NT [TR]] 仅估计实验设计(无输入数据)

[-mask mname] mname = 3d mask dataset名称

[-automask] 从输入数据集自动创建mask(长时间序列的数据集速度将很慢)

[-censor cname] cname = 用于删减异常时间点的censor .1D 时间序列文件名

[-concat rname] rname = 用于连接不同runs的连接时间点列表的文件名

[-nfirst fnum] fnum = 用于反卷积步骤的第1个数据集图像序号(default=max maxlag)

[-nlast lnum] lnum = 用于反卷积步骤的最后1个数据集图像序号(default = last point)

[-polort pnum] pnum = 与零假设(null hypothesis)一致的多项式的次数(可参详3dToutcount)

(default: pnum = 1)

[-legendre] 为零假设使用勒让德多项式

[-nolegendre] 为零假设使用幂函数多项式(power polynomials)(default is -legendre)

[-nodmbase] 不对基线时间序列进行de-mean操作(i.e., polort>1 and -stim_base inputs)

[-dmbase] 对基线时间序列进行de-mean操作(default if polort>0)

[-nocond] 不计算 matrix condition number

[-svd] 用 SVD 代替 Gaussian elimination (default)

[-nosvd] 用 Gaussian elimination替代 SVD

[-rmsmin r] r = 拒绝减化模型的最小均方根误差( rms error )

**** 输入刺激选项

-num_stimts num num = 输入刺激时间序列的数量(num>=0), (default: num = 0)

-stim_file k sname sname = filename of kth time series input stimulus

[-stim_label k slabel] slabel = 第k个输入刺激的label

[-stim_base k] 第k个输入刺激是baseline model的一部分

[-stim_minlag k m] m = 第k个输入刺激的最小time lag(default: m = 0)

[-stim_maxlag k n] n = 第k个输入刺激的最大time lag(default: n = 0)

注:-stim_minlag和stim_maxlag选项中k为第k个输入刺激时间序列,与-stim_file k sname选项对应。其中minlag为最小延迟(通常反应血流动力学反应的潜伏期,默认为0),maxlag为最大延迟,即我们预

期的HRF在刺激后持续的时间(要判断激活—— brief or long) (minlag和maxlag均以TR为单位)

[-stim_nptr k p] p = 每个TR区间刺激函数的时间点数(default: p = 1)

注:该选项需要0层面位移时间(0 slice offset times )

[-stim_times k tname Rmodel] 从文件'tname'中给定的一系列刺激时间生成第k个(k-th)反应模型。该反

应模型由'Rmodel' 参数指定,'Rmodel'可以是下面的一个:

'GAM(p,q)' = 1 parameter gamma variate

'SPMG' = 2 parameter SPM gamma variate + derivative

'POLY(b,c,n)' = n parameter polynomial expansion

'SIN(b,c,n)' = n parameter sine series expansion

'TENT(b,c,n)' = n parameter tent function expansion

'BLOCK(d,p)' = 1 parameter block stimulus of duration 'd'

(can also be called 'IGFUN' which stands for 'incomplete gamma function')

'EXPR(b,c) exp1 ... expn' = n parameter; arbitrary expressions

[-basis_normall a] 为'-stim_times'标准化所有基准函数,以具有幅度'a' (must have a > 0)。每个基准函

数的绝对值的极值将被缩放至'a'。

注:-basis_normall 仅影响-stim_times选项

[-slice_base k sname] 从文件sname中输入第k个刺激时间序列。并指定该回归量(regressor)属于基线,

并指定在输入的3D+time数据集中每个层面的每个回归量不同。Sname文件必须

精确地含有nz列,其中nz=层面数,或者只有1列,这时就等效于'-stim_file k sname'

和 '-stim_base k'选项。

注:* 不能在-stim_minlag or -stim_maxlag or -stim_nptr选项中使用该k值。

* 不能与-input1D或-nodata一起使用该选项。

* 使用该选项的意图是提供层面独立的生理学噪音回归量,e.g. 来自1dCRphase程序的结果

**** 通用线性检验(GLT) 选项:

-num_glt num num = 进行通用线性检验(GLTs)的数目(0 <= num)(default: num = 0)

[-glt s gltname] 执行由文件gltname中包含的矩阵指定的s个同时进行的线性检验

[-glt_label k glabel] glabel = 第k个通用线性检验的label

[-gltsym gltname] 从文件中以符号名读取GLT

[-TR_irc dt] 用'dt'作为-IRC_times选项中积分计算(computation of integrals)的stepsize, 默认作

用'-TR_times'选项中指定的值

**** 输出3D+time 数据集的选项

[-iresp k iprefix] iprefix = 3D+time输出数据集的前缀,该数据集将包含第k个估计冲激响应

[-tshift] 使用立方样条插值函数(cubic spline interpolation) 来对估计的冲激响应函数进行

时间位移,以纠正层面获取时间的差异。注意,该选项只影响-iresp选项产生的

输出3D+time数据集

[-sresp k sprefix] sprefix = 3D+time输出数据集的前缀,该数据集将包含第k个冲激响应函数参数

的标准差。

[-fitts fprefix] fprefix = 3D+time输出数据集的前缀,该数据集将包含(full model)拟合输入数据的

时间序列

[-errts eprefix] eprefix = 3D+time输出数据集的前缀,该数据集将包含(full model)拟合输入数据

的残差(residual error)时间序列

[-TR_times dt] 有'dt' 作为输出-irest和-sresp文件的stepsize。默认与输入3D+time数据集的时间

间隙(spacing)相同,单位为s

**** 控制输出bucket数据集内容的选项:

[-fout –rout –tout –vout –nobout -nocout] 分别表示输出F-statistics, 输出R^2 statistics, 输出t-statistics, 输出sample variance (MSE) map, 抑制输出baseline coefficients(and associated statistics)以及抑制输出regression

coefficients(and associated statistics)

[-full_first] 表示拽定全模型(full model)统计量将出现在输出bucket 数据集的第1个

[-bucket bprefix] 创建1个AFNI 'bucket' 数据集,包含不同的感兴趣的参九,如估计的冲激响应函

数(IRF)、各种系数及全模型拟合统计量。输出bucker数据集前缀 bprefix

[-xsave] 表示保存X矩阵至(仅在同时使用-bucket选项时有效)

[-noxsave] 不保存X矩阵(this is the default)

[-cbucket cprefix] 保存回归系数(regression coefficients),不包括对回归系数的检验统计量,保存的

数据集前缀为'cprefix'。如果可能,该数据集将代替bucket 数据集在-xrestore选项

中被使用。

[-xrestore ] 从以前的run中保存的''文件恢复X矩阵等。当使用-xrestore选项时,大多

数其它命令行选项被忽略。

**** 下面的选项只控制屏幕的输出

[-quiet] 表示抑制大多数屏幕输出

[-xout] 表示输出X和inv(X'X)矩阵至屏幕

[-xjpeg filename] 输出X 矩阵的图形(JPEG 文件)

[-progress n] 每n个体素输出一次统计结果

[-fdisp fval] 输出这些全模型F- statistic > fval 的体素的统计结果

-jobs J 以'J' 个jobs运行程序(子进程)。在多CPU电脑上,该选项将明显地增加程序运

行速度。J的范围为1至CPU的数目。

* 更多的平行处理,参见/afni/doc/misc/afni_parallelize

* 使用 -mask 可以加快速度

**注意事项 **

新版本的3dDeconvolve已经为大多数内部计算使用双精度算法编译。

例1 用反卷积进行简单回归分析:

例2 去除头动伪迹的反卷积分析(在基线中考虑头动参数)

3dDeconvolve -input epi_r1_reg+orig -nfirst 2 -num_stimts 7

-stim_file 1 epi_r1_ideal.1D -stim_label 1 AllStim

-stim_file 2 epi_r1_mot.1D'[0]' -stim_base 2

-stim_file 3 epi_r1_mot.1D'[1]' -stim_base 3

-stim_file 4 epi_r1_mot.1D'[2]' -stim_base 4

-stim_file 5 epi_r1_mot.1D'[3]' -stim_base 5

-stim_file 6 epi_r1_mot.1D'[4]' -stim_base 6

-stim_file 7 epi_r1_mot.1D'[5]' -stim_base 7

-tout -bucket epi_r1_func_mot -fitts epi_r1_fitts_mot

epi_r1_mot.1D为3dvolreg校正时生成的头动参数,0~5分别表示x-, y-, z-方向的位移以及Roll(I-S), Pitch(A-P),

Yaw(R-L)的参数。

去除和不去除头动伪影的激活图效果参见附页 #9。

例3 多元回归的脚本(decon_ht2)

实验介绍:Block Design, Human vs. Tool, and High vs. Low contrast, visual stimulaiton

3dDeconvolve -xout -input rall_vr+orig -num_stimts 4

-stim_file 1 stim_files/scan1to4a_hrf.1D -stim_label 1 Actions *

-stim_file 2 stim_files/scan1to4t_hrf.1D -stim_label 2 Tool

-stim_file 3 stim_files/scan1to4h_hrf.1D -stim_label 3 HighC

-stim_file 4 stim_files/scan1to4l_hrf.1D -stim_label 4 LowC

-concat contrasts/runs.1D **

-glt 1 contrasts/contr_ -glt_label 1 AvsT ***

-glt 1 contrasts/contr_ -glt_label 2 HvsL

-glt 1 contrasts/contr_ -glt_label 3 ATvsHL

-full_first -fout -tout 全模型统计量出现在最前面

-bucket func_ht2

注:键入source decon_ht2(需要几分钟)运行该脚本。

* Stim #1 = 人的运动视觉呈现, Stim #2 = 简单(tool-like)运动视觉呈现,Stims #3和#4 = 高对比和低对比率

** 表示连接成的各个run从哪个时间点开始,在该例中,不[0,108216,324] T

*** 3dDeconvolve提供为每个回归量进行单独地检验,但如果想对每个体素b weights的连合(combination)或对比(contrasts)进行检验,可使用-glt选项。

在上例中,我们有12个b weights:8个基线参数baseline parameters (每个run2个), 它们位于b vector的前面,4个回归幅度regressor magnitudes。(from -stim_file选项)

文件contrasts/contr_ 内容为 0 0 0 0 0 0 0 0 1 -1 0 0 (1行12个数),目标就是对b weights 的线性组合进行检验。例子中的contr_为比较Actions and Tool bs,检验假设bActions– bTool ¹ 0。

文件contrasts/contr_ 内容为 0 0 0 0 0 0 0 0 0 0 1 -1,用来检验bHighC– bLowC ¹ 0

文件contrasts/contr_ 内容为 0 0 0 0 0 0 0 0 1 1 -1 -1,用来检验(bActions+ bTool)– (bHighC+ bLowC) ¹

0

这些统计量显著的区域,在activit viewing tasks和grating viewing tasks将具有不同的BOLD信号改变量。

-glt_label 3 ATvsHL选项用来为结果统计量sub-bricks附1个有意义的标签(label)。

该脚本的结果(bucket dataset)如下图:

Color Overlay is HvsT contrast

a) F-检验来检测一个回归因子(model component)(降低时间序列变异)的显著性。

全模型F-检验:全模型回归量降低数据的变异性-基线回归量降低的变异性

(上图的sub-brick #0)。

对某个回归量的部分F-检验:全模型回归量降低变异性-去除该回归量的全模型降低的变异性

(上图的sub-bricks #19, #22, #25和 #28)

b) Coef sub-bricks为b weights (上图的17, #20, #23, #26)

c) T-statistic sub-brick测量每个系数的作用(检验bj = 0)

3dDeconvolve将会建立1 个新的AFNI功能数据集(名为bucket 数据集)包含有统计图,包括:t-检验,线性模型的β-权重参数,拟合参数等。你现在可以在AFNI中将统计图重叠至原始的EPI数据上进行查看。你应该检查该模型的拟合优度,因为你还会将之用于其它类型数据的线性回归和GLM分析。β-权重参数及其相应的t-test值为我们提供了对每个实验情况反应的功能激活的指标。

被试内应用GLM分析非常灵活,并且也非常复杂。

完整的3dDecnvolve脚本源程序参见附录3。

5.6.b 激活图的clustering, montage, and render

(1) 3d clustering 见5.6.a相关分析部分

(2) 图像的剪辑

完成激活图后,选择适合的视角(Axial, sagittal, coronal image),在图像

窗口下部点击Mont按钮,出现如右侧窗口,可选择横向图像数(Across)、

纵向层面数(Down)、图像层面间隔(Spacing)及边框(Border)等。

(3) 图像的渲染(Render)

RenderPlugins通过Define Datamode New可调出图像演染窗口:

*注:视角控制Roll(about I-S axis), Pitch(about L-R axis, after roll rotation), Yaw(about A-P axis, after roll and pitch)

Volume渲染的概念:

目标是创建由像素(pixel)组成的2D图像。每个像素为视线俯视3D volume时获得。

如果沿从右到左的直线查看,则沿该线条(白色)的所有数据对获得的

像素都有影响。

每个3D 体素包含一个数值。该数值决定该体素的亮度(或色彩)(如果

该体素可见)。

该数值还决定该像素的不透明度(opacity):

Opacity = 0 Þ 透明 (该体素亮度对图像无影响,即该体素不可见)

Opacity = 1 Þ 不透明 (该体素后面的任何体素都不可见)

Opacity = 0.5 Þ 半透明(该体素亮度的50% 作用于图像; 该体素后面的体素贡献其余50%)

Rendering需要较多的CPU和内存资源。可使用3dIntracranial程序对T1加权解剖像进行去除头皮及头皮外组织。有时,该步骤需要对orig数据集操作,然后再保存为Talairach坐标系。

3dIntracranial – 进行颅内区域解剖分割, 从高分辨率的T1加权解剖像去除头皮及其它非脑组织。

Usage: 3dIntracranial

-anat filename => 需进行分割的解剖数据集文件名

[-min_val a] => 最小体素值界限[默认:内部的概率分布函数(PDF)估计下限(?)]

[-max_val b] => 最大体素值界限[默认:内部的概率分布函数(PDF)估计上限(?)]

[-min_conn m] => Minimum voxel connectivity to enter(Default: m=4)

[-max_conn n] => Maximum voxel connectivity to leave(Default: n=2)

[-nosmooth] => 抑制对分割遮罩(segmentation mask)的平滑

[-mask] => 产生功能像遮罩(辅助)[默认:生成解剖像]

[-quiet] => 抑制屏幕输出

-prefix pname => 包含分割后图像数据集的前缀

** 注 **:新版本的3dSkullStrip可能会生成更好的分割结果。

例如:3dIntracranial -anat anat+tlrc -prefix astrip

AFNI可以渲染以任意轴向和体素大小储存的数据集。数据集内在地重定向re-oriente(见3dresample)至轴向层面顺序(axial slice order),所以切除方向具有意义。注意,轴向层面顺序是以+acpc或+tlrc坐标系输出的‘warped’数据集的标准轴向。

OverLay数据集也被重采样,这样其网格间格与UnderLay数据集相匹配。

按Accumulate按钮,接着按DynaDraw按钮,然后Roll按Draw按钮完成第1幅图像Plugin将载入体素值,创建直方图,并准备好进行渲染渲染示例:在Talairach视图下,打开渲染plugin,然后选择astrip+tlrc作为UnderLay数据 将会产生自不同角度的渲染。t几次

如果DynaDraw被关闭(off),必须每次手头按Draw按钮进行新的渲染;

如果Accumulate开启,渲染的图像被保存,可以使用图像查看滑杆进行回顾;如果关闭该功能,下一次生成的渲染图像将清除历史。

Load默认地,当在渲染的各图像移动时,Plugin控制器件(‘widgets’,相当于各种参数)不发生改变;可以选择Script

Widgets将使widgets显示当前显示图像被渲染的设置。

□ 改变从体素值至亮度和透明度的映射:

要想让白质完全白色:拖动#3亮度handle至顶部,超过白质的值;

要想让所有低强度的体素降低不透明度至0:拖动#不透明度handle至底部,超过直方图的波谷值。

□ 切除部分volume,这样可以看见想看见的部分:

每种cutout指定数据集中的空间的一部分将被去除(通过设定这些体素的不透明度为0)

多重cutouts可以两种方式组合:

OR Þ所有cutouts内的体素都被去除;

AND Þ只有包含在所有cutouts内的体素才被去除

(即使在AND组合方式下,也可用Must Do按钮来强迫cutout内的体素被去除)

Right of, left of, anterior to, posterior to, inferior to, superior to分别用来切除指定体素右、左、前、后、下和上面的

体面;

Behind …, front …, above …, below …,分别用来切除指定斜面的体素;

TT Ellipsoid用来切除具有相同Talairach脑比例的椭圆以外的体素;

Expr > 0用通常的数学表达式来指定要去除的脑区域:

表达式的语法与3dcalc相同;

变量‘x’, ‘y’和‘z’可用来代表数据集空间坐标位置(当选择Automate时,变量‘t’也可被使用),使表达式取正值的(x, y, z)位置的体素将被去除

举例:渲染与冠状面(xz-plane)和水平面(xy-plane)之间任意角度的倾斜层面(slab):

Slab内的点的集合可以表示成不等式:

ïy▪cos(α) - z▪sin(α) -sï < w/2

对于给定的角度=a, slab中心位移= s即层面宽度= w,渲染与垂直面呈25º倾斜的厚度为30mm,位移为20mm的冠状斜面,可以使用下面表达式:

abs(y*cosd(25)-z*sind(25)-20)-15

*sind()和cosd()为以度数为单位的三角函数。

通过使用Automate,并依赖于变量‘t’设置角度和位移,可以生成向下旋转和/或向下位移的动态图像。

□ Automate允许一次生成多个渲染图像(自动生成)

大多数(但不是全部)数字输入框具有细微浮雕边框:

当选择Automate时,具有浮雕边框的输入框可以使用含变量‘t’的表达式,e.g.

在Roll输入框输入70+5*t,然后点击Compute按钮,数据集将依次设置‘t’ 为0, 1, 2, 3, 4 来渲染数据集,这样,依次生成Roll 70º,75º,80º,85º,90º的图像。

□ OverLay 面板

点击OverLay按钮打开控制如何生成功能OverLay的面板(与Image Viewer窗口的Devine OverLay面板相似)。下面仅讨论不同的方面:

Color Opacity用来选择着色体素(即超过阈值的体素)的不透明度:

Overlaid体素的不透明度与UnderLay相应位置的不透明度不同,通常设置较高(0.5或更高);该菜单有两个特定的值:

Underlay表示该着色体素的不透明度将由来自Underlay图像的不透明度决定;

ShowThru表示不管Underlay体素的透明度是多少,该着色体素穿过Overlay而显示Underlay体素(透明脑‘glass brain’效果),。

5.7 空间标准化(切换至Talairach坐标系)

一般的fMRI实验都有许多不同的被试参加,以保证实验结果的可靠性和可重复性。对于fMRI数据处理软件来说,它的分析处理过程是基于像素的,所以我们对不同被试的fMRI图像比较时,程序一般都是直接用像素位置来判断解剖位置。问题在于,每个人的脑结构是有差异的,在不同被试的脑图像中,相同的像素位置并不代表相同的解剖结构,因此就需要一个图像处理过程来尽可能地把对应的脑激活区域在不同被试之间对应起来,这个过程就称为“空间归一化(spatial normalization)”。

(1) 各种数据集间的转换

AFNI给出了一些方法,可以用来实现各种数据集之间的相互转换。如上面右图显示。当Talairach标准图谱确定之后,可以将研究对象功能数据集映射到标准图谱,叠加到解剖数据集上,以便确定脑激活区域的相对确定的部位。

a) 父数据集

每个数据集a与三个与之相关的父数据集相互关联:πc(a)表示坐标父数据集,πt(a)表示变换父数据集,πd(a)表示数据父数据集。其中函数πc, πt, πd均以3D数据集作为输入,并以另一个3D数据集作为输出,在实现过程中,是指向3D数据集的指针,在AFNI中是功能数据集和解剖数据集。如上右图,每一列上是以坐标相互匹配的,每一行的左边

的数据集是原始数据集,而右边的数据集是由它们变换而来。

b) 坐标系统

每个数据集可以有两种坐标系统与之相对应:一种是连续坐标系统用r = (x, y, z) T表示,另一种是相应于体素的离散坐标系统,用 q = (i, j, k) T表示。通常在连续坐标系统下进行变换,然后再转化为离散坐标系统,也就是只在需要时才以体素形式表示出来。

c) 数据集之间的几何变换

每个数据集在坐标系统中有3种与之相关的几何变换关系:

i) 将a的变换父数据集πt(a)的连续坐标系统映射到a的连续坐标系统,即πt(a)中坐标为x的点在a中的坐标为

ii) 将a的变换父数据集πd(a)的连续坐标系统映射到a的连续坐标系统,即πd(a)中坐标为x的点在a中的坐标为

iii) Da将a的离散坐标系统映射到a的连续坐标系统,即数组索引为q的体素其连续坐标x = Da(q)。

d) 数据集之间的继承关系

如果数据集b具有坐标父集a ≠ b,并且a具有变换子集,则b也有从a继承而来的变换子集。

在AFNI中,将变换子集a’产生之后,则程序会自动生成b’,如果对a到a’的变换作修正之后,则这种变化自动影响到b’。在空间标准化过程中,通过标定markers对解剖数据集(a)进行变换至Talairach坐标系(a’),则功能数据集(b)可根据需要自动转变至Talairach坐标系(b’)[warped-on-need]。

(*注:实际上,同一目录下所有数据集都进行变换,但都只生成.HEAD文件,而不生成.BRIK文件)

(2) Talairach空间

自1952年,法国神经科学家Jean Talairach等人多次证明:穿过前、后连合(anterior commissures, AC; posterior

commissures, PC)的平面与端脑以及神经核团之间有十分固定的空间位置关系,这就是著名的AC-PC线。

□ Talairach坐标系的定义

4条重要的分界线(基线):

a) AC-PC线:定义了Talairach坐标系统中的水平面,此平面通过前连合的上缘和后连合的下缘,平行于下丘脑的脑沟,恰好把丘脑与下丘脑分开;

b) VCA线:与AC-PC线垂直,并通过前连合后端的平面;

c) VCP线:与AC-PC线垂直,并通过后连合前端的平面;

d) 中线:就是脑中隔(左右脑分界面)

大脑的外边界:

上面:与AC-PC线平行,并与顶叶最高点相切的平面;

下面:与AC-PC线平行,并与颞叶最低点相切的平面;

前面:与VCA线平行,并与额叶最前端相切的平面;

后面:与VCA线平行,并与枕叶最后端相切的平面;

两侧面:与中线平行,并与左、右脑最外缘相切的两个平面。

□ Talairach标准脑划分

a) 沿轴状面方向:AC-PC线以上部分8等分,分别以1,2, …, 8表示;以下部分4等分,分别以9, 10, 11, 12表示。

b) 沿冠状面方向:VCA以前4等分,分别以A, B, C, D表示;VCP以后4等分,分别以F, G, H, I表示;VCA、VCP之间3等分,分别以E1, E2, E3表示。

c) 沿矢状面方向:中线左、右各4等分,从中间至两侧依次以a, b, c, d表示。

(3) 空间标准化(Transforamtion to Talairach Coordinates)

分为三个步骤:

(a) 3D像与解剖像(Se)对齐(Rotation and Shift)

(*注:以Se像为标准,3D像的位置向Se像做对齐,误差在1mm范围内)

(b) 定义Markers来确定AC-PC轴及I-S轴(to +acpc坐标系)

(c) 定义边界点来缩放至Talairach标准脑(to +tlrc坐标系)

a) 对齐:

在AFNI中按New按钮打开另一个AFNI界面,在两个界面中分别切换Underlay(分别选择3D数据集和Se解剖数据集),然后分别点击Axial, Coronal和Sagita Images打开图像,以同一层的同一结构(最好是比较容易定位的结构,如如中脑导水管)进行比较,记录它们在两个数据集中的旋转偏差和位移偏差值。

(注:可选择Define SetLockDataMode All进行图像锁定,有助于进行比较)

然后打开终端窗口,输入命令(有机会试试,应该可以两步合一):

3drotate -prefix 3dro -rotate 2L 3P 4S 3d+orig

3drotate -prefix 3drosh -ashift 8L 0P 16S 3drol+orig

其中,-rotate表示进行旋转,数字表示旋转的度数,R, L, A, P, I, S分别表示right-to-left, left-to-right,

anterior-to-posterior, posterior-to-anterior, inferior-to-superior, superior-to-inferior。

-bshift表示before rotation shift, -ashift表示after rotation shift。

b) 定义Markers来确定AC-PC轴及I-S轴

至少要定义5个markers,分别是

AC superior edge = 前连合上缘

AC posterior margin = 前连合后端

PC inferior edge = 后边合的下缘

First mid-sag point = 中线矢状面上的任意一点

Another mid-sag point = 中线矢状面上的另一点(与上一点之间距离最好>2cm)

进行该操作时应尽可能的精确。

在orig视图下,点击主界面Define Markers(之前先开启See Markers),打开AC-PC GUI窗口如下:

首先,来定义前连合上缘和后端:

选择[AC在矢状面或冠状面中上移,直至水平面中的AC消失,然后再下移1像素将AC位于十字准心焦点的中心(水平和冠状面)水平面:寻找大脑半球之间的纤维连接冠状面:寻找胡须(mustache)样结构(上图红色圈中)矢状面:在胼胝体底部水平,穹窿下寻找前连合 superior 选择[AC将焦点移回AC的中间,向后移动直至在冠状面中AC消失,然后向前移1像素(*通常,在设置好AC上缘后,直接向后向下各移1像素)edge],然后按Set按钮 posterior margin],选择按Set按钮。

其次,再来定义后连合下缘:


本文标签: 数据 时间 进行 体素 使用