【学习笔记】循环组织感知网络《Recurrent Tissue-Aware Network for Deformable Registration of Infant Brain MR Images》

由于婴儿身体快速发育的原因，大脑MRI图像的精确配准是具有挑战的，于是提出了深度循环神经网络。
该作品有以下三个优点：
1).使用脑组织分割图而不是强度图，以应对婴儿的快速发育；
2).单次配准网络一次配准进行训练，再循环多次，逐步逼近正确的形变场；
3).提出自适应平滑层（adaptive smoothing layer）和组织感知反折叠约束（tissue-aware anti-folding constraint）确保不降低配准准确性的同时符合物理规律；
在实现高水平配准的基础上还确保了形变场的平滑。

1.简介

婴儿脑MRI陪准对于研究婴儿时期脑部动态容量和形态变化具有重要意义。一岁的婴儿闹容量会翻倍增长达到成年人的65%，灰质增长108%-149%，白质增长11%左右，皮质形态快速变化，皮质厚度和表面积明显增加。人们对早期大脑发育的纵向轨迹（longitudinal trajectory）和大脑发育成熟的表征感兴趣。正是由于大脑的快速发育才使得婴儿脑配准充满挑战。
传统的配准方法（SyN，diffeomorphic Demons 和NiftyReg）都是最大化F和M间的相似性，同时约束行变场光滑。大多数方法是和一般的问题，并不是全部都适合婴儿的脑图像配准。就是因为婴儿脑结构在发育。在等强度时期（6～8月大）时，灰质和白质的强度相见，难以对齐组织边界。
也有些学者提出一些方法进行婴儿脑配准如：有人使用纵向轨迹和关键点技术进行配准；还有人使用随机森林法预测初始形变场然后使用卷积进行优化。上述的方法需要大量的步骤和超参数，并且难以适应婴儿脑图像强度可变性。
现在深度学习的配准方法可以达到很好的效果，但是仍然需要在准确率和光滑形变场之间权衡（要求光滑即处处可导，有利于反向传播）。有研究同时比较相似性和L2范数（L2-norm）设置损失函数，损失函数的每一项都要预先设置好。通常要求平滑形变场是牺牲相似性得到的；还有学者使用使用缩放和平方积分（scaling and squaring integration）缓解上述问题，缩放和平方积分的步骤数在推理阶段是不同的，在训练阶段固定了步骤数，为了确保图像的相似性需要牺牲步骤数；还有学者使用了一个鉴别器网络测量平滑而非直接进行平滑约束，但是这种方法涉及手动监督在应用中有所限制。
为了提高准确率，一些学者使用级联网络的方法，虽然提升了准确率，却需要浪费大量的内存和训练时间，因此不为临床所用。
故此，提出了循环组织感知网络（recurrent tissue-aware network ，RTA-Net），使用组织感知正则化训练一个变形配准网络（TA-Net），虽然TA-Net网络是一次行训练的，但是经常在推理阶段使用它更新增量形变场（？），虽然每个增量形变场幅度都很小，最后通过整合所有的增量形变场获得最终形变场。对比级联网络本设计在循环推理前仅训练单个网络。

2.相关介绍

A.传统非深度学习配准

在一般的问题尚且不能说最好，在婴儿脑图像上的更是如此。并且迭代优化耗时耗力，深度学习的方法在节省时间的同时还表现出良好的配准性能。

B.深度一次配准网络

深度神经网络是求一对输入图像间的映射关系，输出变形场。大多数是采用单个网络训练常见的有以下三种方法：
1).监督训练——训练时使用金标准图像，有些使用传统的方法获得金标准图像或使用预定义的形变场生成F/M图像对；
2).无监督训练——损失函数由相似性和形变场正则化共同构成。鉴别网络用来检测相似性和平滑形变场；
3).弱监督训练——损失函数中含有人工标注信息，以提高配准目标区域的准确率。标注信息只在训练时使用。
（本人认为可以分为有监督和无监督，弱监督应该包含于有监督，此外有监督还有：双重监督，深度监督等）
虽然许多网络宣称保持了拓扑结构，但是保持解剖结构的拓扑结构并不容易。有研究表明成人大脑MR图像的变形折叠率（不平滑）会随着速度场积分的离散化步骤的增加而减少。不同数据集的最佳积分步骤数可能不同。在SyN（标准对称归一化）中，当积分步骤数到达4以后会有更多负的雅可比行列式（负的导致折叠）。根据形变场振幅的不确定性采用自适应平滑约束，还提出组织感知解决雅可比行列式，解决折叠产生的伪影。

C.级联配准网络

级联配准网络提升了配准的正确性，一些学者使用不同的网络分别或同时训练（保持之前的网络权重）以提升配准准确率。使用这个方法的前提是有足够好的硬件基础，相反，本文提出的方法只需训练一个网络，然后在推理阶段循环调用，减少网络参数，更容易训练和使用。

D.婴儿脑MR图像分割

婴儿脑图像灰质（GM），白质（WM），脑脊液（CSF）的分割对于基于组织分割的配准映射十分重要，但是在等强度时期进行分割很有挑战。先前的方法多使用多地图集融合的方法。现在致力于基于学习的分割方法，有人使用DNN提升多模态输入MRI的分割结果；还有人使用外观特征和估计组织概率映射迭代完善分割结果。

3.方法

RTA-Net结构如上：(a)是整体结构包括TA网络（获得形变场）和循环推理过程；(b)是训练阶段和推理阶段；
RTA-Net包括TA网络和空间变换器，TA网络通过输入F/M图像对生成平滑形变场。然后空间变换器使用暂时得到的变形场对M图像进行重采样，最后输出一个形变场和变形后的M图像。下面详细介绍：

A.一次训练和循环推理

在训练阶段：只使用TA-Net，并且使用自适应组织平滑和组织感知雅可比行列式协助生成平滑且准确的形变场。
在循环推理阶段：TA-Net与空间变换器都进行循环，每次循环TA-Net都会暂时输出一个增量形变场，然后根据最新的形变场对M图像进行变形。循环推理阶段包括以下三个步骤：
1).将{ $I_m,I_f$}图像对输入预先训练好的TA-Net得到行变场${?}_1$和变形后的M图像$I_m({?}_1)$；
2).将新的图像对{ $I_m({?}_1),I_f$}送入TA-Net得到新的增量形变场${?}_2$和$I_m({?}_2?{?}_1)$；
3).进行n次循环后获得行变场$?={?}_n?...{?}_2?{?}_1$和最终扭曲图像$I_m(?)$。
注意：为了避免重复重采样的误差累积，在每次迭代时对运动图像进行变形之前合成增量形变场。(即在每次循环变形前先计算$?={?}_i?...{?}_2?{?}_1$)

B.TA-Net组织感知变形配准

使用TA-Net获得形变场，结构如上图（b）所示由可训练层、采样层、平滑层和积分层（对速度形变场进行积分求得位移形变场）构成，在循环结构中仅训练一个网络结构即可。训练层结构与3D U-Net类似，由速度场进一步生成形变场$?=Id+{\int }_0^{1}vdt$，其中$Id$表示恒等变换，${\int }_0^{1}vdt$是由速度场计算的位移矢量场。
可训练层获得速度场，可以进行有限步数的离散积分（在整合层进行了七步）以获得形变场$?$。整合层通过缩放和平方操作实现，结构图如上图（b）。实验发现速度场的大小在数量上小于形变场的大小，相比于形变场可训练层可以更准确的预测速度场。
整个TA-Net的损失函数由相似度、速度场正则化参数和形变场正则化参数组成：
$$\mathcal{L}=?Sim(I_m(?),I_f)+Reg(v)+Reg(?),?=Id+{\int }_0^{1}vdt\text{\hspace{1em}(1)}$$
$$????????????重点在下边????????????$$
1.$Reg(v)$详细介绍
建模速度场为多元的高斯场$\mathcal{N}({\mu }_v,{\Sigma }_v)$，每个体素都有平均速度估计和方差。**训练阶段**可训练层获得每个体素的${\mu }_v和{\Sigma }_v$，然后根据获得的数据通过重采样实例化速度场，接着对速度场进行平滑并积分求解形变场。**推理阶段**可以通过忽略${\Sigma }_v$将估计的${\mu }_v$值指定为$v$的实例。（我认为作者经过训练得到的${\mu }_v$已经是比较完善的结果，采样层是用来优化训练层的输出结果的，当训练完成标志着已经到达理想结果）以${\Sigma }_v$编码的速度场的不确定性可以帮助自适应平滑见$Reg(?)$。为平滑速度场定义$Reg(v)$如下：
$$Reg(v)=||▽{\mu }_v|{|}^2+||{\Sigma }_v|{|}^2\text{\hspace{1em}(2)}$$
传统的正则化对于避免婴儿MR脑图像配准伪影是不够的，进一步提出了基于不确定性的自适应平滑和组织感知的雅可比行列式正则化。
2.基于不确定性的自适应平滑
形变场的局部应该平滑以确保邻域空间一致性（邻域空间上是连续的），传统的方法使用固定尺寸的高斯滤波器卷积核实现，定义为：$K_{\sigma }=e^{?\frac{1}{2}||\frac{p?p^{\prime }}{\sigma }|{|}^2}$其中p和p’表示体素位置，p’为p的邻域体素。然后将高斯核用于位移场，当$\sigma$取值较大时抑制折叠伪影的效果越好（越平滑），甚至会造成过平滑而降低配准准确率。自适应平滑机制避免同质（homogeneous），调节向同性避免过平滑。
本设计中平滑层使用自适应高斯核，首先计算位移幅度$|d|$的不确定性${\Sigma }_{|d|}$。原则上使用较大的$\sigma$平滑较大的${\Sigma }_{|d|}$，相反如果较小使用较小的值以保证速度评估的准确性。
不确定性计算式如下：$${\Sigma }_{|d|}={\int }_0^{1}3\sqrt{3}?{\Sigma }_v^{max}dt\text{\hspace{1em}(3)}$$
其中，${\Sigma }_v^{max}=max（{\Sigma }_{v,x},{\Sigma }_{v,y},{\Sigma }_{v,z}）$表明最大速度不确定组成，积分层由缩放和平方操作完成。
基于不确定性的自适应平滑
基于不确定性将体素分为三组，每组分别使用相应的核，由小到大依次定义三个核为${\sigma }_1=\frac{max({\Sigma }_{|d|})+2?min(\Sigma |d|)}{3}$、${\sigma }_2=\frac{2?max({\Sigma }_{|d|})+min(\Sigma |d|)}{3}$、${\sigma }_3=max({\Sigma }_{|d|})$，见下图（位移不确定性等级）：

图a是形变场变形的不确定性与运动图像重叠，由蓝到红不确定性依次增加。图b按不确定性将图像分为三部分，右边是不同强度对应的核。
3.组织感知雅可比行列式正则化器
除了上述对速度场正则化和平滑外，还使用雅可比行列式正则化器以获得生理上合理的形变场$J(?)$。提出的组织感知自适应正则化约束区分GM、WM、CSF，定义如下：
$$Reg(?)=\{\begin{array}{llc}{e}^{(|min(J(?))?1|)?1}& ;& 如果是GM或WM\\ e^{(|mean(J(?))?1|)?1}& ;& 其他\end{array}\text{\hspace{1em}(4)}$$
对于GM和Wm雅可比行列式应该是正值，而背景和CSF雅可比行列式应该接近1，因为这些区域不应该有形变场。使用指数运算进行惩罚非常大或负的雅可比行列式的区域。

4.实验与结果

A.数据集

取自47个健康新生儿（2周大、3,6,12个月大）T1w和T2w的纵向驰豫（longitudinal relaxation）是高级儿科脑成像研究（Multi-visit Advanced Pediatric Brain Imaging，MAP）的一部分。T1w MR图像的TR（两次相邻脉冲间隔）为1900ms，TE（射频脉冲与其回波的间隔）为4.38ms，反转角度（flip angle，射频脉冲的参数）为7°，144 张矢状切片，1 mm 各向同性体素分辨率。T2w的参数为，TR=7380ms，TE=119ms，旋转角为150°，64张矢状切片和$1.25\times 1.25\times 1.95m{m}^3$的体素分辨率。
每个个体的图像数量都有所不同，训练集使用27个个体的56张扫描图，另外20个个体的57张扫描图用于推理和验证。为了有效地使用可用数据，我们模拟了形变场来扭曲可用的组织分割图以进行数据增强。

B.预处理

使用UNC婴儿预处理方案（pipeline），包含以下几个步骤：
1).N3 强度不均匀性校正；
2).进行颅骨剥离；
3).使用配准方法移除小脑和脑干；
4).刚性对齐同一个体的纵向驰豫图像和，并刚性对齐T1w相应的T2w图像；
5).使用基于学习的多源集成框架分割组织；
6).使用FLIRT将上述分割图与固定分割图刚性对齐；
7).对组织分割图和强度图进行重采样，得到大小为$256\times 256\times 256的1^{3}mm各向同性体素（分辨率）$；
8).皮层下结构的掩蔽和填充，以及大脑的左/右半球分离；
9).构建拓扑正确和几何准确的白质（WM/GM 界面）和软脑膜（GM/CSF 界面）皮质表面；
10).分割构建的皮质表面，然后执行标签的表面到体积映射来描绘 34 个皮质 ROI（区域）；
标签值通过最近邻域法（最近的顶点）进行填充。

C.评价指标

DSC（Dice similarity coefficient）系数[越大越好]、最大平均豪斯多夫距离（max average Hausdorff distance，MAHD）[越小越好]、折叠点比率（ratio of the folding points，RFP）[越小越好]。

D.对比分析

进行对象自身、对象间和空间配准，并与最先进的一些方法进进行比较。所有的实验使用组织分割图并非基于强度的图像。
1).对象间配准：
随机挑选不同对象的图像作为F/M，共获得1161对训练图像。推理阶段随机选择推理对象的12月大时的图像作为F，其余（2周，3月，6月大的）图像作为M，共获得54对对象间图像用于配准。
使用DSC评估GM与WM，用MAHD评估内/外大脑皮层的配准正确率，使用RFP评估行变场是否平滑。结果如下表1：

进行三组实验2周——12月（2周大的配准到12月大的）、3月——12月、6月——12月。观察表格可知，本论文提出的RAT-Net在三个评价指标上都优于其他方法（除RFP的Demons和NiftyReg）。分析原因上述两种方法可能是由于欠评估（under-estimate）和过平滑（over-smooth）。
配准结果和形变场如下图：

第一行为T1w图像，第二行为组织分割后的图像（可以看出图像分为三部分WM，GM和CSF（与背景色一致）），第三行为形变场（可以看出形变场在较为理想的配准情况下比较平滑（VoxelMorph的形变场就比较紊乱不平滑））。
还评估了34个区域的DSC如下图：

RAT-Net提高了小区域的准确率，图中加*的区域表示RAT-Net效果好的。
2).对象自身的配准：
对象本身的关键是测量婴儿皮层结构随时间的变化，验证了推理数据集所有对象自身的数据对。使用前面的对象间的配准网络，除了1161对对象间训练图像外又加入了64个对象自身的图像对完善网络，仍旧选取12月大的图像作为F，其他时间段为M（64对），共用57个图像对进行整体评估。结果如下表2：
DSC整体效果比较好，但也存在两个较差的情况。MAHD与DSC结果一致，RAT-Net略差于ANTs和NiftyReg相对于深度学习网络本设计的网络结果比较理想，但对比非深度学习结果较差。想对于VoxelMorph的RFP结果较好，形变场平滑。由此可知，RAT-Net对于GM配准和外皮层配准结果理想，WM和内皮层配准与传统配准方法相当。
3).空间配准：
结合前面的配准设置，使用数据集中所有对象的空间配准设置，对于每一个推理（预测）对象空间配准包括两步：（1）.同一对象的所有时间点图像首先与12月大的图像使用上面的对象自身网络配准；（2）.使用对象间网络将12月大的图像作为F，对所有对象进行配准，获得形变场。为了避免插值错误，合成变形场是每个对象的每个时间点都会被扭曲一次，最后与F图像对齐。试验结果如下表3：

RAT-Net的结果就DSC而言最好，RFP比Demons坏些，但是DSC的数据高了4.7%-6.7%。表明在空间配准任务中没有牺牲形变场的平滑来提高配准准确率。
表1与表3对比可知，对于Demons和NiftyReg来说，空间配准比直接进行对象间配准能平稳提升DSC；ANTs两种结果差不多；VoxelMorph的空间配准可以提升DSC，并且有较好的RFP；RAT-Net的结果却是变差了，但是仍好于其他网络结构。
上述实验表明RAT-Net对对象间配准配准好于对象自身的配准，需要注意的是对象间配准的推理时间点的图像与F是相互独立的。相反，在空间配准时，每个时间点都会使用同一对象的12个月大的时间点作为中间目标，然后才变形到固定图像，这可能解释了为什么空间配准的结果差于对象间配准。
在应对更多动态形变时，RAT网络由于其他网络；细微变形时优势则没有那么明显。这表明在空间配准时，对象间的配准（第二步）有利于提升性能。
作者总结了三点：
a.Demons和NiftyReg可能低估了直接对象间变形场，这可以通过空间配准来环节；
b.ANTs的结果相近，但是形变场的折叠情况可能很严重；
c.总体来说RTA网络结果比较好。

E.消融研究（验证网络里的几个模块是否起作用）

本文验证四部分循环推理模块、$Reg(v)$、$Reg(?)$和平滑层。实验结果如下表4：

**严格控制变量（网络参数和数据集）**在不设循环推理时，DSC明显下降，RFP变差；损失函数不包括$Reg(v)$时，结果同上；无$Reg(?)$时同上；无平滑层时，虽然有部分结果好些，其实几乎一样，只是RFP变得差了。因此，循环推理提升配准正确率；$Reg(v)$和平滑层在不牺牲准确率的条件下增加形变场约束；$Reg(?)$主要保持形变场的拓扑结构。

5.讨论

A.循环推理分析

修改循环次数计算GM和WM的DSC和RFP，结果如下图：

依图可见，前四次迭代对DSC和RFP都有较好的结果，但是从第六次开始趋于平缓，说明进行五次迭代效果最佳，也说明了循环推理策略不仅可以完善预估形变场，还可以平滑形变场。本设计还包括组织感知雅可比正则化和自适应平滑层在不牺牲配准率的同时抑制折叠形变点，同时，RAT网络为轻量级网络，训练推理时不需要特别大的GPU空间。

B.基于强度的配准

使用基于强度的婴儿图像进行训练可以吗？
作者使用ANTs进行测试，采用CC和MI作为相似性准则；使用弱监督和原始的VoxelMorph，弱监督的VoxelMorph在训练时不止使用强度图像还使用GM和WM的分割图，但是结果都不太理想，见下图：

C.参数分析

为了解释公式（1）中的各项含义，使用不同权重训练网络，每个对比实验，设置左边权重等于1，右边某一项权重介于0.1-10。结果如下图：

当增加$Reg(v)$，$Reg(?)$时，在降低DSC的同时还提高了RFP；增加相似性权重时，既增加了DSC又减小了RFP，说明损失函数中相似性权重大效果好。

D.缺陷与展望

主要创新点是循环推理策略，在较大变形配准时，提高配准准确率。深度学习的方法很难通过一次推理完成配准，设计的网络可以通过循环网络提升配准性能。主要集中于婴儿脑MRI图像和评估RAT网络对于分割图像的性能的提升，设计的网络对于其他的数据集有好的兼容性。没有分析各种指标对最终性能的影响，虽然那些比较的方法使用不同的度量（CC、MI等）来优化变形场作为它们的推荐配置（使用度量方法反向传播优化形变场），但是作者使用MSE和CC实验的结果却没有太大的改变。

=================================================================================

=================================================================================

补充

1.位移量的不确定性的推导

$v(p)$表示p处的不确定性，是一个具有三个独立分量的多元正态值：
$$v(p)=[v_x(p),v_y(p),v_z(p)]$$
$$|v(p)|=\sqrt{|v_x(p){|}^2+|v_y(p){|}^2+|v_z(p){|}^2}$$
$${\Sigma }_{|v|}(p)=\sqrt{({\Sigma }_v^x(p){)}^2+({\Sigma }_v^y(p){)}^2+({\Sigma }_v^z(p){)}^2}\le \sqrt{3}?{\Sigma }_v^{max}(p)$$
其中，$\Sigma$表示不确定性，上式含义是p点体素的不确定性小于最大不确定性的$\sqrt{3}$倍，由正态分布里的3$\sigma$原则可知，有超过99%体素$v(p)\in [{\mu }_v(p)?3?{\Sigma }_v(p),{\mu }_v(p)+3?{\Sigma }_v(p)]$。故设置最大不确定性等于方差，作为高斯核自适应平滑形变场，
$$\sigma ={\Sigma }_{|d|}(p)={\int }_0^{1}3\sqrt{3}?{\Sigma }_v^{max}(p)dt$$

2.基于不确定性的平滑层

如图，位移不确定性高的地方基本在脑外和靠近皮层的区域，原因可能为：（1）在各向同性区域上很难获得有意义的位移向量，（脑外为黑色区域，缺乏视觉线索信息）因此需要更大的平滑核；（2）由于高度解剖变异（生长），在复杂区域（靠近内外皮层的地方）建立精确关系也较为困难。
RTA网络通过迭代获得很好的效果，如下图：

有位移幅度图（左）可知，经过第一次循环后位移幅度和不确定性大幅下降，可以一次生成一个形变场但是可能具有较高的不确定性，然后通过循环完善不确定性。直接预测位移可以获得较高的DSC结果，但会牺牲形变场的光滑性，本设计使用自适应平滑，根据形变的不确定性选择合适的平滑核（见上$\sigma$），不确定性越高平滑核越大。

3.评价指标

DSC

实质是交并比$$DSC=2?\frac{|V_m\cap V_f|}{|V_m|\cup |V_f|}$$

MAHD

最大平均豪斯多夫距离是计算扭曲运动图像皮层表面$S_m$与固定图像皮层表面$S_f$间的平均最大距离，
d H , m ? a v g ( S m , S f ) = m a x { d ? H , a v g ( S m , S f ) , d ? H , a v g ( S f , S m ) } d_{H,m-avg}(S_m,S_f)=max\{\vec {d}_{H,avg}(S_m,S_f),\vec {d}_{H,avg}(S_f,S_m)\} dH,m?avg?(Sm?,Sf?)=max{ d H,avg?(Sm?,Sf?),d H,avg?(Sf?,Sm?)}
其中，${\vec{d}}_{H,avg}(S_m,S_f)=\frac{1}{|S_m|}{\Sigma }_{v_i\in S_m}d(v_i,S_f)$，$d(v,S)$表示顶点v与曲面S间的最小距离$|S_m|$为顶点数。

RFP

利用折叠点的比值来评价变形场的平滑度，负雅可比矩阵行列式表示折叠和撕裂，大于1时表示形变场膨胀，在0，1之间是表示压缩。

4.增强数据

采用两种策略：
1).从训练集中，选择一个对象按照下面步骤生成99个扭曲的组织分割图像，然后训练RTA-Net-DA1；
2).将训练数据集中所有图像增强9次训练RTA-Net-DA2；
增强方法涉及以下步骤：
a)在10×10×10网格上的每个维度的位移从[?4,4]mm均匀采样；
b)B样条模型用于根据网格点上的位移插入具有所需图像尺寸的密集变形场；
c)组织分割图用生成的变形场和最近邻进行扭曲；（首先形成一个均匀的位移场（阶梯型）大小为10×10×10，采样间隔为0.8mm，将位移场对修改到训练图像的大小）

将不同增强方法获得的网络对比，可以看出主体间配准RTA-Net-DA2的DSC较高，且RFP也较好。使用预定义的形变场增强数据集可以提升DSC的性能。