【作 者】王吉喆，严壮志,，温军玲

1 上海大学通信与信息工程学院，上海市，200444

2 上海大学上海生物医学工程研究所，上海市，200444

随着世界人口老龄化的加剧，阿尔茨海默症（Alzheimer's Disease，AD）患者逐年增加。AD作为一种起病隐匿的进行性发展神经系统退行性疾病，目前尚无特效治疗方法。所以只能通过早诊断，早治疗延缓病症。通过磁共振图像（Magnetic Resonance Imaging，MRI）测量大脑中海马体的萎缩程度，是诊断AD患者或者轻度认知功能障碍患者的最有效的方法之一[1]。所以针对海马体的形态学分析在临床诊断和监控疾病上具有十分重要的意义。海马体的轮廓提取作为形态学分析的重要步骤，在临床上以人工逐层分割为主，完成一组海马图片的手动勾勒的时间在2个小时左右[2]。人工分割具有时间效率低，可重复性差等缺点，而且逐层分割失去了数据的空间连续性等三维空间的性质。直接进行三维分割，又由于海马体复杂的边界形状，其计算代价较高。因此，研究适用于海马体三维分割的高效算法是十分必要的。

近几年，学者们提出很多针对海马体的分割算法。一种是基于图谱的方法，Tong等[3]使用多通道大形变微分同胚度量图在配准模板和待分割图像之间寻找微分同胚对应关系。Kichang等[4]使用单图谱的方法对海马体进行粗分割，并将分割结果作为下一步使用图像切割的方法进行细分割的初始化和先验信息。Wu等[5-6]通过融合不同尺度的图谱图像块，来增加图谱与待分割图像的相似度。Suppa等[7]已经开始尝试将基于图谱的全自动分割、测量方法用于临床辅助诊断，并证明该方法在预测轻度认知功能障碍转化为AD患者的过程中表现良好。另一种方法是基于偏微分方程的形变模型，但由于其仅依赖于从图片中提取的信息，存在着一些不足。而先验知识的引入可以弥补这些不足。最近研究学者们提出了名为边界梯度分配（Gradient Distribution on Boundary，GDB）的方法。通过机器学习建立起一个区域权重图。这张权重图成为GDB[8]，GDB被用作在活动轮廓模型演化时候的先验信息。并且通过文献[9]的扩展，GDB模型成为了基于活动轮廓模型的边界演化的自适应边界梯度分配（Adaptive Gradient Distribution on Boundary，AGDB）。文献[10]构建了三相GDB模型来解决迭代次数与细节分割的矛盾。

基于形变模型的方法具有良好的空间扩展性，可以将成熟的模型直接扩展至三维空间。但求解形变模型中的偏微分方程，是一个十分耗时的过程。格子玻尔兹曼（Lattice Boltzmann，LB）方法作为一种快速求解偏微分方程的方法，其模型经过离散化处理，求解过程简单，易于进行编程实现；本身具有三维模型，可以直接与三维形变模型相结合，加速计算求解；计算求解只与其周围像素有关，可以通过并行化，进一步加速计算，所以便有学者将其应用于图像处理领域[11-12]。

但在处理三维数据时，若仅仅使用中央处理器（Central Processing Unit，CPU）进行串行计算，不能充分发挥LB方法的并行性，所以仍存在计算耗时过长的问题。作为专门为了处理图像显示而设计的图形处理器（Graphics Processing Unit，GPU），由于其设计的初衷就是能够并行处理大量数据，利用GPU进行通用并行计算便应运而生。CUDA（Compute Unified Device Architecture）作为一种GPU并行计算开发环境，以C语言为基础，对C语言进行扩展，使得开发人员不需要去学习特定的显示芯片指令或者特殊的结构就能完成GPU并行程序的编写。

CUDA一经推出便受到了很多研究LB学者的青睐，Balla-Arabé等[13]和Jones等[14]将LB方法与GPU并行计算相结合，加速了图像处理的速度。这充分地显示出LB方法与GPU并行计算的匹配性。所以结合LB方法的算法优势和GPU的硬件平台优势，可以大大加快处理三维数据的速度。本文利用LB方法结合现有的形变模型与并行计算平台，实现对海马体的三维快速分割。

1 基于先验信息的三维LB模型

1.1 三维LB分割模型

LB方法是近二十年来发展起来的一种系统建模和模拟的方法，由于其对偏微分方程的快速求解和图像处理领域相关方法的发展，已在图像去噪、分割、修复、增强等方面取得了显著的成就。本文利用三维LB模型，结合脑部MRI数据，完成对海马体的快速分割。

图1 三维格子玻尔兹曼的空间模型Fig.1 Three dimensions lattice Boltzmann model

三维LB的离散空间模型如图1所示。模型中的每个节点都可以看作是三维数据中的体素，每个节点中的粒子都可沿着矢量方向迁移至相邻节点，即LB模型中的演化过程；随后对各个节点的粒子进行密度重分配，即LB模型中的碰撞过程。

在有源项（即外力项）存在的情况下，LB模型的一般方程为：

其中fa(r+ea△t, t+△t)表示在t+△t时刻在节点r+ea.△t上的粒子分布函数，ea是粒子运动的方向，△t是时间步长，为松弛因子，(r, t)是平衡态分布函数，Fa(r, t) 是有源项。对于三维LB来说，本文设定的各个适量方向上的速度ea为

通过对等式(1)左边使用泰勒展开，等式右边使用Chapmann-Enskog展开，可写作下面的形式：

其中，F为外力项，D为扩散张量

ωα为粒子速度函数，称为平衡态加权系数，当α=0时ω=，当α=1, 2, ...,6时ω=。αα

脑部MRI中，海马体由于其体积较小，且与周围的脑灰质部分边界不是很明显，在分割的时候经常出现过分割的情况。所以本文将边缘信息、局部信息和先验形状信息引入至演化方程的外力项，来弥补分割精度不足的问题。

1.2 基于边缘信息和局部信息的外力项

根据文献[15]，可得D=g(│△I│)，则方程可改写为

F1=αg为分割提供边缘信息，其中g(│△I│)是边缘停止函数，定义为：

其中，│△Gσ*I │是经过高斯平滑的图像梯度，Gσ是高斯核函数，σ是高斯函数的方差。

F2=β(λ1e1+λ2e2)为分割提供了区域信息，λ1和λ2是演化曲线的内部和外部区域信息权重。其中e1，e2定义为：

其中，H(.)是Heaviside函数，本文采用平滑函数Hε近似为：

1.3 基于先验形状信息的外力项

在使用先验形状分割模型时应考虑到个体差异性的存在，需要针对先验信息进行配准和融合，获得一个具有普适性的先验形状信息，随后将先验信息配准至待分割图像，作为分割所使用的外力项之一。本文采用基于最大互信息的配准方法，利用Powell算法进行寻优。将所得20组配准结果进行加和平均，求得融合先验信息，结果如图2所示。

先验形状的描述方法有很多种，其中，距离符号函数能够将形状的轮廓与主动轮廓模型进行融合[16]，在变分法的水平集中有广泛的应用。本文采用欧式距离作为相似度衡量指标，构建的形状先验外力项，即：

图2 配准前后位置的对比Fig.2 Comparison of the position before and after registration

综上，本文使用的LB模型的方程为：

其中α，β，μ是有源项的权重系数，可根据不同的图片进行调整。

2 三维LB分割算法的CUDA实现

基于CUDA的程序设计流程为：首先在CPU端和GPU端进行内存的分配，初始化CPU端相关数据，并将初始化数据复制到GPU端中；设置grid和block的大小，并根据核函数，在GPU端完成并行计算过程；最后，在当计算满足条件的时候，计算结束，并将数据由GPU端复制回CPU端，完成并行计算过程。

由于LB方法的计算只与其周围网格有关，每个格点的计算相对简单，因此影响程序性能的主要原因就是数据的复制。所以，应尽可能地将计算所需数据经由一次复制过程复制到GPU端。其次便是资源的分配。由于并行线程的数量是不可能无限增加的，所以，如何根据具体情况进行线程的分配，也是并行计算的重中之重。

2.1 算法流程

基于上述的模型，本文提出了三维LB分割算法。为了充分发挥LB方法并行化的特性，使用了GPU平台对运算过程进行加速。算法流程图见图3。

步骤1：设置初始化距离函数φ；

步骤2：设置参数α，β，μ，λ1，λ2，根据式(6)计算边缘信息；

步骤3：将边缘信息、距离函数、待分割图像及参数从CPU端复制到GPU端；

图3 三维格子玻尔兹曼分割算法单GPU平台实现流程图Fig.3 Flow chart of three dimensions lattice Boltzmann model on GPU platform

步骤4：根据式(7)计算区域信息，根据式(9)计算先验信息；

步骤5：计算外力项等其他必要参量；

步骤6：根据式(10)执行碰撞迁移过程，并且更新距离函数的值；

步骤7：由CPU判断是否达到迭代步数，若没有达到，则重新启动核函数，返回步骤4；若已经达到，则将计算得到的新的距离函数从GPU端复制回CPU端，进行边界处理，完成分割过程。

GPU特化了计算能力的代价就是做分支判断的能力很低，所以一般情况下，分支判断的部分都交由CPU端运行。由于边界处理涉及到很多判断过程，所以在计算的时候，将边界处理放在了完成迭代之后，交由CPU一并完成。

双GPU平台能够解决因内存固化而导致的显存不足的问题，但是又出现了一个新的问题：数据集的拆分。三维LB算法存在着层与层之间的数据交换，如果只是简单的进行上下拆分，那么在计算的过程中，位于边界上的两层将无法获得来自上面一层（下面一层）扩散来的粒子，造成计算错误。而对于迭代求解的方法，错误是会累积的，所以需要将GPU1中向下扩散的部分搬运至GPU2中，完成GPU1中的最后一层向GPU2中第一层扩散的过程。但是CUDA无法支持GPU对GPU通信，所以需要CPU作为中介，协助两组GPU完成扩散过程。双GPU平台算法流程见图4。

相较于单GPU，双GPU在执行值步骤6时，会分别取GPU1中的最后一层的数据向下扩散的缓存，和GPU2中的第一层数据上扩散的数据缓存进行交换，完成迁移过程。

3 验证实验结果与讨论

3.1 验证实验设计

图4 三维格子玻尔兹曼分割算法双GPU平台实现流程图Fig.4 Flow chart of three dimensions lattice Boltzmann model on dual GPU platform

为了验证三维LB分割算法的有效性，本文将该算法与AGDB方法[11]进行对比。AGDB方法将自适应边界分布图和先验知识引入主动轮廓模型，对海马体进行三维分割。本文的实验硬件平台CPU为Intel®Core™ i7-4790K @ 4.00 GHz，GPU为两块NVIDIA GeForce GTX 660，8 GB内存。软件平台为Microsoft Visual Studio 2012，CUDA版本为6.5。实验图像由 ADNI(The Alzheimer's Disease Neuroimaging Initiative, ADNI)数据库获得，EADC-ADNI 工作组有国际知名的研究海马体勾画的专家团体，该组织定义了一个关于海马体勾画的标准协议[17]，为诊断研究、临床试验和算法验证工作提供了统一的标准。为了保证实验具有统计学意义，本文选取了20组AD患者的脑部MRI图像进行分割，每组图像尺寸均为166×256×30，且每组图像都有相对应的、由专家进行勾画的海马体的金标准。

具体的实验步骤如下：

（1）取20组脑部MRI海马体分割的金标准，进行图像配准及融合，获得一组166×256×30的海马体的先验信息，如图5所示。

图5 融合先验信息的三维重建Fig.5 3-D reconstruction of fuse prior information

（2）使用上述先验信息，对20组脑部MRI图像进行分割，并按照不同情况，设置α，β，μ，λ1，λ2，进行迭代分割。

（3）使用同样的参数，分别在CPU平台、单GPU平台、双GPU平台上进行计算。

（4）针对实验结果进行评估实验。

3.2 评估方法

算法的评估包括对分割准确性的评估和时间效率的评估。为了定量比较算法分割的准确性，本文使用目前在衡量分割精度之中最常采用的Dice相似系数(Dice Similarity Coefficient，DSC)[18]和豪斯多夫距离(Hausdorあ Distance，HD)进行定量比较。时间效率的评估标准为算法达到稳定后所需要的时间和迭代步数。

DSC可以衡量分割结果与金标准之间的重叠率，其计算公式为：

其中：A是分割算法获得的结果，M是金标准，运算符∩计算A与M的重叠的像素数量，│A│+│M│计算A与M中，全部像素的数量。当两个区域完全重叠的时候，DSC=1；当两个区域完全不重叠的时候，DSC=0。

HD计算结果为两个轮廓曲面之间最接近点的最大值，其计算公式为：

其中，A和B表示两个曲面，h(A, B)表示曲面A上的所有点到曲面B上的所有点的距离的最小值，h(B, A)表示曲面B上的所有点到曲面A上的所有点的距离的最小值，这两者的最大值，即为豪斯多夫距离。HD越小，说明曲面A与曲面B越接近。

3.3 实验结果

由于迭代次数会影响到二者的分割精度，在进行三维LB算法分割质量的评估之前，应确定最佳迭代方案。所以本文根据DSC随时间变化的情况，来确定分割的最佳迭代次数。

图6记录了DSC随迭代次数变化的曲线。从图中可以得知，两种方法的DSC均随着迭代次数的增加而增加，并逐渐达到稳定。但是三维LB分割算法在迭代30次左右，DSC就可以达到一个相对稳定的值。相较于AGDB方法的70次迭代，三维LB分割算法能够在迭代次数更少的前提下获得一个较好的分割效果。

图6 三维格子玻尔兹曼分割算法和AGDB方法与金标准的Dice相关系数随迭代次数变化曲线Fig.6 The curve of DSC based on standard for three dimensions lattice Boltzmann method and AGDB

使用得到的最佳迭代次数后，使用本文提出的三维LB分割算法，在20名AD患者的脑部MRI图像上进行分割结果精度的验证。结果见图7。图中白色部分均为金标准。红色部分分别为三维LB分割算法和AGDB方法的分割结果。从图7可以很直观地看出三维LB分割算法的分割精度要高于AGDB方法。而表1记录的两种方法的平均DSC和HD，也从定量的角度上印证了三维LB分割算法分割的准确程度，要高于AGDB方法。

图7 三维格子玻尔兹曼分割算法与AGDB方法分割结果对比Fig.7 Segmentation results between 3-D LB and AGDB

表1 三维格子玻尔兹曼分割算法和AGDB方法的分割结果与金标准的Dice相关系数与豪斯多夫距离Tab.1 The result of DSC and HD for three dimensions lattice Boltzmann method and AGDB with standard

三维LB分割算法分割效率的评估是建立在最优资源分配的前提下，线程块（block）大小对计算效率的影响。线程块是GPU进行计算时的获得硬件资源的最小单位。所以以一组AD患者脑部MRI图像为例，输入实际分割时所使用的参数，进行十次迭代，只计算GPU端核函数的计算时间，以观察block大小对计算效率的影响，结果如表2所示。

表2 不同block大小对同规模并行计算所需时间的影响Tab.2 The eあect of block sizes in parallel computation

从表2中可以看出，在一定范围内，随着block的变大，计算时间在减少，说明硬件资源在逐步饱和，而当block的大小达到256时，硬件资源达到饱和，计算时间不会进一步缩短。这也与多数程序中推荐的block维度相同。

在上述实验的基础上，使用已经达到最优化的block参数，针对同一组脑部MRI进行分割，记录达到最大DSC时，分割所需要的时间，结果如表3所示。

表3 不同平台下三维格子玻尔兹曼方法和AGDB方法完成分割所需时间Tab.3 The time required to complete the segment between AGDB and three dimensions lattice Boltzmann method on diあerent platform

从表3中我们可以看出，在完成分割的时候三维LB分割算法的速度明显优于AGDB方法，而且在使用GPU并行加速后，计算速度又有了进一步的提升。说明三维LB分割算法在时间效率上要优于AGDB方法，并且，在GPU并行计算加入后，计算效率提升了一个数量级，说明了LB方法与GPU并行计算有着很好的契合性。

但是同时发现，在使用双GPU计算的情况下，计算效率反而比单GPU计算时的效率要低。这是因为在使用双GPU在计算的时候，涉及到层间信息交换的过程，这个过程是交由CPU来完成的。数据的GPU端与CPU端互相复制和CPU的串行执行，是导致双GPU计算效率低于单GPU计算效率的主要原因。虽然存在着这样的弊端，但是双GPU计算的优势在于更大的显存和并行计算规模，所以在针对更大数据量的计算时，双GPU可以完成单GPU所无法完成的工作。

3.4 实验结果讨论

本文提出了结合先验信息的LB三维算法，并且证明该模型在分割弱边缘，低对比度的海马体MRI图像时，表现良好。并且通过GPU平台的加入计算，速度得到了明显的提升。原因分析如下：

（1）先验信息作为约束项的加入，使得模型在针对图像边界模糊，待分割区域与背景灰度相似等复杂背景图像进行分割时，有着较强的准确性和鲁棒性。特别是针对海马体体积较小，边缘模糊，组织对比度低，结构复杂的器官，具有较强的分割能力。

（2）LB方法的求解过程只与该像素点周围的区域有关，自身带有可并行化的特点，通过GPU平台的使用，大大减少了计算耗时。

（3）通过对block中线程数目的优化，达到最佳的硬件资源分配；通过对算法流程分析，减少两端复制过程，能够让GPU在不被打断的情况下，完成高密集度的计算过程，进一步提升了计算效率。

4 结论

在脑部MRI中能够快速、准确地分割海马体的结构，对于AD等疾病的临床诊断具有重大意义。本文提出的三维LB分割算法，直接在三维空间中对海马体进行分割，避免了逐层分割带来的层间信息损失和三维重建的过程。并针对脑部MRI中，海马体成像边界模糊、不连续等问题，引入了先验形状等信息作为三维LB分割算法的外力项，使得活动曲面能够快速、准确地演化至海马体的边界，获得海马体的三维轮廓信息。在实验中，通过与AGDB方法的比较，三维LB分割算法在DSC、HD和计算速度等指标上，均优于AGDB方法。而针对三维分割带来的耗时问题，根据LB方法自身可高度并行化的特点，通过对算法流程的优化，在单GPU平台实现了算法的并行计算，进一步提高了计算速度，解决了三维数据带来的计算代价高的问题。为了应对可能出现的板载显存不足的问题，针对双GPU平台计算的特点，重新设计了算法的实现过程，使得LB方法能够在并行平台上进行更大规模的计算。进一步扩展了LB方法在快速图像处理领域中的应用前景。

[1] Dubois B, Feldman H H, Jacova C, et al. Research criteria for the diagnosis of Alzheimer's disease: revising the NINCDS-ADRDA criteria[J]. Lancet Neurol, 2007, 6(8): 734-746.

[2] Dill V, Franco A R, Pinho M S. Automated methods for hippocampus segmentation: the evolution and a review of the state of the art[J]. Neuroinformatics, 2015, 13(2):133-150.

[3] Tong T, Wolz R, Coupé P, et al. Segmentation of MR images via discriminative dictionary learning and sparse coding: application to hippocampus labeling[J]. Neuroimage, 2013, 76(1):11-23.

[4] Kichang K, Uicheul Y, Dong-Kyun L, et al. Fully-automated approach to hippocampus segmentation using a graphcuts algorithm combined with atlas-based segmentation and morphological opening[J]. Magn Reson Imag, 2013, 31(7):1190-1196.

[5] Wu G, Shen D. Hierarchical label fusion with multiscale feature representation and label-specif i c patch partition[C]. Int Conf Med Imag Comput Computer-Assisted Intervent, 2014:299-306.

[6] Wu G, Kim M, Sanroma G, et al. Hierarchical multi-atlas label fusion with multi-scale feature representation and label-specific patch partition[J]. Neuroimage, 2015, 106: 34-46.

[7] Suppa P, Hampel H, Spies L, et al. Fully automated atlas-based hippocampus volumetry for clinical routine: validation in subjects with mild cognitive impairment from the ADNI cohort[J]. J Alzheimers Dis, 2015, 46(1):199-209.

[8] Zarpalas D, Zafeiropoulos A, Daras P, et al. Brain structures segmentation using optimum global and local weights on mixing active contours and neighboring constraints[C]. Int Symp Appl Sci Biomed Commun Tech, ACM, 2011: 127-132

[9] Zarpalas D, Gkontra P, Daras P, et al. Hippocampus segmentation through gradient based reliability maps for local blending of ACM energy terms[C]. Int Symp Biomed Imaging, IEEE, 2013: 53-56.

[10] Zarpalas D, Gkontra P, Daras P, et al. Gradient-Based Reliability Maps for ACM-Based Segmentation of Hippocampus[J]. IEEE Trans Biomed Eng, 2014, 61(4): 1015-1026.

[11] 林笑曼, 严壮志, 蒋皆恢,等. 基于三维格子玻尔兹曼模型的海马结构MRI快速分割[J]. 生物医学工程学进展, 2014, 35(1): 1-7.

[12] Wen J, Jiang J, Yan Z. A new lattice Boltzmann algorithm for assembling local statistical information with MR brain imaging segmentation applications[J]. Multidimen Syst Sign Proc, 2016:1-17.

[13] Balla-Arabé S, Gao X, Wang B. GPU accelerated edge-region based level set evolution constrained by 2D gray-scale histogram[J].IEEE Trans Imag Proc, 2013, 22(7): 2688-2698.

[14] Jones B D, Feng Y T. Effect of image scaling and segmentation in digital rock characterisation[J]. Comput Part Mech, 2016, 3(2):201-213.

[15] Wen J, Yan Z, Jiang J. Novel lattice Boltzmann method based on integrated edge and region information for medical image segmentation[J]. Biomed Mater Eng, 2014, 24(1):1247-1252.

[16] Bresson X, Vandergheynst P, Thiran J P. A Variational model for object segmentation using boundary information and shape prior driven by the Mumford-Shah functional[J]. Int J Comput Vis, 2006,68(2):145-162.

[17] ADNI 数据库[R/OL]. http://www.hippocampal-protocol.net/SOPs/LINK_PAGE/Harmonized Protocol ACPC User Manual_biblio.pdf.

[18] Cabezas M, Oliver A, Lladó X, et al. A review of atlas-based segmentation for magnetic resonance brain images[J]. Comput Meth Program Biomed, 2011, 104(3): 158-177.