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