git.maemo.org Git - opencv/blob - src/cv/cvsnakes.cpp

   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, CvSize win, CvTermCriteria criteria, int scheme )
  80 {
  81     int i, j, k;
  82     int neighbors = win.height * win.width;
  83
  84     int centerx = win.width >> 1;
  85     int centery = win.height >> 1;
  86
  87     float invn;
  88     int iteration = 0;
  89     int converged = 0;
  90
  91
  92     float *Econt;
  93     float *Ecurv;
  94     float *Eimg;
  95     float *E;
  96
  97     float _alpha, _beta, _gamma;
  98
  99     /*#ifdef GRAD_SNAKE */
 100     float *gradient = NULL;
 101     uchar *map = NULL;
 102     int map_width = ((roi.width - 1) >> 3) + 1;
 103     int map_height = ((roi.height - 1) >> 3) + 1;
 104     #define WTILE_SIZE 8
 105     #define TILE_SIZE (WTILE_SIZE + 2)
 106     short dx[TILE_SIZE*TILE_SIZE], dy[TILE_SIZE*TILE_SIZE];
 107     CvMat _dx = cvMat( TILE_SIZE, TILE_SIZE, CV_16SC1, dx );
 108     CvMat _dy = cvMat( TILE_SIZE, TILE_SIZE, CV_16SC1, dy );
 109     CvMat _src = cvMat( roi.height, roi.width, CV_8UC1, src );
 110     cv::Ptr<cv::FilterEngine> pX, pY;
 111
 112     /* inner buffer of convolution process */
 113     //char ConvBuffer[400];
 114
 115     /*#endif */
 116
 117
 118     /* check bad arguments */
 119     if( src == NULL )
 120         return CV_NULLPTR_ERR;
 121     if( (roi.height <= 0) || (roi.width <= 0) )
 122         return CV_BADSIZE_ERR;
 123     if( srcStep < roi.width )
 124         return CV_BADSIZE_ERR;
 125     if( pt == NULL )
 126         return CV_NULLPTR_ERR;
 127     if( n < 3 )
 128         return CV_BADSIZE_ERR;
 129     if( alpha == NULL )
 130         return CV_NULLPTR_ERR;
 131     if( beta == NULL )
 132         return CV_NULLPTR_ERR;
 133     if( gamma == NULL )
 134         return CV_NULLPTR_ERR;
 135     if( coeffUsage != CV_VALUE && coeffUsage != CV_ARRAY )
 136         return CV_BADFLAG_ERR;
 137     if( (win.height <= 0) || (!(win.height & 1)))
 138         return CV_BADSIZE_ERR;
 139     if( (win.width <= 0) || (!(win.width & 1)))
 140         return CV_BADSIZE_ERR;
 141
 142     invn = 1 / ((float) n);
 143
 144     if( scheme == _CV_SNAKE_GRAD )
 145     {
 146         pX = cv::createDerivFilter( CV_8U, CV_16S, 1, 0, 3, cv::BORDER_REPLICATE );
 147         pY = cv::createDerivFilter( CV_8U, CV_16S, 0, 1, 3, cv::BORDER_REPLICATE );
 148         gradient = (float *) cvAlloc( roi.height * roi.width * sizeof( float ));
 149
 150         map = (uchar *) cvAlloc( map_width * map_height );
 151         /* clear map - no gradient computed */
 152         memset( (void *) map, 0, map_width * map_height );
 153     }
 154     Econt = (float *) cvAlloc( neighbors * sizeof( float ));
 155     Ecurv = (float *) cvAlloc( neighbors * sizeof( float ));
 156     Eimg = (float *) cvAlloc( neighbors * sizeof( float ));
 157     E = (float *) cvAlloc( neighbors * sizeof( float ));
 158
 159     while( !converged )
 160     {
 161         float ave_d = 0;
 162         int moved = 0;
 163
 164         converged = 0;
 165         iteration++;
 166         /* compute average distance */
 167         for( i = 1; i < n; i++ )
 168         {
 169             int diffx = pt[i - 1].x - pt[i].x;
 170             int diffy = pt[i - 1].y - pt[i].y;
 171
 172             ave_d += cvSqrt( (float) (diffx * diffx + diffy * diffy) );
 173         }
 174         ave_d += cvSqrt( (float) ((pt[0].x - pt[n - 1].x) *
 175                                   (pt[0].x - pt[n - 1].x) +
 176                                   (pt[0].y - pt[n - 1].y) * (pt[0].y - pt[n - 1].y)));
 177
 178         ave_d *= invn;
 179         /* average distance computed */
 180         for( i = 0; i < n; i++ )
 181         {
 182             /* Calculate Econt */
 183             float maxEcont = 0;
 184             float maxEcurv = 0;
 185             float maxEimg = 0;
 186             float minEcont = _CV_SNAKE_BIG;
 187             float minEcurv = _CV_SNAKE_BIG;
 188             float minEimg = _CV_SNAKE_BIG;
 189             float Emin = _CV_SNAKE_BIG;
 190
 191             int offsetx = 0;
 192             int offsety = 0;
 193             float tmp;
 194
 195             /* compute bounds */
 196             int left = MIN( pt[i].x, win.width >> 1 );
 197             int right = MIN( roi.width - 1 - pt[i].x, win.width >> 1 );
 198             int upper = MIN( pt[i].y, win.height >> 1 );
 199             int bottom = MIN( roi.height - 1 - pt[i].y, win.height >> 1 );
 200
 201             maxEcont = 0;
 202             minEcont = _CV_SNAKE_BIG;
 203             for( j = -upper; j <= bottom; j++ )
 204             {
 205                 for( k = -left; k <= right; k++ )
 206                 {
 207                     int diffx, diffy;
 208                     float energy;
 209
 210                     if( i == 0 )
 211                     {
 212                         diffx = pt[n - 1].x - (pt[i].x + k);
 213                         diffy = pt[n - 1].y - (pt[i].y + j);
 214                     }
 215                     else
 216                     {
 217                         diffx = pt[i - 1].x - (pt[i].x + k);
 218                         diffy = pt[i - 1].y - (pt[i].y + j);
 219                     }
 220                     Econt[(j + centery) * win.width + k + centerx] = energy =
 221                         (float) fabs( ave_d -
 222                                       cvSqrt( (float) (diffx * diffx + diffy * diffy) ));
 223
 224                     maxEcont = MAX( maxEcont, energy );
 225                     minEcont = MIN( minEcont, energy );
 226                 }
 227             }
 228             tmp = maxEcont - minEcont;
 229             tmp = (tmp == 0) ? 0 : (1 / tmp);
 230             for( k = 0; k < neighbors; k++ )
 231             {
 232                 Econt[k] = (Econt[k] - minEcont) * tmp;
 233             }
 234
 235             /*  Calculate Ecurv */
 236             maxEcurv = 0;
 237             minEcurv = _CV_SNAKE_BIG;
 238             for( j = -upper; j <= bottom; j++ )
 239             {
 240                 for( k = -left; k <= right; k++ )
 241                 {
 242                     int tx, ty;
 243                     float energy;
 244
 245                     if( i == 0 )
 246                     {
 247                         tx = pt[n - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;
 248                         ty = pt[n - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;
 249                     }
 250                     else if( i == n - 1 )
 251                     {
 252                         tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[0].x;
 253                         ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[0].y;
 254                     }
 255                     else
 256                     {
 257                         tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;
 258                         ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;
 259                     }
 260                     Ecurv[(j + centery) * win.width + k + centerx] = energy =
 261                         (float) (tx * tx + ty * ty);
 262                     maxEcurv = MAX( maxEcurv, energy );
 263                     minEcurv = MIN( minEcurv, energy );
 264                 }
 265             }
 266             tmp = maxEcurv - minEcurv;
 267             tmp = (tmp == 0) ? 0 : (1 / tmp);
 268             for( k = 0; k < neighbors; k++ )
 269             {
 270                 Ecurv[k] = (Ecurv[k] - minEcurv) * tmp;
 271             }
 272
 273             /* Calculate Eimg */
 274             for( j = -upper; j <= bottom; j++ )
 275             {
 276                 for( k = -left; k <= right; k++ )
 277                 {
 278                     float energy;
 279
 280                     if( scheme == _CV_SNAKE_GRAD )
 281                     {
 282                         /* look at map and check status */
 283                         int x = (pt[i].x + k)/WTILE_SIZE;
 284                         int y = (pt[i].y + j)/WTILE_SIZE;
 285
 286                         if( map[y * map_width + x] == 0 )
 287                         {
 288                             int l, m;
 289
 290                             /* evaluate block location */
 291                             int upshift = y ? 1 : 0;
 292                             int leftshift = x ? 1 : 0;
 293                             int bottomshift = MIN( 1, roi.height - (y + 1)*WTILE_SIZE );
 294                             int rightshift = MIN( 1, roi.width - (x + 1)*WTILE_SIZE );
 295                             CvRect g_roi = { x*WTILE_SIZE - leftshift, y*WTILE_SIZE - upshift,
 296                                 leftshift + WTILE_SIZE + rightshift, upshift + WTILE_SIZE + bottomshift };
 297                             CvMat _src1;
 298                             cvGetSubArr( &_src, &_src1, g_roi );
 299
 300                             cv::Mat _src_ = cv::cvarrToMat(&_src1);
 301                             cv::Mat _dx_ = cv::cvarrToMat(&_dx);
 302                             cv::Mat _dy_ = cv::cvarrToMat(&_dy);
 303
 304                             pX->apply( _src_, _dx_, cv::Rect(0,0,-1,-1), cv::Point(), true );
 305                             pY->apply( _src_, _dy_, cv::Rect(0,0,-1,-1), cv::Point(), true );
 306
 307                             for( l = 0; l < WTILE_SIZE + bottomshift; l++ )
 308                             {
 309                                 for( m = 0; m < WTILE_SIZE + rightshift; m++ )
 310                                 {
 311                                     gradient[(y*WTILE_SIZE + l) * roi.width + x*WTILE_SIZE + m] =
 312                                         (float) (dx[(l + upshift) * TILE_SIZE + m + leftshift] *
 313                                                  dx[(l + upshift) * TILE_SIZE + m + leftshift] +
 314                                                  dy[(l + upshift) * TILE_SIZE + m + leftshift] *
 315                                                  dy[(l + upshift) * TILE_SIZE + m + leftshift]);
 316                                 }
 317                             }
 318                             map[y * map_width + x] = 1;
 319                         }
 320                         Eimg[(j + centery) * win.width + k + centerx] = energy =
 321                             gradient[(pt[i].y + j) * roi.width + pt[i].x + k];
 322                     }
 323                     else
 324                     {
 325                         Eimg[(j + centery) * win.width + k + centerx] = energy =
 326                             src[(pt[i].y + j) * srcStep + pt[i].x + k];
 327                     }
 328
 329                     maxEimg = MAX( maxEimg, energy );
 330                     minEimg = MIN( minEimg, energy );
 331                 }
 332             }
 333
 334             tmp = (maxEimg - minEimg);
 335             tmp = (tmp == 0) ? 0 : (1 / tmp);
 336
 337             for( k = 0; k < neighbors; k++ )
 338             {
 339                 Eimg[k] = (minEimg - Eimg[k]) * tmp;
 340             }
 341
 342             /* locate coefficients */
 343             if( coeffUsage == CV_VALUE)
 344             {
 345                 _alpha = *alpha;
 346                 _beta = *beta;
 347                 _gamma = *gamma;
 348             }
 349             else
 350             {
 351                 _alpha = alpha[i];
 352                 _beta = beta[i];
 353                 _gamma = gamma[i];
 354             }
 355
 356             /* Find Minimize point in the neighbors */
 357             for( k = 0; k < neighbors; k++ )
 358             {
 359                 E[k] = _alpha * Econt[k] + _beta * Ecurv[k] + _gamma * Eimg[k];
 360             }
 361             Emin = _CV_SNAKE_BIG;
 362             for( j = -upper; j <= bottom; j++ )
 363             {
 364                 for( k = -left; k <= right; k++ )
 365                 {
 366
 367                     if( E[(j + centery) * win.width + k + centerx] < Emin )
 368                     {
 369                         Emin = E[(j + centery) * win.width + k + centerx];
 370                         offsetx = k;
 371                         offsety = j;
 372                     }
 373                 }
 374             }
 375
 376             if( offsetx || offsety )
 377             {
 378                 pt[i].x += offsetx;
 379                 pt[i].y += offsety;
 380                 moved++;
 381             }
 382         }
 383         converged = (moved == 0);
 384         if( (criteria.type & CV_TERMCRIT_ITER) && (iteration >= criteria.max_iter) )
 385             converged = 1;
 386         if( (criteria.type & CV_TERMCRIT_EPS) && (moved <= criteria.epsilon) )
 387             converged = 1;
 388     }
 389
 390     cvFree( &Econt );
 391     cvFree( &Ecurv );
 392     cvFree( &Eimg );
 393     cvFree( &E );
 394
 395     if( scheme == _CV_SNAKE_GRAD )
 396     {
 397         cvFree( &gradient );
 398         cvFree( &map );
 399     }
 400     return CV_OK;
 401 }
 402
 403
 404 CV_IMPL void
 405 cvSnakeImage( const IplImage* src, CvPoint* points,
 406               int length, float *alpha,
 407               float *beta, float *gamma,
 408               int coeffUsage, CvSize win,
 409               CvTermCriteria criteria, int calcGradient )
 410 {
 411
 412     CV_FUNCNAME( "cvSnakeImage" );
 413
 414     __BEGIN__;
 415
 416     uchar *data;
 417     CvSize size;
 418     int step;
 419
 420     if( src->nChannels != 1 )
 421         CV_ERROR( CV_BadNumChannels, "input image has more than one channel" );
 422
 423     if( src->depth != IPL_DEPTH_8U )
 424         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
 425
 426     cvGetRawData( src, &data, &step, &size );
 427
 428     IPPI_CALL( icvSnake8uC1R( data, step, size, points, length,
 429                               alpha, beta, gamma, coeffUsage, win, criteria,
 430                               calcGradient ? _CV_SNAKE_GRAD : _CV_SNAKE_IMAGE ));
 431     __END__;
 432 }
 433
 434 /* end of file */