wordpress 大数据平台 tws mirror 多进程 USB串口通信 idea 常用快捷键 golang 哨兵模式 meteor io vue版本 vue引入组件 pmp培训视频 jquery通过class获取元素 jq获取元素宽度 linux获取当前时间 mysql查询结果拼接 maya曲线建模 java使用redis 数据库学习 python自学教程 python函数内定义函数 java中的队列 java实例变量 linux格式化命令 磁盘分区软件 linux操作系统原理 assist是什么意思 qq钱包实名认证 苹果x银色 cmd命令数字雨 怎么看淘龄 quickchm 黑客攻防技术宝典 五子棋大师 键盘指法练习游戏 cad文件 opencv是什么 ae渲染设置 lol游戏环境异常
当前位置: 首页 > 学习教程  > 编程语言

【论文笔记】SYMNet:基于CNN的快速对称微分同胚图像配准

2020/7/24 10:00:55 文章标签:

本文是论文《Fast Symmetric Diffeomorphic Image Registration with Convolutional Neural Networks》的阅读笔记。文章中有很多关于数学的描述,包括微分同胚、雅克比矩阵没有看懂,所以这部分内容有所残缺,建议读者直接阅读原文。

一、摘要

由于之前的方法忽略了形变场的拓扑保持(topology preservation),并且形变场的平滑性只由全局平滑能量函数来约束,文章提出了一个名为SYMNet的无监督对称图像配准模型,该模型同时估计正向和反向的变换。

基于深度学习的配准方法中实质微分同胚性(substantial diffeomorphic properties)没有得到保证,即忽略了变换的拓扑保持和可逆性。模型的输出是一对微分同胚映射,用来将输入图像从两条测地线(geodesic path)映射到两幅图像的中间地带(middle ground)。微分同胚是可微的映射,并且存在可微的逆。

文章的主要贡献有:

  • 提出了一个快速对称微分同胚图像配准方法来保证变换的拓扑保持和可逆性;
  • 提出了新的方向一致性正则化,用负的Jacobian行列式来惩罚局部区域;
  • 文章提出的范式和目标函数很容易在其他地方应用。

二、 可变形图像配准

可变形图像(deformable image)配准的目标是建立起图像对之间的一个非线性关系,并且估计一个对齐两张图像的非线性转换。

可变形图像配准指的是将浮动图像变形以对齐到固定图像的过程,在该过程中两幅图像之间的相似性得到最大化。可变形图像配准是一个非线性配准过程用来建立浮动图像和固定图像之间的密集的体素级的非线性空间关系(dense voxel-wise non-linear spatial correspondence)。

FFMM分别表示固定图像和浮动图像,ϕ\phi来表示形变场,可变形图像配准过程可以表示为:
ϕ=argminϕLsim(F,M(ϕ))+Lreg(ϕ) \phi^{*}=\underset{\phi}{\arg \min } \mathcal{L}_{s i m}(F, M(\phi))+\mathcal{L}_{r e g}(\phi)
其中ϕ\phi^{*}表示最优的形变场,Lsim(,)\mathcal{L}_{sim}(\cdot,\cdot)表示不相似性函数,Lseg()\mathcal{L}_{seg}(\cdot)表示平滑正则函数,M(ϕ)M(\phi)表示变形后的浮动图像。

三、相关工作

1. 微分同胚配准

当前的配准方法通常是将位移场uu参数化,形变场和位移场的关系为:ϕ(x)=x+u(x)\phi(x)=x+u(x),其中xx是恒定变换。这种方式虽然简单,但是对于大的变形来说可能对应的位移场并不存在,所以文章使用的是带有静态速度场的微分同胚变形模型。微分同胚是可微和可逆的,这就保证了平滑的一对一的映射,同时也保留了拓扑性。微分同胚形变场ϕt\phi_t可以由速度场来生成:
dϕtdt=vt(ϕt)=vtϕt \frac{d \phi_{t}}{d t}=\boldsymbol{v}^{t}\left(\phi^{t}\right)=\boldsymbol{v}^{t} \circ \phi^{t}
其中\circ是复合运算符(composition operator),vtv^t表示时间tt时刻的速度场,ϕ0=Id\phi^0=Id是恒等变换。

形变场可以表示为李代数的一员,并将其指数化以产生时间为1的变形ϕ(1)\phi^{(1)},即李群ϕ(1)=exp(v)\phi^{(1)}=\exp(v)的一员。这表示指数化的流场使得映射是微分同胚和可逆的。给定初始形变场ϕ(1/2T)=x+v(x)/2T\phi^{\left(1 / 2^{T}\right)}=x+v(x) / 2^{T},其中TT表示总的时间,那么,ϕ(1/2)\phi^{(1/2)}可以通过重现(recurrence)得到:ϕ(1/2t1)=ϕ(1/2t)ϕ(1/2t)\phi^{\left(1 / 2^{t-1}\right)}=\phi^{\left(1 / 2^{t}\right)} \circ \phi^{\left(1 / 2^{t}\right)}

2. 基于学习的配准

大多数有监督的配准方法采用ground truth形变场或分割图来指导学习过程。文献11(DIRNet)证明了在无监督配准中采用互相关作为相似度评价指标的有效性,文献5(VoxelMorph)通过L2L_2正则损失使得形变场平滑,文献9提出了一个概率微分同胚配准模型。

四、方法

大多数基于学习的配准方法都是只从浮动图像映射到固定图像,而忽略了其逆映射(从固定图像映射到浮动图像)。让XXYY表示两个3D图像,可变形配准问题可以参数化为一个方程:fθ(X,Y)=(ϕXY(1),ϕYX(1))f_{\theta}(X, Y)=\left(\phi_{X Y}^{(1)}, \phi_{Y X}^{(1)}\right),其中θ\theta表示CNN的参数,ϕXY(1)=ϕXY(x,1)\phi_{X Y}^{(1)}=\phi_{X Y}(x, 1)ϕYX(1)=ϕYX(y,1)\phi_{Y X}^{(1)}=\phi_{Y X}(y, 1)分别表示时间为1的从位置xXx\in X变形到yYy\in Y和从位置yYy\in Y变形到xXx\in X的微分同胚形变场。

文章提出分别学习两个时间为0.5的将XXYY变形到它们的平均大小MM的形变场,模型收敛之后,时间为1的将XX变形到YY和将YY变形到XX的形变场就可以通过结合两个估计的时间为0.5的形变场而得到了。从XXYY的变形可以结构为:ϕXY(1)=ϕYX(0.5)(ϕXY(0.5)(x))\phi_{X Y}^{(1)}=\phi_{Y X}^{(-0.5)}\left(\phi_{X Y}^{(0.5)}(x)\right),从YYXX的变形可以结构为:ϕXY(1)=ϕXY(0.5)(ϕYX(0.5)(y))\phi_{XY}^{(1)}=\phi_{X Y}^{(-0.5)}\left(\phi_{Y X}^{(0.5)}(y)\right)。因此方程fθf_\theta可以重写为:fθ(X,Y)=(ϕYX(0.5)(ϕXY(0.5)(x)),ϕXY(0.5)(ϕYX(0.5)(y)))f_{\theta}(X, Y)=\left(\phi_{Y X}^{(-0.5)}\left(\phi_{X Y}^{(0.5)}(x)\right), \phi_{X Y}^{(-0.5)}\left(\phi_{Y X}^{(0.5)}(y)\right)\right)

1. 对称微分同胚神经网络

将方程fθf_\theta参数化为一个全卷积网络(FCN)、缩放层、squaring层和可微的空间变换器(STN),ϕXY(0.5)\phi_{XY}^{(0.5)}ϕYX(0.5)\phi_{YX}^{(0.5)}使用缩放和squaring方法通过估计的速度场vXYv_{XY}vYXv_{YX}来计算。FCN网络的结构类似于UNet,是一个5层的编码器-解码器结构,并且有跳跃连接。FCN以拼接的2通道的图像作为输入,以2个密集的(dense)、非线性的速度场vXYv_{XY}vYXv_{YX}作为输出。在编码器的每一层,都有两个卷积核大小为3×3×33\times3\times3,步长分别为1和2的卷积层。在解码器的每一层采用步长为1的卷积核大小分别为3×3×33\times3\times3的卷积层和步长为1,卷积核为2×2×22\times2\times2的反卷积层。在解码器的最后两层是卷积核大小为5×5×55\times5\times5,步长为1的卷积层,并且跟着一个softsign激活函数(softsign(x)=x1+x\operatorname{softsign}(x)=\frac{x}{1+|x|}),然后乘以常数cc,用来产生两个速度场。乘以常数cc是为了让速度场的范围在[c,c][-c,c]之间。实际操作时c=100c=100。除了最后一个卷积层,FCN的每一个卷积层后面都跟着一个ReLU激活函数。

此外,使用空间变换器(STN)实现缩放层和squaring层,并用它对预测的速度场进行积分。

特别的给定时间步长TT,初始化ϕXY(1/2T)=x+vXY(x)/2T\phi_{X Y}^{\left(1 / 2^{T}\right)}=x+\boldsymbol{v}_{X Y}(x) / 2^{T}ϕYX(1/2T)=x+vYX(x)/2T\phi_{YX}^{\left(1 / 2^{T}\right)}=x+\boldsymbol{v}_{YX}(x) / 2^{T},通过重现ϕ(1/2t1)=ϕ(1/t)ϕ(1/t)\phi\left(1 / 2^{t-1}\right)=\phi^{(1 / t)} \circ \phi^{(1 / t)}计算时间为0.5的形变场,直到t=1t=1

在这里插入图片描述

上图是模型的结构示意图,使用FCN来学习从XYX,Y到他们平均大小MM的对称的时间为0.5的形变场,绿色的路径表示从XXYY的变换,黄色的路径表示从YYXX的变换,为了简介没有在图中画出Lmag\mathcal{L}_{mag}损失。

上图是FCN网络的结构示意图。

2. 对称相似性

Lmean\mathcal{L}_{mean}是对称平均形状相似性损失,Lsim\mathcal{L}_{sim}是图像的相似性损失,有很多相似性度量,比如正则化的互相关NCC、均方误差MSE、距离的平方和SSD和互相关MI,这里选用的是NCC。让IIJJ表示两个输入图像,Iˉ(x)\bar{I}(x)Jˉ(x)\bar{J}(x)分别是IIJJ在位置xx的局部均值,局部均值通过以xx为中心,大小为w3w^3w=7w=7的窗口来计算。NCC的定义如下:
NCC(I,J)=xΩxi(I(xi)Iˉ(x))(J(xi)Jˉ(x))xi(I(xi)Iˉ(x))2xi(J(xi)Jˉ(x))2 \begin{array}{l} N C C(I, J)= \sum_{x \in \Omega} \frac{\sum_{x_{i}}\left(I\left(x_{i}\right)-\bar{I}(x)\right)\left(J\left(x_{i}\right)-\bar{J}(x)\right)}{\sqrt{\sum_{x_{i}}\left(I\left(x_{i}\right)-\bar{I}(x)\right)^{2} \sum_{x_{i}}\left(J\left(x_{i}\right)-\bar{J}(x)\right)^{2}}} \end{array}
其中xix_i表示以xx为中心的窗口内的位置。

Lsim\mathcal{L}_{sim}的表示如下:
Lsim=Lmean+Lpair \mathcal{L}_{s i m}=\mathcal{L}_{\text {mean}}+\mathcal{L}_{\text {pair}}

Lmean=NCC(X(ϕXY(0.5)),Y(ϕYX(0.5))) \mathcal{L}_{\text {mean}}=-N \operatorname{CC}\left(X\left(\phi_{X Y}^{(0.5)}\right), Y\left(\phi_{Y X}^{(0.5)}\right)\right)

Lpair=NCC(X(ϕXY(1)),Y)NCC(Y(ϕYX(1)),X) \mathcal{L}_{\text {pair}}=-N C C\left(X\left(\phi_{X Y}^{(1)}\right), Y\right)-N C C\left(Y\left(\phi_{Y X}^{(1)}\right), X\right)

其中,Lmean\mathcal{L}_{mean}用来衡量XXYY之间的不相似性,它们朝向平均大小MMLpair\mathcal{L}_{pair}衡量变形后的XXYY以及变形后的YYXX之间的不相似性。

3. 局部方向一致性

全局正则化会使配准的精度降低,实际上全局正则化不能保证变换的拓扑保持,为了解决这个问题,提出了可选性雅可比行列式正则,以促进预测形变场的局部的方向一致性。可选性雅克比行列式正则损失为:
LJdet=1NpΩσ(Jϕ(p)) \mathcal{L}_{J d e t}=\frac{1}{N} \sum_{p \in \Omega} \sigma\left(-\left|J_{\phi}(p)\right|\right)
其中NN表示Jϕ|J_\phi|中元素的总数量;σ()\sigma(\cdot)表示激活函数,这里用的是σ()=max(0,)\sigma(\cdot)=\max(0,\cdot),作用等同于ReLU;Jϕ()|J_\phi(\cdot)|表示在位置pp处形变场ϕ\phi的雅克比行列式,其公式如下:
Jϕ(p)=(ϕx(p)xϕx(p)yϕx(p)zϕy(p)xϕy(p)yϕy(p)zϕz(p)xϕz(p)yϕz(p)z) J_{\phi}(p)=\left(\begin{array}{lll} \frac{\partial \phi_{x}(p)}{\partial x} & \frac{\partial \phi_{x}(p)}{\partial y} & \frac{\partial \phi_{x}(p)}{\partial z} \\ \frac{\partial \phi_{y}(p)}{\partial x} & \frac{\partial \phi_{y}(p)}{\partial y} & \frac{\partial \phi_{y}(p)}{\partial z} \\ \frac{\partial \phi_{z}(p)}{\partial x} & \frac{\partial \phi_{z}(p)}{\partial y} & \frac{\partial \phi_{z}(p)}{\partial z} \end{array}\right)
但是并不是用可选性雅克比行列式正则损失代替原来的全局正则化,而是同时使用两者,以产生平滑且拓扑保持的变换。

接着,使用Lreg=pΩ(vXY(p)22+vYX(p)22)\mathcal{L}_{r e g}=\sum_{p \in \Omega}\left(\left\|\nabla \boldsymbol{v}_{X Y}(p)\right\|_{2}^{2}+\left\|\nabla \boldsymbol{v}_{Y X}(p)\right\|_{2}^{2}\right)来加强速度场的平滑,还通过幅值约束(magnitude constraint)Lmag=1N(vXY22vYX22)\mathcal{L}_{m a g}=\frac{1}{N}\left(\left\|\boldsymbol{v}_{X Y}\right\|_{2}^{2}-\right.\left.\left\|\boldsymbol{v}_{Y X}\right\|_{2}^{2}\right)来避免在路径上的偏移。所以总的损失函数为:
L(X,Y)=Lsim+λ1LJdet+λ2Lreg+λ3Lmag \mathcal{L}(X, Y)=\mathcal{L}_{s i m}+\lambda_{1} \mathcal{L}_{J d e t}+\lambda_{2} \mathcal{L}_{r e g}+\lambda_{3} \mathcal{L}_{m a g}

五、实验

在T1权重的脑部MRI数据集OASIS上做了基于图谱的配准实验,先将所有图像重采样到256×256×256256\times256\times256大小,然后使用FreeSurfer做标准的预处理,包括运动校正、去除头骨、仿射空间归一化和皮质下结构分割。然后将图像裁剪为144×192×160144\times192\times160大小。并将数据集划分为255、20、150分别作为训练集、验证集和测试集。随机从测试集选择5个MR图像作为图谱。

使用Dice相似系数(DSC)和雅克比行列式Jϕ|J_\phi|作为衡量标准。DSC的取值范围是[0,1],配准的越好对应的DSC值越大。雅克比行列式可以捕捉形变场的局部行为。选用ANTs包中的SyN、VoxelMorph(VM)和DIF-VM作为baseline。使用SGD优化器,学习率为1e41e^{-4},动量为0.9,λ1=1000λ2=3λ4=0.1\lambda_1=1000,\lambda_2=3,\lambda_4=0.1

在这里插入图片描述
上图是SYMNet与baseline模型的配准结果对比。
在这里插入图片描述
上图是平均DSC值(越高越好)和非正雅克比行列式的平均体素数(越低越好),可以发现SYMNet在DSC上取得最好的效果,在Jϕ|J_\phi|上效果次优。
在这里插入图片描述
上表展示了局部方向一致性损失的影响。

上表对比了各种模型运行所需要的时间,可以发现SYMNet所需时间最短。


本文链接: http://www.dtmao.cc/news_show_50200.shtml

附件下载

相关教程

    暂无相关的数据...

共有条评论 网友评论

验证码: 看不清楚?