相关研究表明猕猴和人类的大脑结构有着一定的相似性，猕猴(具有较高智力的非人灵长类动物)代表了最佳的用于验证假设的可侵入性模型[1]，相对于人脑，其在研究脑结构和认知方面有着无可比拟的优势。由于上述原因，近几年关于猕猴的研究正在逐渐增多[2-3]。
对核磁图像进行分析是研究猕猴大脑的重要手段，作为核磁图像分割领域的重要研究方向，核团分割的准确性及鲁棒性在促进脑研究领域的发展方面具有重要的研究价值和意义。然而对猕猴大脑的分割结果目前大多由研究人员借助人脑分割工具(如FSL)中提供的方法获得，然后再进行手动修正[4]。2017年，BALBASTRE[5]在高斯混合模型(gaussian mixture model，GMM)方法基础上进行改进，应用于猕猴的T2图像进行分割，但其能同时分割的脑组织有限。多图谱分割方法与该方法相比可以引进更多的先验信息，提高分割精度，所以很多研究人员将研究视线投向了多图谱方法。一些研究人员发布猕猴的各种MRI模板、图谱以及数据集[6-8]。MILHAM et al[9]对现有的非人灵长类公开数据集进行了比较。陈伟导等[10]用多图谱分割方法对猕猴脑组织进行了分割。最近，随着深度学习方法的迅速发展，有研究人员将深度学习技术应用于核磁图像分割。HUO et al[11]采用3D-UNET网络对人脑进行三维分割，ZHAO et al[12]利用深度学习对非人灵长类动物的大脑进行提取。但此方法所需的训练数据集较大，还有较大的改进空间。在多图谱方法人脑分割领域，ZHANG et al[13]将稀疏表示的算法引进到多图谱标签融合的领域后，取得了较好的效果。严盟等[14-15]提出的方法在稀疏表示标签融合方面做出了一些努力。夏瑞等[16]提出了一种基于重采样的改进的多图谱分割算法。LIU et al[17]提出了一种基于图集配准和线性化的核稀疏代表分类器的有效分割方法，但算法运行时间较长。
本文选用稀疏表示多图谱方法对猕猴脑组织的MR图像进行分割。为了提高大脑皮下核团的分割准确度，针对稀疏块包含的信息较少和稀疏表示方法会丢失信息的问题，本文提出了一种改进的稀疏表示多图谱分割标签融合算法(multi-atlas segmentation label fusion algorithm based on improved sparse representation，ISRLF).该方法在稀疏表示的过程中引入图谱标签，图像中的标签信息用来增加图像块的先验信息，使计算的稀疏系数更加合理；然后通过改变信息熵的计算方式改进互信息，将其用来衡量线性配准后的图谱与目标图像之间的相似度，将相似度的值与稀疏系数相结合计算出最终分配给图谱的权重，然后进行标签融合；随后通过设计的一个相似性度量指标将非局部块加权融合的结果和稀疏表示融合的结果进行综合，该指标包括Dice系数和余弦距离两部分，最终达到了较好的分割效果。
1 背景描述
稀疏表示最初应用于信号处理领域，其思想是将原信号用一些基本信号的组合代替，这样可以尽可能地减少处理的数据量，如公式(1)所示，T为被表示信号，D为字典矩阵，α为稀疏系数。稀疏表示在图像处理领域已被应用于图像的压缩，超分辨率重建，特征提取表示等方面。稀疏系数α可利用公式(2)求得。式中x为待表示信号。
T=Dα.
(1)
(2)
2 ISRLF方法
2.1 方法框架
本文将稀疏表示多图谱分割方法进行了改进，提出了一种将稀疏表示融合方法与非局部块加权方法综合在一起的新标签融合算法，对猕猴大脑皮下核团进行分割。图1为本方法的框架图。
算法步骤如下：
1) 对图谱灰度图像进行线性配准，得到线性配准后的图像Li.
2) 对线性配准后的图谱灰度图像进行非线性配准，得到配准后的结果Fi.将图谱标签图像非线性配准到目标图像中，得到配准后的结果Ii.
3) 对目标图像采用局部加权融合的方法进行预分割，得到预分割结果。
4) 构建目标图像标签块Tpx和图谱图像标签块Tpi，所有图谱图像标签块cdi组成过完备字典D，即D=[Tp1,Tpi,Tpn]，进行稀疏表示。
5) 引入互信息并对信息熵的计算方法进行改进，然后对目标体素和线性配准后的图谱相对应的平面进行全局的相似性衡量，求平均值，改进稀疏表示计算出的权重。
6) 分别用稀疏表示的标签融合方法和非局部块加权的标签融合方法对目标图像进行分割，得到各自的分割结果。
7) 采用Dice系数和余弦距离相结合的相似性度量标准对分割结果进行度量，将两种融合方法的分割结果再次进行融合，最终得到分割结果。
图1 方法框架图
Fig.1 Method framework
2.2 目标图像标签块和图谱图像标签块的构建
一般的稀疏表示方法其图像块所包含的信息只有目标图像对应位置的灰度值，没有充分利用更多的先验信息。本方法将灰度图像中的灰度信息和标签图像中的标签信息结合，分别引入到目标图像和图谱灰度图像中构成目标图像标签块和图谱图像标签块，使图像块携带的信息更加丰富，从而优化标签融合时各图谱的权重。
本文采用D99-SL、NMT和inia19三个图谱进行分割。为了将图谱中的标签信息引入，需要对目标图像进行预分割，本方法采用局部加权法对目标图像进行预分割。构建初始的仅含有灰度信息的图像块后，将对应的标签值块变成二值化的形式，再与初始图像块结合构成新的图像标签块。二值化的具体方式为以中心体素的标签值为基准，与之标签值一致的标签记为1，否则记为0.
图谱图像标签块构建的具体过程如图2所示，以非线性配准后的图谱灰度图像Fi中位于x处的体素作为中心，以r×r×r为大小构建初始图像块adi，再从非线性配准后的图谱标签图像It中，以对应位置的标签为中心，以r×r×r为大小构建标签块bdi并二值化。将adi与bdi分别变成列向量后进行拼接，组成图谱图像标签块Tpi.
图2 图谱图像标签块的构建
Fig.2 Construction of atlas image label block
目标图像标签块构建的过程与图谱图像标签块的构建过程相似。首先，以目标体素x为中心，r×r×r为大小构建目标图像块adx，然后以预分割结果的对应位置x为中心，r×r×r为大小构建标签块bdx，分别变成列向量，最后将adx与bdx合并，组成目标图像标签块Tpx.
对于目标图像标签块Tpx，选取各图谱图像标签块Tpi组成过完备字典D，即D=[Tp1,Tpi,Tpn].按照公式(2)进行稀疏矩阵的计算，式中的x为目标图像标签块，计算出的稀疏系数作为图谱的初始权重矩阵W1.
2.3 基于互信息的图像相似度计算
稀疏表示的标签融合方法只考虑了以目标体素为中心的图像块中的体素，并没有从全局的角度考虑图谱图像与目标图像的相似度，在标签融合时对整体图像相似性的衡量有助于提高分割效果。互信息是信息论中的一个概念，在图像处理领域常通过衡量两幅图像的信息熵比较两幅图像的相似程度，所以本方法引入互信息度量参与融合的图谱灰度图像的相似度。它的值越大代表两张图片的相似性越高。公式(3)为图像的信息熵计算公式，式中的概率pM由公式(4)求得，公式(4)中通过灰度直方图进行信息熵的计算，hi表示灰度值为i的体素数量，分母表示该图像的体素个数。
(3)
(4)
通常灰度图像之间的互信息是通过灰度直方图进行计算，对于核磁图像来说有不足之处。例如，核磁图像中灰度值为零的点表示的是图像的背景，一般的互信息会将其纳入相似度计算之中，而核磁图像中只关注大脑区域，背景基本上可以忽略不计。假如核磁图像中灰度值为零的点过多，而不为零的点过少，就可能会出现计算出的互信息值大，但真正重要的脑部结构可能相似度并不高的情况。此外，只对灰度图像进行相似度计算引入的先验信息有限，不能很好地衡量两幅图像的相似性。
为了对上述两种情况进行改进，本方法对信息熵的公式进行了修改。首先将灰度值为零的体素剔除，不纳入概率计算中，按照公式(5)计算灰度值的概率pM′.
(5)
其次由于大脑结构中灰白质占有着重要的位置，本方法关注的是皮下核团，所以本方法通过配准后的组织概率模板，按照公式(6)求得整体图像中体素为灰质的概率与体素为灰白质概率的比值，达到衡量脑结构体素相似度的目的，该计算结果作为参数β，引入到灰度图像的信息熵计算之中。
(6)
为了增加先验信息，可以将标签图像引入，将标签图像引进的方法是对标签图像按照核团的标签值进行信息熵计算。标签图像信息熵的计算方式如公式(7).Im表示标签图像，m表示的是图像中核团的数量。由于Dice系数可以表征相同核团的相似度，所以将其作为参数引入到计算过程中以增加先验信息。如公式(8)，求得目标图像分割结果IA与IB关于标签j的Dcie系数，将公式(8)计算出的值作为参数γ引入到改进的信息熵的计算中。
(7)
γ=Dice(IAj,IBj) .
(8)
由此，经过改进的信息熵表达形式为公式(9).对于给出的两幅图像A和B，H1(A)和H1(B)分别为图像的信息熵，计算方式为公式(10)、(11)，H1(A，B)为图像A和B的联合信息熵，计算公式为(12).I1(A,B)为A与B的互信息，计算公式为(13).
H1(M)=βH(M)+γH(IM) .
(9)
(10)
(11)
log2pA′B′(a,b)+γH1(IA,IB) .
(12)
I1(A,B)=H1(A)+H1(B)-H1(A,B) .
(13)
将计算出的矢状面、冠状面、横断面三个不同面的互信息作为S1,S2,S3矩阵，提取矩阵中与目标体素(a,b,c)对应的相似度值，然后按照公式(14)进行计算，得到体素x的相似度值，再按照公式(15)、(16)，得到最终的图谱融合权重。
(14)
W2=W1×S(x) .
(15)
(16)
2.4 标签融合
基于图像块的融合方法有两种，一种是通过一个相似性度量函数度量所有参与标签融合的图像块与目标图像块的相似性的非局部块加权融合方法，一种是稀疏表示方法，将目标图像块由参与融合的图谱图像块稀疏表示。
非局部块加权的标签融合方法有着参与融合的图像块没有被舍弃、冗余信息多的优点，稀疏表示的融合方法与非局部块加权的融合方法相比，该方法会舍弃一些信息。但该融合方法也有着非局部块加权的融合方法难以比拟的优势，其在边界点的分割效果上比非局部块加权的融合方法有着更好的表现。
为了结合这两种方法各自的优点，提高对大脑皮下核团的分割准确度，本方法将非局部块加权方法和稀疏表示的标签融合方法进行了结合，从而提出了一种新的标签融合方法ISRLF.
基于图像块加权融合方法的过程如下：
1) 将各图谱灰度图像和标签图像配准到目标图像T中，得到配准后的灰度图像Fi和标签图像Ii.
2) 提取T中以位置x为中心大小为r×r×r的目标图像块Tpx.(本实验r值为3)然后进行图像块的提取。提取过程为各图谱以当前要融合的体素x为中心，r×r×r作为搜索邻域，提取所有的参与标签融合的图像块Tpi.
3) 采用某种方法计算出各图像块的权重，然后按照公式(17)进行标签值概率值的计算。
4) 按照公式(18)得到标签的结果。
在本方法中提取的图像块为图像标签块，普通的非局部块加权方法采用常用的归一化相关系数函数(normalized correlation coefficient，NCC)衡量目标图像块与提取出的图谱图像块之间的相似度。稀疏表示的标签融合方法是构建目标图像标签块和图谱图像标签块，将图谱图像标签块作为过完备字典D=[Tp1,Tpi,Tpn]，通过公式(2)进行稀疏系数的计算，然后通过公式(14)、(15)、(16)得到最终的各图谱权重。然后按照公式(17)、(18)得到分割结果。
(17)
(18)
为了将两种分割方法进行融合，本方法设计了一个相似性度量函数。该函数包含两个部分，第一部分为Dice系数，第二部分为余弦距离。公式中相似性度量结果与Dice系数成反比，与余弦距离成正比。记目标体素x处预分割结果为R，非局部块加权方法的融合结果为R1，稀疏表示融合方法的融合结果为R2.度量规则表示如下：
cos(p(x),pi(x))i=1,2 .
(19)
式中：Dice(R，Ri)为非局部块加权(或者稀疏表示)融合方法的分割结果和预分割结果的Dice系数，d为常数；添加常数d是为了防止分母为零的情况出现；pi(x)表示非局部块加权的结果或者稀疏表示融合方法的结果以体素x为中心的一个图像块转化成的列向量，cos(p(x),pi(x))为两个列向量的余弦距离；μ和1-μ代表两者的权重；DRi值越大，说明相似度越低。
公式(20)为两种方法的融合公式。当非局部块加权方法的分割结果和稀疏表示的分割结果的标签值相等时，就将该标签值赋予体素x，当两种方法的分割结果不相同时，通过公式(19)，计算相似度结果。若DR1>DR2，说明R1与预分割结果R之间的差距比R2与R之间的差距大，将R2的分割结果赋予体素x，若DR1<DR2，将R1的分割结果赋予体素x.
(20)
3 实验与结果分析
3.1 数据介绍
实验选用的测试数据集为牛津大学发布的猕猴数据集。该数据集用3T的扫描仪采集了20只雄性猕猴的数据。数据集包括T1，静息态的fMRI。本实验使用T1数据。猕猴年龄分布为2.41岁～6.72岁(平均年龄为4.01岁)；重量分布为：4.35 kg～11.7 kg(平均重量6.57 kg)；体素分辨率为0.5 mm×0.5 mm×0.5 mm，TE为4.01 ms，TR为2 500 ms，TI为1 100 ms，翻转角度为8°.
3.2 评价指标
实验选用Dice系数评价分割质量，其计算方法如公式(21)所示。公式中的A表示目标图像的金标准(即专家分割结果)，B表示用分割方法得出的分割结果。‖表示各集合的体素个数，|A∩B|表示集合A与集合B对应体素标签相同的体素个数，|A|与|B|分别表示集合A与集合B中的体素个数。Dice系数取值范围为0到1，越接近1，说明准确率越高。
(21)
3.3 结果分析
本实验选取三个方法分别对海马体、纹状体、带状体这三个核团的分割结果与金标准进行比较。参与对比的方法有多数投票法(majority voting，MV)，联合标签融合法(LW joint)以及本文的方法(ISRLF).
图3(a)-(d)分别为目标图像和三个不同方法的分割结果。其中，3(a)为金标准分割结果，3(b)为MV方法的分割结果，3(c)为LW joint方法的分割结果，3(d)为ISRLF方法的分割结果。从图中可以看出，LW joint方法的分割结果最差，ISRLF方法的效果最好。
图3 三种不同方法的分割结果
Fig.3 Segmentation results of three different methods
表1显示了在牛津大学猕猴数据集上采用上述三种方法对选定核团的分割准确率。图4为3个方法分割结果Dice系数的盒图。从表1与图4可以看出，MV方法虽然对于海马体的分割准确率比LW joint的准确率高，但该方法稳定性很差，准确率波动很大，其海马体的准确率相对于平均值波动达到了10%以上。LW joint方法的准确率相对于平均值波动在10%以内，稳定性较好，但只有带状体的核团分割准确率超过了MV方法。本文提出的方法提高了三个核团的分割准确率，尤其是对纹状体的平均准确率相对于MV方法提高了8%左右。对带状体的分割准确率提高程度有限，相对于准确率的平均值波动在5%以内，稳定性在三个方法中最高，可以看出本文提出的方法与其他方法相比，有着较好的鲁棒性。
表1 三种方法的Dice系数计算结果
Table 1 The results of Dice coefficient calculated
by the three methods
脑部位置MV/%LW joint/%ISRLF/%海马体38.56±12.9431.18±9.6445.78±4.25纹状体67.56±4.0150.46±8.5475.53±3.21带状核33.83±9.6545.25±5.6146.97±3.68
图4 Dice系数
Fig.4 Dice coefficient results
图5为将三个核团的冠状面分割结果单独呈现的效果图。第一行至第三行分别为海马体、纹状体、屏状核。第一列至第四列分别为金标准、MV方法分割结果、LW joint分割结果、ISRLF分割结果。可以看出本方法的分割结果较其他方法有一定的提升。
图5 海马体、纹状体、屏状核横断面分割结果图
Fig.5 Transverseplane segmentation results of
hippocampus, striatum, and claustrum
4 结束语
本文提出了一种改进的稀疏表示多图谱分割方法，然后应用于猕猴大脑皮下核团的分割。该方法以稀疏表示多图谱分割方法为基础，将标签信息二值化后引入到图像块中，以增加图像块的先验信息；然后通过改变信息熵的计算方式改进互信息，将其用来衡量目标图像的目标体素所在平面和线性配准后的图谱中对应平面的全局相似度，使融合时各图谱的权重更加合理。后提出一个相似性指标将非局部块加权融合的分割结果与稀疏表示的融合分割结果进行综合，从而得到最终的分割结果。通过与其他方法得到的结果进行Dice系数的对比，证明本文提出的方法获得了更高的准确率和更好的鲁棒性。
[1] QUALLO M M,PRICE C J,UENO K,et al.Gray and white matter changes associated with tool-use learning in macaque monkeys[J].Proceedings of the National Academy of Sciences,2009,106(43):18379-18384.
[2] BAUMAN M D,SCHUMANN C M.Advances in nonhuman primate models of autism:integrating neuroscience and behavior[J].Experimental Neurology,2018,299:252-265.
[3] WANG P Y,WANG D,ZHANG J,et al.Evaluation of submillimeter diffusion imaging of the macaque brain using readout-segmented EPI at 7 T[J].IEEE Transactions on Biomedical Engineering,2019,66(10):2945-2951.
[4] YOUNG J T,SHI Y,NIETHAMMER M,et al.The UNC-Wisconsin rhesus macaque neurodevelopment database:a structural MRI and DTI database of early postnatal development[J].Frontiers in neuroscience,2017,11:29.DOI:10.3389/fnins.2017.00029.
[5] BALBASTRE Y,RIVIÉRE D,SOUEDET N,et al.Primatologist:a modular segmentation pipeline for macaque brain morphometry[J].Neuroimage,2017,162:306-321.
[6] SEIDLITZ J,SPONHEIM C,GLEN D,et al.A population MRI brain template and analysis tools for the macaque[J].Neuroimage,2018,170(SI):121-131.
[7] REVELEY C,GRUSLYS A,YE F Q,et al.Three-dimensional digital template atlas of the macaque brain[J].Cerebral Cortex,2017,27(9):4463-4477.
[8] ROHLFING T,KROENKE C D,SULLIVAN E V,et al.The INIA19 template and NeuroMaps atlas for primate brain image parcellation and spatial normalization[J].Frontiers in Neuroinformatics,2012,6(8):27.DOI:10.3389/fninf.2012.00027.
[9] MILHAM M P,AI L,KOO B,et al.An open resource for non-human primate imaging[J].Neuron,2018,100(1):61-74.e2.
[10] 陈伟导.基于7T磁共振图像的猕猴脑组织及脑解剖学结构自动分割方法研究[D].杭州：浙江大学，2018.
[11] HUO Y K,XU Z B,ABOUD K,et al.Spatially localized atlas network tiles enables 3D whole brain segmentation from limited data[C]∥International Conference on Medical Image Computing and Computer-assisted Intervention.Springer,Cham,2018,11072:698-705.
[12] ZHAO G Y,LIU F,OLER J A,et al.Bayesian convolutional neural network based MRI brain extraction on nonhuman primates[J].Neuroimage,2018,175:32-44.
[13] ZHANG D Q,GUO Q M,WU G R,et al.Sparse patch-based label fusion for multi-atlas segmentation[C]∥International Workshop on Multimodal Brain Image Analysis.Springer,Berlin,Heidelberg,2012:94-102.
[14] YAN M,LIU H,SONG E M,et al.Sparse patch-based representation with combined information of atlas for multi-atlas label fusion[J].IET Image Processing,2018,12(8):1345-1353.
[15] YAN M,LIU H,XU X Y,et al.An improved label fusion approach with sparse patch-based representation for MRI brain image segmentation[J].International Journal of Imaging Systems and Technology,2017,27(1):23-32.
[16] 夏瑞，马瑜，王文娜,等.基于重采样改进的多图谱分割算法[J].计算机工程与设计,2018,39(8):2587-2592.
XIA R,MA Y,WANG W N,et al.Multi-atlas segmentation algorithm based on resampling[J].Computer Engineering and Design,2018,39(8):2587-2592.
[17] LIU Y,WEI Y,WANG C Y.Subcortical brain segmentation based on atlas registration and linearized kernel sparse representative classifier[J].IEEE Access,2019,7:31547-31557.
