基于深度学习的脑片图像区域划分方法

王松伟，赵秋阳，王宇航，饶小平

（1. 郑州大学电气工程学院，郑州450000； 2. 中国科学院武汉物理与数学研究所，武汉430071）

（∗通信作者电子邮箱wangsongwei_zzu@163.com）

0 引言

动物的大脑由数目庞大、不同形态及功能的细胞组成，通过不同的神经网络主导各种心理、生理活动［1］。研究大脑是治疗脑神经疾病、模拟脑网络功能、优化人工智能的基础。免疫组织化学染色［2］、乙酰胆碱酯酶染色［3-4］和尼氏染色［5-6］等神经组织染色方法的出现和大量应用，促使人们制作完成了多种动物的标准脑图谱［7-9］。标准脑图谱的建立为脑科学研究人员进行脑区研究提供了极大的便利。随着神经环路标记及神经染色技术［10-12］的进一步发展和磁共振成像技术（Magnetic Resonance Imaging，MRI）的应用［13-14］，往往需要对各个脑区的细胞数目、荧光信号强度以及分子表达水平进行量化分析比较，进而对脑网络组成结构及神经细胞分布特性进行研究。

参照标准脑图谱对脑片进行区域划分是脑区研究的基础，直接影响标记神经元及荧光信号在不同脑区分布的统计结果，从而对整体研究结论产生较大影响。自动化区域划分通常需要配准操作。尽管已经开发了许多小鼠脑图谱，但并不存在小鼠荧光染色图谱。因此，采用常见的单模态或多模态配准算法对小鼠脑图谱与小鼠荧光染色脑切片配准时，由于两种图像并不存在映射的对称性，导致配准效果较差。而且为保证小鼠脑切片中各个脑分区标记神经元数量统计的准确性，需要保持脑切片的不变性，将标准脑图谱配准到小鼠脑切片上。目前常见的小鼠脑切片与标准脑图谱配准方法有以下几种:1）参照标准图谱，在脑切片图像上人工识别脑分区，并对脑区轮廓进行手绘，但手绘需要依靠专家经验且标准不一，只适用于少量样本绘制，不适用大量样本规模化处理；2）通过Photoshop 图像处理软件，进行简易脑片图像的半自动区域划分［15］，此方法同样需要专家经验，且费时费力，无法进行规模化的脑区划分。3）采用传统多模态配准方式进行配准［16］，与模态内配准相比，模态间配准是更具挑战性的任务，因为很难定义有效的相似性度量以指导跨模态的局部配准。一般来说，互信息（Mutual Information，MI）［17］及其变体［18］是解决模态间配准问题的常用方法。然而，MI是一种全局相似性度量，其精确进行局部匹配的能力有限，因为局部区域中的体素数不足使得在计算MI时强度分布鲁棒性较差。因此，采用多模态配准算法对本研究任务进行配准时，其配准精度较差，目前仍然缺乏简单高效的脑片自动化区域划分方法。为克服多模态配准相似性度量的缺点，一些学者提出使用结构表示的方法将多模态配准转换为单模态配准。Heinrich 等［19］提出基于局部自相似的模态无关邻域描述符（Modality Independent Neighbourhood Descriptor，MIND）方法进行结构表示，该方法对噪声和强度差异鲁棒性较强。Wachinger等［20］提出基于拉普拉斯的结构表示方法，但该方法对图像噪声敏感且计算成本较高。Oktay 等［21］提出基于结构化决策森林表示的概率边缘图，并用于多模态超声图像配准上。但这些人为设计的低级特征并不能对图像的复杂特征进行很好的表示。

由于传统配准算法的泛化性能差，且传统配准算法采用迭代优化的方式进行参数寻优，处理速度慢，无法做到实时化配准。而深度学习方法具有很好的泛化能力且能够做到实时化处理。因此，随着近年来深度学习的发展，基于深度学习的医学图像配准方法逐渐成为配准热门的研究方向。Simonovsky 等［22］提出使用深度学习网络来估计两个图像的相似性度量以驱动迭代优化策略；Miao 等［23］使用深度学习算法去预测输入图像的配准变换参数；Krebs等［24］使用人工代理优化变形模型参数，该方法对前列腺MRI图像进行配准实验时，表现出了较好的配准结果；Marc-Michel 等［25］提出一种以全卷积神经网络为主要架构的配准网络，并在患者间的心脏图像上进行配准实验，表现出较好的配准精度与配准速度；刘晨等［26］提出一种结合深度学习与尺度不变特征转换（Scale-Invariant Feature Transform，SIFT）的多尺度B 样条配准方法；高宁等［27］提出一种基于级联的卷积神经网络实现渐进分层的人脸配准方法。但这些方法大都需要提供真实的解剖学标志信息。

基于以上脑科学研究需求及配准方法的问题，本文利用ARA（Allen Reference Atlases）［28］中 的 成 年 小 鼠Average Template、Atlas 脑 图 谱 和PCANet-SR（Principal Component Analysis Net-based Structural Representation）、U-net 以及空间变换网络，提出一种无监督多模态脑片图像的区域划分方法。对于非刚性配准任务，与传统的基于优化的配准算法相比，该配准方法可以在测试阶段有效地配准两个图像而无需迭代优化或参数调整。该配准网络架构首先采用PCANet-SR网络提取并融合图像的多级特征，将非同一模态的脑切片与标准脑图谱进行同模态转换，然后送入配准神经网络中进行训练，训练完成后将训练得到的空间变形场作用到Atlas 脑图谱上，最后提取变形后的Atlas 脑图谱的边缘轮廓与原始脑切片融合，完成脑区划分工作。此配准框架较为简单，可以扩展到任何模态间配准任务。本文使用模态内相似性度量局部归一化相关系数（Normalized Correlation Coefficient，NCC）［29］来有效地训练多模态配准网络，而非直接定义不同模态的相似性度量。该方法适用于不同模态的脑片和脑图谱配准且在训练期间不需要解剖标志等监督信息，能够精准快速地进行无监督多模态脑片配准。本文算法的流程如图1所示。

图1 脑片图像区域划分算法流程Fig.1 Flowchart of brain slice image region division algorithm

1 相关知识

1.1 脑片图像获取

小鼠脑片图像获取主要分为三步:嗜神经工具病毒注射、取材及样本制备、免疫组化及共聚焦成像。首先选取8 周龄成年小鼠进行麻醉及消毒等准备工作，然后在小鼠大脑目标区域注入可以表达荧光蛋白的日本脑炎病毒（Japanese Encephalitis Virus，JEV）［30］。待病毒表达7.5 d 后，对小鼠进行过量麻醉，然后进行心脏灌流手术。手术完成后使用多聚甲醛固定液对小鼠脑组织进行固定，再采用震荡切片机对小鼠脑组织进行切片处理并在磷酸盐缓冲液中保存切片。最后使用10%羊血清和抗JEV 病毒对目标切片进行免疫组化处理，使用共聚焦显微镜（Leica SP8）对脑片红绿蓝三通道进行荧光成像，然后导出Tiff 格式图片。小鼠大脑内与目标区域相关的标记神经网络切片如图2所示。

1.2 PCANet-SR特征提取网络

PCANet-SR 网络是一个串级的线性网络，由基本的数据处理组件组成:级联主成分分析（PCA）组件、输出组件、结构表示组件。级联主成分分析组件由两个隐藏层组成:第一个隐藏层通过提取输入图像的主要特征向量产生卷积核，输入图像与卷积核进行二维卷积操作后输出，作为第二个隐藏层的输入；第二个隐藏层的结构与第一个隐藏层类似。输出组件使用sigmoid 函数对第二个隐藏层的输出进行表示。最后使用指数函数将第一个隐藏层的输出和第二个隐藏层的输出进行多特征融合，然后输出PCANet-SR 网络感兴趣的图像特征。

图2 标记小鼠脑神经网络的脑切片和标记神经元局部放大图Fig. 2 Mouse brain neural network labeled brain slice and local enlargement of labeled neurons

1.3 医学图像配准

医学图像配准是指针对一幅图像，寻求一系列空间变换，使之与另一幅图像上的空间对应点达到一致。通常，变换操作将来自目标生理学的同源位置映射到源生理学位置。执行配准的两幅图像中:一幅被称为参考图像，另一幅被称为浮动图像。参考图像由F 表示，浮动图像由M 表示，空间变换由φ表示。配准的目标是估计优化能量函数。能量函数定义如下:

L(F，M，φ) = Lsim(F，M(φ)) + λLR(φ) （1）其中:M(φ)表示被空间变换关系φ 扭曲的浮动图像M；函数Lsim表示对参考图像F和扭曲后的浮动图像M(φ)进行相似性度量；LR表示对空间变换关系φ施加正则化约束，使空间变换更加平滑；λ 表示正则化系数。优化问题会根据相似性度量函数的选择对能量函数进行最大化或最小化。

1.4 配准

配准算法主要涉及三个部分:变形模型、匹配准则和优化方法。常见的变形模型有弹性模型、扩散模型、微分同胚和弹性样条等。匹配准则有基于描述符的匹配、基于属性和强度的匹配及混合方法的匹配。本文的匹配准则即通过将两种模式映射到公共模式，从而将多模态配准转换为单模态配准。优化方法包括梯度下降法、共轭梯度法、拟牛顿法、高斯牛顿法以及贪婪法和进化算法等。

2 本文方法

2.1 脑片图像区域划分

2.1.1 全局粗配准

对图像进行全局刚性粗配准可以校正图像的全局形变，再对图像进行配准时可以取得更好的配准效果。通过嗜神经工具病毒JEV 标记小鼠大脑神经网络得到脑切片图像，脑切片中因个体差异的存在，导致与标准脑图谱形变差异较大，且在脑切片的制作过程中存在大位移误差，这些误差会对配准产生较大影响，进而影响配准精度。通过提取脑切片与脑图谱灰度特征，并对脑切片灰度值进行调整，然后通过高斯滤波平滑处理，再使用优化函数对两幅图像的均方误差进行优化，最后对变换图像进行插值，以达到对标准脑图谱仿射变换的目的。仿射变换后的标准脑图谱如图3所示。

图3 脑切片与仿射变换前后脑图谱Fig.3 Brain slices and the brain maps before and after affine transformation

2.1.2 PCANet-SR模态转换

在该配准框架中，PCANet-SR 方法［31］用于提取图像的多级特征。手工设计的结构表示特征虽然可以对特定任务下的图像进行较好的特征提取，但是这些特征的泛化能力有限，对于新问题往往需要构建新的特征，且对于结构复杂的医学图像特征并不能进行很好地表示，因此面对配准任务并不能得到较高的配准精度。而PCANet-SR可以提取图像的多级特征并进行特征融合，可以将不同模态的图像转换为同一模态，减小模态之间的差异，可以从感兴趣的数据中自动学习特征，且具有优秀的泛化能力，因此PCANet-SR 方法面对配准任务可以得到更好的配准性能。具有两个隐藏层的PCANet-SR模型如图4所示。下文将详细介绍PCANet-SR的结构。

图4 使用PCANet-SR多级特征融合的结构表示Fig.4 Structure representation using PCANet-SR multi-level feature fusion

衰减参数a1和a2可以保证在PCANet-SR 中对于相似性较高的区域有一个较明显的响应。衰减参数对于每个图像像素是局部自适应的。衰减参数与噪声方差成正比［33］，因此衰减参数的选择与图像的梯度相关。

2.1.3 配准网络架构

对配准网络进行建模Reg-Net:(MPSR，FPSR；θ，λ) ⇒φ。对于输入的浮动图像M 和参考图像F，首先使用仿射变换进行预配准，然后使用PCANet-SR进行同模态转换，最后送入无监督深度学习配准框架中进行训练，输出与输入图像相同尺寸的变形场φ，使用空间变换网络Τ以及训练得到的变形场φ对浮动图像进行扭曲以适应参考图像。在测试阶段，只需要输入新的浮动图像和参考图像，通过配准网络Reg-Net即可获得空间变形场φ，配准网络架构如图5所示。

图5 配准网络框架Fig.5 Registration network framework

网络以MPSR和FPSR为输入，配准网络Reg-Net中的参数θ（卷积核参数）在训练时进行学习，根据参数θ 计算空间变形场φ，使用空间变换网络和空间变形场φ 将MPSR扭曲为MPSR(φ)，使得模型能够评估MPSR(φ)与FPSR的相似性并更新θ。

图6 配准网络Reg-NetFig.6 Registration network Reg-Net

空间变换网络通过解码器的输出φ 扭曲浮动图像MPSR，空间变换层使神经网络能够执行全局参数化图像对齐，而无需监督标签［35］。为了使用基于标准梯度的方法，本文构建了基于空间变换器网络的可微分运算来计算MPSR(φ)，因为图像值仅在整数位置处有定义，所以需要在相邻像素处线性插值。

配准网络通过最小化损失函数进行迭代优化。损失函数被定义为:

其中:MPSR(φ)是被MPSR和空间变形场φ 扭曲后的图像；φ 表示从FPSR到MPSR指定的每个像素的向量偏移量；函数Lsim(·)用来度量FPSR和MPSR(φ)的相似性；Lreg(φ)表示在φ 上施加正则化约束；和(φ)分别表示参考图像和扭曲后的浮动图像的去均值后的灰度值。在实验中局部范围取D = 9，使用局部归一化相关系数定义Lsim(·):

其中:p 表示图像中的所有像素值；q 表示局部范围内的像素值；Ω 表示整个图像域；D 表示局部大小为9 × 9 的图像域。局部归一化相关系数越大，说明配准效果越好。通过最小化Lsim(·)来训练配准网络，但是可能产生不连续的变形场φ，因此需要施加正则化约束对φ进行平滑处理:

图7 部分仿射变换及PCANet-SR处理后的训练图像Fig.7 Some training images after affine transformation and PCANet-SR processing

其中:∇表示拉普拉斯变换；λ 表示正则化参数。本文实验中λ取值为0.55。

完整的损失函数定义如下:

2.2 数据集制作及网络训练

用免疫组化及共聚焦成像脑切片与ARA 中的成年小鼠Average Template 脑图谱建立数据集。首先由脑科学专家人工选定与脑切片对应的Average Template 脑图谱，在预处理中，将所有图像裁剪为相同尺寸(512 × 512)，提取图像灰度特征，然后使用仿射变换进行全局粗配准，并使用镜像操作进行数据扩充，最后将扩充后的数据送入PCANet-SR 网络中进行模态转换，生成脑切片FPSR和Average Template 脑图谱MPSR数据集。脑切片FPSR与对应Average Template 脑图谱MPSR为一组图像，共使用64组图像送入配准网络中进行训练，8组图像进行验证，8 组图像进行测试。部分训练图像组如图7所示。

本文方法在以Tensorflow为后端的Keras深度学习框架上实现，Keras 是一个基于Python 编程语言的深度学习框架，它可以对几乎所有类型的深度学习模型方便地定义和训练。Kera 目前有三个后端引擎实现:Tensorflow 后端、Theano 后端和微软认知工具包后端，由于Tensorflow 的应用最为广泛，且可扩展、可用于生产环境，所以本文最终选择Tensorflow 作为Keras 的后端引擎。配准网络在Nvidia RTX2080Ti GPU 上训练。采用自适应矩估计（Adaptive moment estimation，Adam）优化策略，Adam优化器利用梯度的一阶矩估计及二阶矩估计对每个参数的学习率实现动态调整。相对于其他优化器，Adam优化器实现简单、计算高效，适合解决包含大规模数据和参数的优化问题，且超参数具备可解释性，因此本文实验选择Adam 优化器作为优化策略。通过模型训练曲线确定学习率为0.01。每轮训练1 000 次，共训练100 轮，当损失函数没有显著降低时，停止训练。在本文实验中，训练耗时约15 h。

3 实验与结果分析

训练完成后，使用新的图像进行配准效果测试。将测试生成的空间变形场同样作用到与之对应的Atlas 脑图谱上，使之发生同样的形变。

从如图8 所示的配准效果可看出:脑切片与标准脑图谱配准精度较高，基本边缘及内部轮廓重叠程度较好，且变形后的Atlas 脑图谱自动将脑切片图像分为若干亚区，各个亚区之间界限清晰，满足实际脑科学研究需求。

本配准任务的精度通过脑切片与变形前后的Average Template 脑图谱之间的均方根误差（Root Mean Squared Error，RMSE）、相关系数（Correlation Coefficient，CC）以及MI 进行衡量。均方根误差越小，两幅图像的相似度越高；相关系数是广泛使用的图像相似性测度之一，它反映了图像的相关程度，相关系数越大，两幅图像配准效果越好；互信息用于描述两幅图像的统计相关性，两幅图像达到最佳配准时，互信息值达到最大。B 样条配准方法是单模态配准方法，经过PCANet-SR 预处理后可以对多模图像配准。本文实验测量了8 组原始图像仿射变换前后、B 样条配准前后、PCANet-SR+B 样条配准前后以及本文方法配准前后的均方根误差、相关系数以及互信息，最后取8组图像指标的平均值，结果如表1所示。

表1 配准方法性能对比表Tab. 1 Registration method performance comparison

从图9 所示的原始脑切片图像和不同配准方法得到的区域划分结果可看出，本文提出的配准算法边缘轮廓匹配程度优于其他几种对比方法。

由实验结果可知，本配准方法配准结果优于PCANet-SR+B 样条配准方法。配准完成后，RMSE 明显降低，CC 及MI 明显上升，与PCANet-SR+B 样条配准方法相比，本文方法的均方根误差降低了1.6%，相关系数和互信息值分别提高了3.5%、0.78%。PCANet-SR+B 样条配准方法也可以对多模态的图像进行配准，但是与本文方法相比，PCANet-SR+B样条配准方法需要调整内部参数，且配准过程较为耗时，两幅512 ×512 大小的图片进行B 样条配准需要231 s。本文配准方法可以实现无监督多模态配准，只需要训练一次配准网络，即可将配准网络模型泛化到新的配准任务上去，采用训练得到的参数对新的图像对进行配准，只需要2 s 即可完成配准操作，能够实现对多模态图像的快速配准和区域划分。本文区域划分方法的配准指标及实际区域划分效果均优于现有区域划分方法。

图8 部分脑切片配准及区域划分效果Fig.8 Some brain slice registration and region division effects

图9 脑切片图像与不同方法融合结果对比Fig.9 Comparison fusion results of brain slice images with different methods

4 结语

本文方法相对于手工划分区域以及借助PS 软件进行区域划分方法具有以下特点:无需专家经验，只需对大脑的一些标志性结构有所了解，即可规模化对脑区进行划分，且区域划分准确、快速；同时可以进行大规模批量化处理。本文实验展示的是小鼠冠状脑切片的区域划分，但是提出的脑片自动区域划分方法也可应用于小鼠矢状脑切片以及大鼠和其他动物脑切片的区域划分任务上。

本文进行了免疫组化及共聚焦成像脑切片与ARA 中的成年小鼠Average Template 脑图谱的配准及Atlas 脑图谱的区域划分，并且该算法的配准效率和配准精度均优于目前现有的多模态配准算法。此外，提出的配准框架也易于扩展并应用于其他多种形式的配准任务上。未来将侧重于将此配准框架向3D数据库上扩展。