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