基于Laplace方程的脑皮层厚度测量方法

邰滢滢，梁永霞

(辽宁大学 信息学院，辽宁 沈阳 110036)

0 引言

大脑是人体最重要、最复杂的器官，它控制着人体各组织器官功能的正常运作.大脑皮层是大脑最外层的薄皮状结构，是人脑最高级的中枢，因此很多疾病的产生也会通过大脑皮层脑部形态学特征变化体现出来.通过对大脑皮层形态特征的变化规律进行研究，可以高效地诊断和治疗疾病，甚至进一步预防疾病的发生，这对于现代医学发展有着十分重要的意义.人眼视力的好坏也是与大脑皮层结构有着密不可分的关系的.视觉剥夺会导致大脑皮层形态结构的变化，研究发现这些变化既包括局部脑结构的改变，如皮层厚度、灰质体积等，也包括脑结构网络的改变.针对这种情况，如果能够发现近视者早期特征和脑部形态特征变化规律，不仅对视力下降者早期治疗非常重要，也对预防近视起着决定性的作用.先天盲人和正常人在脑皮层形态学上都存在着某些差别.通过研究这两类人群的脑部差异，能够有助于盲人的康复治疗，甚至能够使盲人通过其他的直觉补偿，来得到更好的生活体验.因此，对大脑皮层进行形态学分析有助于发现疾病病理特征和发病规律.大脑皮层形态学指标有很多种，其中脑皮层厚度就是形态学研究中的一个重要指标.通过估计整个大脑皮层的厚度信息，可以发现视力剥夺引起的大脑皮层的形态学变化.

1 相关知识介绍

脑皮层厚度是指灰质皮层外表面到白质皮层外表面之间的距离，皮层厚度的测量方法主要分为基于网格表面测量和基于体数据测量两大类[1].在基于体数据的测量方法中有由Rosenfeld和Pfaltzc[2]提出的基于距离变换的方法.这种方法在进行距离变换时是在二维图像切片上实现的，所以在某些程度上会使得测量结果产生偏差.Lohmann等[3]考虑了二维图像距离变换时会丢失部分三维信息，导致大脑组织结构的特征变化，其算法在三维结构上进行，保证了拓扑不变性，使结果更合理.Jones等[4]提出了一种新的皮层厚度定义方法，对脑皮层构建势能场模型，并且用拉普拉斯方程来描述这个势能场，通过求解这个偏微分方程，获取各个等势面，再根据等势面计算出场线.场线位于皮层内表面和外表面之间，与所有经过的等势面都垂直，这样我们就可以把皮层内过任一点的场线的长度定义为该点的厚度.尽管这种求解二阶偏微分方程的方法计算量很大，但这种方法是目前较为科学的一种测量方法.

基于网格表面的测量方法中有一种比较常见的方法由Fischl等[5]在2000年提出.该方法具有良好的鲁棒性.其原理是基于大脑皮层曲面的网格顶点，计算了一次外表面(灰质)到内表面(白质)顶点的最短距离，一次内表面(白质)到外表面(灰质)网格顶点的最短距离，然后求两次最短距离的平均值作为该点的皮层厚度.这种测量方法的基础在于从分割的MRI图像重建出三维网格皮层表面.因此这类方法测得的结果与网格重建和图像分割等过程紧密相关.例如求图1中A点处的皮层厚度，首先通过三维重建得到网格化的脑皮层表面，再基于网格顶点计算每一点的皮层厚度.具体步骤是先计算A到B的最短欧氏距离，再计算B到C的最短欧氏距离，取AB+BC的平均值就得到A点的皮层厚度值.

图1 脑皮层厚度算法示意图

2 基于网格的脑皮层厚度测量

基于网格的脑皮层厚度测量的基础是从MRI二维图像上重建出网格状的三维皮层表面，再由Free Surfer直接计算出网格上每一点的皮层厚度.这其中需要经过一系列的预处理步骤，如图2所示.

图2 图像预处理过程图

2.1 图像预处理

图像处理包括预处理和迭代计算两部分.在进行分析之前要经过一系列的预处理，图像的预处理包含很多过程，对于不同的后处理其对应的预处理过程也不相同.总的来说，图像预处理大致包括以下几步：去噪、配准、不均匀场校正、组织分割、内外皮层曲面重建等.这些预处理过程都是在医学图像处理中非常热门的研究方向.在本文中主要是大脑皮层的形态学指标的计算和统计，因此，预处理的质量决定着形态学指标测量的准确度.图像的预处理步骤如下：

2.1.1 去噪

图像去噪是指减少数字图像中噪声的过程.现实中的数字图像在数字化和传输过程中常受到成像设备与外部环境噪声干扰等影响，称为含噪图像或噪声图像.图像中噪声的来源有许多种，这些噪声来源于图像采集、传输、压缩等各个方面.在这里的去噪主要是去除核磁共振成像过程中产生的噪声信号和人脑其他组织信号，这些信号对后期进行大脑皮层的处理没有作用，甚至会造成误差影响，因此需要去除掉这些噪声信号.

2.1.2 配准

配准是一个类似于标准化的处理，将个体的核磁共振图像通过平移、旋转、缩放使其与标准模板对齐，并最大化与模板相似的一个过程.通过配准，可以得到一个与模板图像一一对应的个体图像，通过模板又可以得到不同个体图像间的对应关系，从而使得不同空间内的图像之间相互比较成为可能.

2.1.3 不均匀场校正

在核磁共振成像的过程中，磁场强度不能完全保证均匀一致，造成了图像灰度值的不均匀性.同一脑组织(灰质、白质、脑脊液)的灰度值和对比度在不同区域相差很大.这对于主要依靠灰度值或对比度来进行分割算法的正确率有很大的影响，进而直接影响到后续的处理和分析.因此，必须对不均匀场进行校正.

2.1.4 脑组织分割

脑组织的分割是预处理中必要的步骤，对于后续计算和分析有重要作用.脑组织分割一般首先把大脑分割成指定的几类，如灰质、白质和脑脊液三类.首先根据图像的灰度值进行简单的预划分，将体素按照灰度值进行分类.选取出那些最可能属于某一类的体素.对于无法从灰度上区分的剩余体素，再利用大脑几何信息来进行精细分类.

2.1.5 曲面重建

在完成脑组织分割之后，我们就可以根据分割结果对大脑灰质和白质外边界进行三维曲面重建.曲面重建就是要恢复出隐藏在三维体数据中的脑皮层结构的几何信息.因为从一组二维断层图像上理解病人身体内部器官的三维结构是非常困难的，即使是经验丰富的医生也很难准确地描述器官的立体结构和确定病变体的方位.另外，三维曲面重建技术的出现也给医学教学领域带来了很大的帮助，学生可以通过三维重建获得直观立体的人体器官和组织图像，形象地理解和掌握各个器官和组织的构成与功能以及常见病变的发生情况，提高学习的质量和效率.

3 基于体数据的脑皮层厚度测量

基于体特征的脑皮层厚度估计方法是在三维立方体体素网格上进行运算的，并通过求解拉普拉斯方程来得到脑皮层厚度[6].这种方法在计算前也是需要进行预处理的.因为大脑是三维结构的，直接对二维图像进行分析会丢失部分三维信息，所以要预处理转化为三维体素结构来进一步测量.图3为基于体数据的皮层厚度测量流程图.在基于体数据的测量中，本文的核心思想是在尽可能准确分割的基础上，借助二阶偏微分方程测量大脑皮层厚度.

图3 数据处理流程图

3.1 图像分割

由于大脑形态结构极其复杂，仅使用一种分割方法来分割图像并不能得到较好的分割结果，因此本文作者将多种分割方法结合起来，发挥各自的优点，中和单一分割方法带来的不足.首先根据图像的灰度值进行简单的预划分[9]，将体素按照灰度值进行分类.选取出那些最可能属于某一类的体素.对于无法从灰度上区分的剩余体素，再利用大脑几何信息来进行精细分类.本文测量方法的流程图如图3所示.

3.2 脑皮层厚度测量

通过一系列预处理提取出脑皮层区域，并定义出皮层的内外边界，根据数学原理用拉普拉斯方程法测量厚度，主要包括以下几步：在脑皮层边界建立拉普拉斯方程得到势能场；计算单位向量场；计算皮层厚度.传统体数据厚度的测量采用拉普拉斯方程求解方法，在整个三维脑皮层空间内进行迭代求解，迭代公式为

(1)

ψi+1(x，y，z)是第i次迭代时坐标(x，y，z)处像素的势能值.当求解得到ψ之后，可以由式(2)求解得到单位向量场线T，其中∇为梯度算子.

(2)

然后，就可以用偏微分方程计算皮层厚度[11].假设皮层区域内每个点的单位向量场线都是由内表面S指向外表面S′，则-T表示沿同一条场线从外表面S′指向内表面S.这里定义两个长度函数L0和L1，L0表示从内表面到某点的场线长度，L1表示从外表面到某点的场线长度.

根据微分几何，L0和L1都满足下面的一阶线性偏微分方程

(3)

那么任一点x处的脑皮层厚度值即为穿过该点的场线长度L0和L1的和：

L(x)=L0(x)+L1(x)

(4)

其中L(x)就是点x处的皮层厚度值.

4 实验与分析

4.1 实验环境

本实验运行的计算机配置为2.4 GHz Intel Core i3和4 GB RAM.通过上面的一系列预处理步骤，得到一个基于三角网格的大脑皮层表面，实验用FreeSurfer软件来进行MRI二维断层图像的预处理.FreeSurfer软件是用于MRI图像处理和分析的一款开源软件，可以进行多种模态数据的预处理、皮层重建等(http：//surfer.nmr.mgh.harvard.edu/)，该软件只能在Linux系统下运行，所以选择在Windows系统下安装虚拟机，在虚拟机中安装Ubuntu，然后将下载好的FreeSurfer安装包解压并安装.然后计算出网格上的脑皮层顶点的厚度.本文中用于实验的研究对象有先天性盲人组20例，年龄范围为32～45岁，正常组20例，年龄范围为35～49岁，均通过核磁共振扫描得到.选取的两组对象的年龄、性别、智力状态都相差不大，也没有其他的疾病.部分相关信息如表1所示.对这40张二维图像进行预处理得到三维数据.本文在分割步骤中使用SPM5软件进行两次粗分割和细分割，实现对所有实验样本的脑组织进行分割，然后在提取出的大脑皮层上建立拉普拉斯方程，依据3.2节的理论公式用迭代法求解方程，得到脑皮层厚度值.

表1 早期盲人样本资料

(5)

4.2 实验结果分析

经过对两组人群脑皮层厚度的统计发现，早期盲人的视觉区皮层厚度总体上要比正常人视觉区的厚，如图4所示，基于Laplace方程得到的皮层厚度统计，横坐标是大脑皮层不同部位的学术名称，纵坐标是皮层厚度值.从图中可以看出，视力剥夺对于直接处理视觉信息的枕叶脑区影响尤为显著，而且研究发现整个脑区在形态学结构上都与正常组有所不同，这也说明盲人在视觉区受到视力剥夺，为了弥补视力的不足就会在其他功能区域上得到强化[13].

图4 大脑视觉区脑皮层统计图

同时，本文分别使用两种不同的计算方法得到正常组和早期盲人组的皮层厚度信息，表2是两种方法在脑皮层不同部位测得厚度的均值，用同样的检验方法进行了比较，得到对应的差异值q，结果显示，在正常组和盲人组之间的差异中，基于网格的方法对应的q值为0.026 5，基于体数据的方法对应的q值为0.031 8.可以发现在基于体数据的方法中盲人和正常人之间的差异更大，这也说明本文准确分割的体数据方法有更高的测量精度.

表2 两种方法在不同部位的皮层厚度均值

5 结束语

本文比较了基于网格和基于体数据的两种计算脑皮层厚度的方法，并对早期盲人和正常人的皮层厚度进行了统计，结果表明基于体数据的测量方法具有更高的精度，但本文实验所采集的被试样本类别较少，后续可以依据不同的阶段来进行分组，如先天盲、后天盲和正常人，可以对两种方法进行多个对照组的比较，探究出盲人在先天和后天上脑皮层形态的变化.此外，关于脑部结构度量有很多形态学指标，后续我们可以计算更多数据指标，如皮层复杂度、灰质体积、灰质密度、白质体积等，通过分析这些指标可以更好地研究盲人的大脑结构，这对于医学研究将具有十分重要的意义.