首页 > 编程学习 > 主动轮廓线模型Snake模型简介amp;openCV中cvSnakeImage()函数代码分析

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

 

Snake模型称为动态轮廓模型(Active Contour  Model)是Kass与1987年提出的,它对于在噪声和对比度不敏感,能将目标从复杂背景中分割出来,并能有效的跟踪目标的形变和非刚体的复杂运动而被广泛用于图像分割和物体跟踪等图像处理领域。

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

 

Opencv中snake函数解析及例子:

 

[cpp] view plain copy

  1. /*M/// 
  2. // 
  3. //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 
  4. // 
  5. //  By downloading, copying, installing or using the software you agree to this license. 
  6. //  If you do not agree to this license, do not download, install, 
  7. //  copy or use the software. 
  8. // 
  9. // 
  10. //                        Intel License Agreement 
  11. //                For Open Source Computer Vision Library 
  12. // 
  13. // Copyright (C) 2000, Intel Corporation, all rights reserved. 
  14. // Third party copyrights are property of their respective owners. 
  15. // 
  16. // Redistribution and use in source and binary forms, with or without modification, 
  17. // are permitted provided that the following conditions are met: 
  18. // 
  19. //   * Redistribution's of source code must retain the above copyright notice, 
  20. //     this list of conditions and the following disclaimer. 
  21. // 
  22. //   * Redistribution's in binary form must reproduce the above copyright notice, 
  23. //     this list of conditions and the following disclaimer in the documentation 
  24. //     and/or other materials provided with the distribution. 
  25. // 
  26. //   * The name of Intel Corporation may not be used to endorse or promote products 
  27. //     derived from this software without specific prior written permission. 
  28. // 
  29. // This software is provided by the copyright holders and contributors "as is" and 
  30. // any express or implied warranties, including, but not limited to, the implied 
  31. // warranties of merchantability and fitness for a particular purpose are disclaimed. 
  32. // In no event shall the Intel Corporation or contributors be liable for any direct, 
  33. // indirect, incidental, special, exemplary, or consequential damages 
  34. // (including, but not limited to, procurement of substitute goods or services; 
  35. // loss of use, data, or profits; or business interruption) however caused 
  36. // and on any theory of liability, whether in contract, strict liability, 
  37. // or tort (including negligence or otherwise) arising in any way out of 
  38. // the use of this software, even if advised of the possibility of such damage. 
  39. // 
  40. //M*/  
  41. #include "_cv.h"  
  42.    
  43. #define _CV_SNAKE_BIG 2.e+38f  
  44. #define _CV_SNAKE_IMAGE 1  
  45. #define _CV_SNAKE_GRAD  2  
  46.    
  47.    
  48. /*F/// 
  49. //    Name:      icvSnake8uC1R     
  50. //    Purpose:   
  51. //    Context:   
  52. //    Parameters: 
  53. //               src - source image, 
  54. //               srcStep - its step in bytes, 
  55. //               roi - size of ROI, 
  56. //               pt - pointer to snake points array 
  57. //               n - size of points array, 
  58. //               alpha - pointer to coefficient of continuity energy, 
  59. //               beta - pointer to coefficient of curvature energy,  
  60. //               gamma - pointer to coefficient of image energy,  
  61. //               coeffUsage - if CV_VALUE - alpha, beta, gamma point to single value 
  62. //                            if CV_MATAY - point to arrays 
  63. //               criteria - termination criteria. 
  64. //               scheme - image energy scheme 
  65. //                         if _CV_SNAKE_IMAGE - image intensity is energy 
  66. //                         if _CV_SNAKE_GRAD  - magnitude of gradient is energy 
  67. //    Returns:   
  68. //F*/  
  69.    
  70. static CvStatus  
  71. icvSnake8uC1R( unsigned char *src,   //原始图像数据  
  72.                int srcStep,         //每行的字节数  
  73.                CvSize roi,         //图像尺寸  
  74.                CvPoint * pt,       //轮廓点(变形对象)  
  75.                int n,            //轮廓点的个数  
  76.                float *alpha,       //指向α的指针,α可以是单个值,也可以是与轮廓点个数一致的数组  
  77.                float *beta,        //β的值,同α  
  78.                float *gamma,       //γ的值,同α  
  79.                int coeffUsage,   //确定αβγ是用作单个值还是个数组  
  80.         CvSize win,       //每个点用于搜索的最小的领域大小,宽度为奇数  
  81.              CvTermCriteria criteria,   //递归迭代终止的条件准则  
  82. int scheme )         //确定图像能量场的数据选择,1为灰度,2为灰度梯度  
  83. {  
  84.     int i, j, k;  
  85.     int neighbors = win.height * win.width;    //当前点领域中点的个数  
  86.    
  87.    //当前点的位置  
  88.     int centerx = win.width >> 1;            
  89.     int centery = win.height >> 1;           
  90.    
  91.     float invn;        //n 的倒数?  
  92.     int iteration = 0;     //迭代次数  
  93.     int converged = 0;      //收敛标志,0为非收敛  
  94.       
  95.   //能量  
  96.     float *Econt;    //  
  97.     float *Ecurv;   //轮廓曲线能量  
  98.     float *Eimg;    //图像能量  
  99.     float *E;      //  
  100.     
  101.    //αβγ的副本  
  102.     float _alpha, _beta, _gamma;  
  103.    
  104.     /*#ifdef GRAD_SNAKE */  
  105.     float *gradient = NULL;  
  106.     uchar *map = NULL;  
  107.     int map_width = ((roi.width - 1) >> 3) + 1;  
  108.     int map_height = ((roi.height - 1) >> 3) + 1;  
  109.     CvSepFilter pX, pY;  
  110.     #define WTILE_SIZE 8  
  111.     #define TILE_SIZE (WTILE_SIZE + 2)         
  112.     short dx[TILE_SIZE*TILE_SIZE], dy[TILE_SIZE*TILE_SIZE];  
  113.     CvMat _dx = cvMat( TILE_SIZE, TILE_SIZE, CV_16SC1, dx );  
  114.     CvMat _dy = cvMat( TILE_SIZE, TILE_SIZE, CV_16SC1, dy );  
  115.     CvMat _src = cvMat( roi.height, roi.width, CV_8UC1, src );  
  116.    
  117.     /* inner buffer of convolution process */  
  118.     //char ConvBuffer[400];  
  119.    
  120.     /*#endif */  
  121.    
  122.      //检点参数的合理性  
  123.     /* check bad arguments */  
  124.     if( src == NULL )  
  125.         return CV_NULLPTR_ERR;  
  126.     if( (roi.height <= 0) || (roi.width <= 0) )  
  127.         return CV_BADSIZE_ERR;  
  128.     if( srcStep < roi.width )  
  129.         return CV_BADSIZE_ERR;  
  130.     if( pt == NULL )  
  131.         return CV_NULLPTR_ERR;  
  132.     if( n < 3 )                         //轮廓点至少要三个  
  133.         return CV_BADSIZE_ERR;  
  134.     if( alpha == NULL )  
  135.         return CV_NULLPTR_ERR;  
  136.     if( beta == NULL )  
  137.         return CV_NULLPTR_ERR;  
  138.     if( gamma == NULL )  
  139.         return CV_NULLPTR_ERR;  
  140.     if( coeffUsage != CV_VALUE && coeffUsage != CV_ARRAY )  
  141.         return CV_BADFLAG_ERR;  
  142.     if( (win.height <= 0) || (!(win.height & 1)))   //邻域搜索窗口得是奇数  
  143.         return CV_BADSIZE_ERR;  
  144.     if( (win.width <= 0) || (!(win.width & 1)))  
  145.         return CV_BADSIZE_ERR;  
  146.    
  147.     invn = 1 / ((float) n);        //轮廓点数n的倒数,用于求平均?  
  148.    
  149.     if( scheme == _CV_SNAKE_GRAD )  
  150. {  
  151.      //X方向上和Y方向上的Scoble梯度算子,用于求图像的梯度,  
  152. //处理的图像最大尺寸为TILE_SIZE+2,此例为12,算子半长为3即{-3,-2,-1,0,1,2,3}  
  153. //处理后的数据类型为16位符号数,分别存放在_dx,_dy矩阵中,长度为10  
  154.         pX.init_deriv( TILE_SIZE+2, CV_8UC1, CV_16SC1, 1, 0, 3 );  
  155.         pY.init_deriv( TILE_SIZE+2, CV_8UC1, CV_16SC1, 0, 1, 3 );  
  156.        //图像梯度存放缓冲区  
  157.         gradient = (float *) cvAlloc( roi.height * roi.width * sizeoffloat ));  
  158.    
  159.         if( !gradient )  
  160.             return CV_OUTOFMEM_ERR;  
  161.        //map用于标志相应位置的分块的图像能量是否已经求得  
  162.         map = (uchar *) cvAlloc( map_width * map_height );  
  163.         if( !map )  
  164.         {  
  165.             cvFree( &gradient );  
  166.             return CV_OUTOFMEM_ERR;  
  167.         }  
  168.         /* clear map - no gradient computed */  
  169.        //清除map标志  
  170.         memset( (void *) map, 0, map_width * map_height );  
  171. }  
  172. //各种能量的存放处,取每点的邻域的能量  
  173.     Econt = (float *) cvAlloc( neighbors * sizeoffloat ));  
  174.     Ecurv = (float *) cvAlloc( neighbors * sizeoffloat ));  
  175.     Eimg = (float *) cvAlloc( neighbors * sizeoffloat ));  
  176.     E = (float *) cvAlloc( neighbors * sizeoffloat ));  
  177.    //开始迭代  
  178.     while( !converged )    //收敛标志无效时进行  
  179.     {  
  180.         float ave_d = 0;  //轮廓各点的平均距离  
  181.         int moved = 0;      //轮廓变形时,发生移动的数量  
  182.    
  183.         converged = 0;       //标志未收敛  
  184.         iteration++;        //更新迭代次数+1  
  185.    
  186. //计算轮廓中各点的平均距离  
  187.         /* compute average distance */  
  188.       //从点0到点n-1的距离和  
  189.         for( i = 1; i < n; i++ )  
  190.         {  
  191.             int diffx = pt[i - 1].x - pt[i].x;  
  192.             int diffy = pt[i - 1].y - pt[i].y;  
  193.    
  194.             ave_d += cvSqrt( (float) (diffx * diffx + diffy * diffy) );   
  195.         }  
  196.      //再加上从点n-1到点0的距离,形成回路轮廓  
  197.         ave_d += cvSqrt( (float) ((pt[0].x - pt[n - 1].x) *  
  198.                                   (pt[0].x - pt[n - 1].x) +  
  199.                                   (pt[0].y - pt[n - 1].y) * (pt[0].y - pt[n - 1].y)));  
  200.     //求平均,得出平均距离  
  201.         ave_d *= invn;  
  202.         /* average distance computed */  
  203.    
  204.    
  205.       //对于每个轮廓点进行特定循环迭代求解  
  206.         for( i = 0; i < n; i++ )  
  207.         {  
  208.             /* Calculate Econt */  
  209.           //初始化各个能量  
  210.             float maxEcont = 0;  
  211.             float maxEcurv = 0;  
  212.             float maxEimg = 0;  
  213.             float minEcont = _CV_SNAKE_BIG;  
  214.             float minEcurv = _CV_SNAKE_BIG;  
  215.             float minEimg = _CV_SNAKE_BIG;  
  216.             float Emin = _CV_SNAKE_BIG;  
  217.          //初始化变形后轮廓点的偏移量  
  218.             int offsetx = 0;  
  219.             int offsety = 0;  
  220.             float tmp;  
  221.    
  222.         //计算边界  
  223.             /* compute bounds */  
  224.            //计算合理的搜索边界,以防领域搜索超过ROI图像的范围  
  225.             int left = MIN( pt[i].x, win.width >> 1 );  
  226.             int right = MIN( roi.width - 1 - pt[i].x, win.width >> 1 );  
  227.             int upper = MIN( pt[i].y, win.height >> 1 );  
  228.             int bottom = MIN( roi.height - 1 - pt[i].y, win.height >> 1 );  
  229.           //初始化Econt  
  230.             maxEcont = 0;  
  231.             minEcont = _CV_SNAKE_BIG;  
  232.          //在合理的搜索范围内进行Econt的计算  
  233.             for( j = -upper; j <= bottom; j++ )  
  234.             {  
  235.                 for( k = -left; k <= right; k++ )  
  236.                 {  
  237.                     int diffx, diffy;  
  238.                     float energy;  
  239.              //在轮廓点集的首尾相接处作相应处理,求轮廓点差分  
  240.                     if( i == 0 )  
  241.                     {  
  242.                         diffx = pt[n - 1].x - (pt[i].x + k);  
  243.                         diffy = pt[n - 1].y - (pt[i].y + j);  
  244.                     }  
  245.                     else  
  246.              //在其他地方作一般处理  
  247.    
  248.                     {  
  249.                         diffx = pt[i - 1].x - (pt[i].x + k);  
  250.                         diffy = pt[i - 1].y - (pt[i].y + j);  
  251.                     }  
  252.              //将邻域陈列坐标转成Econt数组的下标序号,计算邻域中每点的Econt  
  253.               //Econt的值等于平均距离和此点和上一点的距离的差的绝对值(这是怎么来的?)  
  254.                     Econt[(j + centery) * win.width + k + centerx] = energy =  
  255.                         (float) fabs( ave_d -  
  256.                                       cvSqrt( (float) (diffx * diffx + diffy * diffy) ));  
  257.              //求出所有邻域点中的Econt的最大值和最小值  
  258.                     maxEcont = MAX( maxEcont, energy );  
  259.                     minEcont = MIN( minEcont, energy );  
  260.                 }  
  261.             }  
  262.            //求出邻域点中最大值和最小值之差,并对所有的邻域点的Econt进行标准归一化,若最大值最小  
  263.            //相等,则邻域中的点Econt全相等,Econt归一化束缚为0  
  264.             tmp = maxEcont - minEcont;  
  265.             tmp = (tmp == 0) ? 0 : (1 / tmp);  
  266.             for( k = 0; k < neighbors; k++ )  
  267.             {  
  268.                 Econt[k] = (Econt[k] - minEcont) * tmp;  
  269.             }  
  270.    
  271.    
  272.            //计算每点的Ecurv  
  273.             /*  Calculate Ecurv */  
  274.             maxEcurv = 0;  
  275.             minEcurv = _CV_SNAKE_BIG;  
  276.             for( j = -upper; j <= bottom; j++ )  
  277.             {  
  278.                 for( k = -left; k <= right; k++ )  
  279.                 {  
  280.                     int tx, ty;  
  281.                     float energy;  
  282.                     //第一个点的二阶差分  
  283.                     if( i == 0 )  
  284.                     {  
  285.                         tx = pt[n - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;  
  286.                         ty = pt[n - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;  
  287.                     }  
  288.                    //最后一个点的二阶差分  
  289.                     else if( i == n - 1 )  
  290.                     {  
  291.                         tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[0].x;  
  292.                         ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[0].y;  
  293.                     }  
  294.                    //其余点的二阶差分  
  295.                     else  
  296.                     {  
  297.                         tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;  
  298.                         ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;  
  299.                     }  
  300.                   //转换坐标为数组序号,并求各点的Ecurv的值,二阶差分后取平方  
  301.                     Ecurv[(j + centery) * win.width + k + centerx] = energy =  
  302.                         (float) (tx * tx + ty * ty);  
  303.                   //取最小的Ecurv和最大的Ecurv  
  304.                     maxEcurv = MAX( maxEcurv, energy );  
  305.                     minEcurv = MIN( minEcurv, energy );  
  306.                 }  
  307.             }  
  308.                //对Ecurv进行标准归一化  
  309.             tmp = maxEcurv - minEcurv;  
  310.             tmp = (tmp == 0) ? 0 : (1 / tmp);  
  311.             for( k = 0; k < neighbors; k++ )  
  312.             {  
  313.                 Ecurv[k] = (Ecurv[k] - minEcurv) * tmp;  
  314.             }  
  315.            
  316.            //求Eimg  
  317.             /* Calculate Eimg */  
  318.             for( j = -upper; j <= bottom; j++ )  
  319.             {  
  320.                 for( k = -left; k <= right; k++ )  
  321.                 {  
  322.                     float energy;  
  323.                //若采用灰度梯度数据  
  324.                     if( scheme == _CV_SNAKE_GRAD )  
  325.                     {  
  326.                         /* look at map and check status */  
  327.                         int x = (pt[i].x + k)/WTILE_SIZE;  
  328.                         int y = (pt[i].y + j)/WTILE_SIZE;  
  329.                         //若此处的图像能量还没有获取,则对此处对应的图像分块进行图像能量的求解  
  330.                         if( map[y * map_width + x] == 0 )  
  331.                         {  
  332.                             int l, m;                             
  333.    
  334.                             /* evaluate block location */  
  335.                            //计算要进行梯度算子处理的图像块的位置  
  336.                             int upshift = y ? 1 : 0;  
  337.                             int leftshift = x ? 1 : 0;  
  338.                             int bottomshift = MIN( 1, roi.height - (y + 1)*WTILE_SIZE );  
  339.                             int rightshift = MIN( 1, roi.width - (x + 1)*WTILE_SIZE );  
  340.                           //图像块的位置大小(由于原ROI不一定是8的倍数,所以图像块会大小不一)  
  341.                             CvRect g_roi = { x*WTILE_SIZE - leftshift, y*WTILE_SIZE - upshift,  
  342.                                 leftshift + WTILE_SIZE + rightshift, upshift + WTILE_SIZE + bottomshift };  
  343.                             CvMat _src1;  
  344.                             cvGetSubArr( &_src, &_src1, g_roi );  //得到图像块的数据  
  345.                             //分别对图像的X方向和Y方向进行梯度算子  
  346.                             pX.process( &_src1, &_dx );  
  347.                             pY.process( &_src1, &_dy );  
  348.                          //求分块区域中的每个点的梯度  
  349.                             for( l = 0; l < WTILE_SIZE + bottomshift; l++ )  
  350.                             {  
  351.                                 for( m = 0; m < WTILE_SIZE + rightshift; m++ )  
  352.                                 {  
  353.                                     gradient[(y*WTILE_SIZE + l) * roi.width + x*WTILE_SIZE + m] =  
  354.                                         (float) (dx[(l + upshift) * TILE_SIZE + m + leftshift] *  
  355.                                                  dx[(l + upshift) * TILE_SIZE + m + leftshift] +  
  356.                                                  dy[(l + upshift) * TILE_SIZE + m + leftshift] *  
  357.                                                  dy[(l + upshift) * TILE_SIZE + m + leftshift]);  
  358.                                 }  
  359.                             }  
  360.                             //map相应位置置1表示此处图像能量已经获取  
  361.                             map[y * map_width + x] = 1;  
  362.                         }  
  363.                       //以梯度数据作为图像能量  
  364.                         Eimg[(j + centery) * win.width + k + centerx] = energy =  
  365.                             gradient[(pt[i].y + j) * roi.width + pt[i].x + k];  
  366.                     }  
  367.                     else  
  368.                     {  
  369.                        //以灰度作为图像能量  
  370.                         Eimg[(j + centery) * win.width + k + centerx] = energy =  
  371.                             src[(pt[i].y + j) * srcStep + pt[i].x + k];  
  372.                     }  
  373.                    //获得邻域中最大和最小的图像能量  
  374.                     maxEimg = MAX( maxEimg, energy );  
  375.                     minEimg = MIN( minEimg, energy );  
  376.                 }  
  377.             }  
  378.               //Eimg的标准归一化  
  379.             tmp = (maxEimg - minEimg);  
  380.             tmp = (tmp == 0) ? 0 : (1 / tmp);  
  381.    
  382.             for( k = 0; k < neighbors; k++ )  
  383.             {  
  384.                 Eimg[k] = (minEimg - Eimg[k]) * tmp;  
  385.             }  
  386.             //加入系数  
  387.             /* locate coefficients */  
  388.             if( coeffUsage == CV_VALUE)  
  389.             {  
  390.                 _alpha = *alpha;  
  391.                 _beta = *beta;  
  392.                 _gamma = *gamma;  
  393.             }  
  394.             else  
  395.             {                    
  396.                 _alpha = alpha[i];  
  397.                 _beta = beta[i];  
  398.                 _gamma = gamma[i];  
  399.             }  
  400.    
  401.             /* Find Minimize point in the neighbors */  
  402.             //求得每个邻域点的Snake能量  
  403.             for( k = 0; k < neighbors; k++ )  
  404.             {  
  405.                 E[k] = _alpha * Econt[k] + _beta * Ecurv[k] + _gamma * Eimg[k];  
  406.             }  
  407.             Emin = _CV_SNAKE_BIG;  
  408.         //获取最小的能量,以及对应的邻域中的相对位置  
  409.             for( j = -upper; j <= bottom; j++ )  
  410.             {  
  411.                 for( k = -left; k <= right; k++ )  
  412.                 {  
  413.    
  414.                     if( E[(j + centery) * win.width + k + centerx] < Emin )  
  415.                     {  
  416.                         Emin = E[(j + centery) * win.width + k + centerx];  
  417.                         offsetx = k;  
  418.                         offsety = j;  
  419.                     }  
  420.                 }  
  421.             }  
  422.          //如果轮廓点发生改变,则记得移动次数  
  423.             if( offsetx || offsety )  
  424.             {  
  425.                 pt[i].x += offsetx;  
  426.                 pt[i].y += offsety;  
  427.                 moved++;  
  428.             }  
  429.         }  
  430.    
  431.       //各个轮廓点迭代计算完成后,如果没有移动的点了,则收敛标志位有效,停止迭代  
  432.         converged = (moved == 0);  
  433.      //达到最大迭代次数时,收敛标志位有效,停止迭代  
  434.         if( (criteria.type & CV_TERMCRIT_ITER) && (iteration >= criteria.max_iter) )  
  435.             converged = 1;  
  436.   //到大相应精度时,停止迭代(与第一个条件有相同效果)  
  437.         if( (criteria.type & CV_TERMCRIT_EPS) && (moved <= criteria.epsilon) )  
  438.             converged = 1;  
  439.     }  
  440.    
  441.   //释放各个缓冲区  
  442.     cvFree( &Econt );  
  443.     cvFree( &Ecurv );  
  444.     cvFree( &Eimg );  
  445.     cvFree( &E );  
  446.    
  447.     if( scheme == _CV_SNAKE_GRAD )  
  448.     {  
  449.         cvFree( &gradient );  
  450.         cvFree( &map );  
  451.     }  
  452.     return CV_OK;  
  453. }  
  454.    
  455.    
  456. CV_IMPL void  
  457. cvSnakeImage( const IplImage* src, CvPoint* points,  
  458.               int length, float *alpha,  
  459.               float *beta, float *gamma,  
  460.               int coeffUsage, CvSize win,  
  461.               CvTermCriteria criteria, int calcGradient )  
  462. {  
  463.    
  464.     CV_FUNCNAME( "cvSnakeImage" );  
  465.    
  466.     __BEGIN__;  
  467.    
  468.     uchar *data;  
  469.     CvSize size;  
  470.     int step;  
  471.    
  472.     if( src->nChannels != 1 )  
  473.         CV_ERROR( CV_BadNumChannels, "input image has more than one channel" );  
  474.    
  475.     if( src->depth != IPL_DEPTH_8U )  
  476.         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );  
  477.    
  478.     cvGetRawData( src, &data, &step, &size );  
  479.    
  480.     IPPI_CALL( icvSnake8uC1R( data, step, size, points, length,  
  481.                               alpha, beta, gamma, coeffUsage, win, criteria,  
  482.                               calcGradient ? _CV_SNAKE_GRAD : _CV_SNAKE_IMAGE ));  
  483.     __END__;  
  484. }  
  485.    
  486. /* end of file */  
  487.    
  488.    
  489.    
  490.    
  491.    
  492. 测试应用程序  
  493.    
  494. #include "stdafx.h"  
  495. #include <iostream>  
  496. #include <string.h>  
  497. #include <cxcore.h>  
  498. #include <cv.h>  
  499. #include <highgui.h>  
  500. #include <fstream>  
  501.    
  502.    
  503. IplImage *image = 0 ; //原始图像  
  504. IplImage *image2 = 0 ; //原始图像copy  
  505.    
  506. using namespace std;  
  507. int Thresholdness = 141;  
  508. int ialpha = 20;  
  509. int ibeta=20;  
  510. int igamma=20;  
  511.    
  512. void onChange(int pos)  
  513. {  
  514.      
  515.     if(image2) cvReleaseImage(&image2);  
  516.     if(image) cvReleaseImage(&image);  
  517.    
  518.     image2 = cvLoadImage("grey.bmp",1); //显示图片  
  519.     image= cvLoadImage("grey.bmp",0);  
  520.    
  521.     cvThreshold(image,image,Thresholdness,255,CV_THRESH_BINARY); //分割域值     
  522.    
  523.     CvMemStorage* storage = cvCreateMemStorage(0);  
  524.     CvSeq* contours = 0;  
  525.    
  526.     cvFindContours( image, storage, &contours, sizeof(CvContour), //寻找初始化轮廓  
  527.         CV_RETR_EXTERNAL , CV_CHAIN_APPROX_SIMPLE );  
  528.    
  529.     if(!contours) return ;  
  530.     int length = contours->total;     
  531.     if(length<10) return ;  
  532.     CvPoint* point = new CvPoint[length]; //分配轮廓点  
  533.    
  534.     CvSeqReader reader;  
  535.     CvPoint pt= cvPoint(0,0);;     
  536.     CvSeq *contour2=contours;     
  537.    
  538.     cvStartReadSeq(contour2, &reader);  
  539.     for (int i = 0; i < length; i++)  
  540.     {  
  541.         CV_READ_SEQ_ELEM(pt, reader);  
  542.         point[i]=pt;  
  543.     }  
  544.     cvReleaseMemStorage(&storage);  
  545.    
  546.     //显示轮廓曲线  
  547.     for(int i=0;i<length;i++)  
  548.     {  
  549.         int j = (i+1)%length;  
  550.         cvLine( image2, point[i],point[j],CV_RGB( 0, 0, 255 ),1,8,0 );  
  551.     }  
  552.    
  553.     float alpha=ialpha/100.0f;  
  554.     float beta=ibeta/100.0f;  
  555.     float gamma=igamma/100.0f;  
  556.    
  557.     CvSize size;  
  558.     size.width=3;  
  559.     size.height=3;  
  560.     CvTermCriteria criteria;  
  561.     criteria.type=CV_TERMCRIT_ITER;  
  562.     criteria.max_iter=1000;  
  563.     criteria.epsilon=0.1;  
  564.     cvSnakeImage( image, point,length,&alpha,&beta,&gamma,CV_VALUE,size,criteria,0 );  
  565.    
  566.     //显示曲线  
  567.     for(int i=0;i<length;i++)  
  568.     {  
  569.         int j = (i+1)%length;  
  570.         cvLine( image2, point[i],point[j],CV_RGB( 0, 255, 0 ),1,8,0 );  
  571.     }  
  572.     delete []point;  
  573.    
  574. }  
  575.    
  576. int main(int argc, char* argv[])  
  577. {  
  578.    
  579.      
  580.     cvNamedWindow("win1",0);  
  581.     cvCreateTrackbar("Thd""win1", &Thresholdness, 255, onChange);  
  582.     cvCreateTrackbar("alpha""win1", &ialpha, 100, onChange);  
  583.     cvCreateTrackbar("beta""win1", &ibeta, 100, onChange);  
  584.     cvCreateTrackbar("gamma""win1", &igamma, 100, onChange);  
  585.     cvResizeWindow("win1",300,500);  
  586.     onChange(0);  
  587.    
  588.     for(;;)  
  589.     {  
  590.         if(cvWaitKey(40)==27) break;  
  591.         cvShowImage("win1",image2);  
  592.     }  
  593.      
  594.     return 0;  

 

 

 

本文由以下两篇博客整理得到:

http://blog.csdn.net/lbd2008/article/details/7180506

http://chyyeng.blog.163.com/blog/static/16918230201272724529214/

 

 

 

 

 

 

 

 

 

 

Copyright © 2010-2022 ngui.cc 版权所有 |关于我们| 联系方式| 豫B2-20100000