使用FSL对DTI进行概率纤维追踪（GPU版本）

  FSL是牛津大学开发的十分好用的脑核磁图像处理工具，其中可以利用bedpostx和protrackx对DTI图像进行概率纤维追踪（probabilistic tractography）。但这一过程（尤其是bedpostx），在FSL给出的官方文档中被说明每张图片都需要消耗大约15小时的运算时间，好在FSL也提供了bedpostx和probtrackx2的GPU版本，稍作配置之后便可以通过命令行形式使用，效果拔群。

  FSL官方文档上仅介绍了概率纤维追踪的GUI方法，命令行教程并不清晰，而GUI是无法使用GPU加速的，故本文同时给出了使用部分命令行方法的详细步骤。

配置信息

  OS:Ubuntu 16.04 LTS
  GPU:NVIDIA Tesla V100-PCIE-32GB
  CUDA Runtime version: 9.2
  Cudnn version: 7.4.1
  FSL:6.0

  要注意一下，根据这里，使用probtrackx2_gpu要求你的GPU至少有5GB的显存，否则在运行过程中容易出现out of memory。

配置bedpostx_gpu和probtrackx2_gpu

  在官方网站上下载：bedpostx_gpu和probtrackx2_gpu，注意根据你的CUDA Runtime version选择对应版本。若不支持你的CUDA版本，则需要重新安装CUDA和Cudnn。

• bedpostx_gpu传送门

• probtrackx_gpu传送门

  下载完成后解压缩

• bedpostx_gpu解压缩后有bin和lib两个文件夹，将bin中所有文件cp到$FSLDIR/bin下，将lib文件夹中所有文件cp到$FSLDIR/lib下

• probtrackx_gpu解压缩后仅有一个二进制文件，将其cp到$FSLDIR/bin下即可。

  文件复制完成后，需要将bedpostx_postproc_gpu.sh中第一行的#!/bin/sh改为#!/bin/bash，否则会出现bug导致第13行的if分支无法进入。

bedpostx

  bedpostx采用马尔可夫链蒙特卡罗采样来建立每个体素的扩散参数分布，最终生成概率纤维追踪（probabilistic tractography）所需要的所有参数文件。

数据准备

  使用bedpostx_gpu需要确定你的执行目录下有：

• data: 4D数据文件(.nii或.nii.gz)，由原始数据预处理后得到。

• nodif_brain_mask: 3D二进制掩码（.nii或.nii.gz），对结构像bet后得到。

• bvecs: 3xN的ASCII文本文件，代表DTI图片采集时的磁场参数之一。

• bvals: 1xN的ASCII文本文件，代表DTI图片采集时的磁场参数之一。

  这是我的存放数据的位置，其余文件是我做了DTIFIT的结果。

数据检查

  在执行bedpostx_gpu前要对上述的数据文件进行格式检查，各数据需要满足的格式要求为：

• mask和data文件除mask的第四维为1以外，其它维度大小应对应相等。

• bvecs与bvals是DTI图像采集过程中的磁场参数，每行的字数N需与data文件的第四个维度相同。

  可在命令行中cd到执行文件夹后运行bedpostx_datacheck这一shell脚本来检查各数据文件是否存在且格式是否匹配。

  注意：在笔者的Ubuntu 16.04下运行bedpost_datacheck时，出现了部分语法错误，后经查询后发现是一些符号问题，只需对应更改脚本的103、127、131行，并同样将#!/bin/sh改为#!/bin/bash即可。

执行

  将所有所需文件放好后（我是存放在~/Documents/Probtrack/data/Data，后续以我为例），只需cd到该目录，输入bedpostx_gpu .(点别忘了)即可。

  运行结果如图所示，输出目录为~/Documents/Probtrack/Data.bedpostX。

  输出目录中，以merged开头的文件将是后续概率纤维追踪所需要的，具体含义参考FSL——DTI官方文档。

  相比起官方文档中提到的CPU版本所需要的15小时，GPU版本的bedpostx仅耗时10分钟，可以说运算速度确实是大大提高了。

概率纤维追踪

  FSL使用probtrackx2进行概率纤维追踪，其具体原理参考官方文档。

数据准备

  probtrackx2_gpu需要如下数据：

• bedpostx生成的merged。

• nodif_brain_mask，bedpostx后会自动将其复制到Data.bedpostX文件夹。

• 根据追踪类型不同所需要的seed、target、waypoint

• 根据seed种类和seed space不同，需要的.mat或_wary.nii.gz文件（这一点很重要，稍后会解释）。

执行

  说明执行步骤之前，针对上一部分所重点提到的seed种类和seed space做一点解释。
  概率纤维追踪的最终结果简单来说可显示为由seed开始，经过waypoint，到termnition的纤维通路，其中seed是必须指定的。而seed可以是一个体素（single voxel），也可以是一个或多个mask。
  probtrackx支持任意空间的输入输出（Diffusion、Standard或自定义），但是其过程只能在Diffusion sapce下进行，因此，在使用mask作为seed的情况下若seed space不是diff空间下，则需要提供将其线性配准到diff空间的变换矩阵.mat文件，或非线性配准到diff的正向和反向waryfiles。
  具体步骤我们以这里提供的数据为例。
  文件解压后，其中subj1_model2_3fibres.bedpostX文件夹为bedpostx后生成的文件夹。

Single Voxel

  Single模式下，追踪计算量较小，因此probtrackx2_gpu没有无法使用–simple参数选项。若需要命令行模式运行，则使用probtrackx2即可，命令如下：

probtrackx2 --simple --seedref=nodif_brain_mask -o sigle_voxel_result -x coordinates.txt -l --onewaycondition --forcedir --opd -s merged -m nodif_brain_mask --dir=single_voxel_result

  其中需要说明的是-x参数后的.txt文件是我自己创建的，内容为68 55 31，代表voxel的坐标。

  用fsleyes查看追踪结果：

Single Mask

  现在我们使用subj1_model2_3fibres.bedpostX中丘脑外侧膝状核掩码（SEED_LGN_LEFT.nii.gz）作为seed。

  注意，该mask是在标准空间内，因此需要提供将其配准到diff空间的变换文件，此例中，已提供所需的.mat或waryfiles在xfms子目录下。

  命令如下（使用非线性配准到diff）：

probtrackx2_gpu -x SEED_LGN_LEFT.nii.gz -l --onewaycondition  --xfm=xfms/standard2diff_warp.nii.gz --invxfm=xfms/diff2standard_warp.nii.gz --forcedir --opd -s merged -m nodif_brain_mask --dir=single_mask_result

  可以看出，GPU版本下，除过13s的参数加载过程，只需2s就可完成追踪。

  用fsleyes查看追踪结果：

  可以看出与LGN有连通的纤维被追踪到。

Including Target

  probtrackx_gpu可使用额外的参数设置纤维追踪的途经点、终点等，以此来进行更具针对性的追踪，这些高级追踪模式包括：Waypoint mask、Exclusion mask、Termination mask、Classification targets。本文只简单介绍Waypoint和Termination，详细介绍请读者在官方文档中自行查阅。

  使用probtrackx_gpu的--waypoints和--stop参数可设定纤维追踪的途经点和终点，用这种方法可以更加针对性地开展纤维追踪工作。

  此例中，我们使用LGN作为seed，初级视皮层（primary visual cortex）掩码TARGET_V1_LEFT.nii.gz同时作为waypiont和termination，以此只追踪LGN到V1的纤维通路。

  命令如下：

probtrackx2_gpu -x SEED_LGN_LEFT.nii.gz -l --onewaycondition --xfm=xfms/standard2diff_warp.nii.gz --invxfm=xfms/diff2standard_warp.nii.gz --forcedir --opd -s merged -m nodif_brain_mask --dir=waypoint_result --stop=TARGET_V1_LEFT.nii.gz --waypoints=TARGET_V1_LEFT.nii.gz --waycond=AND 

  使用fsleyes查看结果：

  可明显看出，与不指定waypoint和termination相比，这种模式下可以更具针对性地追踪到seed到ROI的纤维通路。

后记

  从接手服务器到配置环境再到尝试用FSL进行概率纤维追踪，整个过程还是踩了不少坑，好在官方文档非常详细，网络资源也很丰富，总的来说主要还是有这样额外几点总结：

• 开展工作前尽量认真地阅读官方文档。一开始我使用protrackx2_gpu时，总是报非法内存访问和out of memory（32G显存你告诉我OTM？？？），换了三个CUDA版本后还是不好使，后来仔细阅读了一下官方文档后发现原来是我的数据不在diff空间，我也没指定xfm和invxfm。

• 不明白命令行参数怎么设置的话，使用probtrackx2_gpu --help可以查看参数解释，实在不行就跟着官方的GUI教程跑一下CPU，直接kill掉后在输出文件夹里有log文件可以查看命令。

• 若要使用自己的seed mask、waypoint mask等，可以在脑网络组图谱中下载相应的图谱后使用MATLAB spm12工具得到特定区域的mask，只是还是要注意配准的问题。

  由于笔者也是最近一段时间开始DTI图像处理相关的学习工作，故本文存在疏漏在所难免，还望大家不吝赐教，若有部分内容需要指教或讨论，请邮箱联系matrixliu@csu.edu.cn。若有转载需要，请注明出处即可。

Reference

  FSL Wiki
  Example Box: Tractography
  Bedpostx GPU
  Probtrackx GPU