三维大脑磁共振图像大脑皮层表面最大主方向场弥散方法
技术领域
本发明涉及一种三维大脑磁共振图像大脑皮层表面最大主方向场弥散方法，利用alpha-expansion图割(graph cuts)方法的三维大脑磁共振图像的大脑皮层表面上最大主方向场弥散方法，属于医学图像处理领域。
背景技术
从三维大脑磁共振图像重构出来的大脑皮层表面是嵌入在三维空间的、极其卷曲和褶皱的二维流型。基于大脑皮层表面的分析非常适合于刻画和分析大脑皮层的本质结构，例如曲率、测地距离、皮层厚度等。大脑皮层表面上的方向场包含重要的信息，例如，最大主方向场和最小主方向场分别指示主曲率达到最大和最小的方向。但是，可靠计算大脑皮层表面上的主方向场并不容易，例如，在平坦的皮层表面区域所计算的主方向场通常噪声较大和不可靠，因为在平坦区域曲率在各个方向上的变化不大，微小的噪声就可能严重影响所计算的主方向场。因此，为了利用主方向场对大脑皮层表面进行分析，必须实施主方向场弥散，其目的是平滑主方向场中的噪声，同时保持主方向场中主要的几何结构和不连续性。弥散后的主方向场可以用于基于流场跟踪的脑沟或脑回盆地分割和方向导向的皮层褶皱形态分析等。
目前已有的大脑皮层表面上的最大主方向场弥散方法将该问题看做一个能量最小化问题，该能量函数包含一个数据项和一个平滑项。数据项用于使弥散后的最大主方向场与原始计算的最大主方向场尽可能相似，平滑项用于使弥散后的最大主方向场平滑地变化。该能量函数通过变分方法或梯度下降方式求解。已有的大脑皮层表面上的最大主方向场弥散方法具有以下的主要缺陷：其一、不同皮层表面区域的权重设为相同，导致弥散的最大主方向场在脑沟底部和脑回顶部过于平滑。其二、目前的能量最小化方法容易陷入到低劣的局部最优。
发明内容
要解决的技术问题
为了避免现有方法的不足之处，本发明提出一种三维大脑磁共振图像大脑皮层表面最大主方向场弥散方法，适用于人类三维大脑磁共振图像重构出的、由三角形表示的大脑皮层表面上最大主方向场弥散。
技术方案
本发明的基本思想是：将大脑皮层表面上的最大主方向场弥散问题公式化为一个能量最小化问题，进而转化为一个离散标记问题，即给大脑皮层表面上的每一个顶点赋予一个标记来表示弥散的最大主方向，该离散标记问题可以利用alpha-expansion图割方法有效的求解，因为图割方法可以保证得到某些类型的能量函数的强局部最优。
一种三维大脑磁共振图像的大脑皮层表面上最大主方向场弥散方法，其特征在于步骤如下：
步骤1对三维大脑磁共振图像进行预处理和大脑皮层表面重建：利用变形模型方法去除脑壳，利用配准方法去除非大脑组织，利用高斯混合模型和马尔科夫随机场方法对图像进行组织分割，得到白质、灰质和脑脊髓液三种组织类型表示的图像，利用基于图的方法对白质图像进行拓扑校正，利用Marching Cubes方法重建由三角形表示的大脑皮层表面；
步骤2计算大脑皮层表面上顶点的主曲率、主方向以及主曲率的导数：利用有限差分方法计算大脑皮层表面上顶点的主曲率、主方向以及主曲率的导数，当每个顶点的最大主曲率在最大主方向的导数为正值时，将该顶点的最大主方向改为相反方向；
步骤3对最大主方向场进行弥散：利用alpha-expansion图割方法最小化能量函数
E=Σ(x,y)∈NVx,y(lx,ly)+Σx∈SDx(lx)]]>
其中，![]()
![]()
x和y为大脑皮层表面上的顶点坐标，p(x)为顶点x处原始的最大主方向；(x，y)∈N表示x和y为相邻顶点的集合，S为大脑皮层表面上的所有顶点的集合；||·||为2范数；![]()
为顶点x弥散的最大主方向，lx为顶点x弥散的最大主方向在解空间的标识且lx∈L；![]()
为顶点y弥散的最大主方向，ly为顶点y弥散的最大主方向在解空间的标识且ly∈L；
所述的L＝{l1，l2，...，ln}为离散的弥散的最大主方向场的解空间Θ＝{v1，v2，...，vn}的标识，其中n＝nθ×(nφ-2)+2；nθ为x-y平面中离散角度的数目为12～36；nφ为z轴方向中离散角度的数目为9～18；
所述解空间Θ中v1＝(0，0，1)，vn＝(0，0，-1)；
vi＝(cos(((i-2)％nθ)·2π/nθ)sin((1+(i-2)/nθ)·π/nφ)，sin(((i-2)％nθ)·2π/nθ)sin((1+(i-2)/nθ)·π/nφ)，cos((1+(i-2)/nθ)·π/nφ))，，1＜i＜n；
所述g(x)＝exp(-λ·|c(x)|)，所述g(y)＝exp(-λ·|c(y)|)，所述h(x)＝1-g(x)，其中λ为权重参数为5.0～10.0，c(·)为最大主曲率；
利用alpha-expansion图割方法求解，得到弥散的最大主方向场；
步骤4将弥散的最大主方向场投影到每个顶点的切平面内。
有益效果
本发明提出的三维大脑磁共振图像大脑皮层表面最大主方向场弥散方法，使得三维大脑磁共振图像中大脑皮层表面上最大主方向场弥散方法的可行性体现在：首先，随着磁共振成像设备的精度不断提高和三维大脑磁共振图像的预处理方法进一步成熟，获取几何结构精确、拓扑结构正确的大脑皮层表面相对容易；同时，将最大主方向场弥散看做一个能量最小化问题，并转换为一个离散标记问题，利用图割方法有效求解该离散标记问题是可行的。
本发明相对于其它方法具有以下优点：1、该方法在大脑皮层表面上不同区域的平滑项和数据项的权重设为不同，这样可以更好的保持最大主方向场中固有的几何结构和不连续性；2、将最大主方向场弥散看做一个能量最小化问题，并转换为一个离散标记问题，利用图割方法可以有效求解该能量函数的强局部最优。
附图说明
图1：本发明实施例的三维大脑磁共振图像大脑皮层表面上最大主方向场弥散，
(a)显示了实施例图像中大脑皮层表面上的最大主曲率，
(b)显示了图1(a)中矩形区域中大脑皮层表面上的最大主方向场，
(c)显示了图1(b)中对应的大脑皮层表面上的弥散的最大主方向场；
图2：仿真表面上的最大主方向场弥散的一个例子，
(a)显示了三维大脑磁共振图像仿真表面上的理想最大主方向场，
(b)显示了三维大脑磁共振图像仿真噪声的最大主方向场，
(c)显示了三维大脑磁共振图像仿真弥散的最大主方向；
图3：三维大脑磁共振图像大脑皮层表面上的最大主方向场弥散前和弥散后的一致性比较实例，
(a)显示实施例中三维大脑磁共振图像大脑皮层表面上的最大主方向场的一致性
(b)显示实施例中三维大脑磁共振图像大脑皮层表面上弥散的最大主方向场的一致性；
图4：本实例大脑皮层表面上的最大主方向场弥散前和弥散后的一致性的分布直方图；
图5：12个大脑皮层表面上的最大主方向场弥散前和弥散后的平均一致性的比较；
具体实施方式
现结合实施例、附图对本发明作进一步描述：
根据本发明提出的基于图割方法的大脑皮层表面上的最大主方向场弥散方法，我们用C++语言实现了一个该方法的原型系统。图像数据的来源是：实际中正常人的三维大脑磁共振图像。
具体的实施步骤如下：
1.预处理和大脑皮层表面重建：
利用变形模型方法去除脑壳，利用配准方法去除非大脑组织，利用高斯混合模型和马尔科夫随机场方法对图像进行组织分割，得到白质、灰质和脑脊髓液三种组织类型表示的图像，利用基于图的方法对白质图像进行拓扑校正，利用Marching Cubes方法重建由三角形表示的大脑皮层表面；
2.计算大脑皮层表面上顶点的主曲率、主方向以及主曲率的导数：
利用有限差分方法计算大脑皮层表面上顶点的主曲率、主方向以及主曲率的导数，检验每个顶点的最大主曲率在最大主方向的导数是否为正值，如果为正值，则翻转该顶点的最大主方向为相反方向；图1(a)显示了一例大脑皮层表面上的最大主曲率；图1(b)显示了图1(a)中矩形区域中大脑皮层表面上的最大主方向场。
3.对最大主方向场进行弥散：
大脑皮层表面上最大主方向场的弥散可以表示为如下的能量函数
E=Σ(x,y)∈N(g(x)+g(y))||v(x)-v(y)||+Σx∈Sh(x)||v(x)-p(x)||]]>
其中，v(x)＝(u(x)，v(x)，w(x))是顶点x处弥散的最大主方向，v(y)是顶点y处弥散的最大主方向，p(x)是顶点x处原始的最大主方向。(x，y)∈N表示x和y为相邻顶点的集合，S表示大脑皮层表面上的所有顶点的集合。||·||表示2范数。g(x)＝exp(-λ·|c(x)|)，g(y)＝exp(-λ·|c(y)|)，h(x)＝1-g(x)，其中λ为权重参数为8.0，c(·)为最大主曲率。在大脑皮层表面上强烈弯曲的区域，对应于脑沟底部和脑回冠部区域，该处c(·)的绝对值较大，最大主方向场可以相对可靠的计算；而在平坦的区域所计算的最大主方向场的噪声较大，比较不可靠。因此g(·)设为在脑沟底部和脑回冠部的值较小，而h(·)在脑沟底部和脑回冠部的值较大。![]()
为数据项，用于使弥散的最大主方向场和原始的最大主方向场尽可能相似。![]()
为平滑项，用于使弥散的最大主方向场平滑变化。可以看出，在脑沟底部和脑回冠部，数据项的权重较大，迫使弥散的最大主方向场与原始最大主方向场相似；在平坦的皮层区域，平滑项的权重较大，迫使弥散的最大主方向场平滑的变化。为了将最大主方问场弥散转化为一个离散标记问题，将v(x)表示为v(x)＝(cosθsinφ，sinθsinφ，cosφ)，其中θ∈[0，2π)，φ∈[0，π]。令θ＝i·2π/nθ，i∈[0，nθ-1]，nθ表示θ中离散角度的数目，设为12。令φ＝j·π/nφ，j∈[0，nφ]，nφ表示φ中离散角度的数目，设为9。令L＝{l1,l2，...，ln}表示离散的弥散的最大主方向场的解空间Θ＝{v1，v2，...，vn}，其中n＝nθ×(nφ-2)+2。以上的能量函数可以转化为：
E=Σ(x,y)∈NVx,y(lx,ly)+Σx∈SDx(lx)]]>
其中，![]()
![]()
利用alpha-expansion图割方法有效的求解该能量函数。
步骤4将弥散的最大主方向场投影到切平面内。
图1(c)显示了图1(b)中对应的大脑皮层表面上的弥散的最大主方向场。
为了定量测试大脑皮层表面上最大主方向场弥散方法的效果，将该方法应用到一个仿真的表面上。在仿真的表面上可以精确解析计算每个顶点的理想最大主方向，然后对每个顶点的位置添加随机噪声，并利用有限差分方法计算添加噪声后每个顶点的最大主方向。图2(a)、(b)和(c)分别显示了仿真表面上的理想最大主方向场、噪声的最大主方向场和弥散的最大主方向。通过以下公式计算弥散后的最大主方向场、添加噪声的最大主方向场与理想最大主方向场的平均角度差：
AngleDiff=1Σx∈S1Σx∈Sarccos⟨p(x)·v(x)⟩]]>
其中，p(x)表示理想最大主方向，v(x)表示弥散的最大主方向或添加噪声的最大主方向。添加噪声的最大主方向场与理想最大主方向场的平均角度差为67.7度，而弥散后的最大主方向场与理想最大主方向场的平均角度差为6.3度，表明该方法达到了很好的弥散结果。
将该方法应用到12个正常人的大脑皮层表面上。定义一致性来评价弥散结果：
Coh(x)=1Σ(x,y)∈N1Σ(x,y)∈Nv(x)·v(y)]]>
其中，v(x)表示顶点x的方向，(x，y)∈N表示相邻顶点集合。大的一致性值表明方向场更加平滑和一致。因此，该值可以在某种程度上反映弥散结果。图3(a)和(b)分别显示一例大脑皮层表面上的最大主方向场和弥散的最大主方向场的一致性。在理想情况下，在弥散的最大主方向场中，一致性在平坦的皮层区域应该接近于1.0，一致性在脑回冠部和脑沟底部接近于0。可以明显的看出，和原始的最大主方向场相比，弥散的最大主方向场中，一致性在平坦的皮层区域增加了，而在脑沟底部和脑回冠部区域很好的保持。同时计算了大脑皮层表面上一致性的分布直方图，如图4所示，可以看出相比于原始的最大主方向场，弥散的最大主方向场中一致性明显增大。图5显示了在12个大脑皮层表面上，平均一致性在最大主方向场弥散前和弥散后的变化，可以看出平均一致性明显的增大，表明了本方法的有效性。