FSL 的核磁数据配准方法

本文对 FSL 软件的配准方式进行说明，并提供一个实用脚本，用于将数据转换到标准空间

FSL 的核磁数据配准方法

头动校正和时间校正

在核磁实验中，被试需要在扫描仪内接受较长时间的数据采集，在这个过程中，被试需要保持头部不动，才能保证最优的数据采集效果，但由于这个要求过于严格，被试往往无法达到。因此，研究一般采用刚体平移变换和旋转变换进行头动校正，头动校正时采用的变换为 6 个自由度的刚体变换，它们是三个方向上的平移和旋转。

%E6%A0%B8%E7%A3%81%E6%95%B0%E6%8D%AE%E9%85%8D%E5%87%86%2005e7b9c249914c81a44bd40829d8e5a7/image1.png

除此之外，由于核磁数据是分时和分层采集的，虽然从逻辑上不同层属于同一张立体图像，但它们采集的时刻却不同，是均匀地分布在一个 TR 之内。然而，在进行分析时，我们从逻辑上认为整个立体图像是在同一个时刻采集的，各个层的采集时刻都是在 TR 开始的时间点上。我们采用时间校正的方法解决这种时间上的不一致。从计算方面，我们假定信号在小时间范围内（一个 TR 范围内）具有二阶连续性，因此采用二阶多项式插值的方式对某层的时间序列进行插值，插值的目标是通过实际采集时间的数值，估计出 TR 开始的时间点时的信号数值。

%E6%A0%B8%E7%A3%81%E6%95%B0%E6%8D%AE%E9%85%8D%E5%87%86%2005e7b9c249914c81a44bd40829d8e5a7/image2.png

标准化映射

接下来，需要将被试个体空间的数据映射到标准空间上。由于被试个体空间中的大脑形态和位置与标准空间差异较大，因此需要 12 个自由度的仿射变换来进行配准，采用的计算方法是 12 个自由度的仿射变换，这是因为在欧氏空间中有 x，y，z 三个坐标轴，每个坐标轴对应的有平移、旋转、缩放、倾斜四种变换，共 12 种即 12 个自由度。在下图中，Rigid 对应平移和旋转变换，Affine 对应缩放和倾斜变换。由于受到计算资源、技术能力和使用场景的限制，目前的核磁数据分析一般不包含 Deformable 效应的校正。

%E6%A0%B8%E7%A3%81%E6%95%B0%E6%8D%AE%E9%85%8D%E5%87%86%2005e7b9c249914c81a44bd40829d8e5a7/image3.jpeg

虽然变换的各类繁多，但总可以用代数方程来表示

其中，_X_, Y 代表增广的原始和变换后图像，它们前三个元素是欧氏空间中的实际位置，它们的最后一个元素为 1，用于实现缩放；M 代表变换矩阵，它具有如下形式，其中mij 代表 12 个自由度。

%E6%A0%B8%E7%A3%81%E6%95%B0%E6%8D%AE%E9%85%8D%E5%87%86%2005e7b9c249914c81a44bd40829d8e5a7/image4.png

采用二步映射的方法来实施这个过程，第一步是将被试的核磁数据（EPI）向被试个体空间的结构像（Structural）进行映射；第二步是再将映射结果向标准空间模板进行映射（Standard），本研究采用的模板是 MNI 模板。

%E6%A0%B8%E7%A3%81%E6%95%B0%E6%8D%AE%E9%85%8D%E5%87%86%2005e7b9c249914c81a44bd40829d8e5a7/image5.png

而在数值计算过程中，由于变换是以矩阵相乘的方式实现的，因此我们总可以将二步映射表示成

而由简单的推导可知，全部的实矩阵M构成群，该群对乘法封闭。

∀M1, M2 ∈ M, ∃M1 ⋅ M2 ∈ M

这说明经过两次变换后，我们仍然可以用 12 个自由度的变换矩阵来表示变换结果。另外，由于矩阵不具有对称性，说明它们不构成可交换群。因此，二步映射的顺序不可改变。

在第一步中，对被试的功能像核磁数据进行时间平均，将平均图像映射到被试的结构像中；在第二步操作中，再将第一步的映射结果映射到标准空间中。经过两步映射后，我们认为被试的核磁数据的全部体素，与标准空间中的脑区具有确定的对应关系。

%E6%A0%B8%E7%A3%81%E6%95%B0%E6%8D%AE%E9%85%8D%E5%87%86%2005e7b9c249914c81a44bd40829d8e5a7/image6.png

%E6%A0%B8%E7%A3%81%E6%95%B0%E6%8D%AE%E9%85%8D%E5%87%86%2005e7b9c249914c81a44bd40829d8e5a7/image7.png

下面提供一个脚本，该脚本用于在 FSL 软件中，对数据向标准模板进行映射。

FSL: convert fMRI data into standard MNI space

Precondition

The subject has been analyzed using FSL, and the .feat folder has been generated.

Aim

The aim of the script is to convert the filtered_func_data.nii.gz into standard MIN space.

All the requirements of the conversion have been saved in the reg folder.

Operation

The operations will be noted if necessary. Beware that, no words will be written if the code finishes itself.

The FSL uses the flirt to estimate the affine transformation parameters. It is operated by FSL with a large time and CPU costing.

FLIRT version 6.0

Usage: flirt [options] -in <inputvol> -ref <refvol> -out <outputvol>
       flirt [options] -in <inputvol> -ref <refvol> -omat <outputmatrix>
       flirt [options] -in <inputvol> -ref <refvol> -applyxfm -init <matrix> -out <outputvol>

Available options are:
        -in  <inputvol>                    (no default)
        -ref <refvol>                      (no default)
        -init <matrix-filname>             (input 4x4 affine matrix)
        -omat <matrix-filename>            (output in 4x4 ascii format)
        -out, -o <outputvol>               (default is none)
        -datatype {char,short,int,float,double}                    (force output data type)
        -cost {mutualinfo,corratio,normcorr,normmi,leastsq,labeldiff,bbr}        (default is corratio)
        -searchcost {mutualinfo,corratio,normcorr,normmi,leastsq,labeldiff,bbr}  (default is corratio)
        -usesqform                         (initialise using appropriate sform or qform)
        -displayinit                       (display initial matrix)
        -anglerep {quaternion,euler}       (default is euler)
        -interp {trilinear,nearestneighbour,sinc,spline}  (final interpolation: def - trilinear)
        -sincwidth <full-width in voxels>  (default is 7)
        -sincwindow {rectangular,hanning,blackman}
        -bins <number of histogram bins>   (default is 256)
        -dof  <number of transform dofs>   (default is 12)
        -noresample                        (do not change input sampling)
        -forcescaling                      (force rescaling even for low-res images)
        -minsampling <vox_dim>             (set minimum voxel dimension for sampling (in mm))
        -applyxfm                          (applies transform (no optimisation) - requires -init)
        -applyisoxfm <scale>               (as applyxfm but forces isotropic resampling)
        -paddingsize <number of voxels>    (for applyxfm: interpolates outside image by size)
        -searchrx <min_angle> <max_angle>  (angles in degrees: default is -90 90)
        -searchry <min_angle> <max_angle>  (angles in degrees: default is -90 90)
        -searchrz <min_angle> <max_angle>  (angles in degrees: default is -90 90)
        -nosearch                          (sets all angular search ranges to 0 0)
        -coarsesearch <delta_angle>        (angle in degrees: default is 60)
        -finesearch <delta_angle>          (angle in degrees: default is 18)
        -schedule <schedule-file>          (replaces default schedule)
        -refweight <volume>                (use weights for reference volume)
        -inweight <volume>                 (use weights for input volume)
        -wmseg <volume>                    (white matter segmentation volume needed by BBR cost function)
        -wmcoords <text matrix>            (white matter boundary coordinates for BBR cost function)
        -wmnorms <text matrix>             (white matter boundary normals for BBR cost function)
        -fieldmap <volume>                 (fieldmap image in rads/s - must be already registered to the reference image)
        -fieldmapmask <volume>             (mask for fieldmap image)
        -pedir <index>                     (phase encode direction of EPI - 1/2/3=x/y/z & -1/-2/-3=-x/-y/-z)
        -echospacing <value>               (value of EPI echo spacing - units of seconds)
        -bbrtype <value>                   (type of bbr cost function: signed [default], global_abs, local_abs)
        -bbrslope <value>                  (value of bbr slope)
        -setbackground <value>             (use specified background value for points outside FOV)
        -noclamp                           (do not use intensity clamping)
        -noresampblur                      (do not use blurring on downsampling)
        -2D                                (use 2D rigid body mode - ignores dof)
        -verbose <num>                     (0 is least and default)
        -v                                 (same as -verbose 1)
        -i                                 (pauses at each stage: default is off)
        -version                           (prints version number)
        -help

Here, after the FSL computation is done, the fsl_rigid_register rewrites the header of the .nii file, like filtered_func_data.nii.gz

USAGE: fsl_rigid_register

Required Arguments:
   -r refvol    : reference/target volume
   -i inputvol  : input/moveable volume
   -o outputvol : input resampled to reference

Optional Arguments
   -fslmat fsmatfile  : spec explicitly
   -regmat regmatfile : get reg matrix as register.dat file
   -xfmmat xfmmatfile :  get reg matrix as MNI xfm file
   -ltamat ltamatfile :  get reg matrix as MGH lta file
   -noinitgeom : do not initialize matrix based on geometry
   -applyxfm xfmfile : do not reg, just apply xfm to input
   -applyinitxfm     : do not reg, just apply init xfm to input
   -initxfm  xfmfile : use this as an initial matrix (instead of geom)
   -maxangle maxangle : only search over +/- maxangle degrees
   -interp method : <trilinear>, nearestneighbour, sinc
   -dof dof : use dof instead of 6
   -bins bins : number of bins to use (default 256)
   -cost cost : objective function (default corratio)
      valid costs are mutualinfo corratio normcorr normmi leastsq
   -tmp  tmpdir (default is ). Implies -nocleanup
   -tmpdir  tmpdir : same as -tmp
   -nocleanup  : do not delete temporary files
   -cleanup  : delete temporary files (default)
   -subject subject  : only puts it in the register.dat file

   -version : print version and exit
   -help    : print help and exit

Since my aim is not to estimate the affine parameters again, the -applyxfm option is used to prevent reg again. It directly injects the parameters into the header of the -o outputvol.

As a result, the final script is

# Write the FSL pre-computed affine transformation parameters,
# into -o outputvol.

fsl_rigid_register -r task.feat/reg/example_func2standard.nii.gz -i task.feat/filtered_func_data.nii.gz -applyxfm task.feat/reg/example_func2standard.mat -o standard_filtered_func_data.nii.gz

The command costs several minutes in one CPU core.

The one more problem of the transformation is the size of the file increases,


# The transformed file size is 753M
$ ll standard_filtered_func_data.nii.gz
-rw-rw-r-- 1 zcc zcc 753M Sep 22 16:14 standard_filtered_func_data.nii.gz

# The raw file size is 446M
$ ll task.feat/filtered_func_data.nii.gz
-rw-rw-r-- 1 zcc zcc 446M Sep 22 15:16 task.feat/filtered_func_data.nii.gz