写在前面

可能比较科普而且风格还很意识流，毕竟学着学着我已经感到……数学是我一辈子都追不到的伟大男神。然后Diffusion真的挺复杂的，刚刚起步，还在路上。
这篇科普会分成下面几个部分：

1. DWI信号的基本处理
   -DTI
   -SD&CSD
2. 纤维追踪
   -追得更合理：ACT
   -断舍离：SIFT
3. 纤维的标量化
   -TDI
4. connectome的建立和分析

我比较懒，怕自己鸽，所以一次填坑一点点。会尽力写的，有新的知识也会更改一下前面的:)

1. DWI信号的基本处理

1.1 从DTI说起

其实好像之前写过DTI，思路是用一个tensor模型去定义每个voxel中的各个方向上的diffusion情况。大概长这样：

图1
关于这个tensor怎么得到的后面会再提，现在只需要知道这个tensor是对称的，所以虽然看似有9个值，实际上只有6个。其中Dxx，Dyy，Dzz表示这个voxel内x，y，z方向上的ADC(apparent diffusion coefficient)值，反映了组织自身的扩散能力（在某个方向上），扩散方面的知识可以戳这里；而Dxy，Dxz，Dyz不太容易说清楚，比较类似一种共变性的量，但需要注意不能把它想成是（1，1，0），（1，0，1），（0，1，1）这种向量方向上的ADC值。或许可以用一个柿子来理解。 式1

用这个tensor来表征voxel内各个方向上的扩散有什么意义呢？最直观的作用就是让扩散情况简化了，去除了过多的噪音干扰。用一个图来表示这种变化。

图2
左图把DWI的值（S）直接放在球面极坐标下，这时候我们可以看到整个球面并不光滑，这主要是噪声造成的。中图是根据原始DWI值计算出的ADC，右图是根据Tensor算出的ADC值球面表征，这样就平滑多了。
虽然看起来光滑了，但DTI模型有一些局限性，特别是针对于一个voxel内部有两束甚至多束不同朝向的纤维。
这要涉及到这个tensor的得出，刚才留了一嘴，现在给它讲明白。
式2
式2其实是完整的方程式，结合式1，我们简单的概括出，Tensor中的6个参数的求解过程其实就是通过远大于6个方程（一个磁场梯度方向就能成为一个方程）来求解6个未知数的过程。
简单的说，一个椭球函数包括6个参数，而我们在这一步实际做的就是把DWI的数据整合到一个椭球里面。
说清楚这点我们就基本说清楚DTI的局限性了。刚刚说了这种椭球的性质不利于区分同一个voxel中的两束不同朝向的纤维。就像标准的椭圆中最大的peak方向只有一个，椭球也会体现出这种特性（有两个的时候，这个椭球更像一个飞盘）。

1.2 SD、CSD和MSMT CSD

1.2.1 SD

那么现在目标很明确了，就是要用更多的参数去拟合函数，简化/平滑化我们的DWI信号。
由于DWI信号是一个球面函数，很容易想到用球面谐波函数的卷积去完成这件事情。（才怪，在此之前我听都没听说过什么球面谐波函数）

SD基于一个假设：streamline刨除朝向问题，其实扩散情况都一毛一样。它的逻辑非常简单，一张图就可以概括（图3）。

图3
下面来解释一下这张图：
左边两个“眼镜”表示两束不同朝向的纤维的信号延迟S1分别与其体积占比的乘积，表征了生理基础上，voxel内diffusion信号延迟。
他们的和就是这个voxel内总体的信号情况（假设这个voxel只有两束纤维）。
这可以抽象成两个球面谐波函数的乘积，R(theta)是一束指定了方向纤维的扩散模式（也就是前面说的那个“所有纤维都一样”的扩散模式）。F(theta,phi)就是纤维的朝向信息。

先定性地说一下（喂喂喂你真的有本事定量讲么？）：
在这个过程中我们完成了一次信息的重分配，把纤维束1号、纤维束2号的信息重新分配给了扩散模式函数和纤维朝向函数。
不用说，大家都知道主要想得到的是纤维朝向函数。现在图3中间的函数我们是已知的，要求出F(theta,phi)，要先知道R(theta)。而这一般是靠FA阈值之类的人工标准去选出voxel，这些voxel被认为是装满了朝向一致的纤维。把它们摆正到每个voxel里的纤维都往一个方向，我们就可以算出R(theta)了。
根据这个就可以求出F(theta,phi)。

下面试着定量，这样就可以讲清楚为什么SD能实现之前说的目标了（心虚：
为了说清楚参数什么的，首先要说什么是球谐展开，球谐函数具有正交完备性，可以作为一组正交完备基，展开任意“性质良好”的函数，不知道为什么不妨碍接下来的理解，至于想知道为什么的朋友，请看某位大佬的博客，反正我是看不懂，如果能简单解释一下欢迎戳我。
总之是说：定义在球面上的连续函数可以展开成一组球谐函数的基函数的加权和。

式3
式3所示，柿子右边的f是每一个基函数前面的加权，而Y就是基函数。到了这一步，我们要决定的东西只有：需要展开到什么程度。
展开这个东西有点像本科时候困扰我们的泰勒展开，我们可以决定到底要展开多少层，当然了，你展开的层数越多（项数越多），系数（f的量）就越多，整个就拟合就更精确。
（不过并不是阶数越多越好，这与你原始数据中梯度方向的数量有关，因为如果对噪音也进行了拟合，可能会带来过拟合问题）

项数实际上由阶数l决定，l越大，基函数就越多，因为diffusion信号中朝向信息具有对称性，我们只采用偶次项，所以n阶为最大阶数时，会有(n+1)(n+2)/2个基函数。而在采集R(theta)信息的过程中，只采用0次项。所以在计算R(theta)和F(theta,phi)乘积的时候，为了避免矩阵不齐的问题，我们手动把R(theta)中缺的基函数补上，让同阶非零的基函数权重和其0次项的一致。
这里附上喂鸡百科里的阶数基函数表，可以参考一下。

用一个栗子来结合定性和定量的部分。

图4

这里一共获得了28个参数去解决问题，而且考虑了用多个纤维朝向去拟合，就比DTI更符合生理意义了~

1.2.2 CSD

但是SD还是有一个问题：计算出来的FOD会存在负值。
我们知道FOD求出的权重系数是不同方向上的扩散，这里的负值来自噪音干扰，不存在生理意义，无法解释。
所以进一步地，对求出来的FD权重都施加一个惩罚项。自此，损失函数变成了这样。

图5

1.2.3 MSMT_CSD

从生理上来说，只有白质（WM）中存在水能单向扩散的神经纤维，那么，我们在具体求解的过程中应该尽量避免灰质（GM）和脑脊液（CSF）的影响。因为在这些组织里的扩散情况不太一样。

图6

如图6所示，CSF的流动性更强，所以扩散的幅度相比GM和WM更大；而GM和WM在扩散中的各向异性差异较大，也可以区别开。

区分这些组织的方法有多种，使用分辨率高的T1像当然是其中之一。在我们没有得到T1像时，使用多种b值的DWI图像也可以达到目的。
也就是利用图6中不同组织的扩散峰值随着b值变化的速率（图像斜率）不同，达到区分组织的目的。具体算法如下柿，按照此法，可以算出不同组织的R(theta)：

式4
同理，利用这个R(theta)去继续求解FOD，不过按照划分组织的假设，CSF和GM中的扩散都具有各向同性，其FOD值各维度为1。
其实在求解过程中，唯一的变化就是维度升高了一维。其他部分和基本的CSD没啥差异。

用一张3D一点的图庆贺一下，终于写完了。

撒花

参考文献：
Jeurissen B , Tournier J D , Dhollander T , et al. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data[J]. NeuroImage, 2014, 103:411-426.
Tournier, Jacques-Donald & Calamante, Fernando & G Gadian, David & Connelly, Alan. (2004). Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage. 23. 1176-85. 10.1016/j.neuroimage.2004.07.037.
Van Hecke, Wim & Emsell, Louise & Sunaert, Stefan. (2015). Diffusion Tensor Imaging: A Practical Handbook. 10.1007/978-1-4939-3118-7.