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
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)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
cvSnakeImage(const IplImage * src,CvPoint * points,int length,float * alpha,float * beta,float * gamma,int coeffUsage,CvSize win,CvTermCriteria criteria,int calcGradient)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 */
439