MRI图像神经纤维束的追踪

大脑的内部有灰质和白质，白质的主要成分就是神经，水分子的弥散方向倾向于沿着神经细胞轴突的走向，也就是沿着神经的方向，通过测定的弥散加权像(DWI)就可以推断出大脑内部白质神经纤维束的走向

 需要用到的工具有：MRtrix3、FSL、freesurfer(可选)、ANTs，这里MRtrix3有图形界面mrview，用于查看nii.gz文件

数据来源openneuro，这里我们下载了三组弥散项数据以及一组结构项数据，包括第一个被试sub-cIIIA01的第一个session ses-s1Bx1的T1w像，和dwi下面的run103/104/311的所有文件,包括：.nii.gz/.bval/.bvec/.json

b值为弥散梯度，b值的增加，弥散权重会增加，会提升不同组织之间的对比度，但会影响图像的信噪比，目前倾向于收集多个b值的数据进行综合分析。我们先简单看一下这些数据的信息。

(zf) hzhb@hzhb:~/Fengz/demo$ mrinfo *311*.nii.gz
************************************************
Image name:          "sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n40r21x21x22peAPP_run-311_dwi.nii.gz"
************************************************
  Dimensions:        112 x 112 x 54 x 47
  Voxel size:        2.14286 x 2.14286 x 2.2 x 2.9
  Data strides:      [ -1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         signed 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 416.72052001953125
  Transform:               0.9948    -0.08041    -0.06207      -102.7
                          0.07542       0.994    -0.07894        -116
                          0.06804     0.07385      0.9949      -28.04
  comments:          TE=79;Time=154737.100

 从这里我们可以得到文件311包含47个全脑，以及体素的大小，另外三个是文本文件，bval记录了每一个全脑对应的弥散梯度。文件名中的b1000n40代表有40个弥散梯度为1000的全脑

(zf) hzhb@hzhb:~/Fengz/demo$ cat *311*.bval
0.1 1000 1000 1000 1000 1000 1000 1000 0.1 1000 1000 1000 1000 1000 1000 1000 1000 0.1 1000 1000 1000 1000 1000 1000 1000 1000 0.1 1000 1000 1000 1000 1000 1000 1000 1000 0.1 1000 1000 1000 1000 1000 1000 1000 1000 0.1 1000 0

bvec记录了弥散梯度的方向，一共3行47列，每列分别记录了47个全脑对应的弥散梯度方向

(zf) hzhb@hzhb:~/Fengz/demo$ cat *311*.bvec
0.113788 0.113788 -0.968727 -0.216716 -0.623099 -0.715945 -0.52519 -0.362848 -0.00638426 -0.00638426 0.871859 0.19485 0.767846 0.664229 0.440518 0.461325 0.506283 -0.775606 -0.775606 0.317664 0.553614 0.969743 0.111843 0.0990063 0.903755 -0.409427 0.0486998 0.0486998 -0.832826 -0.306651 0.362747 -0.660762 0.655123 0.828161 -0.18888 0.0873047 0.0873047 0.981366 -0.376661 -0.671461 0.285902 -0.182966 -0.721001 0.631451 -0.810531 -0.810531 0
-0.926411 -0.926411 -0.183337 0.330202 -0.446984 0.675006 0.0759306 -0.850491 0.472956 0.472956 -0.154863 -0.875436 -0.463984 0.741069 0.452956 0.0602404 -0.83168 -0.469402 -0.469402 0.818146 -0.537712 -0.239773 -0.0112147 0.994176 0.317305 0.568948 0.655811 0.655811 0.00446179 -0.540238 -0.923775 0.176767 0.720473 -0.548681 -0.300922 -0.758077 -0.758077 0.0220674 -0.916785 -0.226687 -0.670548 -0.185768 -0.587721 -0.340454 0.45739 0.45739 0
0.358908 0.358908 -0.167201 0.918695 0.641836 0.178296 -0.847591 -0.380797 -0.881063 -0.881063 -0.464628 -0.442317 0.441736 -0.0980633 0.775097 -0.885184 0.228003 -0.422016 -0.422016 -0.479298 -0.635907 0.0459107 0.993663 0.0425783 -0.287306 -0.713209 0.753353 0.753353 -0.553517 0.783651 -0.122701 0.729484 0.227449 0.114445 -0.934757 0.646295 0.646295 -0.190877 0.132784 -0.705517 -0.684563 0.965409 0.367072 0.696678 0.365833 0.365833 0

json文件记录了相位编码方向（从A到P或者从P到A），这个从文件名也可看出，j代表PA方向

(zf) hzhb@hzhb:~/Fengz/demo$ cat *311*.json
{"PhaseEncodingDirection": "j"}

可以用FSLeyes查看nii.gz文件

下面进行数据预处理。

先把数据从nii.gz格式转化为mif格式，只需要一个文件就可以储存四个文件的信息，为方便操作，进行如下操作

(zf) hzhb@hzhb:~/Fengz/demo$ b1000=sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n40r21x21x22peAPP_run-311_dwi
(zf) hzhb@hzhb:~/Fengz/demo$ b2000=sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b2000n56r21x21x22peAPP_run-103_dwi
(zf) hzhb@hzhb:~/Fengz/demo$ b1000AP=sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n3r21x21x22peAPA_run-104_dwi
(zf) hzhb@hzhb:~/Fengz/demo$ T1w=sub-cIIIsA01_ses-s1Bx1_anat_sub-cIIIs01_ses-s1Bx1_acq-r10x10x10_T1w

用mrconvert转换格式，这里没有将json信息放进去

(zf) hzhb@hzhb:~/Fengz/demo$ mrconvert -fslgrad "$b1000".bvec "$b1000".bval "$b1000".nii.gz raw1000.mif -force
mrconvert: [WARNING] existing output files will be overwritten
mrconvert: [100%] uncompressing image "sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n40r21x21x22peAPP_run-311_dwi.nii.gz"
mrconvert: [100%] copying from "sub-cIIIsA...x22peAPP_run-311_dwi.nii.gz" to "raw1000.mif"

得到了mif文件，我们来看它的信息，包含了维度，体素大小，以及弥散梯度大小和方向的信息

(zf) hzhb@hzhb:~/Fengz/demo$ mrinfo raw1000.mif 
************************************************
Image name:          "raw1000.mif"
************************************************
  Dimensions:        112 x 112 x 54 x 47
  Voxel size:        2.14286 x 2.14286 x 2.2 x 2.9
  Data strides:      [ -1 2 3 4 ]
  Format:            MRtrix
  Data type:         signed 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 416.721
  Transform:               0.9948    -0.08041    -0.06207      -102.7
                          0.07542       0.994    -0.07894        -116
                          0.06804     0.07385      0.9949      -28.04
  command_history:   mrconvert -fslgrad sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n40r21x21x22peAPP_run-311_dwi.bvec sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n40r21x21x22peAPP_run-311_dwi.bval sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n40r21x21x22peAPP_run-311_dwi.nii.gz raw1000.mif -force  (version=3.0.3)
  comments:          TE=79;Time=154737.100
  dw_scheme:         -0.06098603046,-0.9577865493,0.2809370571,0.1000000002
  [47 entries]       -0.06098603046,-0.9577865493,0.2809370571,1000.000002
                     ...
                     0.7468558701,0.486905443,0.4529121316,999.9998979
                     -0,0,0,0
  mrtrix_version:    3.0.3

b2000，b1000AP用同样的命令转化即可

(zf) hzhb@hzhb:~/Fengz/demo$ mrconvert -fslgrad "$b2000".bvec "$b2000".bval "$b2000".nii.gz raw2000.mif -force
mrconvert: [WARNING] existing output files will be overwritten
mrconvert: [100%] uncompressing image "sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b2000n56r21x21x22peAPP_run-103_dwi.nii.gz"
mrconvert: [100%] copying from "sub-cIIIsA...x22peAPP_run-103_dwi.nii.gz" to "raw2000.mif"
(zf) hzhb@hzhb:~/Fengz/demo$ mrconvert -fslgrad "$b1000AP".bvec "$b1000AP".bval "$b1000AP".nii.gz raw1000AP.mif -force
mrconvert: [WARNING] existing output files will be overwritten
mrconvert: [100%] uncompressing image "sub-cIIIsA01_ses-s1Bx1_dwi_sub-cIIIs01_ses-s1Bx1_acq-b1000n3r21x21x22peAPA_run-104_dwi.nii.gz"
mrconvert: [100%] copying from "sub-cIIIsA...x22peAPA_run-104_dwi.nii.gz" to "raw1000AP.mif"

下面要做的是把1000和2000对应的数据合并起来，我们用mrcat命令，生成名为dwi的数据

(zf) hzhb@hzhb:~/Fengz/demo$ mrcat raw1000.mif raw2000.mif dwi.mif -force
mrcat: [WARNING] existing output files will be overwritten
mrcat: [100%] concatenating "raw1000.mif"
mrcat: [100%] concatenating "raw2000.mif"

接下来对这个dwi数据进行一系列预处理，首先用dwi2mask创建mask，然后用maskfilter将其变大，生成文件preproce_mask.mif

(zf) hzhb@hzhb:~/Fengz/demo$ dwi2mask dwi.mif - | maskfilter - dilate preproce_mask.mif -npass 3 -force
maskfilter: [WARNING] existing output files will be overwritten
dwi2mask: [100%] preloading data for "dwi.mif"
dwi2mask: [done] computing dwi brain mask
dwi2mask: [done] applying mask cleaning filter
maskfilter: [100%] applying dilate filter to image -

下面进行空间的降噪，物理机制参考文献，noiselevel.mif可以查看降噪降了些什么，添加mask文件限制降噪区域

(zf) hzhb@hzhb:~/Fengz/demo$ dwidenoise dwi.mif denoise.mif -noise noiselevel.mif -mask preproce_mask.mif -force
dwidenoise: [WARNING] existing output files will be overwritten
dwidenoise: [100%] preloading data for "dwi.mif"
dwidenoise: [100%] running MP-PCA denoising

降噪完成后，下面进行吉布斯环伪影去除，参考文献

(zf) hzhb@hzhb:~/Fengz/demo$ mrdegibbs denoise.mif degibbs.mif -force
mrdegibbs: [WARNING] existing output files will be overwritten
mrdegibbs: [100%] performing Gibbs ringing removal

接着进行头动校正和畸变校正，参考文献，需要用到两个编码方向相反的像，先生成一个b0_pair.mif文件

(zf) hzhb@hzhb:~/Fengz/demo$ dwiextract raw1000AP.mif b0_AP.mif -bzero -force
dwiextract: [WARNING] existing output files will be overwritten
dwiextract: [100%] extracting volumes
(zf) hzhb@hzhb:~/Fengz/demo$ dwiextract raw1000.mif - -bzero | mrconvert - -coord 3 0 b0_PA.mif -force
mrconvert: [WARNING] existing output files will be overwritten
dwiextract: [100%] extracting volumes
mrconvert: [100%] copying from "/tmp/mrtrix-tmp-O72Lwd.mif" to "b0_PA.mif"
(zf) hzhb@hzhb:~/Fengz/demo$ mrcat b0_AP.mif b0_PA.mif b0_pair.mif -force
mrcat: [WARNING] existing output files will be overwritten
mrcat: [100%] concatenating "b0_AP.mif"
mrcat: [100%] concatenating "b0_PA.mif"

下面正式开始校正，大概需要几十分钟，得到最终的geomcorr.mif(未成功)

(zf) hzhb@hzhb:~/Fengz/demo$ dwifslpreproc degibbs.mif geomcorr.mif -pe_dir PA -rpe_pair -se_epi b0_pair.mif -eddy_options " --data_is_shelled --slm=linear --niter=5 " -force

下面进行bias field correction，能得到更准确的mask(非必要)，这里需要安装ants

dwibiascorrect ants geomcorr.mif biascorr.mif -bias biasfield.mif

下面将最终的弥散像配准到结构像T1w

选取b0的数据配准，因为b0数据配准效果好（由于前两步校正未成功，我们的操作对象仍是degibbs.mif）

(zf) hzhb@hzhb:~/Fengz/demo$ dwiextract degibbs.mif -bzero - | mrmath -axis 3 - mean b0.nii -force
mrmath: [WARNING] existing output files will be overwritten
dwiextract: [100%] extracting volumes
mrmath: [100%] preloading data for "/tmp/mrtrix-tmp-EyAt89.mif"
mrmath: [100%] computing mean along axis 3...

为什么要生成nii格式文件呢？因为配准软件为fsl，之后用如下命令配准，得到b0数据和T1w数据之间的转换矩阵，将其应用到degibbs.mif上，得到最终的align.mif，可以用mrview查看配准结果，实际上每一步预处理结束之后都得看处理质量怎么样，有没有达到预期效果

(zf) hzhb@hzhb:~/Fengz/demo$ flirt.fsl -dof 6 -cost normmi -ref "$T1w".nii.gz -in b0.nii -omat T_fsl.txt
(zf) hzhb@hzhb:~/Fengz/demo$ transformconvert T_fsl.txt b0.nii "$T1w".nii.gz flirt_import T_DWItoT1.txt
(zf) hzhb@hzhb:~/Fengz/demo$ mrtransform -linear T_DWItoT1.txt degibbs.mif align.mif -force
mrtransform: [WARNING] existing output files will be overwritten
mrtransform: [100%] copying from "degibbs.mif" to "align.mif"...

下一步可以获取大脑的分区，运行需要六个小时，要用到freesurfer，得到一个文件夹

(zf) hzhb@hzhb:~/Fengz/demo$ recon-all -i "$T1w".nii.gz -subjid sub-cIIIsA01 -sd . -all

(zf) hzhb@hzhb:~/Fengz/demo$ cd sub-cIIIsA01/
(zf) hzhb@hzhb:~/Fengz/demo/sub-cIIIsA01$ ls
label  mri  scripts  stats  surf  tmp  touch  trash

我们主要用到的是mri中的mgz格式的文件，有很多经典的分区

预处理的最后一步，将大脑分成5种组织，皮层灰质，皮层下灰质，白质，脑脊液，病理组织，我们可以获得白质和灰质的交界，用于提取神经纤维

首先将一个分区文件转换成MRtrix可以读的nii.gz文件

(zf) hzhb@hzhb:~/Fengz/demo/sub-cIIIsA01$ mrconvert mri/aparc.a2009s+aseg.mgz aparc.a2009s+aseg.nii.gz
mrconvert: [100%] uncompressing image "mri/aparc.a2009s+aseg.mgz"
mrconvert: [100%] copying from "mri/aparc.a2009s+aseg.mgz" to "aparc.a2009s+aseg.nii.gz"
mrconvert: [100%] compressing image "aparc.a2009s+aseg.nii.gz"

用5ttgen将上一步的分区文件分为5个区，得到5ttseg.mif分区结果

(zf) hzhb@hzhb:~/Fengz/demo/sub-cIIIsA01$ 5ttgen freesurfer aparc.a2009s+aseg.nii.gz 5ttseg.mif -force
5ttgen: 
5ttgen: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
5ttgen: 
5ttgen: Generated scratch directory: /home/hzhb/Fengz/demo/sub-cIIIsA01/5ttgen-tmp-IHZTDB/
Command:  mrconvert /home/hzhb/Fengz/demo/sub-cIIIsA01/aparc.a2009s+aseg.nii.gz /home/hzhb/Fengz/demo/sub-cIIIsA01/5ttgen-tmp-IHZTDB/input.mif
5ttgen: Changing to scratch directory (/home/hzhb/Fengz/demo/sub-cIIIsA01/5ttgen-tmp-IHZTDB/)
Command:  labelconvert input.mif /home/hzhb/Fengz/AD/freesurfer/freesurfer/FreeSurferColorLUT.txt /home/hzhb/anaconda3/envs/zf/share/mrtrix3/_5ttgen/FreeSurfer2ACT.txt indices.mif
Command:  mrthreshold indices.mif - -abs 0.5 | mrgrid indices.mif crop indices_cropped.mif -mask -
Command:  mrcalc indices_cropped.mif 1 -eq cgm.mif
Command:  mrcalc indices_cropped.mif 2 -eq sgm.mif
Command:  mrcalc indices_cropped.mif 3 -eq  wm.mif
Command:  mrcalc indices_cropped.mif 4 -eq csf.mif
Command:  mrcalc indices_cropped.mif 5 -eq path.mif
Command:  mrcat cgm.mif sgm.mif wm.mif csf.mif path.mif - -axis 3 | mrconvert - result.mif -datatype float32
Command:  mrconvert result.mif /home/hzhb/Fengz/demo/sub-cIIIsA01/5ttseg.mif
Command:  5ttcheck result.mif
5ttgen: Changing back to original directory (/home/hzhb/Fengz/demo/sub-cIIIsA01)
5ttgen: Deleting scratch directory (/home/hzhb/Fengz/demo/sub-cIIIsA01/5ttgen-tmp-IHZTDB/)

获取5ttseg.mif中白质和灰质的交界处，可以用mrview查看生成的5ttseg.mif(5个分区)，5tt_gmwmi.mif(灰质和白质的交界处)

(zf) hzhb@hzhb:~/Fengz/demo/sub-cIIIsA01$ 5tt2gmwmi 5ttseg.mif 5tt_gmwmi.mif -force
5tt2gmwmi: [WARNING] existing output files will be overwritten
5tt2gmwmi: [100%] Generating GMWMI seed mask

至此，预处理阶段结束，下面开始白质纤维束的提取（神经纤维束追踪）

先生成一个弥散张量图像DTI

对于我们刚才配准得到的文件align.mif，首先生成一个mask

(zf) hzhb@hzhb:~/Fengz/demo$ dwi2mask align.mif - | maskfilter - dilate dwi_mask.mif -force
maskfilter: [WARNING] existing output files will be overwritten
dwi2mask: [100%] preloading data for "align.mif"
dwi2mask: [done] computing dwi brain mask
dwi2mask: [done] applying mask cleaning filter
maskfilter: [100%] applying dilate filter to image -

生成弥散张量dt.mif，可用mrview查看，每个体素有六个值(描述一个椭圆体)，标志着水分子的弥散方向，用椭球方向标志纤维的走向，用FA值衡量各向异性，各向异性越大，FA值越大

(zf) hzhb@hzhb:~/Fengz/demo$ dwi2tensor -mask dwi_mask.mif align.mif dt.mif -force
dwi2tensor: [WARNING] existing output files will be overwritten
dwi2tensor: [100%] computing tensors

 还可以用该弥散张量文件生成多种特征图（fa、ad、rd等）

(zf) hzhb@hzhb:~/Fengz/demo$ tensor2metric dt.mif -fa dt_fa.mif -ad dt_ad.mif -force
tensor2metric: [WARNING] existing output files will be overwritten
tensor2metric: [100%] computing metrics

由于弥散像体素比纤维大得多，无法反映神经纤维束的真实走向，这里我们不用dt.mif估算神经的方向，因此使用constrained spherical deconvolution（CSD）算法，计算各个方向纤维束的贡献，这里的kernel称为反应函数

先对灰质、白质、脑脊液分别进行反应函数kernel的估计，得到三个反应函数，并且可以选择输出估计时用到的体素，可以使用shview查看反应函数(3个b值0/1000/2000的都可以看)

(zf) hzhb@hzhb:~/Fengz/demo$ dwi2response msmt_5tt align.mif sub-cIIIsA01/5ttseg.mif ms_5tt_wm.txt ms_5tt_gm.txt ms_5tt_csf.txt -voxels ms_5tt_voxels.mif
dwi2response: 
dwi2response: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
dwi2response: 
dwi2response: Generated scratch directory: /home/hzhb/Fengz/demo/dwi2response-tmp-MWHD2Y/
dwi2response: Importing DWI data (/home/hzhb/Fengz/demo/align.mif)...
Command:  mrconvert /home/hzhb/Fengz/demo/sub-cIIIsA01/5ttseg.mif /home/hzhb/Fengz/demo/dwi2response-tmp-MWHD2Y/5tt.mif
dwi2response: Changing to scratch directory (/home/hzhb/Fengz/demo/dwi2response-tmp-MWHD2Y/)
dwi2response: Computing brain mask (dwi2mask)...
Command:  5ttcheck 5tt.mif
Command:  dwi2tensor dwi.mif - -mask mask.mif | tensor2metric - -fa fa.mif -vector vector.mif
Function: shutil.copy('vector.mif', 'dirs.mif')
Command:  mrtransform 5tt.mif 5tt_regrid.mif -template fa.mif -interp linear
Command:  mrconvert 5tt_regrid.mif - -coord 3 2 -axes 0,1,2 | mrcalc - 0.95 -gt mask.mif -mult wm_mask.mif
Command:  mrconvert 5tt_regrid.mif - -coord 3 0 -axes 0,1,2 | mrcalc - 0.95 -gt fa.mif 0.2 -lt -mult mask.mif -mult gm_mask.mif
Command:  mrconvert 5tt_regrid.mif - -coord 3 3 -axes 0,1,2 | mrcalc - 0.95 -gt fa.mif 0.2 -lt -mult mask.mif -mult csf_mask.mif
dwi2response: Selecting WM single-fibre voxels using 'tournier' algorithm
Command:  dwi2response tournier dwi.mif wm_ss_response.txt -mask wm_mask.mif -voxels wm_sf_mask.mif -scratch /home/hzhb/Fengz/demo/dwi2response-tmp-MWHD2Y/
Command:  amp2response dwi.mif wm_sf_mask.mif dirs.mif wm.txt -shells 0,1000,2000
Command:  amp2response dwi.mif gm_mask.mif dirs.mif gm.txt -shells 0,1000,2000 -isotropic
Command:  amp2response dwi.mif csf_mask.mif dirs.mif csf.txt -shells 0,1000,2000 -isotropic
Function: shutil.copyfile('wm.txt', '/home/hzhb/Fengz/demo/ms_5tt_wm.txt')
Function: shutil.copyfile('gm.txt', '/home/hzhb/Fengz/demo/ms_5tt_gm.txt')
Function: shutil.copyfile('csf.txt', '/home/hzhb/Fengz/demo/ms_5tt_csf.txt')
Command:  mrcat csf_mask.mif gm_mask.mif wm_sf_mask.mif voxels.mif -axis 3
Command:  mrconvert voxels.mif /home/hzhb/Fengz/demo/ms_5tt_voxels.mif
dwi2response: Changing back to original directory (/home/hzhb/Fengz/demo)
dwi2response: Deleting scratch directory (/home/hzhb/Fengz/demo/dwi2response-tmp-MWHD2Y/)

 得到反应函数后，我们估算FOD，该步骤大概需要几分钟，可以用mrview查看得到的文件dwi_gmCsd.mif等，得到每个体素的不同纤维走向的贡献程度(可能性大小)

(zf) hzhb@hzhb:~/Fengz/demo$ dwi2fod msmt_csd align.mif ms_5tt_wm.txt dwi_wmCsd.mif ms_5tt_gm.txt dwi_gmCsd.mif ms_5tt_csf.txt dwi_csfCsd.mif
dwi2fod: [100%] preloading data for "align.mif"
dwi2fod: [100%] performing MSMT CSD (3 shells, 3 tissues)

 下面用FOD进行神经纤维束的追踪，把一个个相连体素的方向画线，这就完成了纤维束的追踪。可以选择全脑的追踪和兴趣区域的追踪。

具体实现如下，我们先做一个全脑的神经纤维追踪，由于是在白质中提取纤维素，我们先生成一个白质的mask文件dwi_wmMask.mif

(zf) hzhb@hzhb:~/Fengz/demo$ mrthreshold -abs 0.2 dt_fa.mif - | mrcalc - dwi_mask.mif -mult dwi_wmMask.mif -force
mrcalc: [WARNING] existing output files will be overwritten
mrcalc: [100%] computing: (/tmp/mrtrix-tmp-nWeFsV.mif * dwi_mask.mif)

设定好算法和参数，开始追踪

(zf) hzhb@hzhb:~/Fengz/demo$ tckgen -algo iFOD2 -act sub-cIIIsA01/5ttseg.mif -backtrack -crop_at_gmwmi -cutoff 0.05 -angle 45 -minlength 20 -maxlength 200 -seed_image dwi_wmMask.mif -select 200k dwi_wmCsd.mif fibs_200k_angle45_maxlen200_act.tck
tckgen: [100%]   367099 seeds,   259465 streamlines,   200000 selected

将得到的追踪像 .tck文件和结构像T1w放在mrview中打开，选择Tractography选项，就得到了对比图如下，可以选择不同的切面

基于兴趣区的追踪与之类似，改变某些参数即可，下图为视觉纤维束

下一步计划：找到感兴趣的纤维束，并且对其中每一条纤维采样一部分点，获取这些点的坐标，然后将这些坐标应用到DTI中找到这些点对应的特征，就得到了患者在一个时间点的数据，