基于三维纹理特征仿真和脑图谱的小鼠脑磁共振显微成像分割算法

付振荣，林岚，王婧璇，邬雪涛，顾克楠，吴水才

北京工业大学 生命科学与生物工程学院，北京 100124

引言

利用高分辨率的磁共振显微成像（Magnetic Resonance Microscopy，MRM）进行小鼠大脑的神经影像学分析已经成为脑科学的研究热点[1-3]。和人脑MRI图像处理类似，非均匀场校正[4]、图像配准[5]以及脑区分割[6]等都是小鼠脑MRM图像处理的关键步骤，其准确性将直接影响后续分析结果。图像分割就是把图像分成若干个特定的、具有独特性质的区域并提取感兴趣目标的技术和过程。小鼠的脑图像分割就是按照小鼠大脑的解剖结构，将小鼠大脑分为若干个感兴趣区的过程。在小鼠脑的形态学分析中，往往对分割的精度有较高的要求[7]。

目前，小鼠脑解剖结构的“金标准”，是基于光学组织切片划分的Paxions-Franklin脑图谱（P-F Atlas）[8]，它一般可将小鼠大脑划分为上百个感兴趣区。与光学组织切片相比，活体采集的小鼠MRM脑图像分辨率和对比度都比较低，感兴趣区之间的边缘比较模糊，不但脑区划分有限[5]，往往数据量也偏少[9]。基于P-F Atlas来实现小鼠MRM感兴趣区的手工划分是非常复杂且枯燥无味的，需要耗费大量的时间。随着研究的不断深入，传统的手工划分耗时耗力且精细程度远不如解剖图谱[2,10]，对研究造成了一定的障碍。如果能够在P-F Atlas与活体MRM脑图像之间建立一一对应的空间关系，就可以将P-F Atlas中脑区映射到活体脑图像中，满足研究中更精细的分割需求。

马尔科夫随机场模型是目前最主要的二维纹理合成模型之一，其原理为待填充点的纹理只由其一定范围内的邻域决定。按照规定好的顺序合成输出图像中的某个待合成点时，该点的值由已合成的点的一定范围内的“邻域”来决定[11]。Efros等[12]提出了一个基于马尔科夫随机场模型的基于点的合成算法，首先从样本图像中选择一些纹理小块作为种子点；然后，对于待合成的像素点，通过计算样本图中所有像素与其邻域的关系，找到满足条件的点进行合成。但是，其每合成一个像素都要对样本图中像素点进行穷尽计算，计算量非常大。Wei[13]对该算法进行了改进，采用“L”形邻域以及通过多分辨率和树形结构矢量量化进行加速，在一定程度上提高了计算的效率，但计算量还是比较大。基于块合成的方法可以大大提高合成的效率，Efros等[14]提出了一种二维纹理块状合成算法，首先随机选择样图中的一块放到输出图像的左上角，然后在输入样图中搜索下一个纹理块，两个纹理块需达到“最佳匹配”，但较多的重复纹理会造成不匹配的情况。Liang[15]也提出了一种基于块的二维纹理合成算法，该种算法通过计算已合成区域与样图中纹理块重叠区域的相关性，在样图中找到相关性高于设定阈值的纹理块进行合成。我们这里将Liang提出的算法扩充到三维，通过计算三维立体的图像块的边界空间重叠区域的相关性进行合成。

本文主要分为3部分：首先提出一种基于P-F Atlas的三维块合成的小鼠MRM脑图像的纹理生成算法；然后，评估图像分割的准确性；最后，对该算法及流程的不足之处和未来研究方向进行讨论。

1　材料与方法

1.1　数据采集

本次研究所用到的小鼠MRM脑图像数据为新加坡国立大学计算功能解剖实验室（Computational Functional Anatomy Lab，CFAL）在网上提供的共享数据（http://www.bioeng.nus.edu.sg/cfa/mouse_atlas.html）。该数据包含5只C57BL/6雄性小鼠T2加权MRM图像（TR=2000 ms，TE=46 ms，视场FOV=9 mm×13 mm×25 mm，矩阵=88×140×256），采集设备为布鲁克7-T/20-cm ClinScan MRI系统，图像的分辨率为100 μm×98 μm×98 μm。这些图像都通过非参数不均匀强度归一化（Nonparametric Non-Uniform Intensity Normalization，N3）非均匀场校正[16]进行了预处理，然后刚体配准到P-F立体定位空间进行划分，基于组织切片图谱将每个鼠脑划分为相互独立的39个解剖结构（表1）。研究中使用的电子版小鼠解剖脑图谱为P-F Atlas[8]，小鼠的种类为C57BL/6J，跟活体MRM数据相同。

1.2　解剖图谱的划分

首先，我们将P-F Atlas划分为与活体MRM数据相对应的39个脑区，该部分工作全部在MRIcro软件上完成。由于P-F Atlas在小脑的划分上没有对白质和灰质进行区分，而MRM数据中的小脑进行了白质和灰质的区分，所以我们在研究中去除了分析中不常用的小脑部分。其次，由于海马体的划分原则比较模糊，因此我们将4个属于海马体的脑区合并为一个进。最终，我们从图谱中分割出的脑区为34个。图谱划分的流程图，见图1。手工分割出的脑区以及其与个体脑图谱的对应编号，见表2，合并脑区以对应的第一个编号为准。

图1　解剖图谱手工分割流程

表1　划分的39个脑区的名字及其相应的编号

表2　手工分割解剖图谱与个体图谱对应的脑区编号和名字

1.3　基于立体块的小鼠MRM三维纹理合成算法和仿真

基于立体块的小鼠MRM三维纹理合成算法的主要步骤为：① 将图像反转作为样本图像输入；② 在样本图像中随机选择一个体积一定的立方体放在待合成空白图像的左上角，作为第一个纹理合成的块；③ 在样图中搜索出所有与上述体积相同的立方体，将其分别放置在待合成位置，在边缘上与已合成区域有一定的体积重合（根据待合成立方体的位置，重合的区域有7种情况），然后计算这些重合区域的互相关性，将高于某一阈值的立方体放在一个集合中，最后，从集合中随机选择一个立体块作为待合成区域的纹理进行合成；④ 重复步骤③直到输出图像合成完毕；⑤ 将图像再次反转，输出纹理图像。

立体块大小的选择可以根据实际情况进行调整，如果所选立体块太小则会增加计算时间，而太大则可能会忽略局部纹理特征。本文选用的立体块大小为6×6×6个体素。实验证明，该算法可以很好地合成小鼠较大脑区的纹理特征。本文中，在解剖图谱中划分出的21个较大脑区由我们提出的这种三维纹理合成算法进行合成。因为小鼠大脑的体积非常小，且拥有较多体积微小的脑区，包括一些白质脑区，如胼胝体等，其形状都是细长怪异的，故这些脑区无法利用块状算法进行有效合成。此外，由于脑脊液的纹理特征不明显，我们也无法对脑室进行纹理合成。因此，剩余13个脑区利用蒙特卡洛法和高斯平滑完成合成。表2中，第一列为可进行纹理填充的21个脑区的编号和名字，第二列为剩余的13个脑区的编号与名字。

基于该纹理合成算法的仿真流程，见图2。该整体流程分为3个部分，分别为纹理合成、解剖图谱处理以及全脑仿真运算。纹理合成部分：首先，将个体图谱中的39个脑区进行合并与删除，生成与解剖图谱划分的34个脑区相对应的脑区；然后，利用我们提出的三维纹理合成算法生成21个较大脑区的三维纹理；最后，利用蒙特卡洛法生成剩余13个脑区的灰度特征高斯分布图像，然后利用高斯平滑（FWHM=[222]）对高斯分布图像进行平滑处理。

图2　小鼠MRM脑图像全脑仿真流程图

解剖图谱处理部分：首先，将划分的34个脑区二值化；然后，我们需要对划分的解剖图谱进行两次重采样，第一次是由于解剖图谱和原图像的维度不同，我们需要对解剖图谱在z轴方向上进行插值，在x轴和y轴方向上进行下采样，以得到与研究数据中MRM图像一样的物理维度。全脑仿真部分：根据纹理法生成的纹理图像和蒙特卡洛法生成的随机图像对图谱中各脑区进行填充，得到初步的仿真全脑图像；最后，对初步仿真全脑中的34个脑区的边缘进行三维平滑，得到最终的仿真小鼠大脑图。整体仿真流程的图像示意图，见图3～4，其中箭头的颜色与图2中对应。

图3　MRM原图像纹理与灰度分布合成流程图

图4　小鼠MRM全脑仿真图像

1.4　分割结果的评估

小鼠的全脑MRM仿真图像创建成功后，我们需要对仿真的准确性进行量化评估。因此，我们将仿真图像与原始图像进行配准，通过表面一致性（Regional Surface Concordance Ratio，SCR）以及曲面平均距离（Surface Average Distance，SAD）方法对仿真图像进行量化评估[17]。在进行配准之前，我们对表2中的图谱进行编号赋值，对应关系和表2中相同。在脑图像的配准算法中，我们选择af fi ne算法和geodesic-SyN算法。geodesic-SyN的参数设置与之前的研究有所改动[17]。由于仿真图像在配准过程中更容易出现振荡，我们将全局平滑参数由2设置为3。由于仿真图像和真实活体图像有着较大的差异，为了更好地进行配准，我们将金字塔层数由3层改为4层，具体参数设置为400×200×100×30。

2　结果

仿真图像与原始MRM图像使用不同配准算法的SCR（City Block 距离为2），见图5。af fi ne和geodesic-SyN两种配准算法的SCR均值和标准差分别为0.7449±0.1324和0.822±0.1253。在af fi ne算法中，SCR值大于85%的脑区个数占所有的34个脑区的11.77%，而geodesic-SyN配准算法中该比例为41.18%。SCR值小于60%的脑区，在全部34个脑区中所占的比例为14.71%（af fi ne）和8.82%（geodesic- SyN）。通过单因素方差（置信区间为95%）对两种算法的SCR值进行统计性分析，两种配准算法的SCR值存在显著差异（F=9.717，P=0.003）。

图5　仿真图像与原始图像在不同配准算法下的表面一致性展示图

仿真图像与真实图像使用两种不同配准算法的SAD示意图，见图6。SAD即所得两个感兴趣区表面之间的平均最短欧几里得距离。两种配准算法的SAD的均值和标准差分别为1.6685±0.6627（af fi ne） 和1.4120±0.7086（geodesic-SyN）。在af fi ne算法中，SAD大于2个体素的脑区占全部34个脑区的25.59%，geodesic-SyN算法中该比例为8.82%。在af fi ne算法中，SAD小于1个体素的脑区占全部34个脑区的2.94%，geodesic-SyN算法中该比例为32.35%。通过单因素方差（置信区间为95%）对两种算法的SAD值进行统计性分析，结果表明，两种配准算法的SAD均值在统计学上不存在显著差异（F=2.668，P=0.107）。

最后，我们将MRM原图像通过逆变换到解剖图谱的空间，来证明我们的纹理合成算法以及流程能够实际应用于小鼠脑组织的细分割。将变换到解剖图谱空间的MRM真实图像与解剖图谱进行融合，见图7。从该图中可以看出，我们的算法及流程能够较好的对活体小鼠MRM进行脑组织结构的细分割。

图6　仿真图像与原始图像在不同配准算法下的曲面平均距离展示图

图7　小鼠MRM真实图像逆变换与解剖图谱融合图

3　讨论与总结

本文将二维的纹理合成算法优化改进为三维小鼠MRM纹理合成算法并提出一套完整的小鼠大脑MRM仿真流程，在解剖图谱和活体MRM之间建立了良好的一一对应关系，对感兴趣区的进一步划分具有重要的意义。

从结果可以看出，我们对虚拟大脑的仿真比较成功，和真实大脑有较高的相似度，非线性配准算法相对于线性配准算法对配准结果有一定的改进。但是，和我们之前的研究[17]相比，在不同程度上仍存在一定的差异。我们认为，造成这种差异的原因主要分为两大部分：一部分是由小鼠MRM全脑仿真造成的误差，另外一部分是由配准造成的误差。

由小鼠MRM全脑仿真造成的误差主要分为以下5个因素：① 数据的个体图谱拥有39个脑区，而解剖图谱拥有上百个脑区，这就涉及到合并的问题，由于数据的划分原则没有细致到具体说明这一百多个脑区分别属于哪个脑区，所以划分原则上可能是存在一定的差异的；② MRM图像时按照等间隔采样的，而解剖图谱是按照光学组织切片的厚度不等间隔采样的，图像的插值会带来一定误差；③ 13个纹理不明显的脑区我们采用了蒙特卡洛法和高斯平滑的方式进行了灰度特征分布合成，与实际的纹理相比还是存在一定的差异；④ 解剖图谱是利用小鼠大脑光学组织切片划分的，但是脑脊液流失导致的脑室崩塌、人工固定以及切片过程会使大脑产生一定的形变，导致解剖结构和活体采集方式采集的图像在局部形态上存在一定的差异；⑤ 原始图像个体图谱的划分也存在一定的误差。配准造成的误差主要是由于geodesic-SyN算法在配准过程中，图像对不同的参数设置比较敏感，而我们使用的参数为针对MRM图像配准时的优化参数。但仿真图像配准更加容易出现由振荡所造成的收敛。虽然我们通过提高全局的平滑系数来防止这种情况的出现，但是增大平滑系数就会导致变形域精细度下降。我们使用的geodesic-SyN的参数设置对于仿真图像的配准来说不一定是全局最优，后期还需要进一步优化。

[参考文献]

[1] Lin L,Fu Z,Xu X,et al.Mouse brain magnetic resonance microscopy: Applications in Alzheimer disease[J].Microsc Res Techniq,2015,78(5):416-424.

[2] Badea A,Alisharief AA,Johnson GA.Morphometric analysis of the C57BL/6J mouse brain[J].Neuroimage,2007,37(3):683-693.

[3] 付振荣,林岚,张柏雯,等.基于MRM的小鼠脑模板创建的研究进展[J].中国医疗设备,2016,31(2):25-30.

[4] Lin L,Wu S,Bin G,et al.Intensity inhomogeneity correction using N3 on mouse brain magnetic resonance microscopy[J].J Neuroimaging,2013,23(4):502-507.

[5] Bai J,Trinh TL,Chuang KH,et al.Atlas-based automatic mouse brain image segmentation revisited: model complexity vs.image registration[J].Magn Reson Imaging,2012,30(6):789-798.

[6] Wu T,Min HB,Zhang M,et al.A prior feature SVMMRF based method for mouse brain segmentation[J].Neuroimage,2012,59(3):2298-306.

[7] Pagani M,Bifone A,Gozzi A.Structural covariance networks in the mouse brain[J].Neuroimage,2016,129:55.

[8] Franklin KBJ,Paxinos G.The mouse brain: in stereotaxic coordinates[A].The Mouse Brain in Stereotaxic Coordinates[C].Pittsburgh:Academic Press,2001:6.

[9] Ma Y,Smith D,Hof PR,et al.In vivo 3D digital atlas database of the adult C57BL/6J mouse brain by magnetic resonance microscopy[J].Front Neuroanat,2008,2(2):1.

[10] Lin L,Chen K,Valla J,et al.MRI Template and atlas toolbox for the C57BL/6J mouse brain[A].International IEEE Embs Conference on Neural Engineering,2005.Conference Proceedings[C].New York:IEEE Xplore,2005:6-8.

[11] 岳晓菊.基于样图的纹理合成算法的研究及应用[D].西安:西北大学,2010.

[12] Efros AA,Leung TK.Texture synthesis by non-parametric sampling[A].The Proceedings of the Seventh IEEE International Conference on Computer Vision[C].New York:IEEE,1999:1033.

[13] Wei LY.Fast texture synthesis using tree-structured vector quantization[A].Conference on Computer Graphics and Interactive Techniques[C].Boston:ACM Press/Addison-Wesley Publishing Co.,2000:479-488.

[14] Efros AA,Freeman WT.Image quilting for texture synthesis and transfer[A].Conference on Computer Graphics and Interactive Techniques[C].New York:ACM,2001:341-346.

[15] Liang L.Real-time texture synthesis by patch-based sampling[J].ACM T Graphic,2001,20(3):127-150.

[16] Sled JG,Zijdenbos AP,Evans AC.A nonparametric method for automatic correction of intensity nonuniformity in MRI data[J].IEEE T Med Imaging,1998,17(1):87-97.

[17] Fu Z,Lin L,Tian M,et al.Evaluation of fi ve diffeomorphic image registration algorithms for mouse brain magnetic resonance microscopy[J].J Microsc,2017,268(2):141-154.