主动轮廓的snake模型

#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;  
} 

基于能量泛函的切割方法:
该类方法主要指的是活动轮廓模型(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模型的基本思想非常easy,它以构成一定形状的一些控制点为模板(轮廓线),通过模板自身的弹性形变,与图像局部特征相匹配达到调和,即某种能量函数极小化,完毕对图像的切割。再通过对模板的进一步分析而实现图像的理解和识别。
简单的来讲,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就停在图像的边缘位置了,也就完毕了切割。那么这个的前提就是目标在图像中的边缘比較明显了,否则非常easy就越过边缘了。)
弹性能量和弯曲能量合称内部能量(内部力),用于控制轮廓线的弹性形变,起到保持轮廓连续性和平滑性的作用。而第三项代表外部能量,也被称为图像能量,表示变形曲线与图像局部特征吻合的情况。内部能量只跟snake的形状有关,而跟图像数据无关。而外部能量只跟图像数据有关。在某一点的α和β的值决定曲线能够在这一点伸展和弯曲的程度。
终于对图像的切割转化为求解能量函数Etotal(v)极小化(最小化轮廓的能量)。在能量函数极小化过程中,弹性能量迅速把轮廓线压缩成一个光滑的圆,弯曲能量驱使轮廓线成为光滑曲线或直线,而图像力则使轮廓线向图像的高梯度位置靠拢。基本Snakes模型就是在这3个力的联合作用下工作的。
由于图像上的点都是离散的,所以我们用来优化能量函数的算法都必须在离散域里定义。所以求解能量函数Etotal(v)极小化是一个典型的变分问题(微分运算中,自变量通常是坐标等变量,因变量是函数;变分运算中,自变量是函数,因变量是函数的函数,即数学上所谓的泛函。对泛函求极值的问题,数学上称之为变分法)。
在离散化条件(数字图像)下,由欧拉方程可知终于问题的答案等价于求解一组差分方程:(欧拉方程是泛函极值条件的微分表达式,求解泛函的欧拉方程,就可以得到使泛函取极值的驻函数,将变分问题转化为微分问题。)

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

 以上仅仅是对snake简单的理解,如要深入,请參考其它很多其它专业文献。水平有限,错误在所难免,还望不吝指正

猜你喜欢

转载自blog.csdn.net/qq_25254777/article/details/80165730