构建fODF群组模板，改进结果比较

构建fODF是后续分析白质信息的重要一步，目前在MRtrix3 工具包中，提供了population_template 工具实现其功能。

population_template的代码非常长，而且缺少说明文档。如果要改进，十分不便。通过反复阅读代码，整理代码思路，以备使用。

构造初始化模板

Line 940 'Generating initial template'

1. 利用mraverageheader构造一个标准的空间，确定一个模板grid size。

在code中，把这个标准空间称之为average space header

mraverageheader -fill [input_c0/*.mif (9 items)] - | mrgrid - regrid -voxel 1.5,1.5,1.5 average_header.mif
shutil.move('average_header.mif', 'average_header_c0.mif')

在Line967 n_contrasts == 1

2. 利用mrtransform 使得所有的fODF具有和初始模板相同的grid size，保存在inp.ims_transformed变量中，但此时只有空间还没有模板。

在Line987中，initial_alignmen Not = ‘none’

mrregister input_c0/1_32.mif average_header_c0.mif -rigid_scale 1  -rigid_niter 0  -type rigid  -rigid_lmax 0  -rigid_init_translation mass  -datatype float32 -rigid linear_transforms_initial/1_32.txt
mrtransform input_c0/1_32.mif -reorient_fod yes  -linear linear_transforms_initial/1_32.txt input_transformed_c0/1_32.mif_translated.mif -datatype float32

生成了配准到初始模版average_header_c0.mif的个体文件，以备平均，生成更新后的模版。

3. 利用aggregate平均所有的inp.ims_transformed模板。

mraverageheader input_transformed_c0/1_32.mif_translated.mif input_transformed_c0/2_33.mif_translated.mif input_transformed_c0/3_42.mif_translated.mif input_transformed_c0/4_68.mif_translated.mif input_transformed_c0/5_48.mif_translated.mif input_transformed_c0/6_33.mif_translated.mif input_transformed_c0/7_59.mif_translated.mif input_transformed_c0/8_46.mif_translated.mif input_transformed_c0/9_23.mif_translated.mif average_header_tight.mif
mrgrid average_header_tight.mif pad -uniform 10 - | mrgrid - regrid -voxel 1.5,1.5,1.5 average_header.mif
posix.remove('average_header_tight.mif')

Reclice 所得到的fODF数据，使它们和average_header.mif具有相同的grid size。

mrtransform  -reorient_fod yes input_transformed_c0/1_32.mif_translated.mif input_transformed_c0/1_32.mif  -interp linear -template average_header.mif -nan -datatype float32
posix.remove('input_transformed_c0/1_32.mif_translated.mif')

后面跟着多个python函数，包括了inplace_nan_mask、calculate_isfinite、和aggregate等。inplace_nan_mask没有被执行，而calculate_isfinite实际上执行等是下面的两段代码，
nanmask_input: False ' ' leave_one_out: True

* sub-num
mrconvert input_transformed_c0/1_32.mif -coord 3 0 - | mrcalc - -finite isfinite_c0/1_32.mif

mrmath [isfinite_c0/*.mif (9 items)] sum isfinite_c0.mif

而aggregate执行的是下面的codes，生成最终的template。

mrmath [input_transformed_c0/*.mif (9 items)] mean -keep_unary_axes initial_template_c0.mif

线性配准过程

在codes中说的是优化过程，会有一个总数为niter，的循环优化过程。主要分为两大步骤，第一步骤是将各个subject配准到模板上；第二个步骤是将所得到的数据给平均起来，生成新一轮的模板。

Line 1149 'Optimising template with linear registration (stage {0} of {1}; {2})'

步骤一、个体配准

1. 根据预先设定的配准配置信息，如rigid/affine， scale，lmax等，生成配准方法的参数；

2. 根据isfinite_count的信息，用mrcalc修正template？；

接上面的calculate_isfinite的codes，

mrmath [isfinite_c0/*.mif (9 items)] sum isfinite_c0.mif

它把所有的模版区域加在了一起，这样就可以定义当前odf的区域，用来做LOO。

mrcalc isfinite_c0.mif isfinite_c0/1_32.mif -sub - | mrcalc initial_template_c0.mif isfinite_c0.mif -mult input_transformed_c0/1_32.mif 1 -mult  -sub - -div loo_initial_template_c0.mif

上面这一步不是很懂是在为什么要这么干，感觉像是在对average的odf进行某种re-weighting操作，把当前的odf剔除出去。Not sure？

mrregister input_c0/1_32.mif loo_initial_template_c0.mif -rigid_scale 0.3 -rigid_niter 100 -rigid_lmax 2 -type rigid -datatype float32 -rigid linear_transforms_00/1_32.txt

3. check_linear_transformation。

拆分transformation matrix 成为三个分量，translation，rotaion and stretch，和 shear。分别保存在scaling、shear、angles、angle_axis、translation、R、S的key中。

transformcalc linear_transforms_00/1_32.txt decompose linear_transforms_00/1_32.txtdecomp
posix.remove('linear_transforms_00/1_32.txtdecomp')
posix.remove('loo_initial_template_c0.mif')

如果出现错误，程序将提示，并终止执行。这里检查了上面标红的三个keys。生成了linear_transforms_00/1_32.txtdecomp之后，如果没有问题，程序继续执行，删掉生成的文件。

https://github.com/MRtrix3/mrtrix3/tree/master/lib/mrtrix3
matrix.dot Dot product

Here we ensure the template doesn’t drift or scale
TODO matrix avarage might produce a large FOV for large rotations # pylint: disable=fixme

步骤二、模板生成

1. 使用transformcalc修正变换矩阵；

为了避免模版drift 或者scale

transformcalc linear_transforms_00/1_32.txt linear_transforms_00/2_33.txt linear_transforms_00/3_42.txt linear_transforms_00/4_68.txt linear_transforms_00/5_48.txt linear_transforms_00/6_33.txt linear_transforms_00/7_59.txt linear_transforms_00/8_46.txt linear_transforms_00/9_23.txt average linear_transform_average.txt -quiet
* transformcalc linear_transform_average.txt rigid linear_transform_average.txt -quiet
transformcalc linear_transform_average.txt invert linear_transform_average_inv.txt -quiet

这里的linear_transform_average.txt 用来生成 linear_transform_average_inv.txt文件，然后用matrix的python接口进行运算，更新transformation的值。But Why does this work？ seems not work even make it worse.

-----------------------round 1-------------------------------------- round3---------------

这里自己用shell写了一个matrix.dot的矩阵点乘函数，代替python中的代码。然后通过计算发现，到了第三轮Rigid配准，模板就飞了。于是，删掉这样的操作，模板似乎质量更好了。不过，在目前的实验中，并没有采用LOO的方式，而是把所有的数据都用上，以生成模板。

mrtransform  -reorient_fod yes input_c0/1_32.mif -template initial_template_c0.mif -linear linear_transforms_00/1_32.txt input_transformed_c0/1_32.mif -nan -datatype float32

更新之后的文件，用来做进一步的线性配准，执行如上一节类似的计算。具体操作如下，
2. calculate_isfinite 计算finite影像
3. aggregate平均生成模板。

最后一个迭代流程中，为非线性配准生成初始模版。

mrmath [isfinite_c0/*.mif (9 items)] sum isfinite_c0.mif
mrmath [input_transformed_c0/*.mif (9 items)] mean -keep_unary_axes nl_template00_c0.mif

非线性配准过程

执行多轮的线性配准过后（default n=10），开始非线性配准。

Line 1287 'Optimising template with non-linear registration (stage {0} of {1})'

经过短短几行的settings，从Line 1316 开始，首先如linear部分，定义当前odf的区域，用来做LOO。根据如上Linear部分生成的odf模版， nl_template00_c0.mif，开始非线性配准。

mrcalc isfinite_c0.mif isfinite_c0/1_32.mif -sub - | mrcalc nl_template00_c0.mif isfinite_c0.mif -mult input_transformed_c0/1_32.mif 1 -mult  -sub - -div loo_nl_template00_c0.mif

然后开始第一轮的非线性配准。

mrregister input_c0/1_32.mif loo_nl_template00_c0.mif -type nonlinear -nl_niter 5 -nl_warp_full warps_01/1_32.mif -transformed input_transformed_c0/1_32.mif  -nl_update_smooth 2.0 -nl_disp_smooth 1.0 -nl_grad_step 0.5 -nl_init warps_00/1_32.mif -datatype float32 -nan -nl_lmax 2
posix.remove('loo_nl_template00_c0.mif')
posix.remove('warps_00/1_32.mif')

根据配准生成的形变场，把所有subject的数据mapping到标准空间，也就是-transformed input_transformed_c0/1_32.mif 跟的文件名 。然后更新模版。

mrmath [input_transformed_c0/*.mif (9 items)] mean -keep_unary_axes nl_template01_c0.mif

紧接着，会对变形场进行进一步缩放的操作。

mrgrid warps_01/3_42.mif regrid -scale 1.250000 tmp.mif
shutil.move('tmp.mif', 'warps_01/3_42.mif')

然后就行进下一轮的配准

mrcalc isfinite_c0.mif isfinite_c0/1_32.mif -sub - | mrcalc nl_template01_c0.mif isfinite_c0.mif -mult input_transformed_c0/1_32.mif 1 -mult  -sub - -div loo_nl_template01_c0.mif
mrregister input_c0/1_32.mif loo_nl_template01_c0.mif -type nonlinear -nl_niter 5 -nl_warp_full warps_02/1_32.mif -transformed input_transformed_c0/1_32.mif  -nl_update_smooth 2.0 -nl_disp_smooth 1.0 -nl_grad_step 0.5 -nl_init warps_01/1_32.mif -datatype float32 -nan -nl_lmax 2
posix.remove('loo_nl_template01_c0.mif')
posix.remove('warps_01/1_32.mif')

缩放六轮之后，也不再缩放。

第一轮配准的结果

配准后结果

改进后atlas 比较

找了最接近的图片，哪一个效果更好呢？