图像形状特征（五）--自由式变形模板

变形模板分自由式和参数式，这里先说自由式，典型的自由式变形模板，就是本文的snake模型。

原文转自：http://blog.csdn.net/zouxy09/article/details/8712287 

基于能量泛函的分割方法：

       该类方法主要指的是活动轮廓模型（active contour model）以及在其基础上发展出来的算法，其基本思想是使用连续曲线来表达目标边缘，并定义一个能量泛函使得其自变量包括边缘曲线，因此分割过程就转变为求解能量泛函的最小值的过程，一般可通过求解函数对应的欧拉(Euler．Lagrange)方程来实现，能量达到最小时的曲线位置就是目标的轮廓所在。

       主动轮廓线模型是一个自顶向下定位图像特征的机制，用户或其他自动处理过程通过事先在感兴趣目标附近放置一个初始轮廓线，在内部能量（内力）和外部能量（外力）的作用下变形外部能量吸引活动轮廓朝物体边缘运动，而内部能量保持活动轮廓的光滑性和拓扑性，当能量达到最小时，活动轮廓收敛到所要检测的物体边缘。

一、曲线演化理论

       曲线演化理论在水平集中运用到，但我感觉在主动轮廓线模型的分割方法中，这个知识是公用的，所以这里我们简单了解下。

       曲线可以简单的分为几种：

 

       曲线存在曲率，曲率有正有负，于是在法向曲率力的推动下，曲线的运动方向之间有所不同：有些部分朝外扩展，而有些部分则朝内运动。这种情形如下图所示。图中蓝色箭头处的曲率为负，而绿色箭头处的曲率为正。

       简单曲线在曲率力（也就是曲线的二次导数）的驱动下演化所具有的一种非常特殊的数学性质是：一切简单曲线，无论被扭曲得多么严重，只要还是一种简单曲线，那么在曲率力的推动下最终将退化成一个圆，然后消逝（可以想象下，圆的所有点的曲率力都向着圆心，所以它将慢慢缩小，以致最后消逝）。

       描述曲线几何特征的两个重要参数是单位法矢和曲率，单位法矢描述曲线的方向，曲率则表述曲线弯曲的程度。曲线演化理论就是仅利用曲线的单位法矢和曲率等几何参数来研究曲线随时间的变形。曲线的演变过程可以认为是表示曲线在作用力 F的驱动下，朝法线方向 N以速度 v演化。而速度是有正负之分的，所以就有如果速度 v的符号为负，表示活动轮廓演化过程是朝外部方向的，如为正，则表示朝内部方向演化，活动曲线是单方向演化的，不可能同时往两个方向演化。

       所以曲线的演变过程，就是不同力在曲线上的作用过程，力也可以表达为能量。世界万物都趋向于能量最小而存在。因为此时它是最平衡的，消耗最小的（不知理解对不？）。那么在图像分割里面，我们目标是把目标的轮廓找到，那么在目标的轮廓这个地方，整个轮廓的能量是最小的，那么曲线在图像任何一个地方，都可以因为力朝着这个能量最小的轮廓演变，当演变到目标的轮廓的时候，因为能量最小，力平衡了，速度为0了，也就不动了，这时候目标就被我们分割出来了。

        那现在关键就在于：1）这个轮廓我们怎么表示；2）这些力怎么构造，构造哪些力才可以让目标轮廓这个地方的能量最小？

       这两个问题的描述和解决就衍生出了很多的基于主动轮廓线模型的分割方法。第一个问题的回答，就形成了两大流派：如果这个轮廓是参数表示的，那么就是参数活动轮廓模型（parametric active contour model），典型为snake模型，如果这个轮廓是几何表示的，那么就是几何活动轮廓模型（geometric active contour model），即水平集方法（Level Set），它是把二维的轮廓嵌入到三维的曲面的零水平面来表达的（可以理解为一座山峰的等高线，某个等高线把山峰切了，这个高度山峰的水平形状就出来了，也就是轮廓了），所以低维的演化曲线或曲面，表达为高维函数曲面的零水平集的间接表达形式（这个轮廓的变化，直观上我们就可以调整山峰的形状或者调整登高线的高度来得到）。

       那对于第二个问题，是两大流派都遇到的问题，是他们都需要解决的最关键的问题。哪些力才可以达到分割的目标呢？这将在后面聊到。

二、Snakes模型

       自1987年Kass提出Snakes模型以来，各种基于主动轮廓线的图像分割理解和识别方法如雨后春笋般蓬勃发展起来。Snakes模型的基本思想很简单，它以构成一定形状的一些控制点为模板（轮廓线），通过模板自身的弹性形变，与图像局部特征相匹配达到调和，即某种能量函数极小化，完成对图像的分割。再通过对模板的进一步分析而实现图像的理解和识别。

        简单的来讲，SNAKE模型就是一条可变形的参数曲线及相应的能量函数，以最小化能量目标函数为目标，控制参数曲线变形，具有最小能量的闭合曲线就是目标轮廓。

       构造Snakes模型的目的是为了调和上层知识和底层图像特征这一对矛盾。无论是亮度、梯度、角点、纹理还是光流，所有的图像特征都是局部的。所谓局部性就是指图像上某一点的特征只取决于这一点所在的邻域，而与物体的形状无关。但是人们对物体的认识主要是来自于其外形轮廓。如何将两者有效地融合在一起正是Snakes模型的长处。Snakes模型的轮廓线承载了上层知识，而轮廓线与图像的匹配又融合了底层特征。这两项分别表示为Snakes模型中能量函数的内部力和图像力。

       模型的形变受到同时作用在模型上的许多不同的力所控制，每一种力所产生一部分能量，这部分能量表示为活动轮廓模型的能量函数的一个独立的能量项。

         Snake模型首先需要在感兴趣区域的附近给出一条初始曲线，接下来最小化能量泛函，让曲线在图像中发生变形并不断逼近目标轮廓。

        Kass等提出的原始Snakes模型由一组控制点：v(s)=[x(s), y(s)]   s∈[0, 1]组成，这些点首尾以直线相连构成轮廓线。其中x(s)和y(s)分别表示每个控制点在图像中的坐标位置。 s 是以傅立叶变换形式描述边界的自变量。在Snakes的控制点上定义能量函数（反映能量与轮廓之间的关系）：

      其中第1项称为弹性能量是v的一阶导数的模，第2项称为弯曲能量，是v的二阶导数的模，第3项是外部能量（外部力），在基本Snakes模型中一般只取控制点或连线所在位置的图像局部特征例如梯度：

也称图像力。（当轮廓C靠近目标图像边缘，那么C的灰度的梯度将会增大，那么上式的能量最小，由曲线演变公式知道该点的速度将变为0，也就是停止运动了。这样，C就停在图像的边缘位置了，也就完成了分割。那么这个的前提就是目标在图像中的边缘比较明显了，否则很容易就越过边缘了。）

        弹性能量和弯曲能量合称内部能量（内部力），用于控制轮廓线的弹性形变，起到保持轮廓连续性和平滑性的作用。而第三项代表外部能量，也被称为图像能量，表示变形曲线与图像局部特征吻合的情况。内部能量仅仅跟snake的形状有关，而跟图像数据无关。而外部能量仅仅跟图像数据有关。在某一点的α和β的值决定曲线可以在这一点伸展和弯曲的程度。

       最终对图像的分割转化为求解能量函数Etotal(v)极小化（最小化轮廓的能量）。在能量函数极小化过程中，弹性能量迅速把轮廓线压缩成一个光滑的圆，弯曲能量驱使轮廓线成为光滑曲线或直线，而图像力则使轮廓线向图像的高梯度位置靠拢。基本Snakes模型就是在这3个力的联合作用下工作的。

        因为图像上的点都是离散的，所以我们用来优化能量函数的算法都必须在离散域里定义。所以求解能量函数Etotal(v)极小化是一个典型的变分问题（微分运算中，自变量一般是坐标等变量，因变量是函数；变分运算中，自变量是函数，因变量是函数的函数，即数学上所谓的泛函。对泛函求极值的问题，数学上称之为变分法）。

        在离散化条件（数字图像）下，由欧拉方程可知最终问题的答案等价于求解一组差分方程：（欧拉方程是泛函极值条件的微分表达式，求解泛函的欧拉方程，即可得到使泛函取极值的驻函数，将变分问题转化为微分问题。）

       记外部力 F = −∇ P， Kass等将上式离散化后，对x(s)和y(s)分别构造两个五对角阵的线性方程组，通过迭代计算进行求解。在实际应用中一般先在物体周围手动点出控制点作为Snakes模型的起始位置，然后对能量函数迭代求解。

 

总结：

       主动轮廓线模型又称为SNAKE模型，自Kass于1987年提出以来，已广泛应用于数字图像分析和计算机视觉领域。由于SNAKE模型具有良好的提取和跟踪特定区域内目标轮廓的能力，因此非常适合于医学图像如CT和MR图像的处理，以获得特定器官及组织的轮廓。简单的来讲，SNAKE模型就是一条可变形的参数曲线及相应的能量函数，以最小化能量目标函数为目标，控制参数曲线变形，具有最小能量的闭合曲线就是目标轮廓。SNAKE模型具有一些经典方法所无法比拟的优点：图像数据、初始估计、目标轮廓及基于知识的约束统一于一个过程中；经适当的初始化后，它能自主地收敛于能量极小值状态；尺度空间中由初到精地极小化能量可以极大地扩展捕获区域和降低复杂性。同时，SNAKE模型也有其自身的缺点：对初始位置敏感，需要依赖其他机制将SNAKE放置在感兴趣的图像特征附近；由于SNAKE模型的非凸性，它有可能收敛到局部极值点，甚至发散。

Snake模型发展10多年来，许多学者对于经典的snake模型做了改进，提出各种改进的snake模型，其中梯度矢量流（Gradient Vector  Flow,GVF）模型扩大了经典snake的外力作用范围，加强了对目标凹轮廓边缘的吸引力，提高了传统的snake模型。

 

Snake模型主要研究的方面：

1.表示内部能量的曲线演化    2.外力    3.能量最小化

Snake模型初始轮廓的选择

由于snake模型对于初始位置比较敏感，因此要求初始轮廓尽可能的靠近真实轮廓，而当图像边缘模糊，目标比较复杂或与其他的物体靠的比较近时，其初始轮廓更不易确定。

现有的初始轮廓确定的方法有以下几种：1.人工勾勒图像的边缘    2.序列图像差分边界    3.基于序列图像的前一帧图像边界的预测  4.基于传统图像分割结果进行边界选取

分水岭算法

分水岭算法是由S.Beucher  F.Meyer最早引入图像分割领域，它的基本思想是来源于测地学上的侧线重构，其内容是把图像看做是测地学上的拓扑地貌。进行分水岭模型计算的比较经典的算法是L  Vincent提出的，在该算法中首先是对每个像素的灰度级进行从低到高排序，然后用等级对垒模拟淹没，初始时，等级队列中为淹没的初始点，在从低到高实现淹没的过程中，对每一个局部极小值在H阶高度的影响域采用先进先出（FIFO）结构进行判断及标注，直到最后一个值被淹没，从而正确划分各个区域。

整个洪水淹没的循环迭代过程可以通过以下两个步骤表示：

分水岭算法的优点：

1.分水岭算法对于图像中由于像素差别较小而产生微弱边缘具有良好的响应，可以得到封闭连续的边缘，而且可以保证在照明，阴影等影响下分割边缘的封闭性和连续性

分水岭算法对于目标物体之间或者是目标物体同背景物体之间粘连的情况有较好的处理效果。能够较好的分割这类目标物体。

3.图像内部的阴暗变化对于分水岭算法影响较小，可以在一定程度上减小由于阴暗便哈带来的图像分割影响

与其他边缘分割算子比较：

Canny算子可以很好的勾勒出物体的轮廓，过分的强调轮廓的特性，而没有强调物体的轮廓必须是封闭的，在图像中显示的轮廓是不封闭的，物体内部阴暗变化也被当做边界检测出来，形成大量的伪边缘。

分水岭算法分割得到的轮廓曲线时连续封闭的，图像内部的阴暗变化没有生成独立的轮廓线。

Snake模型的缺陷：

对初始位置敏感，易陷入局部极值，无法收敛到轮廓深度凹陷部分，不具备自动拓扑变换功能等。

Snake模型的改进算法：

1.Cohen提出的气球（balloon）理论模型:应用压力和高斯能力一起增大吸引范围的方法，该压力可使模型扩大或缩小，因此不再要求将模型初始化在所期望的对象边界附近。在图像的梯度力场上叠加气球里，以使轮廓线作为一个整体进行膨胀或收缩，从而扩大了模型寻找图像特证的范围。

优势：对初始边界不敏感            

存在的缺点：存在弱边界，漏出边界间隙等问题。

2.Xu提出梯度矢量流（GVF）概念，用GVF场代替经典外力场，GVF场可以看做是对图像梯度场得逼近，这不仅使模型捕捉的范围得到了提高，而且能使活动轮廓进入凹陷区。

优势：有良好的收敛性，深入目标边缘的凹陷区域           

存在的缺点：仍不能解决曲线的拓扑变化问题

局部优化算法：

1.Amini提出基于动态规划的snake算法。 2.变分法  3.贪婪算法  4.有限差分法   5.有限元法

全局优化算法：

1.模拟退火     2.遗传算法    3.神经网络

Snake模型的蚁群算法（Ant Colony Optimization）模型

蚁群算法是最近几年有意大利学者M.Dorigo等人首次提出的一种新型的模拟进化算法，称为蚁群系统，蚁群算法通过候选解组成的群体的进化过程来寻求最优解，该过程包括两个基本阶段：适应阶段和协同工作阶段，算法本身采用正反馈原理，加快了进化过程，不易陷入局部最优解，而且个体之间不断进行信息交流和传递，有利于对解空间的进一步探索，因此有很强的发展解的能力。。

 

 

Snake的进化模型

1.McInerney 提出一种拓扑自适应snake模型（Topology Adaptive  Snake,T-Snake）

该算法基于仿射细胞图像分解（Affine Cell Image  Decomposition，ACID）先在待分割图像上加以三角网格，然后在图像区域的适当位置做一条初始曲线，最后取曲线与网格的交点作为snake的初始离散点，其第i个snake的离散点的坐标为其中，相邻两点，之间由一条弹性样条连接而成

由于T-Snake模型可借助三角形网格和网格点的特征函数来确定边界三角形，可促使snake模型演化过程中的分裂和合并，从而保证了其具有能够处理拓扑结果复杂图像的能力，因此能够很好的满足医学图像拓扑结果复杂的特点。此算法用于脑部MR切片有良好的性能。

2.双T-Snake模型

双T-Snake模型（Dual-T-Snakes）是在T-Snake模型的基础上产生的，其主要思想是采用内外两个初始轮廓，其中一个轮廓从目标外向内收缩和分裂，另一个轮廓从目标内部向外膨胀，两个初始轮廓可以离目标边界较远，迭代的过程中对能量较大的轮廓增加驱动力，使其靠近与之相对应的轮廓，直到连个轮廓收敛到同一个为止

3.Loop  Snake 模型

Loop  Snake模型是一种加强了拓扑控制的T-Snake模型，这种方法的关键集中在曲线的每一步进化中都要形成循环，其基本思想是，确保图像轮廓曲线精确地线性地映射到适当的分类中，然后在额外的记过loop-Tree的帮助下，尽可能少的时间内运用已经被snake探究的循环来决定是否进行区域划分，这种模型的实质是对T-Snake模型的一种改进。由于加强了拓扑控制，使得Loop Snake模型既可以忽略背景中强噪声又可以在演化过程中进行多次分裂。

4.连续snake模型

在Snake模型中，轮廓曲线由一条给定容许误差范围的光滑曲线组成，相对于离散snake来说，连续snake模型所需要的控制点少，比离散的更具优越性。

5.B-Snake模型

B-Snake模型是通过B样条曲线来定义的，其轮廓曲线由各曲线段光滑相连而成，每一个曲线段都是由一个给定次数多项式表示，这种多项式是B样条曲度函数的一种线性组合，并以控制点为系数。在有些B-Snake模型中并没有明确应用内部能量，这是因为B样条本身就含有内部能量，snake轮廓曲线只受外力影响着图像边缘移动。可用于对图像切片分割区域的描述与跟踪而用于器官的三维重建。

 

应用snake的优势：由于生物或人体组织解剖结构的复杂性，以及软组织形状的易变性，那些仅依赖于图像本身的灰度，纹理属性等低层次视觉属性来进行分割的图像分割方法难以获得理想的分割效果，因此医学图像分割迫切需要有一种灵活的框架，能将基于图像本身低层次视觉属性（边缘，纹理，灰度，色彩）和人们对于待分割目标的知识经验，如目标形状的描述，亮度，色彩的经验统计，医生的经验等，可以一种有机的方式整合起来，得到待分割区域的完整表达。

 

Open源码分析：

 

/*M///
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                        Intel License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of Intel Corporation may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "_cv.h"
 
#define _CV_SNAKE_BIG 2.e+38f
#define _CV_SNAKE_IMAGE 1
#define _CV_SNAKE_GRAD  2
 
 
/*F///
//    Name:      icvSnake8uC1R    
//    Purpose:  
//    Context:  
//    Parameters:
//               src - source image,
//               srcStep - its step in bytes,
//               roi - size of ROI,
//               pt - pointer to snake points array
//               n - size of points array,
//               alpha - pointer to coefficient of continuity energy,
//               beta - pointer to coefficient of curvature energy, 
//               gamma - pointer to coefficient of image energy, 
//               coeffUsage - if CV_VALUE - alpha, beta, gamma point to single value
//                            if CV_MATAY - point to arrays
//               criteria - termination criteria.
//               scheme - image energy scheme
//                         if _CV_SNAKE_IMAGE - image intensity is energy
//                         if _CV_SNAKE_GRAD  - magnitude of gradient is energy
//    Returns:  
//F*/
 
static CvStatus
icvSnake8uC1R( unsigned char *src,   //原始图像数据
               int srcStep,         //每行的字节数
               CvSize roi,         //图像尺寸
               CvPoint * pt,       //轮廓点(变形对象)
               int n,            //轮廓点的个数
               float *alpha,       //指向α的指针，α可以是单个值，也可以是与轮廓点个数一致的数组
               float *beta,        //β的值，同α
               float *gamma,       //γ的值，同α
               int coeffUsage,   //确定αβγ是用作单个值还是个数组
        CvSize win,       //每个点用于搜索的最小的领域大小,宽度为奇数
             CvTermCriteria criteria,   //递归迭代终止的条件准则
int scheme )         //确定图像能量场的数据选择，1为灰度，2为灰度梯度
{
    int i, j, k;
    int neighbors = win.height * win.width;    //当前点领域中点的个数
 
   //当前点的位置
    int centerx = win.width >> 1;          
    int centery = win.height >> 1;         
 
    float invn;        //n 的倒数？
    int iteration = 0;     //迭代次数
    int converged = 0;      //收敛标志，0为非收敛
    
  //能量
    float *Econt;    //
    float *Ecurv;   //轮廓曲线能量
    float *Eimg;    //图像能量
    float *E;      //
  
   //αβγ的副本
    float _alpha, _beta, _gamma;
 
    /*#ifdef GRAD_SNAKE */
    float *gradient = NULL;
    uchar *map = NULL;
    int map_width = ((roi.width - 1) >> 3) + 1;
    int map_height = ((roi.height - 1) >> 3) + 1;
    CvSepFilter pX, pY;
    #define WTILE_SIZE 8
    #define TILE_SIZE (WTILE_SIZE + 2)       
    short dx[TILE_SIZE*TILE_SIZE], dy[TILE_SIZE*TILE_SIZE];
    CvMat _dx = cvMat( TILE_SIZE, TILE_SIZE, CV_16SC1, dx );
    CvMat _dy = cvMat( TILE_SIZE, TILE_SIZE, CV_16SC1, dy );
    CvMat _src = cvMat( roi.height, roi.width, CV_8UC1, src );
 
    /* inner buffer of convolution process */
    //char ConvBuffer[400];
 
    /*#endif */
 
     //检点参数的合理性
    /* check bad arguments */
    if( src == NULL )
        return CV_NULLPTR_ERR;
    if( (roi.height <= 0) || (roi.width <= 0) )
        return CV_BADSIZE_ERR;
    if( srcStep < roi.width )
        return CV_BADSIZE_ERR;
    if( pt == NULL )
        return CV_NULLPTR_ERR;
    if( n < 3 )                         //轮廓点至少要三个
        return CV_BADSIZE_ERR;
    if( alpha == NULL )
        return CV_NULLPTR_ERR;
    if( beta == NULL )
        return CV_NULLPTR_ERR;
    if( gamma == NULL )
        return CV_NULLPTR_ERR;
    if( coeffUsage != CV_VALUE && coeffUsage != CV_ARRAY )
        return CV_BADFLAG_ERR;
    if( (win.height <= 0) || (!(win.height & 1)))   //邻域搜索窗口得是奇数
        return CV_BADSIZE_ERR;
    if( (win.width <= 0) || (!(win.width & 1)))
        return CV_BADSIZE_ERR;
 
    invn = 1 / ((float) n);        //轮廓点数n的倒数，用于求平均？
 
    if( scheme == _CV_SNAKE_GRAD )
{
     //X方向上和Y方向上的Scoble梯度算子，用于求图像的梯度,
//处理的图像最大尺寸为TILE_SIZE+2，此例为12，算子半长为3即{-3，-2，-1，0，1，2，3}
//处理后的数据类型为16位符号数,分别存放在_dx,_dy矩阵中,长度为10
        pX.init_deriv( TILE_SIZE+2, CV_8UC1, CV_16SC1, 1, 0, 3 );
        pY.init_deriv( TILE_SIZE+2, CV_8UC1, CV_16SC1, 0, 1, 3 );
       //图像梯度存放缓冲区
        gradient = (float *) cvAlloc( roi.height * roi.width * sizeof( float ));
 
        if( !gradient )
            return CV_OUTOFMEM_ERR;
       //map用于标志相应位置的分块的图像能量是否已经求得
        map = (uchar *) cvAlloc( map_width * map_height );
        if( !map )
        {
            cvFree( &gradient );
            return CV_OUTOFMEM_ERR;
        }
        /* clear map - no gradient computed */
       //清除map标志
        memset( (void *) map, 0, map_width * map_height );
}
//各种能量的存放处，取每点的邻域的能量
    Econt = (float *) cvAlloc( neighbors * sizeof( float ));
    Ecurv = (float *) cvAlloc( neighbors * sizeof( float ));
    Eimg = (float *) cvAlloc( neighbors * sizeof( float ));
    E = (float *) cvAlloc( neighbors * sizeof( float ));
   //开始迭代
    while( !converged )    //收敛标志无效时进行
    {
        float ave_d = 0;  //轮廓各点的平均距离
        int moved = 0;      //轮廓变形时，发生移动的数量
 
        converged = 0;       //标志未收敛
        iteration++;        //更新迭代次数+1
 
//计算轮廓中各点的平均距离
        /* compute average distance */
      //从点0到点n-1的距离和
        for( i = 1; i < n; i++ )
        {
            int diffx = pt[i - 1].x - pt[i].x;
            int diffy = pt[i - 1].y - pt[i].y;
 
            ave_d += cvSqrt( (float) (diffx * diffx + diffy * diffy) ); 
        }
     //再加上从点n-1到点0的距离,形成回路轮廓
        ave_d += cvSqrt( (float) ((pt[0].x - pt[n - 1].x) *
                                  (pt[0].x - pt[n - 1].x) +
                                  (pt[0].y - pt[n - 1].y) * (pt[0].y - pt[n - 1].y)));
    //求平均，得出平均距离
        ave_d *= invn;
        /* average distance computed */
 
 
      //对于每个轮廓点进行特定循环迭代求解
        for( i = 0; i < n; i++ )
        {
            /* Calculate Econt */
          //初始化各个能量
            float maxEcont = 0;
            float maxEcurv = 0;
            float maxEimg = 0;
            float minEcont = _CV_SNAKE_BIG;
            float minEcurv = _CV_SNAKE_BIG;
            float minEimg = _CV_SNAKE_BIG;
            float Emin = _CV_SNAKE_BIG;
         //初始化变形后轮廓点的偏移量
            int offsetx = 0;
            int offsety = 0;
            float tmp;
 
        //计算边界
            /* compute bounds */
           //计算合理的搜索边界，以防领域搜索超过ROI图像的范围
            int left = MIN( pt[i].x, win.width >> 1 );
            int right = MIN( roi.width - 1 - pt[i].x, win.width >> 1 );
            int upper = MIN( pt[i].y, win.height >> 1 );
            int bottom = MIN( roi.height - 1 - pt[i].y, win.height >> 1 );
          //初始化Econt
            maxEcont = 0;
            minEcont = _CV_SNAKE_BIG;
         //在合理的搜索范围内进行Econt的计算
            for( j = -upper; j <= bottom; j++ )
            {
                for( k = -left; k <= right; k++ )
                {
                    int diffx, diffy;
                    float energy;
             //在轮廓点集的首尾相接处作相应处理,求轮廓点差分
                    if( i == 0 )
                    {
                        diffx = pt[n - 1].x - (pt[i].x + k);
                        diffy = pt[n - 1].y - (pt[i].y + j);
                    }
                    else
             //在其他地方作一般处理
 
                    {
                        diffx = pt[i - 1].x - (pt[i].x + k);
                        diffy = pt[i - 1].y - (pt[i].y + j);
                    }
             //将邻域陈列坐标转成Econt数组的下标序号，计算邻域中每点的Econt
              //Econt的值等于平均距离和此点和上一点的距离的差的绝对值（这是怎么来的？）
                    Econt[(j + centery) * win.width + k + centerx] = energy =
                        (float) fabs( ave_d -
                                      cvSqrt( (float) (diffx * diffx + diffy * diffy) ));
             //求出所有邻域点中的Econt的最大值和最小值
                    maxEcont = MAX( maxEcont, energy );
                    minEcont = MIN( minEcont, energy );
                }
            }
           //求出邻域点中最大值和最小值之差,并对所有的邻域点的Econt进行标准归一化,若最大值最小
           //相等，则邻域中的点Econt全相等，Econt归一化束缚为0
            tmp = maxEcont - minEcont;
            tmp = (tmp == 0) ? 0 : (1 / tmp);
            for( k = 0; k < neighbors; k++ )
            {
                Econt[k] = (Econt[k] - minEcont) * tmp;
            }
 
 
           //计算每点的Ecurv
            /*  Calculate Ecurv */
            maxEcurv = 0;
            minEcurv = _CV_SNAKE_BIG;
            for( j = -upper; j <= bottom; j++ )
            {
                for( k = -left; k <= right; k++ )
                {
                    int tx, ty;
                    float energy;
                    //第一个点的二阶差分
                    if( i == 0 )
                    {
                        tx = pt[n - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;
                        ty = pt[n - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;
                    }
                   //最后一个点的二阶差分
                    else if( i == n - 1 )
                    {
                        tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[0].x;
                        ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[0].y;
                    }
                   //其余点的二阶差分
                    else
                    {
                        tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;
                        ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;
                    }
                  //转换坐标为数组序号，并求各点的Ecurv的值,二阶差分后取平方
                    Ecurv[(j + centery) * win.width + k + centerx] = energy =
                        (float) (tx * tx + ty * ty);
                  //取最小的Ecurv和最大的Ecurv
                    maxEcurv = MAX( maxEcurv, energy );
                    minEcurv = MIN( minEcurv, energy );
                }
            }
               //对Ecurv进行标准归一化
            tmp = maxEcurv - minEcurv;
            tmp = (tmp == 0) ? 0 : (1 / tmp);
            for( k = 0; k < neighbors; k++ )
            {
                Ecurv[k] = (Ecurv[k] - minEcurv) * tmp;
            }
         
           //求Eimg
            /* Calculate Eimg */
            for( j = -upper; j <= bottom; j++ )
            {
                for( k = -left; k <= right; k++ )
                {
                    float energy;
               //若采用灰度梯度数据
                    if( scheme == _CV_SNAKE_GRAD )
                    {
                        /* look at map and check status */
                        int x = (pt[i].x + k)/WTILE_SIZE;
                        int y = (pt[i].y + j)/WTILE_SIZE;
                        //若此处的图像能量还没有获取，则对此处对应的图像分块进行图像能量的求解
                        if( map[y * map_width + x] == 0 )
                        {
                            int l, m;                           
 
                            /* evaluate block location */
                           //计算要进行梯度算子处理的图像块的位置
                            int upshift = y ? 1 : 0;
                            int leftshift = x ? 1 : 0;
                            int bottomshift = MIN( 1, roi.height - (y + 1)*WTILE_SIZE );
                            int rightshift = MIN( 1, roi.width - (x + 1)*WTILE_SIZE );
                          //图像块的位置大小（由于原ROI不一定是8的倍数，所以图像块会大小不一）
                            CvRect g_roi = { x*WTILE_SIZE - leftshift, y*WTILE_SIZE - upshift,
                                leftshift + WTILE_SIZE + rightshift, upshift + WTILE_SIZE + bottomshift };
                            CvMat _src1;
                            cvGetSubArr( &_src, &_src1, g_roi );  //得到图像块的数据
                            //分别对图像的X方向和Y方向进行梯度算子
                            pX.process( &_src1, &_dx );
                            pY.process( &_src1, &_dy );
                         //求分块区域中的每个点的梯度
                            for( l = 0; l < WTILE_SIZE + bottomshift; l++ )
                            {
                                for( m = 0; m < WTILE_SIZE + rightshift; m++ )
                                {
                                    gradient[(y*WTILE_SIZE + l) * roi.width + x*WTILE_SIZE + m] =
                                        (float) (dx[(l + upshift) * TILE_SIZE + m + leftshift] *
                                                 dx[(l + upshift) * TILE_SIZE + m + leftshift] +
                                                 dy[(l + upshift) * TILE_SIZE + m + leftshift] *
                                                 dy[(l + upshift) * TILE_SIZE + m + leftshift]);
                                }
                            }
                            //map相应位置置1表示此处图像能量已经获取
                            map[y * map_width + x] = 1;
                        }
                      //以梯度数据作为图像能量
                        Eimg[(j + centery) * win.width + k + centerx] = energy =
                            gradient[(pt[i].y + j) * roi.width + pt[i].x + k];
                    }
                    else
                    {
                       //以灰度作为图像能量
                        Eimg[(j + centery) * win.width + k + centerx] = energy =
                            src[(pt[i].y + j) * srcStep + pt[i].x + k];
                    }
                   //获得邻域中最大和最小的图像能量
                    maxEimg = MAX( maxEimg, energy );
                    minEimg = MIN( minEimg, energy );
                }
            }
              //Eimg的标准归一化
            tmp = (maxEimg - minEimg);
            tmp = (tmp == 0) ? 0 : (1 / tmp);
 
            for( k = 0; k < neighbors; k++ )
            {
                Eimg[k] = (minEimg - Eimg[k]) * tmp;
            }
            //加入系数
            /* locate coefficients */
            if( coeffUsage == CV_VALUE)
            {
                _alpha = *alpha;
                _beta = *beta;
                _gamma = *gamma;
            }
            else
            {                  
                _alpha = alpha[i];
                _beta = beta[i];
                _gamma = gamma[i];
            }
 
            /* Find Minimize point in the neighbors */
            //求得每个邻域点的Snake能量
            for( k = 0; k < neighbors; k++ )
            {
                E[k] = _alpha * Econt[k] + _beta * Ecurv[k] + _gamma * Eimg[k];
            }
            Emin = _CV_SNAKE_BIG;
        //获取最小的能量，以及对应的邻域中的相对位置
            for( j = -upper; j <= bottom; j++ )
            {
                for( k = -left; k <= right; k++ )
                {
 
                    if( E[(j + centery) * win.width + k + centerx] < Emin )
                    {
                        Emin = E[(j + centery) * win.width + k + centerx];
                        offsetx = k;
                        offsety = j;
                    }
                }
            }
         //如果轮廓点发生改变，则记得移动次数
            if( offsetx || offsety )
            {
                pt[i].x += offsetx;
                pt[i].y += offsety;
                moved++;
            }
        }
 
      //各个轮廓点迭代计算完成后，如果没有移动的点了，则收敛标志位有效，停止迭代
        converged = (moved == 0);
     //达到最大迭代次数时，收敛标志位有效，停止迭代
        if( (criteria.type & CV_TERMCRIT_ITER) && (iteration >= criteria.max_iter) )
            converged = 1;
  //到大相应精度时，停止迭代（与第一个条件有相同效果）
        if( (criteria.type & CV_TERMCRIT_EPS) && (moved <= criteria.epsilon) )
            converged = 1;
    }
 
  //释放各个缓冲区
    cvFree( &Econt );
    cvFree( &Ecurv );
    cvFree( &Eimg );
    cvFree( &E );
 
    if( scheme == _CV_SNAKE_GRAD )
    {
        cvFree( &gradient );
        cvFree( &map );
    }
    return CV_OK;
}
 
 
CV_IMPL void
cvSnakeImage( const IplImage* src, CvPoint* points,
              int length, float *alpha,
              float *beta, float *gamma,
              int coeffUsage, CvSize win,
              CvTermCriteria criteria, int calcGradient )
{
 
    CV_FUNCNAME( "cvSnakeImage" );
 
    __BEGIN__;
 
    uchar *data;
    CvSize size;
    int step;
 
    if( src->nChannels != 1 )
        CV_ERROR( CV_BadNumChannels, "input image has more than one channel" );
 
    if( src->depth != IPL_DEPTH_8U )
        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
 
    cvGetRawData( src, &data, &step, &size );
 
    IPPI_CALL( icvSnake8uC1R( data, step, size, points, length,
                              alpha, beta, gamma, coeffUsage, win, criteria,
                              calcGradient ? _CV_SNAKE_GRAD : _CV_SNAKE_IMAGE ));
    __END__;
}
 
/* end of file */
 
 
 
 
 
测试应用程序
 
#include "stdafx.h"
#include <iostream>
#include <string.h>
#include <cxcore.h>
#include <cv.h>
#include <highgui.h>
#include <fstream>
 
 
IplImage *image = 0 ; //原始图像
IplImage *image2 = 0 ; //原始图像copy
 
using namespace std;
int Thresholdness = 141;
int ialpha = 20;
int ibeta=20;
int igamma=20;
 
void onChange(int pos)
{
   
    if(image2) cvReleaseImage(&image2);
    if(image) cvReleaseImage(&image);
 
    image2 = cvLoadImage("grey.bmp",1); //显示图片
    image= cvLoadImage("grey.bmp",0);
 
    cvThreshold(image,image,Thresholdness,255,CV_THRESH_BINARY); //分割域值   
 
    CvMemStorage* storage = cvCreateMemStorage(0);
    CvSeq* contours = 0;
 
    cvFindContours( image, storage, &contours, sizeof(CvContour), //寻找初始化轮廓
        CV_RETR_EXTERNAL , CV_CHAIN_APPROX_SIMPLE );
 
    if(!contours) return ;
    int length = contours->total;   
    if(length<10) return ;
    CvPoint* point = new CvPoint[length]; //分配轮廓点
 
    CvSeqReader reader;
    CvPoint pt= cvPoint(0,0);;   
    CvSeq *contour2=contours;   
 
    cvStartReadSeq(contour2, &reader);
    for (int i = 0; i < length; i++)
    {
        CV_READ_SEQ_ELEM(pt, reader);
        point[i]=pt;
    }
    cvReleaseMemStorage(&storage);
 
    //显示轮廓曲线
    for(int i=0;i<length;i++)
    {
        int j = (i+1)%length;
        cvLine( image2, point[i],point[j],CV_RGB( 0, 0, 255 ),1,8,0 );
    }
 
    float alpha=ialpha/100.0f;
    float beta=ibeta/100.0f;
    float gamma=igamma/100.0f;
 
    CvSize size;
    size.width=3;
    size.height=3;
    CvTermCriteria criteria;
    criteria.type=CV_TERMCRIT_ITER;
    criteria.max_iter=1000;
    criteria.epsilon=0.1;
    cvSnakeImage( image, point,length,&alpha,&beta,&gamma,CV_VALUE,size,criteria,0 );
 
    //显示曲线
    for(int i=0;i<length;i++)
    {
        int j = (i+1)%length;
        cvLine( image2, point[i],point[j],CV_RGB( 0, 255, 0 ),1,8,0 );
    }
    delete []point;
 
}
 
int main(int argc, char* argv[])
{
 
   
    cvNamedWindow("win1",0);
    cvCreateTrackbar("Thd", "win1", &Thresholdness, 255, onChange);
    cvCreateTrackbar("alpha", "win1", &ialpha, 100, onChange);
    cvCreateTrackbar("beta", "win1", &ibeta, 100, onChange);
    cvCreateTrackbar("gamma", "win1", &igamma, 100, onChange);
    cvResizeWindow("win1",300,500);
    onChange(0);
 
    for(;;)
    {
        if(cvWaitKey(40)==27) break;
        cvShowImage("win1",image2);
    }
   
    return 0;
}

 

 

 

 

Reference：

http://blog.csdn.net/hongxingabc/article/details/51606520

李天庆等，Snake模型综述，计算机工程，2005，第31卷 第9期