作者: Guha Balakrishnan, Amy Zhao, Mert R. Sabuncu, John Guttag, and Adrian V. Dalca
期刊/会议: CVPR
年份: 2019
标签: registration, machine learning, convolutional neural networks
核心贡献
- 根据“摊销优化”的设计哲学，设计了一种根据输入的两张图片输出配准变形场的网络，代替了传统的需要逐个图片配准的算法，实现了更快的配准。网络通过大量数据学习了两张图片之间的配准函数，把计算耗时从测试阶段转移到了训练阶段。
方法概述
整体方法是无监督的范式，不需要标注好的变形场作为label。而是通过两个无监督的损失函数约束参数更新，相似度损失$L_{sim}$，用于计算目标图像$f$和配准图像$m$的相似度，$L_{sim}$有多种选择，可以根据不同图像之间对比度是否一致来选择不同的指标作为损失。
例如配准后图像$m \circ\phi$和目标图像$f$计算mse，适合两个图像之间对比度变化不大的情况： $$ MSE(f,m\circ\boldsymbol{\phi})=\frac{1}{|\Omega|}\sum_{p\in\Omega}\left[f(\mathbf{p})-[m\circ\boldsymbol{\phi}](\mathbf{p})\right]^2. $$ 或者用互信息损失，计算两张图像的相关系数，由于计算的是相关性不是直接计算误差，更适合两张图像之间对比度有明显变化的情况： $$ \begin{aligned}&CC(f,m\circ\phi)=\\&\sum_{\mathbf{p}\in\Omega}\frac{\left(\sum_{\mathbf{p}_i}(f(\mathbf{p}_i)-\hat{f}(\mathbf{p}))([m\circ\boldsymbol{\phi}](\mathbf{p}_i)-[\hat{m}\circ\boldsymbol{\phi}](\mathbf{p}))\right)^2}{\left(\sum_{\mathbf{p}_i}(f(\mathbf{p}_i)-\hat{f}(\mathbf{p}))^2\right)\left(\sum_{\mathbf{p}_i}([m\circ\boldsymbol{\phi}](\mathbf{p}_i)-[\hat{m}\circ\boldsymbol{\phi}](\mathbf{p}))^2\right)}.\end{aligned} $$
复习一下皮尔森相关系数的公式 $$ > r=\frac{\sum_{i=1}^n(X_i-\overline{X})(Y_i-\overline{Y})}{\sqrt{\sum_{i=1}^n(X_i-\overline{X})^2}\sqrt{\sum_{i=1}^n(Y_i-\overline{Y})^2}} > $$
平滑损失$L_{smooth}$，用于约束变形场平滑，不产生突变，约束变形场符合物理规律： $$ \mathcal{L}_{smooth}(\phi)=\sum_{\mathbf{p}\in\Omega}\lVert\nabla\mathbf{u}(\mathbf{p})\rVert^2, $$ 这里计算的是变形场的微分矩阵的F范数平方，但是由于变形场在离散的网格上而不是连续的，所以用有限差分（相邻体素的差值）来近似，计算成本也比较低
此外，作者还测试了把分割标签作为附加约束$L_{seg}$，这里计算的是配准后分割区域的dice系数 $$ \mathcal{L}_{seg}(s_f,s_m\circ\phi)=-\frac{1}{K}\sum_{k=1}^K\mathrm{Dice}(s_f^k,s_m^k\circ\boldsymbol{\phi}). $$ 通过计算分割区域形变后的dice系数来判断配准的准确性。
网络结构是UNet，但是对于本工作来说网络结构的设计不是重点，重点是上文提到的损失函数，输入是两通道的3D图像，将$m$和$f$拼接在channel维作为输入。网络输出变形场$\phi$后通过一个Spatial Transform进行$m \circ\phi$的配准操作。
Spatial Transform公式如下： $$ m\circ\boldsymbol{\phi}(\mathbf{p})=\sum_{\mathbf{q}\in\mathcal{Z}(\mathbf{p}^{\prime})}m(\mathbf{q})\prod_{d\in\{x,y,z\}}(1-|\mathbf{p}_d^{\prime}-\mathbf{q}_d|), $$ 其中$\mathbf{p}$表示一个voxel的坐标，$\mathbf{p}'=\mathbf{p}+\mathbf{u}(\mathbf{p})$中，$\mathbf{u}$表示变形操作。$\mathcal{Z}(\mathbf{p}^{\prime})$表示$\mathbf{p}^{\prime}$的8个邻居坐标（3维空间中有8个邻居，如果是2维就只有4个）。这里要注意的是，坐标$\mathbf{p}$和$\mathcal{Z}(\mathbf{p}^{\prime})$都是整数，但是$\mathbf{p}'$是经过变形得到的，通常不是整数，因此需要通过对$\mathbf{p}'$的邻居进行插值得到$\mathbf{p}'$处的强度值。另外一个要注意的点是值的流向是$\mathbf{p}'$到$\mathbf{p}$，也就是通过插值得到$\mathbf{p}'$位置的值后将其赋值给坐标$\mathbf{p}$，如下图所示：
实验设置
- 数据集：使用了多个公开数据集OASIS,ABIDE,ADHD200,MCIC,PPMI,HABS,Harvard GSP,FreeSurfer Buckner40
-
对比方法：
ANTs工具包和NiftyReg工具包
-
评价指标：
Dice系数，计算时间，以及变形场雅各比矩阵非正的数量$|J_{\boldsymbol{\phi}}|\leq0$
$|J_{\boldsymbol{\phi}}|\leq0$表示这一点发生了局部的折叠，翻转或者局部体积被压缩到0，也就是表示不合理的变形，不符合微分同胚的性质或物理上不合理，这个点的多少可以反应变形的合理程度
结果分析
在dice系数上所提方法和两个对比方法差别不大，显著优于只有刚性配准的结果，但是在计算时间时间上所提方法更快，这也是这个方法的出发点，通过摊销优化减少推理时间，所以其他指标接近即可。
互信息CC相对于L2作为相似度度量的模型鲁棒性更好。
instance-specific的方式能取得比大数据集上训练更好的结果，但是时间更长。
个人思考
- 可借鉴的点：首先摊销优化这种设计哲学可以作为一个网络设计的出发点。通过大量数据训练的网络，代替原本需要的逐subject计算。相当于用训练时间代替测试时间。其次是这个网络结构可以直接作为配准的baseline使用。
- 存在的问题：
- 与我研究的关联：对比方法
- 我的理解程度：
⭐️⭐️⭐️⭐️