Snake,_Shape,_and_Gradient_Vector_Flow(报告).docx

Snake, Shape, and Gradient Vector Flow 姓名：崔光民 学号：I201122047专业：计算机应用技术 Snake, Shape, and Gradient Vector Flow Chenyang Xu, Student Member IEEE, and Jerry L. Prince, Senior Member, IEEE摘要 Snake，即活动轮廓，广泛地用于计算机视觉和图像处理应用， 特别对于配置对象的轮廓。 但是对于凹轮廓怎么进行初始化、怎么向它收敛， 这两个问题仍然影响这个方法的应用。本文提出了把新的外部力（external force）引进到活动的轮廓， 以解决上述两个问题。这个外部力叫梯度向量流（gradient vector flow-GVF），以灰度的梯度向量扩散或者图像的二值边缘映射。 这种外部力跟传统的snake外部力根本不一样，因为传统的snake外部力不能写成某种保守函数的负梯度。除此之外，与这种外部力对应的snake， 不是由别的公式，而是由平衡力条件直接构成。 把这个方法用于二维（2D）对象和三维（3D）对象，这实验证明了GVF方法的捕捉范围比较大，而且让snake 变形到凹轮廓。索引活动轮廓模型（Active contour models），可变型的曲面模型（deformable surface models），梯度向量流(gradient vector flow), 图像分割（image processing），外形表现和向量（shape representation and vector）， snake。1. 前言Snake1，即活动的轮廓，是在图像领域中定义的曲线，它受到内部力（internal force）和外部力（external force）的影响。内部力是由曲线本身引起的，而外部力是用图像资料进行计算的。 内部力和外部力使得snake适应对象轮廓或者图像的特征。 Snake用于很多方面：例如，边缘检测1、形态模型化2，3、分割4，5、运动跟踪4，6。 按文献，把活动轮廓可以分成两个类：参数活动轮廓（parametric active contours）1和几何活动轮廓（geometric active contours）7-9。本文仅关注参数活动轮廓。 参数活动轮廓在图像内部产生一种参数曲线，并且使它向想要的特征，例如边缘，变形。 一般，这种变形是由保守力引起的。 保守力是势函数的负梯度。 除此以外，类似于压力的其他力跟保守力在一起构成外部力。 以外，还有内部力保持曲线。 对参数活动轮廓算法，存在两个难题。 第一，是要在边缘附近进行初始化，否则会收敛到错结果。 有若干方法，它包括多重分辨率方法11、压力方法10和距离势力12。 这思路用于增加外部力的捕捉域并使得曲线向想要的轮廓变形。 第二，是很难把活动轮廓向凹地方拉开。 有若干方法提出来，例如压力方法（pressure forces）10、控制点方法（control points）13、领域适应（domain-adaptivity）15、指向性吸引方法（dirctional attraction）14 和 solenidal 方法，但是不能解决这个问题。 但是，大多数方法仅解决了某一个问题，却造成新的问题。例如，多重分辨率方法涉及捕捉域的问题，但是没提出怎么向别的分辨率过去。 至于压力方法，它把曲线向凹轮廓推进，但是不能太强推，否则“弱”边缘会断裂17。 还有，为了推出或者推进，要进行初始化，这时，初始化条件很严格。 本文，对上述两个问题，提出了活动轮廓模型中的新形式的外部力。它叫做梯度向量流（GVF - gradient vector flow），是把某种能量函数最小化的密度大的向量场。 用一双可分离的篇微方程式求这个最小化问题的解，它使得灰度或者由图像算出来的二值化边缘映射的梯度向量扩散。以GVF作为外部力的活动轮廓叫做GVF snake。 这种外部力跟传统的snake外部力根本不一样，因为传统的snake外部力不能写成某种保守函数的负梯度。因此，它不能用标准能量最小化方式表示，但是能用平衡力条件直接表示。 GVF的一个特点是对初始化不太敏感，并且稳定地收敛于凹轮廓。 从本文容易看出，在对象内部初始化也可以，在外面也可以，穿过对象也都可以。 和压力方法不同，GVF snake不需要关于变形方向的先验（prior）信息， 也就是说膨胀还是收缩。 还有，GVF snake的捕捉域很大。 这个区域是用扩散处理求得的，不使得边缘模糊，所以不需要多重分辨率方法。 与GVF的思路最接近的外部力模型是Cohen和Cohen12的距离保守力。 与GVF一样，采用图像的边缘映射并提供很大的捕捉域。但是，和 GVF不一样，把snake不能推进凹轮廓去。 可见，采用与GVF类似的非保守外部力，是活动轮廓模型的以后研究方向。 注意，本文的部分内容出现在会议论文18中。2. 背景1）参数snake模型（Parametric Snake Model） 传统的snake以曲线 xs=xs, ys, s0, 1 表示，这个模型在图像的空间域最小化如下的函数 E=0112xs2+xs2+Eext(xs)ds （1） 这里， 和 是分别控制snake张度和刚性， xs 和 xs 分别表示关于 s 的一阶、二阶偏导数。 外部能量函数 Eext 由图像导出，以使得它在边缘类似的特征有更小的值。 灰度图像 Ix, y 是连续位置变量 x, y 的函数，并且传统的外部能如下 Eext1x, y=-Ix, y2 （2） Eext2x, y=-Gx, y*Ix, y2 （3） 让活动轮廓推向阶梯边缘 1，这里，Gx, y 是带有标准偏差 的二维Gaussian函数， 是梯度算子。 如果图像只包括线图形（白色上的黑色）， 适当的外部力会包括下列的公式 10 Eext3x, y=Ix, y （4） Eext4x, y=Gx, y*Ix, y （5） 从这些定义，容易看出， 越大，轮廓也越模糊。但是这么大的 也会用于增长活动轮廓的捕捉范围。 使得函数E最小化的snake满足Euler方程式 xs-xs-Eext=0。 （6） 可以把它看成一种力平衡方程式 Fint+ Fext(p)=0 (7)这里，Fint=xs-xs ， Fext(p)=-Eext。 内部力 Fint 阻止拉伸和弯曲，而外部保守力（external potential force）Fext(p) 把snake推向理想的图像边缘。 为了得到公式（6）的解，引入除了s以外，还依赖时间t的函数 x ，也就是 xs, t，以让snake动态。下次，函数x的偏导数由公式（6）算得，如下： Xts, t=xs-xs-Eext。 （8） 当解 Xs, t 稳定的时候，Xts, t项消失，可以得到方程式（6）的解。把方程式离散化，并把它的离散系迭代解下去，可以得到方程式（8）的数值解（cf.1）。 一般来说，确定snake模型的时候， 大多采用下列的两个参数中的一个； 一个是乘以 Xt，以控制时空间步长大小（temporal step-size）的，另一个是乘以 Eext， 以允许单独控制外部力的强度。本文中，我们把外部力规则化，以使它的最大值为一，并且在所有的是实验，采用单位的时空间步长大小。2）传统snake的情况 图1给出了传统snake的例子。 图1(a)给出了U-型对象的6464像素的线图形，它上面有凹轮廓。 图片上还有表示传统snake（=0.6, =0.0）的迭代过程的曲线序列，它的初始化在对象的外面而在保守力场的捕捉区域内进行。 保守力场为 Fext(p)=-Eext4，=1.0 像素时的结果如图1(b)。 用snake公式的Euler方程式求得图1(a)的最后解，但凹区域仍保持分离。 Snake粗收敛的原因在图1(c)中说明了，图片上给出的是凹轮廓里面的外部力场。虽然外部力正确指向对象的边缘，但是在凹轮廓的内部，却指向水平的相反方向。 因此，活动轮廓向U-型对象的“手指”部分拉开，但是不能向凹里的部分前进。这里没有参数 和 能解决这个问题。 关于传统snake公式重要的另一个问题是捕捉区域的限制，通过图1(c), 很容易看出这点。 图中，我们能看到外部力的大小在边缘附近很快消失。 在方程式（5）中， 增大的时候，区域也增长，但是轮廓的细节越来越模糊和不清楚，最后对于太大的 ，凹部分本身最终被抹去。 Cohen和Cohen 12 提出了一种外部力模型，这个模型明显地增加了传统snake的捕捉区域。这些外部力是某种势函数的负梯度，这种势函数由Euclidean（或者chamfer）距离映射计算。 这些力叫距离保守力，以区别于第二节A中说明的传统保守力区分。 图2给出了基于保守力的snake工作性能。 图2(a)给出了U-型对象，同时给出了一系列的轮廓，显示 snake从离对象很远的初始到最后状态怎么变化。 图2(b)上，在离对象很远的地方，距离保守力的向量的大小仍然很大， 这说明了基于这种外部力模型的捕捉区域为什么那么大。 从图2(a)，容易看出，这种snake也仍然不能收敛到轮廓的凹部分。把距离保守力放大到图2(c)，仔细观察，可以说明这个原因。 可以看到与传统保守力一样，这些力仍然指向水平的相反方向，也就是说能把snake拉开，而不能拉到凹轮廓的里面。可见，Cohen和Cohen的修正方法是对距离映射进行非线性变换12， 只能变更力的大小，而不能变更它的方向。 因此，关于收敛到凹轮廓的问题，用距离保守力也不能解决。 3）一般化的平衡力方程式 图1(a)和2(a)上的两种snake解分别满足关于该能量模型的Euler方程式（6）。因此，收敛到对象函数(1)的局部最小会产生很粗的结果。 有学者直接用平衡力方程式定义snake，得到了这个问题的解， 这里把标准外部力 Fext(p) 换成更一般化的外部力Fext(g)： Fint+ Fext(g)=0。 （9） Fext(g) 的选择对snake的实行和行为有很大影响。 一般来说，外部力 Fext(g) 可以分成两个部分：静态和动态。静态力由图像数据计算，所以在snake过程中保持不变。 标准snake保守力是静态外部力。 动态力在snake过程中会发生变化。 本文提出了新形式的静态外部力， 它不随时间的变化，也不依赖snake本身的位置。Helmholtz理论主张把最一般的静态向量场分成两个分量：无旋（curl-free）分量和无源（divergence-free）分量。 静态无旋场由势函数的梯度表示，所以外部保守力在平衡力方程式以静态无旋的形式出现。 因此，允许包含无旋和无源这两个分量，就可以得到更一般的静态场。 本文提出了一种更合适地设计外部力场的方法，这个力场包含期望的两个特点，即更大的捕捉域，以及存在指向凹轮廓的力。 推到的公式产生外部力场， 这个力场会包含无旋和无源这两个分量。3. 梯度向量流Snake上述方法以snake设计的出发点，采用了平衡力方程式（7）。 下面，提出了新的静态外部力场 Fextg=vx, y，这个力场叫gradient vector flow（GVF）场。 由式（8）得到下面的方程式 Xts, t=xs-xs-v (10) 满足上面的方程式的参数曲线，叫GVF snake。 与传统snake一样，经过离散化和迭代，可以求数值解。1） 边缘映射首先定义由图像 Ix, y 导出的边缘映射 f(x, y) （edge map f(x, y)）， 它的特点是在边缘附近更大。 可以采用在图像处理文献提到的任何一种灰度映射或者二值边缘映射（cf.,20）；例如，可以采用 fx, y=-Eextix, y （11） 这里，i=1, 2, 3, 4。 本文中，边缘映射的三个特点非常重要。 第一，边缘映射的梯度 f 带有指向边缘的向量，在边缘上垂直于边缘。 第二，这些向量只在离边缘很近的地方有较大的大小。 第三，在各向同性的区域，就是说 Ix, y 近于常数的区域内，f 也近于零。 考虑以边缘映射的梯度作为外部力，对传统的snake，会产生什么样的影响？ 因为第一个特点，如果在边缘附近进行初始化，那条snake就在离边缘很近的地方稳定。这是想要的性质。 但是，一般因为第二个特点，捕捉区域很小。因为第三个特点，在各向同性的区域，不存在外部力。 后面的两个性质不是期望的。本文提出到的方法，保持了边缘附近的梯度特点，同时，通过计算扩散过程，将梯度映射从边缘拓展到各向同性的区域。一个很大的好处是这种扩散处理产生指向凹轮廓的向量。 2） 梯度向量流本文提出到的梯度向量流场是向量场 vx, y=ux, y, v(x, y)，它把能量函数最小化 =(ux2+uy2+vx2+vx2)+f2v-f2dxdy。 （12）特别，在 f 很小的时候，能量由向量场的偏导数的平方和决定，导致缓慢变化的场。 在 f 很大的时候，重积分中第两个项起决定作用，所以在 v=f 的时候， 最小。 这得到期望的结果，也就是说在边缘映射的梯度很大的时侯，让 v 等于它，而在各向同性的区域，让力慢慢地变化。 参数 是规则化参数，控制重积分中第一项和第二项之间的关系。 根据图像内存在的噪音量，可以确定这参数（噪音越大， 越大）。 第一项也代表关于向量旋度和散度的评价函数（equal pernalty）22。 所以，把式（12）最小化的向量场会是完全无旋无源的。 可以用下面的Euler方程式，求GVF场 2u-u-fxfx2+fy2=0 （13a） 2v-v-fyfx2+fy2=0 （13b） 这里，2 是Laplacian算子。在各向同性的区域 Ix, y 是常数，因为 fx, y 的梯度等于零，每个方程式的第二个项等于零。 所以在这样的区域内，u 和 v 分别由Laplace方程式确定。3） 数值计算为了求式（13a）和（13b）的解，把 u 和 v 看作时间的函数，就得到如下的方程式 utx,y,t=2ux,y,t-ux,y,t-fxx,y fxx,y2+fyx,y2 （14a） vtx,y,t=2vx,y,t-vx,y,t-fyx,y fx(x,y)2+fy(x,y)2 （14b） （14）的方程式称为一般化扩散方程式 24。 为了方便处理，把式（14）改写成如下： utx,y,t=2ux,y,t-bx,yux,y,t+c1x,y （15a） utx,y,t=2vx,y,t-bx,yvx,y,t+c2x,y （15b） 这里， bx,y=fx(x,y)2+fy(x,y)2 c1x,y=bx,yfxx,y c2x,y= bx,yfyx,y。为了计算 fx 和 fy， 可以采用任何数字图像梯度算子（cf.，20）。 在本文的例子中，采用了简单的中心差（difference）。 然后，计算系数 bx,y， c1x,y 和 c2x,y，它们在整个迭代过程中保持固定。 为了处理迭代解，用系数 i, j 和 n 分别代替 x, y 和 t。 像素之间的空间是 x 和 y，并且每个迭代的时间步长是 t。 则需要的篇导数近似如下 ut= 1t (ui,jn+1-ui,jn) vt= 1t (vi,jn+1-vi,jn) 2u=1xy(ui+1,j+ui,j+1+ui-1,j+ui,j-1-4ui,j) 2v=1xy(vi+1,j+vi,j+1+vi-1,j+vi,j-1-4vi,j)如果把这些近似式代入式（15），就得到GVF的迭代解： ui,jn+1=1-bi,jtui,jn+rui+1,jn+ ui,j+1n+ui-1,jn+ui,j-1n-4ui,jn +ci,j1t （16a） vi,jn+1=1-bi,jtvi,jn+rvi+1,jn+ vi,j+1n+vi-1,jn+vi,j-1n-4vi,jn +ci,j2t （16b）这里， r=txy （17） 如果 b，c1 和 c2带有一定的限制，而且Courant-Friedrichs-Lewy步长限制满足r1/4，解（16）是稳定的。 为了保证GVF的收敛，时间步长 t 要满足下列的限制： txy4 （18）从这个条件可以看到：首先，在粗图像，也就是说在 x 和 y 比较大的时候，收敛速度很快。其次， 越大，且GVF越平滑，收敛速度会越慢（因为 t 要小）。 本文中，用MATLAB代码，进行了2-D GVF的计算。 对SGI Indego-2 上 N=256 256 像素的图像，进行了实验。 对传统的势力，典型的运算时间为8秒（以C编程）。对距离势力，为155秒（Euclidean距离映射，以MATLAB编程）。对GVF力，为420秒（以MATLAB编程，用 N 迭代）。 用以C或者FORTRAN编程的优化代码，可以本质上减少运算时间。 例如，在以C实现3-D GVF的时候，对256 256 60 立体素（voxel）图像，以150迭代完成运算要31分的时间。 对2-D（）像素的图像，大概为53秒。以后，类似于多重网（multigrid）的方法可以更大地改善上述运算。4. GVF场和Snake：论证在本节中，对简单的对象，进行GVF场的运算，并且论证GVF snake的几个基本性质。 对所有的snake，采用了 =0.6 和 =0.0 ，并且对GVF，采用了 =0.2 。 所有在GVF运算中所用的边缘映射都归一化到 0, 1 的区域。流线（Streamlines