本文档的操作流程基于 Ububtu 18.04LTS 以及 Freesurfer 6.0，安装过程记录在另一篇文档 在 Ubuntu18.04LTS 中安装 FreeSurfer 中。本文档操作全部在终端中进行，Freesurfer 的 GUI 仅用作查看操作结果。

创建于2020-03-07

修改1：2020-03-08

​ 修正SUBJECTS_DIR的设置方法与单subject多图像的文件夹结构

​ 修正recon-all在单subject多图像情况下的命令

修改2：2020-03-12

​ 删除修改1中增加的单subject多图像情况

​ 新增多图像的多线程处理内容

修改3：2020-04-15

​ 新增DICOM转nii内容

修改4：2020-04-17

​ 新增批处理提取统计学结果的内容

修改5：2020-04-27

​ 新增PET处理内容

修改6：2020-05-12

​ 加入PET后的批处理流程

修改7：2020-07-23

​ 对PET分割限定线程数避免报错 补充了一些后续步骤

修改8：2020-11-02 新增海马子区域分割流程

SUBJECTS_DIR

首先需要用户自定义 SUBJECTS_DIR，不推荐使用放在 Freesurfer 安装目录下的默认文件夹，自定义 SUBJECTS_DIR 的命令如下：操作中严格避免中文路径的出现！！！

# 设置到桌面上
export SUBJECTS_DIR=/home/用户名/Desktop/Subjects
# 快速设置为当前文件夹
export SUBJECTS_DIR=$(pwd)

# 查看当前SUBJECTS_DIR
echo $SUBJECTS_DIR

在 SUBJECTS_DIR 下，可以直接放置 sample-001.mgz, sample-002.mgz 等文件。下方文件夹结构为分割完成后的文件夹结构：

SUBJECTS_DIR
├── sample-001
│   ├── label
│   ├── mri
│   │   ├── orig
│   │   └── transforms
│   │       └── bak
│   ├── scripts
│   ├── stats
│   ├── surf
│   ├── tmp
│   ├── touch
│   └── trash
└── sample-002
    ├──...

recon-all

recon-all 为 Freesurfer 的核心命令，使用命令 recon-all --help 可以查看详细帮助。在操作时可以使用以下命令切换到 SUBJECTS_DIR 下并执行 recon-all，如果需要同时处理多个图像，可以使用后面提到的 Parallel；如果需要处理文件夹形式的 DICOM 图像，可先用 dcm2niix 转换成 nii 格式。

cd $SUBJECTS_DIR
# bert 是可自己设置的 subject 名字
recon-all -s bert -i sample-001.mgz -all

上面的命令成功执行后会持续数个小时，处理结果会分类存放在上一部分列出的文件夹结构中。

freeview

在终端中输入 freeview 进入 GUI，然后通过 File - Load Volume 加载 mri -> T1.mgz / aseg.mgz / wm.mgz / brainmask.mgz，aseg.mgz是灰质白质等解构的分割结果，可以选择 Color map 为 Lookup Table 并调低透明度进行观察。

recon-all的结果中还包含了曲面信息，使用 File - Load Surface 加载 surf -> lh.pial / lh.white / lh.inflated 以及对应的右脑文件 surf -> rh.pial / rh.white / rh.inflated。

下面的命令行操作与上面的结果完全相同，无需重复操作：

$> freeview -v \
    bert/mri/T1.mgz \
    bert/mri/wm.mgz \
    bert/mri/brainmask.mgz \
    bert/mri/aseg.mgz:colormap=lut:opacity=0.2 \
    -f \
    bert/surf/lh.white:edgecolor=blue \
    bert/surf/lh.pial:edgecolor=red \
    bert/surf/rh.white:edgecolor=blue \
    bert/surf/rh.pial:edgecolor=red

stats

Freesurfer的 recon-all 命令已经处理的MR数据的体积以及皮层厚度信息，使用 asegstats2table, aparcstats2table 两个命令可以分别提取出。体积与皮层厚度结果均是按照已经分割好的脑区进行统计，皮层厚度只能对左脑和右脑分别输出。另外，这两个命令仅能输出txt格式文档，查看时需要使用excel进行转换。使用命令如下：

# https://surfer.nmr.mgh.harvard.edu/fswiki/asegstats2table
asegstats2table --subjects bert --common-segs --meas volume --stats=aseg.stats --table=segststs.txt

# https://surfer.nmr.mgh.harvard.edu/fswiki/aparcstats2table
# aparcstats2table 需要分别提取左右脑的结果
aparcstats2table --subjects bert --hemi lh --meas thickness --parc=aparc --tablefile=aparc.txt
aparcstats2table --subjects bert --hemi rh --meas thickness --parc=aparc --tablefile=aparc.txt

Parallel

Freesurfer 本身无法同时对多个图像（多个 subjects ）进行计算，为了充分利用电脑的多核性能节省计算时间，需借助一个多线程工具，名为 GNU Parallel。Parallel 在 Ubuntu下可以直接使用 sudo apt-get install parallel 进行安装。多线程进行 recon-all 的命令如下，使用此方法会将每个图像的文件名作为分割该图像时的subject名字。执行此命令时可能没有信息在终端输出，执行完成后会一次性输出，如需查看进行状态，可以到 ./scripts/recon-all.log 查看 recon-all 日志.

# mgz 为后缀，需要指定参与 CPU 参与计算的核心数则在 parallel 后加：--jobs 核心数
ls *.mgz | parallel recon-all –s {.} –i {} –all

# 统计数据提取也可以使用 Parallel
ls *.nii | parallel asegstats2table --subjects {.} --common-segs --meas volume --stats=aseg.stats --table=segststs_{.}.txt
ls *.nii | parallel aparcstats2table --subjects {.} --hemi lh --meas thickness --parc=aparc --tablefile=aparc_{.}_lh.txt
ls *.nii | parallel aparcstats2table --subjects {.} --hemi rh --meas thickness --parc=aparc --tablefile=aparc_{.}_rh.txt

对于上述 Parallel 导出的一系列统计txt文件，我编写了 python 程序转换成更易读的 csv 文件，将下面的文件保存为 txt2csv.py 然后放在当前的 SUBJECTS_DIR 下，在终端中执行 python txt2csv.py 即可进行转换。

import csv
import os


def txt2csv(inputfile, outputfile):
    datacsv = open(outputfile, 'w')
    csvwriter = csv.writer(datacsv, dialect=("excel"))
    mainfileH = open(inputfile, 'rb')
    for line in mainfileH.readlines():
        # print "Debug: " + line.replace('\n', '')
        csvwriter.writerow([a for a in line.decode().split('\t')])
    datacsv.close()
    mainfileH.close()


if __name__ == "__main__":
    txt_list = [
        x for x in os.listdir('.')
        if os.path.isfile(x) and os.path.splitext(x)[1] == '.txt'
    ]
    for txt in txt_list:
        outputfile = txt.split('txt')[0] + 'csv'
        txt2csv(txt, outputfile)

dcm2niix

recon-all 命令无法对文件夹形式的 DICOM 图像直接进行处理，借用 dcm2niix 库以及 parallel 可以进行批量转换。dcm2niix 默认将转换结果放到原 DICOM 文件夹内，不方便后续处理，所以用设置 -o 参数指定输出文件夹。

# 安装
sudo apt-get install pigz dcm2niix

# 使用
# 将多个 DICOM 文件夹放到同一目录下，并避免有其他非 DICOM 文件
ls | parallel dcm2niix -o 目标输出文件夹 {}

PETSurfer

PETSurfer 是 Freesurfer 6 中引入的新功能，包含部分容积校正以及力学建模等功能。

  1. 创建一个分割

    # 替换 SUBJECT
gtmseg --s SUBJECT

    此步操作会使用生成一个高精度的分割结果 SUBJECT/mri/gtmseg.mgz ，使用到原来的 ageg.mgz 提供皮层下结构，使用 ?h.aparc.annot 提供皮层结构，还将估计一些额外的脑结构。

  2. PET 配准到 MR

    # 替换 PET_IMAGE, template名字也可以修改
# PET_IMAGE 为PET路径/PET文件名
mri_coreg --s subject --mov PET_IMAGE --reg template.reg.lta

    检查配准结果使用以下命令：

    tkregisterfv --mov PET_IMAGE --reg template.reg.lta --surfs

  3. 部分容积校正

    # 其中 PET_IMAGE FWHM 均需要替换
# PET_IMAGE 为PET路径/PET文件名
# 如果不需要参考脑桥代谢值，可以再命令中补充 --no-rescale
# FWHM 取决于设备，医大一院用的2.8
mri_gtmpvc --i PET_IMAGE --reg template.reg.lta --psf 2.8 --seg gtmseg.mgz --default-seg-merge  --auto-mask PSF .01 --mgx .01 --o gtmpvc.output 

   
官网对各参数的解释如下：
   
   >--psf FWHM is the full-width/half-max of the the point-spread function (PSF) of the scanner as measured in image space (also known as the burring function). The blurring function depends on the scanner and reconstruction method; here it is modeled as an isotropic Gaussian filter of the given FWHM. Eg, an HR+ is typically about 6mm. This parameter has nothing to do with applying smoothing. Set this to 0 to turn off PVC. If you don't know what to set this to, ask your PET physicist.
   >--seg gtmseg.mgz is the segmentation created with gtmseg
   >--default-seg-merge will merge several segmentations, eg, all the ventricular CSF segments will be merged into one ROI
   >--auto-mask FWHM .01 will create a mask to exclude voxels from the analysis (saves time). Output volumes will be reduced to a bounding box around this mask (saves space)
   >--mgx .01 Run Muller-Gartner analysis. 01 is the GM threshold. Only necessary if you want to do a voxel-wise analysis.
   >--o gtmpvc.output This is the output folder.

此步生成了一个文件 ./gtmpvc.output/gtm.stats ，里面存放了各个 ROI 的代谢信息，阅读参考如下内容：

9 17 Left-Hippocampus subcort_gm 473 174.083 1.406 0.1216

9 = ninth row 17 = index for ROI Left-Hippocampus = name of ROI subcort_gm = tissue class 473 = number of PET voxels in the ROI 174 = variance reduction factor for ROI (based on GLM/SGTM) 1.406 = PVC uptake of ROI relative to Pons 0.1216 = resdiual varaince across voxels in the ROI

引入PET后的批处理

简单的PET批处理方法需要从转DICOM开始进行操作。首先切换到在目标目录打开终端并设置为SUBJECTS_DIR ，并建立两个存放DICOM的文件夹，然后分别讲MR和PET图像存进去：

export SUBJECTS_DIR=$(pwd)
# 分别为PET图像和MR图像建立一个名为PET和MR的文件夹，分隔开
mkdir MR
mkdir PET

下面开始转换格式，转换的nii使用人名命名，如果有空格或者重复人名需要手动修正。

# 切换到MR文件夹，输出文件直接放到 SUBJECTS_DIR
cd MR
ls | parallel dcm2niix -o $SUBJECTS_DIR -f %n -b n {}

# 切换到PET文件夹
cd ../PET
ls | parallel dcm2niix -o $(pwd) -f pet_%n -b n {}

# 返回 SUBJECTS_DIR
cd ..

然后对MR进行分割

#此时目录应该是SUBJECTS_DIR
ls *.nii | parallel recon-all –s {.} –i {} –all

然后处理PET，各命令含义见上一节：

ls *.nii | parallel --jobs 4 gtmseg --s {.}

ls *.nii | parallel mri_coreg --s {.} --mov PET/pet_{} --reg {.}.reg.lta
# 配准结果检查不使用批处理

# 如果不需要参考脑桥代谢值，可以再命令中补充 --no-rescale
# FWHM 需要指定一个值
ls *.nii | parallel --jobs 4 mri_gtmpvc --i PET/pet_{} --reg {.}.reg.lta --psf FWHM --seg gtmseg.mgz --default-seg-merge  --auto-mask PSF .01 --mgx .01 --o {.}_gtmpvc.output

提取recon-all的统计结果

mkdir asegStats
ls *.nii | parallel asegstats2table --subjects {.} --common-segs --meas volume --stats=aseg.stats --table=./asegStats/{.}_segststs.txt

# https://surfer.nmr.mgh.harvard.edu/fswiki/aparcstats2table
# aparcstats2table 需要分别提取左右脑的结果
mkdir aparcStats
ls *.nii | parallel aparcstats2table --subjects {.} --hemi lh --meas thickness --parc=aparc --tablefile=./aparcStats/{.}_aparc_lh.txt
ls *.nii | parallel aparcstats2table --subjects {.} --hemi rh --meas thickness --parc=aparc --tablefile=./aparcStats/{.}_aparc_rh.txt

提取 gtmseg 及其统计结果

mkdir gtmseg
ls *.nii | parallel cp ./{.}/mri/gtmseg.mgz ./gtmseg/{.}_gtmseg.mgz

mkdir gtmVolume
ls *.nii | parallel asegstats2table --subjects {.} --common-segs --meas volume --stats=gtmseg.stats --table=./gtmVolume/{.}_gtmVolume.txt

提取PET处理结果 gtm.stats.dat

mkdir gtmStats
ls *.nii | parallel cp ./{.}_gtmpvc.output/gtm.stats.dat ./gtmStats/{.}.gtm.stats.dat

海马子区域分割

进行此步骤的前提是正确安装了matlab runtime，然后只需要执行 recon-all 命令

如果是一个新的待分割数据，需要重新开始 recon-all 过程

 # bert可以自定义
 recon-all -all -s bert -hippocampal-subfields-T1

如果之前进行过 recon-all 分割，只需要执行

# bert为原来的subject名称
recon-all -s bert -hippocampal-subfields-T1

分割结果存放在文件夹 $SUBJECTS_DIR/bert/mri/中，可以使用freeview查看分割结果

cd ./bert/mri
freeview -v nu.mgz -v lh.hippoSfLabels-T1.v10.mgz:colormap=lut -v rh.hippoSfLabels-T1.v10.mgz:colormap=lut

有关分割结果的介绍，fs文档原文如下

The output will consist of six files (three for each hemisphere) that can be found under the subject’s mri directory (in this example, $SUBJECTS_DIR/bert/mri/):

[lr]h.hippoSfLabels-T1.v10.mgz: they store the discrete segmentation volume (lh for the left hemisphere, rh for the right) at 0.333 mm resolution.

[lr]h.hippoSfLabels-T1.v10.FSvoxelSpace.mgz: they store the discrete segmentation volume in the FreeSurfer voxel space (normally 1mm isotropic, unless higher resolution data was used in recon-all with the flag -cm).

[lr]h.hippoSfVolumes-T1.v10.txt: these text files store the estimated volumes of the subfields and of the whole hippocampi.

Note that [lr]h.hippoSfLabels-T1.v10.mgz covers only a patch around the hippocampus, at a higher resolution than the input image. The segmentation and the image are defined in the same physical coordinates, so you can visualize them simultaneously with (run from the subject’s mri directory):

On the other hand [lr]h.hippoSfLabels-T1.v10.FSvoxelSpace.mgz lives in the same voxel space as the other FreeSurfer volumes (e.g., orig.mgz, nu.mgz, aseg.mgz), so you can use it directly to produce masks for further analyses, but its resolution is lower than that of [lr]h.hippoSfLabels-T1.v10.mgz.