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 70static CvStatus 71icvSnake8uC1R( 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 408CV_IMPL void 409cvSnakeImage( 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 */ 439