CN104102839A - 一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法 - Google Patents

Abstract

本发明公开了一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，步骤如下：确定阿尔茨海默病患者(AD)和正常人(NC)两组样本，并等比例分成样本集和测试集；依据样本的脑部磁共振(MRI)图像提取多尺度网格曲面；针对各顶点计算局部点面距离(LVPD)和平均曲率；以平滑过的LVPD和平均曲率为观测指标，提取具有显著统计学差异的区域，并筛选出两种指标意义下的种子点；对训练集每个样本提取特征行向量，构成特征矩阵，用降维后的特征矩阵结合相应样本类别来训练分类器；用测试集样本测试分类器性能。本发明克服了现有技术易受皮层分割误差影响和可能遗漏某个尺度差异的缺点，实现了依据脑皮层多尺度形状特征为依据的两组样本分类。

Description

一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法

技术领域

本发明属于医学图像处理技术领域，尤其是涉及一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，该方法可以利用脑皮层的形状特征检测阿尔茨海默病病例，具有临床辅助诊断的作用。

背景技术

随着人口老龄化问题的日益突出，关注老年人生存质量，关注包括阿尔茨海默病(Alzheimer’s Disease，AD)在内的老年疾病，探索认识、预防和治疗此类疾病的新方法和新手段对提高社会的整体生活质量有实际价值。

AD的成因复杂，发展过程是渐进的，并且临床表现各异，特别是早期并没有明显的临床症状，在成像技术成熟以前，AD只有等到患者去世后对其进行脑部解剖才能确诊。脑部磁共振成像(Magnetic Resonance Imaging，MRI)技术已经可以获得较好空间分辨率(典型值1×1×1毫米)的三维大脑内部图像，从而为非侵入式地研究AD患者的大脑形态与疾病之间的关系提供了数据基础。因而采用MRI进行脑部扫描，提取具有预测或者指示性的特征，进而自动地完成诊断分类，已成为重要的辅助诊断措施。

现有的基于MRI的AD病例分类方法主要包括以下两类：

1、基于全脑的分类方法：典型代表为基于体素(voxel-based)或皮层厚度(corticalthickness)的分类方法，然而在基于体素的框架下，选择更多的特征并不能提升分类性能但却明显地增加了计算时间，而且年龄较大的NC组病例(正常老化)和较年轻的AD病人很容易被混淆和误判。而皮层厚度方法不可避免地涉及到组织分割或者功能结构分割，由于存在较大的个体差异，该方法对分割算法的精度和操作者的经验和耐心都是极大的挑战，并且不具备多尺度特性，可能存在遗漏某个尺度上的差异；

2、基于感兴趣区域的分类方法：比如基于海马体和颞叶局部区域的分类方法，涉及到组织分割和功能结构分割，因此与皮层厚度方法类似，该方法同样容易受到组织分割精度和绝对观测值误差的影响。

发明内容

本发明的目的是针对上述已有技术的不足，提出一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，以克服常用AD脑皮层分类方法的易受皮质分割误差影响和可能存在某个尺度上的差异遗漏的缺点，实现多尺度的阿尔茨海默病样本分类。

为解决上述技术问题，本发明的技术方案是：

一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，从形态学的角度，在多个尺度上分析正常人(NC)和阿尔茨海默病(AD)的脑皮层形状特征具有显著统计学差异的区域，采用基于种子点选择和区域生长准则的方法对每个样本进行稳定特征提取，使用训练集数据训练分类器，从而实现两组样本的分类，包括如下步骤：

(1)按临床医学诊断标准，分别采集正常人对照组(NC)和阿尔茨海默病组(AD)的大脑磁共振图像，两组样本的数目相等或接近，并从两组样本中等比例地选择一部分样本用作训练集，剩余的作为测试集；

(2)针对所有样本的大脑核磁共振图像，经过图像强度矫正、配准至标准脑、灰质/白质/脑脊液分割、脑皮层表面提取与配准、顶点采样等过程得到各样本脑皮层的各尺度三角网格曲面G_j(1≤j≤L)，j称为网格曲面G_j的尺度，且G_L为最精细网格曲面，G₀为最粗糙网格曲面；

(3)记第j级三角网格曲面G_j中任意顶点为p，定义顶点p的1-环邻域内的顶点构成的集合为E_p，并用集合E_p内的所有顶点构造一个平面f_p，记顶点p到该平面f_p的距离为顶点p的局部点面距离d_p；

(4)计算各尺度网格上顶点p的局部点面距离d_p和该顶点的平均曲率c_p，并用半高全宽半径为10毫米的扩散核函数对d_p和c_p分别进行平滑，得到平滑后的局部点面距离d′_p和平均曲率c′_p；

(5)结合样本人口学信息，建立关于观测值与性别、年龄和分组类别的多元线性回归模型，以各尺度网格上顶点的局部点面距离d′_p为观测值，通过双样本T检验对AD和NC两组的训练集进行组间比较，得到局部点面距离意义下各尺度网格上顶点的统计学差异显著性水平值；以各尺度网格上顶点的平均曲率c′_p为观测值，通过双样本T检验对AD和NC两组的训练集进行组间比较，得到平均曲率意义下各尺度网格上顶点的统计学差异显著性水平值；

(6)针对训练集样本，以各尺度网格上顶点的统计学差异显著性水平值为标准，初步筛选出局部点面距离意义下的种子点p_di，i为筛选出的第i个种子点，1≤i≤N₀，N₀为符合初步筛选条件的所有尺度上种子点的总个数，d表示以局部点面距离为指标；初步筛选出平均曲率意义下的种子点p_cq，q为筛选出的第q个种子点，1≤q≤M₀，M₀为符合初步筛选条件的所有尺度上种子点的总个数，c表示以平均曲率为指标；

(7)针对训练集样本，对每个初步筛选出的种子点p_di进行区域生长，得到每个种子点对应的差异区域顶点集合J_di，其中1≤i≤N，N是经过二次筛选后的局部点面距离意义下的所有尺度上种子点的总个数；对每个初步筛选出的种子点p_cq进行区域生长得到每个种子点对应的差异区域顶点集合J_cq，其中1≤q≤M，M是经过二次筛选后的平均曲率意义下的所有尺度上种子点的总个数；

(8)针对训练集每个样本，计算各个差异区域顶点集合J_di(1≤i≤N)内所有顶点的局部点面距离均值u_di(1≤i≤N)，以种子点序号从小到大的顺序排列u_di，得到表征每个样本各尺度上局部点面距离大小的行向量U_d＝[u_di，i＝1...N]；

(9)针对训练集每个样本，计算每一个差异区域顶点集合J_cq(1≤q≤M)内所有顶点的平均曲率均值u_cq(1≤q≤M)，以种子点序号从小到大的顺序排列u_cq，得到表征每个样本各尺度上平均曲率大小的行向量U_c＝[u_cq，q＝1...M]；

(10)针对训练集每个样本，将行向量U_d＝[u_di，i＝1...N]和行向量U_c＝[u_cq，q＝1...M]合并成为得到表征每个样本所有尺度上形状特征大小的行向量U＝[u_di，u_cq]，i＝1...N，q＝1...M，将样本集内所有样本的行向量U构成每行元素个数为N+M的特征矩阵T；

(11)对特征矩阵T进行归一化，并使用主分量分析进行降维处理，选取每一行的前K个元素，K＜N+M，构成训练集的降维特征矩阵

(12)对测试集的所有样本采用与步骤(6)到步骤(11)相同的方法，计算测试集的降维特征矩阵

(13)将训练集的降维特征矩阵和各样本的分组类别输入支撑向量机完成分类器的训练，输出经过训练得到的分类器的各项系数；

(14)为了测试分类器的性能，将测试集的降维特征矩阵输入相应的分类器，并将分类器输出结果与样本相应的分类组别相比较，输出分类的准确率、敏感性和特异性的数值。

进一步地，所述步骤(3)中用集合E_p内的所有顶点构造得到一个平面f_p，按如下步骤进行：

(3a)设平面f_p的方程为z＝a₀x+a₁y+a₂，其中a₀，a₁，a₂为平面待定系数，a₀为x的待定系数，a₁为y的待定系数，a₂为常数项的待定系数，x，y为自变量，z为因变量；

(3b)为确定平面待定系数a₀，a₁，a₂，构造集合E_p内的所有顶点到平面f_p的偏移平方和函数S：

$S=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(a_0{x}_i+a_1{y}_i+a_2-z_i)}^2---1)$

其中(x_i，y_i，z_i)为集合E_p中第i个顶点的坐标，i＝1，2...n，n为集合E_p中顶点的个数；

(3c)根据极小值的判定方法，利用函数S分别对平面待定系数a₀，a₁，a₂求偏导，得到如下方程组：

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)x_i=0\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)y_i=0\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)=0\end{array}\right.$

其中，(x_i，y_i，z_i)为集合E_p中第i个顶点的坐标，i＝1，2...n；

(3d)解上面的方程组得到平面待定系数a₀，a₁，a₂为：

$\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)={\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x_i}^2& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y_i}^2& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i& n\end{array}\right]}^{-1}\left(\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{z}_i\end{array}\right)---2)$

其中[·]^-1表示矩阵[·]的逆阵；

(3e)将平面待定系数a₀，a₁，a₂代入平面f_p的方程z中，即得到了所要构造的平面f_p，此时集合E_p内的所有顶点到平面f_p的偏移平方和函数S为极小值。

进一步地，步骤(5)所述建立多元线性回归模型如下：

Y～b₀+b₁A+b₂S+b₃P    3)

其中，Y为观测值，在本发明中具体指各顶点处的局部点面距离或平均曲率，A是每个样本对应的年龄，S是每个样本对应的性别，P是每个样本对应组别(即AD或NC)，b₀，b₁，b₂，b₃是待定系数。

进一步地，步骤(6)所述的初步筛选，通过如下方法进行：

(6a)把各尺度上所选观测值意义下存在显著差异的簇包含的所有顶点，将其对应到解剖结构模板上，得到各尺度上具有显著差异的解剖学结构；

(6b)针对各尺度上所选观测值意义下每一个具有显著差异解剖学结构内所有顶点，检测该解剖学结构的显著性水平值的极小值，每个极小值顶点被筛选作该尺度上所选观测值意义下的一个种子点；

进一步地，步骤(7)所述的区域生长，通过如下方法进行：

(7a)设定阈值ξ，0＜ξ＜0.25，将种子点所在尺度的脑皮层网格曲面所有顶点的显著性水平值从小到大排列，选择前N_j×ξ个顶点的平均显著性水平值作为该尺度上种子点区域生长的阈值δ，其中j为种子点所在尺度，N_j为j尺度上网格曲面G_j的顶点数，1≤j≤L；

(7b)从每个种子点开始，在该种子点所在尺度的网格曲面上逐渐向外围扩展，并依次累加扩展过程中碰到的顶点的显著性水平值，直到累加值达到阈值δ或扩散到其它解剖结构时停止，此时扩展过程涉及到的顶点的集合即为该种子点对应的差异区域顶点集合。

进一步地，步骤(7)所述的二次筛选，通过如下方法进行：若种子点经过步骤(7)所述的区域生长以后，得到对应的差异区域顶点集合中顶点的个数小于该种子点所在尺度上相应解剖结构内的总顶点数的5％，则舍弃该种子点，同时也舍弃其生长的区域顶点集合。

采用上述技术方案，本发明的有益效果为：

(1)由于皮质分割误差并不会改变大脑皮层形状整体走势和变化规律，本发明可以避免全脑法和感兴趣区域法中分割精度造成的不利影响。

(2)本发明通过比较AD和NC样本的脑皮层曲面的形状特征的组间差异，能够有效地检测出阿尔茨海默病人与正常人之间的脑皮层形状差异。

(3)本发明采用选择种子点并对种子点进行区域生长的方法进行特征提取，避免了特征区域过大或过小的问题并且能够确保每个特征区域内的顶点集中在一个的解剖学结构区域，提高了特征的可靠性和抵抗随机噪声的能力。

(4)本发明通过在多个尺度上提取AD和NC两组样本的脑皮层形状特征差异信息，获得了可靠完备脑皮层形状差异信息，避免了遗漏某个尺度上的差异。

附图说明

图1是本发明的三角网格剖分示意图。

图2是本发明计算局部点面距离d_p的示意图。

图3是本发明计算平均曲率c_p的示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

实施例

本发明实施例中的基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，实现步骤如下：

第一步，将所有样本按临床医学诊断标准分成两组，分别称为对照组NC和阿尔茨海默病AD。

表1给出了一个建议的诊断标准，但在实际操作之中仍然可以做出调整或者更换，只需要满足特定的临床特征即可。在本步骤中，样本的大脑核磁共振图像(MRI)数据和人口学数据完备可靠，NC和AD两组样本的数目相等或接近，这是为了确保采用统计学方法分析得到的结果的准确性和可靠性。

表1 样本的医学诊断标准

注1：MMSE的全称是简易智力状态检查(Mini-Mental State Examination，MMSE)，它采用问卷量表的形式对受试者进行智力和精神状态的检查，满分是30分。临床认为，MMSE＞＝25正常；21＜＝MMSE＜＝24轻度痴呆；10＜＝MMSE＜＝20中度痴呆；MMSE＜＝9重度痴呆。

注2：CDR的全称是临床痴呆分级(Clinical Dementia Rating，CDR)，通过医学评判和诊断将痴呆分为五个等级，即CDR＝0，正常；CDR＝0.5，很轻微痴呆；CDR＝1，轻度痴呆；CDR＝2，中度痴呆；CDR＝3，重度痴呆。

第二步，针对所有样本的大脑核磁共振图像，经过图像预处理后，采用三角形网格逐级剖分法，得到各样本合格的脑皮层各尺度网格曲面G_j(1≤j≤L)，j称为网格曲面G_j的尺度，且G_L为最精细网格曲面，G₀为最粗糙网格曲面。

(2a)样本的MRI图像由于存在各种环境和扫描因素的差异，为了保证皮质表面处理结果的可靠性和准确性，对MRI图像进行强度矫正、配准至标准脑(比如ICBM152模板)、灰质(GM)/白质(WM)/脑脊液(CSF)分割等预处理过程，得到脑皮层图像；

(2b)对配准的脑皮层图像采用三角形网格逐级剖分法提取网格曲面，参照附图1，具体实现如下：本发明选取的初始网格为十二个顶点的二十面体G₀，则G₀为最粗糙网格，各尺度网格曲面G_j(1≤j≤L)是由最粗糙网格G₀经过j级剖分后得到，其剖分过程是在第j-1级网格中的某两个相邻顶点之间新增加一个j级网格顶点，新增网格顶点P_j要尽可能地使剖分后的网格能够更精确地逼近脑皮层外形轮廓，将新增网格顶点P_j与已有顶点合并，得到第j级三角形网格G_j，则有：

G_j＝G_j-1∪P_j，1≤j≤L    1)

其中∪表示集合的并运算，j称为网格G_j的尺度，很显然在经过L次剖分以后得到的G_L为最精细网格。在这个步骤中，除了顶点的合并以外，还需要按照空间位置关系重新构造出网格上的三角形，第j尺度下网格曲面的总顶点数为10×4^j+2，三角形个数为20×4^j，0≤j≤L。

第三步，记第j级三角网格G_j中任意顶点为p，定义顶点p的1-环邻域内的顶点构成的集合为E_p。

第四步，用集合E_p内的所有顶点构造得到一个平面f_p。

(4a)设平面f_p的方程为z＝a₀x+a₁y+a₂，其中a₀，a₁，a₂为平面待定系数，a₀为x的待定系数，a₁为y的待定系数，a₂为常数项的待定系数，x，y为自变量，z为因变量；

(4b)为确定平面待定系数a₀，a₁，a₂，构造集合E_p内的所有顶点到平面f_p的偏移平方和函数S：

$S=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(a_0{x}_i+a_1{y}_i+a_2-z_i)}^2---2)$

其中(x_i，y_i，z_i)为集合E_p中第i个顶点的坐标，i＝1，2...n，n为集合E_p中顶点的个数；

(4c)根据极小值的判定方法，利用函数S分别对平面待定系数a₀，a₁，a₂求偏导，得到如下方程组：

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)x_i=0\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)y_i=0\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)=0\end{array}\right.$

其中，(x_i，y_i，z_i)为集合E_p中第i个顶点的坐标，i＝1，2...n；

(4d)解上面的方程组得到平面待定系数a₀，a₁，a₂为：

$\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)={\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x_i}^2& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y_i}^2& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i& n\end{array}\right]}^{-1}\left(\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{z}_i\end{array}\right)---3)$

其中[·]^-1表示矩阵[·]的逆阵；

(4e)将平面待定系数a₀，a₁，a₂代入平面f_p的方程z中，即得到了所要构造的平面f_p，此时集合E_p内的所有顶点到平面f_p的偏移平方和函数S为极小值。

第五步，计算顶点p到平面f_p的距离d_p，并用半高全宽半径为10毫米的扩散核函数对d_p进行平滑，得到平滑后的局部点面距离d′_p。

由第四步得到了平面f_p的方程z＝a₀x+a₁y+a₂，参照附图2，任取平面内一点q，设其坐标为(x′，y′，z′)，则有z′＝a₀x′+a₁y′+a₂；设顶点p的坐标为(x_p，y_p，z_p)，连接顶点p和q点，并过q点做平面f_p的单位法向量所求顶点p到平面f_p的距离d_p为向量与的内积数值大小，即：

$\begin{array}{c}{d}_p=|<\stackrel{\&RightArrow;}{\mathrm{pq}},\stackrel{\&RightArrow;}{n}>|=\frac{|a_0(x_p-x^{\&prime;})+a_1(y_p-y^{\&prime;})-(z_p-z^{\&prime;})|}{\sqrt{{a_0}^2+{a_1}^2+1}}\\ =\frac{|a_0{x}_p+a_1{y}_p-z_p-(a_0{x}^{\&prime;}+a_1{y}^{\&prime;}-z^{\&prime;})|}{\sqrt{{a_0}^2+{a_1}^2+1}}\end{array}---4)$

第六步，计算顶点p处的平均曲率c_p，并用半高全宽半径为10毫米的扩散核函数对c_p进行平滑，得到平滑后的平均曲率c′_p。

参照附图3，α_i和a_i分别是点p的1-环邻域内每个三角形的夹角和面积，β_j为相邻三角形所在平面的法向量和之间的夹角，e_i为1-环邻域各边的长度，e_j为点p与1-环邻域各点连线的长度，则点p的平均曲率为：

$c_p=\frac{\frac{1}{4}{\mathrm{\&Sigma;}}_j{\mathrm{\&beta;}}_j\&CenterDot;e_j}{\frac{1}{2}{\mathrm{\&Sigma;}}_i{a}_i-\frac{1}{8}{\mathrm{\&Sigma;}}_i{e}_i^2\cot \left({\mathrm{\&alpha;}}_i\right)}---5)$

第七步，从两组样本集中等比例(比如60％)的选择AD样本和NC样本作为训练集，剩余的作为测试集。

第八步，结合样本人口学信息，建立关于观测值与性别、年龄和分组类别的多元线性回归模型，分别以局部点面距离和平均曲率为观测指标，通过统计分析获取各尺度上脑皮层网格上具有显著统计学差异的区域。

(8a)结合样本人口学信息，建立关于观测值与性别、年龄和分组类别的多元线性回归模型：

Y～b₀+b₁A+b₂S+b₃P，    6)

其中，Y为观测值，在本发明中具体指各顶点的局部点面距离或平均曲率，A是每个样本对应的年龄，S是每个样本对应的性别，P是每个样本对应组别(即AD或NC)，b₀，b₁，b₂，b₃是待定系数；

(8b)以脑皮层各尺度网格上各顶点的局部点面距离为观测值，通过双样本T检验对AD和NC两组的训练集进行组间比较，得到局部点面距离意义下各尺度上各顶点的统计学差异显著性水平值，从而获取局部点面距离意义下各尺度上具有显著统计学差异的簇；

(8c)以脑皮层各尺度网格上各顶点的平均曲率为观测值，通过双样本T检验对AD和NC两组的训练集进行组间比较，得到平均曲率意义下各尺度上各顶点的统计学差异显著性水平值，从而获取平均曲率意义下各尺度上具有显著统计学差异的簇。

第九步，针对训练集每个样本，初步筛选出局部点面距离意义下的种子点p_di(i为筛选出的第i个种子点，1≤i≤N₀，N₀为符合初步筛选条件的所有尺度上种子点的总个数，d表示以局部点面距离为指标)和平均曲率意义下的种子点p_cq(1≤q≤M₀，M₀为符合初步筛选条件的所有尺度上种子点的总个数，c表示以平均曲率为指标)。

(9a)针对局部点面距离意义下各尺度上存在显著差异的簇包含的所有顶点，将其对应到解剖结构模板上，得到局部点面距离意义下各尺度上具有显著差异的解剖学结构；

(9b)针对局部点面距离意义下各尺度上每一个存在显著差异解剖学结构内所有顶点，检测该解剖学结构局部点面距离意义下显著性水平值的极小值，每个极小值顶点被筛选作该尺度上局部点面距离意义下一个种子点p_di；

(9c)针对平均曲率意义下各尺度上存在显著差异的簇包含的所有顶点，将其对应到解剖结构模板上，得到平均曲率意义下各尺度上具有显著差异的解剖学结构；

(9d)针对平均曲率意义下各尺度上每一个存在显著差异解剖学结构内所有顶点，检测该解剖学结构平均曲率意义下显著性水平值的极小值，每个极小值顶点被筛选作该尺度上平均曲率意义下的一个种子点p_cq。

第十步，针对训练集每个样本，按照一定的区域生长准则对各尺度上每个初步筛选出的种子点进行区域生长，得到各尺度上每个种子点对应的差异区域顶点集合。

(10a)设定阈值ξ，0＜ξ＜0.25，将局部点面距离意义下的种子点所在尺度的脑皮层网格曲面所有顶点的局部点面距离意义下的显著性水平值从小到大排列，选择前N_j×ξ个顶点的平均显著性水平值作为该尺度上局部点面距离意义下的种子点区域生长的阈值δ，其中j为种子点所在尺度，N_j为j尺度上网格曲面G_j的顶点数，1≤j≤L；

(10b)从每个局部点面距离意义下的种子点开始，在该种子点所在尺度的网格曲面上逐渐向外围扩展，并依次累加扩展过程中碰到的顶点的局部点面距离意义下的显著性水平值，直到累加值达到阈值δ或扩散到其它解剖结构时停止，此时扩展过程涉及到的顶点的集合即为局部点面距离意义下该种子点对应的差异区域顶点集合；

(10c)针对生长后的局部点面距离意义下的种子点进行二次筛选，若种子点对应的差异区域顶点集合中顶点的个数小于该种子点所在尺度上相应解剖结构内的总顶点数的5％，则舍弃该种子点，同时也舍弃其生长的区域顶点集合，从而得到局部点面距离意义下种子点二次筛选后的差异区域顶点集合J_di，其中1≤i≤N，N是经过二次筛选后的局部点面距离意义下所有尺度上种子点的总个数；

(10d)设定阈值ξ，0＜ξ＜0.25，将平均曲率意义下的种子点所在尺度的脑皮层网格曲面所有顶点的平均曲率意义下的显著性水平值从小到大排列，选择前N_j×ξ个顶点的平均显著性水平值作为该尺度上平均曲率意义下的种子点区域生长的阈值δ，其中j为种子点所在尺度，N_j为j尺度上网格曲面G_j的顶点数，1≤j≤L；

(10e)从每个平均曲率意义下的种子点开始，在该种子点所在尺度的网格曲面上逐渐向外围扩展，并依次累加扩展过程中碰到的顶点的平均曲率意义下的显著性水平值，直到累加值达到阈值δ或扩散到其它解剖结构时停止，此时扩展过程涉及到的顶点的集合即为平均曲率意义下该种子点对应的差异区域顶点集合；

(10f)针对生长后的平均曲率意义下的种子点进行二次筛选，若种子点对应的差异区域顶点集合中顶点的个数小于该种子点所在尺度上相应解剖结构内的总顶点数的5％，则舍弃该种子点，同时也舍弃其生长的区域顶点集合，从而得到平均曲率意义下种子点二次筛选后的差异区域顶点集合J_cq，其中1≤q≤M，是经过二次筛选后的平均曲率意义下所有尺度上种子点的总个数。

第十一步，针对训练集每个样本，在各尺度上对每个种子点对应的差异区域顶点集合进行特征提取得到表征脑皮层所有尺度上形状特征的行向量U。

(11a)针对每个样本，计算所有尺度上每个差异区域顶点集合J_di(1≤i≤N)内所有顶点的局部点面距离均值u_di(1≤i≤N)，以种子点序号从小到大的顺序排列u_di，得到表征每个样本所有尺度上局部点面距离大小的行向量U_d＝[u_di，i＝1...N]；

(11b)针对每个样本，计算各尺度上每个差异区域顶点集合J_cq(1≤q≤M)内所有顶点的平均曲率均值u_cq(1≤q≤M)，以种子点序号从小到大的顺序排列u_cq，得到表征每个样本所有尺度上平均曲率大小的行向量U_c＝[u_cq，q＝1...M]；

(11c)将行向量U_d＝[u_di，i＝1...N]和行向量U_c＝[u_cq，q＝1...M]合并得到表征每个样本形状特征大小的行向量U＝[u_di，u_cq]，i＝1...N，q＝1...M。

第十二步，训练样本集内所有样本的行向量U构成每行元素个数为N+M的训练集特征矩阵T。

第十三步，针对训练集样本的特征矩阵T，进行归一化，并进行主分量分析降维处理，选取每一行的前K个元素，K＜N+M，构成训练集的降维特征矩阵

第十四步，计算测试集的降维特征矩阵

在本步骤中，对测试集样本采用与第九步到第十三步针对训练集处理的相同方法，得到测试集的降维特征矩阵

第十五步，训练集特征矩阵和各样本相应的分组类别训练分类器。

在本步骤中，将训练集的降维特征矩阵和各样本相应的分组类别输入支撑向量机完成分类器的训练，输出经过训练得到的分类器的各项系数。

第十六步，测试集特征矩阵和各样本相应的分组类别测试分类器的性能。

(16a)将测试集的降维特征矩阵输入分类器，输出每个测试集样本的分类结果；

(16b)将测试集样本的分类结果与样本相应的分组类别相比较，输出分类的准确率、敏感性和特异性的数值。

本发明不局限于上述具体的实施方式，本领域的普通技术人员从上述构思出发，不经过创造性的劳动，所做出的种种变换，均落在本发明的保护范围之内。

Claims (6)

1.一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，其特征在于，包括如下步骤：

(1)按临床医学诊断标准，分别采集正常人对照组(NC)和阿尔茨海默病组(AD)的大脑磁共振图像，两组样本的数目相等或接近，并从两组样本中等比例地选择一部分样本用作训练集，剩余的作为测试集；

(2)针对所有样本的大脑核磁共振图像，经过图像强度矫正、配准至标准脑、灰质/白质/脑脊液分割、脑皮层表面提取与配准、顶点采样等过程得到各样本脑皮层的各尺度三角网格曲面G_j(1≤j≤L)，j称为网格曲面G_j的尺度，且G_L为最精细网格曲面，G₀为最粗糙网格曲面；

(3)记第j级三角网格曲面G_j中任意顶点为p，定义顶点p的1-环邻域内的顶点构成的集合为E_p，并用集合E_p内的所有顶点构造一个平面f_p，记顶点p到该平面f_p的距离为顶点p的局部点面距离d_p；

(4)计算各尺度网格上顶点p的局部点面距离d_p和该顶点的平均曲率c_p，并用半高全宽半径为10毫米的扩散核函数对d_p和c_p分别进行平滑，得到平滑后的局部点面距离d_p和平均曲率c′_p；

(5)结合样本人口学信息，建立关于观测值与性别、年龄和分组类别的多元线性回归模型，以各尺度网格上顶点的局部点面距离d′_p为观测值，通过双样本T检验对AD和NC两组的训练集进行组间比较，得到局部点面距离意义下各尺度网格上顶点的统计学差异显著性水平值；以各尺度网格上顶点的平均曲率c′_p为观测值，通过双样本T检验对AD和NC两组的训练集进行组间比较，得到平均曲率意义下各尺度网格上顶点的统计学差异显著性水平值；

(6)针对训练集样本，以各尺度网格上顶点的统计学差异显著性水平值为标准，初步筛选出局部点面距离意义下的种子点p_di，i为筛选出的第i个种子点，1≤i≤N₀，N₀为符合初步筛选条件的所有尺度上种子点的总个数，d表示以局部点面距离为指标；初步筛选出平均曲率意义下的种子点p_cq，q为筛选出的第q个种子点，1≤q≤M₀，M₀为符合初步筛选条件的所有尺度上种子点的总个数，c表示以平均曲率为指标；

(7)针对训练集样本，对每个初步筛选出的种子点p_di进行区域生长，得到每个种子点对应的差异区域顶点集合J_di，其中1≤i≤N，N是经过二次筛选后的局部点面距离意义下的所有尺度上种子点的总个数；对每个初步筛选出的种子点p_cq进行区域生长得到每个种子点对应的差异区域顶点集合J_cq，其中1≤q≤M，M是经过二次筛选后的平均曲率意义下的所有尺度上种子点的总个数；

(8)针对训练集每个样本，计算各个差异区域顶点集合J_di(1≤i≤N)内所有顶点的局部点面距离均值u_di(1≤i≤N)，以种子点序号从小到大的顺序排列u_di，得到表征每个样本所有尺度上局部点面距离大小的行向量U_d＝[u_di，i＝1...N]；

(9)针对训练集每个样本，计算每一个差异区域顶点集合J_cq(1≤q≤M)内所有顶点的平均曲率均值u_cq(1≤q≤M)，以种子点序号从小到大的顺序排列μ_cq，得到表征每个样本所有尺度上平均曲率大小的行向量U_c＝[u_cq，q＝1...M]；

(10)针对训练集每个样本，将行向量U_d＝[u_di，i＝1...N]和行向量U_c＝[u_cq，q＝1...M]合并成为得到表征每个样本所有尺度上形状特征大小的行向量U＝[u_di，u_cq]，i＝1...N，q＝1...M，将训练集所有样本行向量U构成每行元素个数为N+M的训练集特征矩阵T；

(11)对训练集特征矩阵T进行归一化，并使用主分量分析进行降维处理，选取每一行的前K个元素，K＜N+M，构成训练集的降维特征矩阵

(12)对测试集所有样本采用与步骤(6)到步骤(11)相同的方法，计算测试集的降维特征矩阵

(13)将训练集的降维特征矩阵和各样本的分组类别输入支撑向量机完成分类器的训练，输出经过训练得到的分类器的各项系数；

(14)为了测试分类器的性能，将测试集的降维特征矩阵输入相应的分类器，并将分类器输出结果与样本相应的分类组别相比较，输出分类的准确率、敏感性和特异性的数值。

2.根据权利要求1所述的一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，其特征在于，其中所述步骤(3)中用集合E_p内的所有顶点构造得到一个平面f_p，按如下步骤进行：

(3a)设平面f_p的方程为z＝a₀x+a₁y+a₂，其中a₀，a₁，a₂为平面待定系数，a₀为x的待定系数，a₁为y的待定系数，a₂为常数项的待定系数，x，y为自变量，z为因变量；

(3b)为确定平面待定系数a₀，a₁，a₂，构造集合E_p内的所有顶点到平面f_p的偏移平方和函数S：

$S=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(a_0{x}_i+a_1{y}_i+a_2-z_i)}^2---1)$ 其中(x_i，y_i，z_i)为集合E_p中第i个顶点的坐标，i＝1，2...n，n为集合E_p中顶点的个数；

(3c)根据极小值的判定方法，利用函数S分别对平面待定系数a₀，a₀，a₂求偏导，得到如下方程组：

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)x_i=0\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)y_i=0\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}2(a_0{x}_i+a_1{y}_i+a_2-z_i)=0\end{array}\right.$

其中，(x_i，y_i，z_i)为集合E_p中第i个顶点的坐标，i＝1，2...n；

(3d)解上面的方程组得到平面待定系数a₀，a₁，a₂为：

$\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)={\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x_i}^2& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y_i}^2& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i& n\end{array}\right]}^{-1}\left(\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{z}_i\end{array}\right)---2)$ 其中[·]^-1表示矩阵[·]的逆阵；

(3e)将平面待定系数a₀，a₁，a₂代入平面f_p的方程z中，即得到了所要构造的平面f_p，此时集合E_p内的所有顶点到平面f_p的偏移平方和函数S为极小值。

3.根据权利要求1所述的一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，其特征在于，其中步骤(5)所述建立多元线性回归模型如下：

Y～b₀+b₁A+b₂S+b₃P    3)

其中，Y为观测值，在本发明中具体指各顶点处的局部点面距离或平均曲率，A是每个样本对应的年龄，S是每个样本对应的性别，P是每个样本对应组别(即AD或NC)，b₀，b₁，b_2，b₃是待定系数。

4.根据权利要求1所述的一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，其特征在于，其中步骤(6)所述的初步筛选，通过如下方法进行：

(6a)把各尺度上所选观测值意义下存在显著差异的簇包含的所有顶点，将其对应到解剖结构模板上，得到各尺度上具有显著差异的解剖学结构；

(6b)针对各尺度上所选观测值意义下每一个具有显著差异解剖学结构内所有顶点，检测该解剖学结构的显著性水平值的极小值，每个极小值顶点被筛选作该尺度上所选观测值意义下的一个种子点。

5.根据权利要求1所述的一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，其特征在于，其中步骤(7)所述的区域生长，通过如下方法进行：

(7a)设定阈值ξ，0＜ξ＜0.25，将种子点所在尺度的脑皮层网格曲面所有顶点的显著性水平值从小到大排列，选择前N_j×ξ个顶点的平均显著性水平值作为该尺度上种子点区域生长的阈值δ，其中j为种子点所在尺度，N_j为j尺度上网格曲面G_j的顶点数，1≤j≤L；

(7b)从每个种子点开始，在该种子点所在尺度的网格曲面上逐渐向外围扩展，并依次累加扩展过程中碰到的顶点的显著性水平值，直到累加值达到阈值δ或扩散到其它解剖结构时停止，此时扩展过程涉及到的顶点的集合即为该种子点对应的差异区域顶点集合。

6.根据权利要求1所述的一种基于多尺度网格曲面形状特征的阿尔茨海默病脑皮层自动分类方法，其特征在于，其中步骤(7)所述的二次筛选，通过如下方法进行：

若种子点经过步骤(7)所述的区域生长以后，得到对应的差异区域顶点集合中顶点的个数小于该种子点所在尺度上相应解剖结构内的总顶点数的5％，则舍弃该种子点，同时也舍弃其生长的区域顶点集合。