专利名称：三维大脑磁共振图像的大脑皮层表面上脑沟区域分割方法
技术领域：
本发明涉及一种三维大脑磁共振图像的大脑皮层表面上脑沟区域分割方法，属于 医学图像处理，计算神经解剖学等领域。适用于人类三维大脑核磁共振图像重构出的 三角化的大脑皮层表面上脑沟区域的分割。
背景技术：
人类大脑皮层是一个极其复杂巻曲的解剖结构，主要由脑沟和脑回构成，分别对 应于大脑皮层上的谷和脊。尽管不同人之间的脑沟和脑回的精确几何模式变化很大， 但是最主要的几条脑沟和脑回是大脑皮层上的共有的解剖界标。因此，主要的脑沟和 脑回已经被广泛的用于辅助非线性大脑核磁共振图像配准，分析正常人大脑的解剖结 构变化规律，以及用来区分正常人和疾病患者。但是由于手工分割和标定脑沟及其费 时，并且容易受到外界主观影响。
近些年来脑沟区域的自动分割已经成为研究的热点课题。各种各样脑沟区域分割 方法已经被提出来，但仍存在很多问题。方法1:基于大脑皮层表面上脑沟深度的分 水岭方法，该方法首先利用主动表面方法找到脑回区域，然后利用分水岭方法提取脑 沟区域，最后利用启发式规则合并过分割的脑沟区域。其缺点在于，脑沟深度不能够 很好的区分脑沟区域和脑沟中深埋的脑回区域，因为脑沟区域和脑沟中深埋的脑回区 域的深度都很大，分水岭方法很容易产生过分割现象，将一个脑沟区域分割为多个脑 沟区域，而且启发式规则合并过分割脑沟区域较难控制；方法2:基于大脑皮层表面 上平均曲率的图剪切方法，该方法将大脑皮层表面看作一个连通图结构，根据平均曲 率信息，利用图剪切方法将大脑皮层表面分割为脑沟区域和脑回区域。特点在于图剪 切方法效率很高，但平均曲率不能很好的区分脑沟和脑回区域；方法3:基于平均曲
4率和脑沟深度的贝叶斯方法，该方法首先利用贝叶斯方法将大脑皮层表面分割为脑沟 和脑回区域，然后利用分水岭区域增长方法提取单个脑沟区域，最后利用启发式规则 合并过分割的脑沟区域。缺点在于该方法所定义的欧氏脑沟深度并非真正的脑沟深 度，同时平均曲率不能很好的区分脑沟和脑回区域。
目前己有的大脑皮层表面上脑沟区域分割方法具有以下三个主要的缺陷其一、 利用平均主曲率来区分脑沟和脑回区域，但是平均主曲率是最大主曲率和最小主曲率 的平均值，而脑沟和脑回区域的最小主曲率都非常小，因此平均曲率不能很好的区分 脑沟和脑回区域。其二、利用脑沟深度来区分脑沟和脑回区域，但是脑沟区域和脑沟 中深埋的脑回区域的深度都很大，因此脑沟深度也不能很好的区分脑沟和脑回区域。 其三、利用分水岭方法和启发式合并规则提取单个脑沟区域，但是分水岭方法很容易 将一个脑沟区域分割为多个脑沟区域，而且启发式合并规则较难控制。
发明内容
要解决的技术问题
为了避免现有方法的不足之处，本发明提出一种三维大脑磁共振图像的大脑皮层 表面上脑沟区域分割方法，可以获取几何结构精确和拓扑结构正确的大脑皮层表面图 像。
技术方案
本发明的基本思想是大脑皮层表面顶点的最大主曲率在脑沟区域和脑回区域分 别为负值和正值，我们利用双混合高斯模型来建模最大主曲率的直方图分布，为了同 时将统计信息和空间邻域信息合并到高斯混合模型中，利用隐马尔可夫随机场期望最 大化框架来进行脑沟区域的分割。
本发明的技术特征具体步骤如下步骤1对三维大脑磁共振图像进行预处理利用可变性模型方法去除脑壳,利用配 准方法去除非大脑组织，利用高斯混合模型方法对大脑图像进行组织分割，得到白质， 灰质和脑脊髓液三种组织类型表示的图像；
步骤2大脑皮层表面重建利用Marching Cubes方法从上述步骤得到的组织分割 后的大脑图像中重构三角化的大脑皮层表面；
步骤3利用有限差分方法估计得到大脑皮层表面上每个顶点的最大主曲率
步骤4统计得到最大主曲率分布直方图，通过隐马尔可夫随机场期望最大化模型 分割脑沟区域和脑回区域-
包括以下步骤-
a、 根据最大主曲率分布直方图，利用0tsu自适应阈值方法得到大脑皮层表面上 初始的脑沟区域和脑回区域分割xw，其中义-(x,,...,、)，"是大脑皮层表面上总顶 点个数，e{0, 1}类标记，0表示脑沟区域，l表示脑回区域，估计初始高斯混合模型 参数^。所述初始模型参数为均值/^和方差^w，其中/表示脑沟区域或脑回区域；
b、 利用迭代条件模ICM，求解MRF-MAP，
X"^argmax(log尸(71义,eW) + logP(X))，将每个顶点划分为脑沟或脑回，其中 y = Ov..，>0， 乂是顶点z'的最大主曲率；
C 、根据当前脑沟和脑回分割结果，采用A'+')和
(= W)",-A)2 ，更新高斯混合模型参数,其中s表示大脑皮层表面上
所有的顶点，后验分布为尸(')(/1 乂) = g(')a',)(/ K')，!)是x,=《和^ {A，g/}
力,) '
的马尔科夫局部相关概率，p(乂)是先验概率，g(y,;《)为高斯函数d、 重复步骤b c直到相邻迭代结果的变化小于总顶点数目的万分之一；
e、 将所有脑回标记为一种颜色，利用连通成分分析，将同一个连通的脑沟区域标 记为同一个颜色，得到大脑皮层表面上的脑回区域和一系列脑沟区域。
有益效果
本发明提出的三维大脑磁共振图像的大脑皮层表面上脑沟区域分割方法，首先， 随着磁共振成像设备的精度不断提高和三维大脑磁共振图像的预处理方法进一步成 熟，获取几何结构精确和拓扑结构正确的大脑皮层表面相对容易；同时，最大主曲率 在脑沟区域为负值和在脑回区域为正值，利用最大主曲率来分割脑沟区域和脑回区域 是可行的。
本发明相对于其它方法具有以下优点1、利用最大主曲率能很好的区分脑沟区域 和脑回区域，同时最大主曲率能够区别脑沟区域和深埋在脑沟中的脑回区域；2、利用 隐马尔可夫随机场期望最大化方法同时考虑了表面上顶点的最大主曲率的统计信息和 邻域信息。
图l:本发明方法的基本流程图。
图2: 1个大脑皮层左半球内表面上最大主曲率分布图。
图3: 12个真实正常人大脑皮层左半球上脑沟区域分割结果。
图4: 1个脑沟中存在深埋的脑回的分割结果。
具体实施例方式
现结合实施例、附图对本发明作进一步描述
根据本发明提出的基于最大主曲率和隐马尔可夫随机场期望最大化模型的大脑皮层表面上脑沟区域的分割方法，我们用0++语言实现了一个脑沟区域分割原型系统。
图像数据的来源是实际中正常人的三维大脑核磁共振图像。
首先对三维大脑核磁共振图像进行预处理和大脑皮层表面重建包括去除脑壳和 非大脑组织，对大脑图像进行脑组织分割(分割为白质，灰质和脑脊髓液三种类型)， 从组织分割后的大脑图像中重构几何结果精确，拓扑结构正确的大脑皮层表面，该皮 层表面由一系列顶点和三角形表示。然后，估计大脑皮层表面上顶点的最大主曲率。 最后，根据最大主曲率和隐马尔可夫随机场期望最大化框架将大脑皮层表面上分割为 脑沟区域和脑回区域，并利用连通性分析标记每一个脑沟区域。
所述的大脑皮层表面上顶点的最大主曲率的估计方法是首先计算每个三角形面 的Weingarten矩阵，然后每个顶点的Weingarten矩阵计算为该顶点周围一圈三角形面 Weingarten矩阵的加权平均，最后计算个顶点的Weingarten矩阵的特征值和特征向量， 其中绝对值最大的特征值即为最大主曲率。
基于隐马尔可夫随机场期望最大化框架的脑沟和脑回区域分割方法的原理是首 先，根据最大主曲率分布直方图，通过Otsu自适应阈值方法得到大脑皮层表面上脑沟 和脑回区域的初始分割结果和混合高斯模型的参数，包括均值和方差；然后，通过隐 马尔可夫随机场期望最大化模型估计每个顶点的类标记(脑沟或脑回)；接着，根据当 前的脑沟和脑回分割结果更新混合高斯模型的参数；通过迭代执行第二步和第三步得 到最终的脑沟和脑回分割结果。最后，利用连通性分析将每个连通脑沟区域标记为一 个值。
本发明整个流程可以参考附图1，具体的实施步骤如下 1.预处理和大脑皮层表面重建
对三维大脑核磁共振图像进行去除头骨和非大脑组织，大脑组织分割和重建几何
8结果精确，拓扑结构正确的大脑皮层表面。
2. 大脑皮层表面上最大主曲率估计
首先计算每个三角形面的Weingarten矩阵，然后加权平均每个顶点周围一圈三角 形面Weingarten矩阵来计算每个顶点的Weingarten矩阵，最后计算每个顶点的 Weingarten矩阵的特征值和特征向量，其中绝对值最大的特征值即为最大主曲率。附 图显示了 1个大脑皮层左半球内表面(大脑白质和灰质的交界面)上估计的最大主曲率。
3. 基于隐马尔可夫随机场期望最大化框架的脑沟和脑回区域分割 假设在顶点,'，x是最大主曲率值，x,.e(O, 1}类标记，其中0表示脑沟区域，l表
示脑回区域。脑沟区域和脑回区域分割问题可以公式化为最大化一个后验概率分布 尸(xir)，这里义(X,...,;O， y = Ov..，>o，"是重建的三角化的大脑皮层表面上总
顶点个数。该分割问题可以重新公式化为寻找真正的类标记i = 使其满足:
f = argmax(尸(y |义)尸(JT)〉
该问题具体求解如下
1) 根据最大主曲率分布直方图，通过Otsu自适应阈值化方法，得到大脑皮层表 面上脑沟和脑回区域的初始分割结果^fW，并且估计脑沟和脑回区域各自高斯模型的 参数e""，包括均值和方差；
2) 通过马尔可夫随机场模型估计每个顶点的类标记
X(D = argmax(logP(J |义,<9(')) + logP(I)}
3) 根据当前的类标记，估计脑沟和脑回区域各自高斯模型的参数"'+1): 这里S表示大脑皮层表面上所有的顶点，后验分布计算为g(乂;《)是一个高斯函数
，)
4)通过迭代执行第2步和第3步得到最终的脑沟和脑回区域分割结果。
最后，利用连通性分析将每个连通脑沟区域标记为一个值。
为了测试该脑沟区域分割方法的准确度，我们将该方法用于12个真实正常人的大 脑核磁共振图像所重建出的大脑皮层表面。图3显示12个真实正常人大脑皮层左半球 上脑沟区域分割结果，其中每个脑沟区域被标为一个颜色。为了定量的评价该方法， 我们采用过分割和欠分割两种度量方法。过分割是指和专家的目测相比将本来一个脑 沟区域分割为多个脑沟区域。欠分割是指没有将相邻的多个脑沟区域合适的分割开。 我们利用利用三个主要的脑沟，包括中央沟，中央后沟和颞上沟来验证。在以上12 个正常人大脑左半球上，所有的中央沟和颞上沟都被正确分割出来，因此没有过分割 和欠分割错误。中央后沟中存在一例过分割错误和没有欠分割错误。从实验结果看出 我们的脑沟区域分割方法具有很好的表现。该分割方法可以很好的处理脑回深埋在脑 沟中情况，附图4显示了一个例子 一个脑回深埋在脑沟中的例子的分割结果。(a) 和(b)是图(c)中红色矩形区域所框定区域的最大主曲率图和脑沟区域分割结果的放大 图。(c)是整个大脑皮层表面半球上脑沟区域分割结果，其中每个颜色表示一个连通的 脑沟区域。
权利要求
1.三维大脑磁共振图像的大脑皮层表面上脑沟区域分割方法，其特征在于步骤1对三维大脑磁共振图像进行预处理利用可变性模型方法去除脑壳，利用配准方法去除非大脑组织，利用高斯混合模型方法对大脑图像进行组织分割，得到白质，灰质和脑脊髓液三种组织类型表示的图像；步骤2大脑皮层表面重建利用Marching Cubes方法从上述步骤得到的组织分割后的大脑图像中重构三角化的大脑皮层表面；步骤3利用有限差分方法估计得到大脑皮层表面上每个顶点的最大主曲率步骤4统计得到最大主曲率分布直方图，通过隐马尔可夫随机场期望最大化模型分割脑沟区域和脑回区域包括以下步骤a、根据最大主曲率分布直方图，利用Otsu自适应阈值方法得到大脑皮层表面上初始的脑沟区域和脑回区域分割X(0)，其中X＝(x1，...，xn)，n是大脑皮层表面上总顶点个数，xi∈{0，1}类标记，0表示脑沟区域，1表示脑回区域，估计初始高斯混合模型参数θ(0)所述初始模型参数为均值μl0)和方差σl(0)，其中l表示脑沟区域或脑回区域；b、利用迭代条件模ICM，求解MRF-MAP，<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msup> <mi>X</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><mo>=</mo><munder> <mrow><mi>arg</mi><mi>max</mi> </mrow> <mi>x</mi></munder><mo>{</mo><mi>log</mi><mi>P</mi><mrow> <mo>(</mo> <mi>Y</mi> <mo>|</mo> <mi>X</mi> <mo>,</mo> <msup><mi>θ</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo></mrow> </msup> <mo>)</mo></mrow><mo>+</mo><mi>log</mi><mi>P</mi><mrow> <mo>(</mo> <mi>X</mi> <mo>)</mo></mrow><mo>}</mo><mo>,</mo> </mrow>]]></math> id="icf0001" file="A2009100217850002C1.tif" wi="74" he="7" top= "193" left = "38" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>将每个顶点划分为脑沟或脑回，其中Y＝(y1，...，yn)，yi是顶点i的最大主曲率；C、根据当前脑沟和脑回分割结果，采用<maths id="math0002" num="0002" ><math><![CDATA[ <mrow><msubsup> <mi>μ</mi> <mi>l</mi> <mrow><mo>(</mo><mi>t</mi><mo>+</mo><mn>1</mn><mo>)</mo> </mrow></msubsup><mo>=</mo><mfrac> <mrow><msub> <mi>Σ</mi> <mrow><mi>i</mi><mo>∈</mo><mi>S</mi> </mrow></msub><msup> <mi>P</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><mrow> <mo>(</mo> <mi>l</mi> <mo>|</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow><msub> <mi>y</mi> <mi>i</mi></msub> </mrow> <mrow><msub> <mi>Σ</mi> <mrow><mi>i</mi><mo>∈</mo><mi>S</mi> </mrow></msub><msup> <mi>P</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><mrow> <mo>(</mo> <mi>l</mi> <mo>|</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow> </mrow></mfrac> </mrow>]]></math> id="icf0002" file="A2009100217850002C2.tif" wi="35" he="10" top= "215" left = "143" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>和<maths id="math0003" num="0003" ><math><![CDATA[ <mrow><msup> <mrow><mo>(</mo><msubsup> <mi>σ</mi> <mi>l</mi> <mrow><mo>(</mo><mi>t</mi><mo>+</mo><mn>1</mn><mo>)</mo> </mrow></msubsup><mo>)</mo> </mrow> <mn>2</mn></msup><mo>=</mo><mfrac> <mrow><msub> <mi>Σ</mi> <mrow><mi>i</mi><mo>∈</mo><mi>S</mi> </mrow></msub><msup> <mi>P</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><msup> <mrow><mrow> <mo>(</mo> <mi>l</mi> <mo>|</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow><mrow> <mo>(</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>-</mo> <msub><mi>μ</mi><mi>l</mi> </msub> <mo>)</mo></mrow> </mrow> <mn>2</mn></msup> </mrow> <mrow><msub> <mi>Σ</mi> <mrow><mi>i</mi><mo>∈</mo><mi>S</mi> </mrow></msub><msup> <mi>P</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><mrow> <mo>(</mo> <mi>l</mi> <mo>|</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow> </mrow></mfrac><mo>,</mo> </mrow>]]></math> id="icf0003" file="A2009100217850002C3.tif" wi="52" he="10" top= "230" left = "30" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>更新高斯混合模型参数θ(t+1)，其中S表示大脑皮层表面上所有的顶点，后验分布为<maths id="math0004" num="0004" ><math><![CDATA[ <mrow><msup> <mi>P</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><mrow> <mo>(</mo> <mi>l</mi> <mo>|</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow><mo>=</mo><mfrac> <mrow><msup> <mi>g</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><mrow> <mo>(</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>;</mo> <msub><mi>θ</mi><mi>l</mi> </msub> <mo>)</mo></mrow><mo>·</mo><msup> <mi>P</mi> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow></msup><mrow> <mo>(</mo> <mi>l</mi> <mo>|</mo> <msub><mi>x</mi><msub> <mi>N</mi> <mi>l</mi></msub> </msub> <mo>)</mo></mrow> </mrow> <mrow><mi>p</mi><mrow> <mo>(</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow> </mrow></mfrac><mo>,</mo> </mrow>]]></math> id="icf0004" file="A2009100217850002C4.tif" wi="58" he="10" top= "246" left = "87" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths> id="icf0005" file="A2009100217850002C5.tif" wi="17" he="6" top= "249" left = "148" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>是xi＝l和θ＝{μl，σl}的马尔科夫局部相关概率，p(yi)是先验概率，g(yi；θl)为高斯函数<maths id="math0005" num="0005" ><math><![CDATA[ <mrow><mi>g</mi><mrow> <mo>(</mo> <mi>y</mi> <mo>;</mo> <msub><mi>θ</mi><mi>l</mi> </msub> <mo>)</mo></mrow><mo>=</mo><mfrac> <mn>1</mn> <msqrt><mn>2</mn><mi>π</mi><msubsup> <mi>σ</mi> <mi>l</mi> <mn>2</mn></msubsup> </msqrt></mfrac><mi>exp</mi><mrow> <mo>(</mo> <mo>-</mo> <mfrac><msup> <mrow><mo>(</mo><mi>y</mi><mo>-</mo><msub> <mi>μ</mi> <mi>l</mi></msub><mo>)</mo> </mrow> <mn>2</mn></msup><mrow> <mn>2</mn> <msubsup><mi>σ</mi><mi>l</mi><mn>2</mn> </msubsup></mrow> </mfrac> <mo>)</mo></mrow><mo>;</mo> </mrow>]]></math></maths>d、重复步骤b～c直到相邻迭代结果的变化小于总顶点数目的万分之一；e、将所有脑回标记为一种颜色，利用连通成分分析，将同一个连通的脑沟区域标记为同一个颜色，得到大脑皮层表面上的脑回区域和一系列脑沟区域。
全文摘要
本发明涉及一种三维大脑磁共振图像的大脑皮层表面上脑沟区域分割方法，技术特征在于首先对3维大脑核磁共振图像进行预处理和大脑皮层表面重建包括去除脑壳和非大脑组织，对大脑图像进行脑组织分割，从组织分割后的大脑图像中重构几何上精确的，拓扑结构正确的大脑皮层表面，该皮层表面是由一系列顶点和三角形表示。然后，估计大脑皮层表面的各顶点的最大主曲率和最小主曲率。最后，根据最大主曲率和利用隐马尔可夫随机场期望最大化框架将大脑皮层表面上分割为脑沟和脑回区域，利用连通性分析标记每一个脑沟区域。本发明相对其它方法具有算法简单有效，分割准确度高的优点。
文档编号G06T7/00GK101515367SQ20091002178
公开日2009年8月26日 申请日期2009年4月1日 优先权日2009年4月1日
发明者刘天明, 刚 李, 聂晶鑫, 雷 郭 申请人:西北工业大学