CN101515368A - 三维大脑磁共振图像的大脑皮层表面上脑沟盆地分割方法 - Google Patents

Abstract

本发明涉及一种三维大脑磁共振图像的大脑皮层表面上脑沟盆地分割方法，技术特征在于：用于在三角化的大脑皮层表面上弥散原始杂乱的主方向场，以生成一个平滑的最大主方向流场。最大主方向弥散的核心思想是：通过最小化一个能量函数，在脑沟和脑回区域，这里最大主曲率的绝对值很大，产生的最大主方向流场应该接近于原始的主方向场；在其他平坦的大脑皮层区域，产生的最大主方向流场应该平滑的变化。本发明相对于其它方法的主要优点是：最大主方向流场跟踪方法，充分利用了大脑皮层表面的几何结构信息，能够自然的形成脑沟盆地间的精确边界，尽管脑沟盆地有可能是非对称结构。

Description

三维大脑磁共振图像的大脑皮层表面上脑沟盆地分割方法

技术领域

本发明涉及一种三维大脑磁共振图像的大脑皮层表面上脑沟盆地分割方法，属于医学图像处理，计算神经解剖学等领域。适用于人类三维大脑磁共振图像重构出的三角化的大脑皮层表面上的脑沟盆地的分割。

背景技术

人类大脑皮层是一个极其复杂卷曲的解剖结构，主要由脑沟和脑回构成，分别对应于大脑皮层上的谷和脊。尽管不同人之间脑沟和脑回的精确几何模式变化很大，但是最主要的几条脑沟和脑回是大脑皮层上共有的解剖界标。因此，主要的脑沟和脑回已经被广泛的用于辅助非线性大脑磁共振图像配准，分析正常人大脑的解剖结构变化规律，以及用来区分正常人和疾病患者。但是手工分割和标定脑沟及其费时，并且容易受到外界主观因素的影响。

近些年来脑沟盆地的自动分割已经成为研究的热点课题。脑沟盆地是脑回冠部所界定的区域，并且相邻的脑沟盆地在大脑皮层表面上的脑回冠部曲线相会。将大脑皮层表面分割为不同的脑沟盆地产生了该大脑皮层表面的一个完整分割。几种脑沟盆地的分割方法已经被提出来，但仍存在很多问题。方法1：在三维大脑磁共振图像上利用区域增长方法提取脑沟盆地，该方法首先对大脑图像进行分割得到大脑白质图像，然后对该图像进行形态学闭操作，将该结果减去大脑白质图像得到脑沟内部区域，接着计算脑沟内部区域的距离变换作为脑沟深度，最后从脑沟最深处开始进行区域增长并进行区域合并得到脑沟盆地。其缺点是区域增长和区域合并中参数较难控制，而且由于脑沟盆地可能是非对称结构，因此区域增长得到的脑沟盆地的边界并非一定位于真正的脑回冠部区域；方法2：基于大脑皮层表面上脑沟深度的分水岭方法，该方法首先利用主动表面方法找到脑回区域，然后利用分水岭方法提取脑沟区域，接着利用启发式规则合并过分割的脑沟区域，最后通过计算脑沟盆地的沿着表面的影响区域作为脑沟盆地。其缺点在于，分水岭方法很容易产生过分割现象，将一个脑沟区域分割为多个脑沟区域，而且启发式规则合并过分割脑沟区域较难控制，同上由于脑沟盆地可能是非对称结构，因此脑沟区域的沿着表面的影响区域的边界并不一定位于真正的脑回冠部区域。

目前已有的大脑皮层脑沟盆地分割方法具有以下两个主要的缺陷：其一、由于脑沟盆地是非对称结构，因此区域增长或分水岭方法得到的脑沟盆地的边界并非真正的脑回冠部区域，因此所分割的脑沟盆地的边界不精确。其二、利用区域增长或分水岭方法进行脑沟盆地分割容易产生过分割，需要后处理进行区域合并，但是区域合并的规则和参数较难控制。

发明内容

要解决的技术问题

为了避免现有分割方法中脑沟盆地的边界不精确的不足之处，本发明提出一种三维大脑磁共振图像的大脑皮层表面上脑沟盆地分割方法。

技术方案

本发明的基本思想是：大脑皮层表面上最大主方向指向最大主曲率最陡减小的方法，同时脑回冠部区域和脑沟底部区域分别具有大的正的和负的最大主曲率，我们可以从脑回冠部区域开始沿着最大主方向一直到达脑沟底部区域。所有流到同一个脑沟底部区域的顶点自然的划分为一个脑沟盆地，这样我们就能将精确的将大脑皮层表面分割为不同的脑沟盆地。但是实际中计算得到最大主方向场中噪声较大，尤其在平坦的皮层表面区域最大主方向场杂乱无章，因为该区域迪两个主曲率都很小，微小的结构变化将导致最大主方向的剧烈变化。为了处理这个问题，受到灰度图像中梯度向量场弥散方法的启发，我们设计了一个新的方法，用于在三角化的大脑皮层表面上弥散原始杂乱的主方向场，以生成一个平滑的最大主方向流场。最大主方向弥散的核心思想是：通过最小化一个能量函数，在脑沟和脑回区域，这里最大主曲率的绝对值很大，产生的最大主方向流场应该接近于原始的主方向场；在其他平坦的大脑皮层区域，产生的最大主方向流场应该平滑的变化。

本发明的技术特征在于步骤如下：

步骤1对三维大脑核磁共振图像进行预处理：利用可变性模型方法去除脑壳，利用配准方法去除非大脑组织，利用高斯混合模型方法对大脑图像进行组织分割，得到白质，灰质和脑脊髓液三种组织类型表示的图像；

步骤2大脑皮层表面重建：利用Marching Cubes方法从组织分割后的大脑图像中重构三角化的大脑皮层表面；

步骤3利用有限差分方法估计三角化的大脑皮层表面上每个顶点的最大主曲率和最大主方向；

步骤4：在满足v(x)·n(x)＝0的条件下，最小化能量函数

$\mathrm{\ϵ}=\underset{x\∈S}{\∫}\mathrm{\λ}{|\▿v\left(x\right)|}^2+f\left(x\right){|v\left(x\right)-p\left(x\right)|}^2\mathrm{dx},$

得到最大主方向流场v(x)＝(u(x)，v(x)，w(x))：其中λ是权重参数为0.05～0.15，

是梯度算子，p(x)是顶点x的最大主方向，f(x)是顶点x的最大主曲率强度的递减函数，f(x)＝|c(x)|，其中c(x)是顶点x的最大主曲率；n(x)是顶点x的法向量，S表示大脑皮层表面上所有顶点的集合；所述的v(x)·n(x)＝0限定最大主方向流场在大脑皮层表面上的切平面内；

步骤5利用大脑皮层表面上最大主方向流场跟踪方法，得到脑沟盆地分割：具体步骤如下：在最大主方向流场中，计算大脑皮层表面上任意一个顶点x所流经的下一个顶点 $x^{\&prime;}=\underset{x_i}{\min }(\arccos <v\left(x\right)\&CenterDot;\stackrel{\&RightArrow;}{{\mathrm{xx}}_i}/|\left|\stackrel{\&RightArrow;}{{\mathrm{xx}}_i}\right||>),$ 其中x_i是顶点x的第一圈邻域；如果两个连续流向量的夹角小于π/2，利用以上公式计算下一个流经的顶点，流场跟踪步骤继续执行；否则，流场跟踪步骤就停止，到达了一个脑沟底部区域，将流场跟踪停止顶点标记为脑沟底部顶点，同时将该流场跟踪停止顶点记录为流路径所经过的所有顶点的对应脑沟底部顶点；利用连通成分分析方法，将所有同一连通的脑沟底部区域标记为同一颜色，将所有对应于同一个脑沟底部顶点的顶点标记为与其对应脑沟底部顶点相同的颜色，得到大脑皮层表面上脑沟盆地的分割。

所述的计算最小化能量函数通过将v看作是时间的方程，计算得到最大主方向流场，具体步骤如下：

步骤1：根据v(x，t)(t＝0，1，2...)，将方程 $v_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}v(x,t)-(v(x,t)-p\left(x\right))f\left(x\right)$ 分解为以下方程 $\left\{\begin{array}{c}{u}_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}u(x,t)-(u(x,t)-p_u\left(x\right))f\left(x\right)\\ v_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}v(x,t)-(v(x,t)-p_v\left(x\right))f\left(x\right)\\ w_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}w(x,t)-(w(x,t)-p_w\left(x\right))f\left(x\right)\end{array}\right.$

得到弥散的最大主方向场v(x，t+1)，上述公式中，v_t(x，t)表示v(x，t)对于时间t的偏导数，拉普拉斯项 ${\&dtri;}^{2}u\left(x\right)=\underset{x_i}{\mathrm{\&Sigma;}}{w}_{x,x_i}(u\left(x\right)-u\left(x_i\right)),$ 其中x_i是顶点x周围一圈的邻域，

是一个权重系数计算为顶点x和x_i之间距离的倒数。同理，可以估计

和

步骤2：投影弥散的最大主方向场v(x，t+1)到切平面内；

步骤3：归一化投影的弥散的最大主方向场v(x，t+1)；

步骤4：将v(x，t)替换为投影的弥散的归一化的v(x，t+1)，重复步骤1至3共30～50次迭代停止。

有益效果

本发明提出的大脑皮层表面上脑沟盆地分割方法的可行性体现在，首先，随着磁共振成像设备的精度不断提高和三维大脑磁共振图像的预处理方法进一步成熟，获取几何结构精确，拓扑结构正确的大脑皮层表面相对容易；同时，大脑皮层表面上最大主方向指向最大主曲率最陡减小的方法，同时最大主方向场弥散进一步平滑了该向量场，将脑沟底部和脑回冠部的最大主方向传播到其他平坦的皮层表面区域，可以通过流场跟踪方法到达脑沟底部区域，所有流到同一个脑沟底部区域的顶点划分为同一个脑沟盆地，因此利用最大主方向流场跟踪来分割脑沟盆地是可行的。

本发明相对于其它方法的主要优点是：最大主方向流场跟踪方法，充分利用了大脑皮层表面的几何结构信息，能够自然的形成脑沟盆地间的精确边界，尽管脑沟盆地有可能是非对称结构。

附图说明

图1：本发明方法的基本流程图

图2：1个大脑皮层左半球内表面上最大主曲率和最大主方向分布图

图3：12个真实正常人大脑皮层左半球上内表面上脑沟盆地分割结果，其中每个颜色表示一个脑沟区域。

图4：两个专家手工标定的和自动分割的中央沟盆地的一个比较。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

首先对三维大脑磁共振图像进行预处理和大脑皮层表面重建：包括去除脑壳和非大脑组织，对大脑图像进行脑组织分割(分割为白质，灰质和脑脊髓液三种类型)，从组织分割后的大脑图像中重构几何结果精确，拓扑结构正确的大脑皮层表面，该皮层表面由一系列顶点和三角形表示。然后，估计大脑皮层表面上顶点的最大主曲率和最大主方向。接着，通过最小化一个能量函数，弥散大脑皮层表面上最大主方向场来产生一个平滑的切向最大主方向流场。最后，根据弥散得到的最大主方向流场，利用流场跟踪方法将大脑皮层表面分割为不同的脑沟盆地。

大脑皮层表面上最大主方向场弥散的原理是：通过最小化一个能量函数，在脑沟和脑回区域，这里最大主曲率的绝对值很大，产生的最大主方向流场应该接近于原始的主方向场；在其他平坦的皮层区域，产生的最大主方向流场应该平滑的变化。为了最小化该能量函数，利用计算变分将该问题转为求解一个偏微分方程，然后通过将偏微分方程看作是时间的方程来求解得到最终的最大主方向流场。

大脑皮层表面上最大主方向流场跟踪的原理是：大脑皮层表面上最大主方向指向最大主曲率最陡减小的方法，同时脑回冠部区域和脑沟底部区域分别具有大的正的和负的最大主曲率，可以从脑回冠部区域沿着最大主方向一直到达脑沟底部区域。如果流路径中两个连续流向量的夹角小于π/2，流场跟踪步骤继续执行。否者，两个连续流向量将会互相指向对方，流场跟踪步骤就停止。这时候就到达了一个脑沟底部区域。所有流到同一个脑沟底部区域的顶点自然的划分为一个脑沟盆地，这样我们就能将大脑皮层表面分割为不同的脑沟盆地。

根据本发明提出的基于最大主方向流场跟踪的大脑皮层表面上脑沟盆地的分割方法，我们用C++语言实现了一个脑沟盆地分割的原型系统。图像数据的来源是：实际中正常人的三维大脑磁共振图像。

本发明整个流程可以参考附图1，具体的实施步骤如下：

1.预处理和大脑皮层表面重建：

对三维大脑磁共振图像进行去除头骨和非大脑组织，大脑组织分割，重建几何结果精确，拓扑结构正确的大脑皮层三角化表面。

2.大脑皮层表面上最大主曲率和最大主方向估计：

利用有限差分方法估计大脑皮层表面上每个顶点的最大主曲率，最大主方向，以及最大主曲率在最大主方向上的导数。如果某个顶点最上大主曲率在最大主方向上的导数为正，将最大主方向翻转为其相反方向，这样所有的最大主方向将均匀的指向脑沟底部区域。附图2显示了1个大脑皮层左半球内表面(白质和灰质交界面)上最大主曲率和最大主方向分布图，一个大脑皮层内表面(大脑中白质和灰质的交界面)上估计的最大主曲率和最大主方向。(a)是最大主曲率图；(b)是图(a)中矩形框定区域的最大主方向图。颜色条显示在图上方。

3.大脑皮层表面上最大主方向场弥散：

给定一个重建的大脑皮层表面和估计的最大主曲率和最大主方向，在顶点x，最大主方向流场v(x)＝(u(x)，v(x)，w(x))在世界坐标系下定义为以下能量函数的解：

$\mathrm{\&epsiv;}=\underset{x\&Element;S}{\&Integral;}\mathrm{\&lambda;}{|\&dtri;v\left(x\right)|}^2+f\left(x\right){|v\left(x\right)-p\left(x\right)|}^2\mathrm{dx}$

在满足v(x)·n(x)＝0的条件下，λ是一个权重参数，

是梯度算子。p(x)是顶点x的最大主方向，f(x)是顶点x的最大主曲率强度的递减函数，设定为：f(x)＝|c(x)|，其中c(x)是顶点x的最大主曲率。n(x)是顶点x的法向量，S表示大脑皮层表面上所有顶点的集合。公式v(x)·n(x)＝0限定生成的主方向流场在大脑皮层表面上的切平面内。根据以上的能量函数，在平坦或均匀的表面区域，这里最大主曲率的绝对值很小，能量函数由第一个偏导数项所主导，确保流场平滑变化。在最大主曲率绝对值大的表面区域，对应于脑沟底部和脑回冠部区域，该区域估计的出的主方向比较可靠，能量函数由第二项所主导，确保流场尽可能接近于原始的主方向场。参数λ确定第一个平滑项和第二个保真项的折中关系，设定为0.1。

为了最小化公式以上的能量函数，利用计算变分得到以下的偏微分方程：

$\mathrm{\&lambda;}{\&dtri;}^{2}v\left(x\right)-(v\left(x\right)-p\left(x\right))f\left(x\right)=0$

以上偏微分方程可以通过将v看作是时间的方程来求解：

$\left\{\begin{array}{c}{v}_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}v(x,t)-(v(x,t)-p\left(x\right))f\left(x\right)\\ v(x,0)=p\left(x\right)\end{array}\right.$

这里v_t(x，t)表示v(x，t)对于时间t的偏导数。上面的方程可以进一步分解为：

$u_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}u(x,t)-(u(x,t)-p_u\left(x\right))f\left(x\right)$

$v_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}v(x,t)-(v(x,t)-p_v\left(x\right))f\left(x\right)$

$w_t(x,t)=\mathrm{\&lambda;}{\&dtri;}^{2}w(x,t)-(w(x,t)-p_w\left(x\right))f\left(x\right)$

表面上的拉普拉斯项估计为：

${\&dtri;}^{2}u\left(x\right)=\underset{x_i}{\mathrm{\&Sigma;}}{w}_{x,x_i}(u\left(x\right)-u\left(x_i\right))$

这里x_i是顶点x周围一圈的邻域，

是一个权重系数计算为顶点x和x_i之间距离的倒数。利用相似方法估计

和

最大主方向场弥散方法总结如下：

1)求解分解的方程得到弥散的最大主方向场；

2)将弥散的最大主方向场投影到切平面内；

3)归一化弥散的最大主方向场；

4)重复步骤1-3，直到迭代停止。

4.大脑皮层表面上最大主方向流场跟踪：

给定一个顶点x和给顶点的主方向流向量v(x)，在主方向流场中x所流经的下一个顶点x′计算为：

$x^{\&prime;}=\underset{x_i}{\min }(\arccos <v\left(x\right)\&CenterDot;\frac{\stackrel{\&RightArrow;}{{\mathrm{xx}}_i}}{\left|\right|\stackrel{\&RightArrow;}{{\mathrm{xx}}_i}\left|\right|}>)$

这里x_i是顶点x的第一圈邻域。在流路径中两个连续流向量的夹角计算为：

θ＝arccos<v(x)·v(x′)>

当两个连续流向量的夹角小于π/2，或者说两个连续流向量的内积为正，流场跟踪步骤继续执行。否者，两个连续流向量将会互相指向对方，即内积为负，流场跟踪步骤就停止，这时候就到达了一个脑沟底部区域。

大脑皮层表面上最大主方向流场跟踪方法总结如下：

1)选取一个不在任何处理过顶点的流轨迹上的顶点作为初始顶点x⁰；

2)根据xⁿ计算xⁿ⁺¹(n＝0，1，2...)；

3)计算xⁿ⁺¹和xⁿ的流向量的夹角θ_n。如果θ_n大于π/2，停止；否则，如果xⁿ⁺¹已经在处理过顶点的流轨迹上，将从顶点x⁰到xⁿ与xⁿ⁺¹对应的流轨迹中止点相关联，停止；

4)将xⁿ替换为xⁿ⁺¹，返回到第2步。

利用这种方式，在流轨迹的顶点就定义了一条平滑的路径通向脑沟底部区域。所有流到同一个脑沟底部区域的顶点自然的产生该脑沟区域对应的脑沟盆地，通过连通性分析来判断流场跟踪终点是否为在同一个脑沟底部区域。最后，如果两个相邻脑沟盆地的边界处的平均最大主曲率为负，就将其合并为一个脑沟盆地，因为通常脑沟盆地又脑回冠部界定，而该区域的最大主曲率为正。

为了测试该脑沟盆地分割方法的准确度，我们将该方法用于12个真实正常人的大脑磁共振图像所重建出的大脑皮层表面。图3显示12个正常人大脑皮层表面左半球上脑沟盆地分割结果，其中每个脑沟区域被标为一个颜色。为了定量的评价该方法，我们以两个专家手工标定的大脑左半球内表面上的中央沟盆地为标准。利用区域重叠率来比较自动分割的和专家手工标定的脑沟盆地。区域重叠率定义为：

$O(R_a,R_m)=\frac{S(R_a\&cap;R_m)}{(S\left(R_a\right)+S\left(R_m\right))/2}$

这里R_a是自动提取的脑沟盆地，R_m是手工标定的脑沟盆地。∩是两个区域的交算子。S(.)是区域的面积。表1显示了12个正常人大脑左半球内表面上中央沟盆地的区域重叠率情况。12个人中手工和自动结果的平均区域重叠率高于0.96，表明该方法能够精确的分割出脑沟盆地。表1中也显示了12个人大脑左半球内表面上所提取出的中央沟盆地的面积。该面积是通过求和所有位于脑沟盆地内的三角形面积。图4显示了左半球内表面上自动分割的和专家手工分割的中央沟盆地的一个比较。(a)和(b)是两个专家手工标定的脑沟盆地；(c)是自动分割的脑沟盆地。(d)是将自动分割结果叠加在大脑皮层表面上。从实验结果看出：我们的脑沟区域分割方法具有很好的表现。

表1  12个正常人大脑左半球内表面上中央沟盆地的面积和区域重叠率

Claims (2)

1.三维大脑磁共振图像的大脑皮层表面上脑沟盆地分割方法，其特征在于：

步骤1对三维大脑核磁共振图像进行预处理：利用可变性模型方法去除脑壳，利用配准方法去除非大脑组织，利用高斯混合模型方法对大脑图像进行组织分割，得到白质，灰质和脑脊髓液三种组织类型表示的图像；

步骤2大脑皮层表面重建：利用Marching Cubes方法从组织分割后的大脑图像中重构三角化的大脑皮层表面；

步骤3利用有限差分方法估计三角化的大脑皮层表面上每个顶点的最大主曲率和最大主方向；

步骤4：在满足v(x)·n(x)＝0的条件下，最小化能量函数

$\mathrm{\&epsiv;}=\underset{x\&Element;S}{\&Integral;}\mathrm{\&lambda;}{|\&dtri;v\left(x\right)|}^2+f\left(x\right)|v\left(x\right)-p\left(x\right){|}^2\mathrm{dx},$

得到最大主方向流场v(x)＝(u(x)，v(x)，w(x))：其中λ是权重参数为0.05～0.15，

是梯度算子，p(x)是顶点x的最大主方向，f(x)是顶点x的最大主曲率强度的递减函数，f(x)＝|c(x)|，其中c(x)是顶点x的最大主曲率；n(x)是顶点x的法向量，S表示大脑皮层表面上所有顶点的集合；所述的v(x)·n(x)＝0限定最大主方向流场在大脑皮层表面上的切平面内；

步骤5利用大脑皮层表面上最大主方向流场跟踪方法，得到脑沟盆地分割：具体步骤如下：在最大主方向流场中，计算大脑皮层表面上任意一个顶点x所流经的下一个顶点 $x^{\&prime;}=\underset{x_i}{\min }(\arccos <v\left(x\right)\&CenterDot;\stackrel{\&RightArrow;}{{\mathrm{xx}}_i}/|\left|\stackrel{\&RightArrow;}{{\mathrm{xx}}_i}\right||>),$ 其中x_i是顶点x的第一圈邻域；如果两个连续流向量的夹角小手π/2，利用以上公式计算下一个流经的顶点，流场跟踪步骤继续执行；否则，流场跟踪步骤就停止，到达了一个脑沟底部区域，将流场跟踪停止顶点标记为脑沟底部顶点，同时将该流场跟踪停止顶点记录为流路径所经过的所有顶点的对应脑沟底部顶点；利用连通成分分析方法，将所有同一连通的脑沟底部区域标记为同一颜色，将所有对应于同一个脑沟底部顶点的顶点标记为与其对应脑沟底部顶点相同的颜色，得到大脑皮层表面上脑沟盆地的分割。

2、根据权利要求1所述的三维大脑磁共振图像的大脑皮层表面上脑沟盆地分割方法，其特征在于：所述的计算最小化能量函数通过将v看作是时间的方程，计算得到最大主方向流场，具体步骤如下：

步骤1：根据v(x，t)(t＝0，1，2...)，将方程 $v_t(x,t)={\mathrm{\&lambda;}\&dtri;}^{2}v(x,t)-(v(x,t)-p\left(x\right))f\left(x\right)$

分解为以下方程 $\left\{\begin{array}{c}{u}_t(x,t)={\mathrm{\&lambda;}\&dtri;}^{2}u(x,t)-(u(x,t)-p_u\left(x\right))f\left(x\right)\\ v_t(x,t)={\mathrm{\&lambda;}\&dtri;}^{2}v(x,t)-(v(x,t)-p_v\left(x\right))f\left(x\right)\\ w_t(x,t)={\mathrm{\&lambda;}\&dtri;}^{2}w(x,t)-(w(x,t)-p_w\left(x\right))f\left(x\right)\end{array}\right.$

得到弥散的最大主方向场v(x，t+1)，上述公式中，v_t(x，t)表示v(x，t)对于时间t的偏导数，拉普拉斯项 ${\&dtri;}^{2}u\left(x\right)=\underset{x_i}{\mathrm{\&Sigma;}}{w}_{x,x_i}(u\left(x\right)-u\left(x_i\right)),$ 其中x_i是顶点x周围一圈的邻域，

是一个权重系数计算为顶点x和x_i之间距离的倒数。同理，可以估计

和

步骤2：投影弥散的最大主方向场v(x，t+1)到切平面内；

步骤3：归一化投影的弥散的最大主方向场v(x，t+1)；

步骤4：将v(x，t)替换为投影的弥散的归一化的v(x，t+1)，重复步骤1至3共30～50次迭代停止。

Images