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

   1 /*M///////////////////////////////////////////////////////////////////////////////////////\r
   2 //\r
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
   4 //\r
   5 //  By downloading, copying, installing or using the software you agree to this license.\r
   6 //  If you do not agree to this license, do not download, install,\r
   7 //  copy or use the software.\r
   8 //\r
   9 //\r
  10 //                        Intel License Agreement\r
  11 //                For Open Source Computer Vision Library\r
  12 //\r
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.\r
  14 // Third party copyrights are property of their respective owners.\r
  15 //\r
  16 // Redistribution and use in source and binary forms, with or without modification,\r
  17 // are permitted provided that the following conditions are met:\r
  18 //\r
  19 //   * Redistribution's of source code must retain the above copyright notice,\r
  20 //     this list of conditions and the following disclaimer.\r
  21 //\r
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,\r
  23 //     this list of conditions and the following disclaimer in the documentation\r
  24 //     and/or other materials provided with the distribution.\r
  25 //\r
  26 //   * The name of Intel Corporation may not be used to endorse or promote products\r
  27 //     derived from this software without specific prior written permission.\r
  28 //\r
  29 // This software is provided by the copyright holders and contributors "as is" and\r
  30 // any express or implied warranties, including, but not limited to, the implied\r
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.\r
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,\r
  33 // indirect, incidental, special, exemplary, or consequential damages\r
  34 // (including, but not limited to, procurement of substitute goods or services;\r
  35 // loss of use, data, or profits; or business interruption) however caused\r
  36 // and on any theory of liability, whether in contract, strict liability,\r
  37 // or tort (including negligence or otherwise) arising in any way out of\r
  38 // the use of this software, even if advised of the possibility of such damage.\r
  39 //\r
  40 //M*/\r
  41 \r
  42 #include "_cv.h"\r
  43 #include "_cvlist.h"\r
  44 \r
  45 #define halfPi ((float)(CV_PI*0.5))\r
  46 #define Pi     ((float)CV_PI)\r
  47 #define a0  0 /*-4.172325e-7f*/   /*(-(float)0x7)/((float)0x1000000); */\r
  48 #define a1 1.000025f        /*((float)0x1922253)/((float)0x1000000)*2/Pi; */\r
  49 #define a2 -2.652905e-4f    /*(-(float)0x2ae6)/((float)0x1000000)*4/(Pi*Pi); */\r
  50 #define a3 -0.165624f       /*(-(float)0xa45511)/((float)0x1000000)*8/(Pi*Pi*Pi); */\r
  51 #define a4 -1.964532e-3f    /*(-(float)0x30fd3)/((float)0x1000000)*16/(Pi*Pi*Pi*Pi); */\r
  52 #define a5 1.02575e-2f      /*((float)0x191cac)/((float)0x1000000)*32/(Pi*Pi*Pi*Pi*Pi); */\r
  53 #define a6 -9.580378e-4f    /*(-(float)0x3af27)/((float)0x1000000)*64/(Pi*Pi*Pi*Pi*Pi*Pi); */\r
  54 \r
  55 #define _sin(x) ((((((a6*(x) + a5)*(x) + a4)*(x) + a3)*(x) + a2)*(x) + a1)*(x) + a0)\r
  56 #define _cos(x) _sin(halfPi - (x))\r
  57 \r
  58 /****************************************************************************************\\r
  59 *                               Classical Hough Transform                                *\r
  60 \****************************************************************************************/\r
  61 \r
  62 typedef struct CvLinePolar\r
  63 {\r
  64     float rho;\r
  65     float angle;\r
  66 }\r
  67 CvLinePolar;\r
  68 \r
  69 /*=====================================================================================*/\r
  70 \r
  71 #define hough_cmp_gt(l1,l2) (aux[l1] > aux[l2])\r
  72 \r
  73 static CV_IMPLEMENT_QSORT_EX( icvHoughSortDescent32s, int, hough_cmp_gt, const int* )\r
  74 \r
  75 /*\r
  76 Here image is an input raster;\r
  77 step is it's step; size characterizes it's ROI;\r
  78 rho and theta are discretization steps (in pixels and radians correspondingly).\r
  79 threshold is the minimum number of pixels in the feature for it\r
  80 to be a candidate for line. lines is the output\r
  81 array of (rho, theta) pairs. linesMax is the buffer size (number of pairs).\r
  82 Functions return the actual number of found lines.\r
  83 */\r
  84 static void\r
  85 icvHoughLinesStandard( const CvMat* img, float rho, float theta,\r
  86                        int threshold, CvSeq *lines, int linesMax )\r
  87 {\r
  88     int *accum = 0;\r
  89     int *sort_buf=0;\r
  90     float *tabSin = 0;\r
  91     float *tabCos = 0;\r
  92 \r
  93     CV_FUNCNAME( "icvHoughLinesStandard" );\r
  94 \r
  95     __BEGIN__;\r
  96 \r
  97     const uchar* image;\r
  98     int step, width, height;\r
  99     int numangle, numrho;\r
 100     int total = 0;\r
 101     float ang;\r
 102     int r, n;\r
 103     int i, j;\r
 104     float irho = 1 / rho;\r
 105     double scale;\r
 106 \r
 107     CV_ASSERT( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1 );\r
 108 \r
 109     image = img->data.ptr;\r
 110     step = img->step;\r
 111     width = img->cols;\r
 112     height = img->rows;\r
 113 \r
 114     numangle = cvRound(CV_PI / theta);\r
 115     numrho = cvRound(((width + height) * 2 + 1) / rho);\r
 116 \r
 117     CV_CALL( accum = (int*)cvAlloc( sizeof(accum[0]) * (numangle+2) * (numrho+2) ));\r
 118     CV_CALL( sort_buf = (int*)cvAlloc( sizeof(accum[0]) * numangle * numrho ));\r
 119     CV_CALL( tabSin = (float*)cvAlloc( sizeof(tabSin[0]) * numangle ));\r
 120     CV_CALL( tabCos = (float*)cvAlloc( sizeof(tabCos[0]) * numangle ));\r
 121     memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );\r
 122 \r
 123     for( ang = 0, n = 0; n < numangle; ang += theta, n++ )\r
 124     {\r
 125         tabSin[n] = (float)(sin(ang) * irho);\r
 126         tabCos[n] = (float)(cos(ang) * irho);\r
 127     }\r
 128 \r
 129     // stage 1. fill accumulator\r
 130     for( i = 0; i < height; i++ )\r
 131         for( j = 0; j < width; j++ )\r
 132         {\r
 133             if( image[i * step + j] != 0 )\r
 134                 for( n = 0; n < numangle; n++ )\r
 135                 {\r
 136                     r = cvRound( j * tabCos[n] + i * tabSin[n] );\r
 137                     r += (numrho - 1) / 2;\r
 138                     accum[(n+1) * (numrho+2) + r+1]++;\r
 139                 }\r
 140         }\r
 141 \r
 142     // stage 2. find local maximums\r
 143     for( r = 0; r < numrho; r++ )\r
 144         for( n = 0; n < numangle; n++ )\r
 145         {\r
 146             int base = (n+1) * (numrho+2) + r+1;\r
 147             if( accum[base] > threshold &&\r
 148                 accum[base] > accum[base - 1] && accum[base] >= accum[base + 1] &&\r
 149                 accum[base] > accum[base - numrho - 2] && accum[base] >= accum[base + numrho + 2] )\r
 150                 sort_buf[total++] = base;\r
 151         }\r
 152 \r
 153     // stage 3. sort the detected lines by accumulator value\r
 154     icvHoughSortDescent32s( sort_buf, total, accum );\r
 155 \r
 156     // stage 4. store the first min(total,linesMax) lines to the output buffer\r
 157     linesMax = MIN(linesMax, total);\r
 158     scale = 1./(numrho+2);\r
 159     for( i = 0; i < linesMax; i++ )\r
 160     {\r
 161         CvLinePolar line;\r
 162         int idx = sort_buf[i];\r
 163         int n = cvFloor(idx*scale) - 1;\r
 164         int r = idx - (n+1)*(numrho+2) - 1;\r
 165         line.rho = (r - (numrho - 1)*0.5f) * rho;\r
 166         line.angle = n * theta;\r
 167         cvSeqPush( lines, &line );\r
 168     }\r
 169 \r
 170     __END__;\r
 171 \r
 172     cvFree( &sort_buf );\r
 173     cvFree( &tabSin );\r
 174     cvFree( &tabCos );\r
 175     cvFree( &accum );\r
 176 }\r
 177 \r
 178 \r
 179 /****************************************************************************************\\r
 180 *                     Multi-Scale variant of Classical Hough Transform                   *\r
 181 \****************************************************************************************/\r
 182 \r
 183 #if defined _MSC_VER && _MSC_VER >= 1200\r
 184 #pragma warning( disable: 4714 )\r
 185 #endif\r
 186 \r
 187 //DECLARE_AND_IMPLEMENT_LIST( _index, h_ );\r
 188 IMPLEMENT_LIST( _index, h_ )\r
 189 \r
 190 static void\r
 191 icvHoughLinesSDiv( const CvMat* img,\r
 192                    float rho, float theta, int threshold,\r
 193                    int srn, int stn,\r
 194                    CvSeq* lines, int linesMax )\r
 195 {\r
 196     uchar *caccum = 0;\r
 197     uchar *buffer = 0;\r
 198     float *sinTable = 0;\r
 199     int *x = 0;\r
 200     int *y = 0;\r
 201     _CVLIST *list = 0;\r
 202 \r
 203     CV_FUNCNAME( "icvHoughLinesSDiv" );\r
 204 \r
 205     __BEGIN__;\r
 206 \r
 207 #define _POINT(row, column)\\r
 208     (image_src[(row)*step+(column)])\r
 209 \r
 210     uchar *mcaccum = 0;\r
 211     int rn, tn;                 /* number of rho and theta discrete values */\r
 212     int index, i;\r
 213     int ri, ti, ti1, ti0;\r
 214     int row, col;\r
 215     float r, t;                 /* Current rho and theta */\r
 216     float rv;                   /* Some temporary rho value */\r
 217     float irho;\r
 218     float itheta;\r
 219     float srho, stheta;\r
 220     float isrho, istheta;\r
 221 \r
 222     const uchar* image_src;\r
 223     int w, h, step;\r
 224     int fn = 0;\r
 225     float xc, yc;\r
 226 \r
 227     const float d2r = (float)(Pi / 180);\r
 228     int sfn = srn * stn;\r
 229     int fi;\r
 230     int count;\r
 231     int cmax = 0;\r
 232 \r
 233     CVPOS pos;\r
 234     _index *pindex;\r
 235     _index vi;\r
 236 \r
 237     CV_ASSERT( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1 );\r
 238     CV_ASSERT( linesMax > 0 && rho > 0 && theta > 0 );\r
 239 \r
 240     threshold = MIN( threshold, 255 );\r
 241 \r
 242     image_src = img->data.ptr;\r
 243     step = img->step;\r
 244     w = img->cols;\r
 245     h = img->rows;\r
 246 \r
 247     irho = 1 / rho;\r
 248     itheta = 1 / theta;\r
 249     srho = rho / srn;\r
 250     stheta = theta / stn;\r
 251     isrho = 1 / srho;\r
 252     istheta = 1 / stheta;\r
 253 \r
 254     rn = cvFloor( sqrt( (double)w * w + (double)h * h ) * irho );\r
 255     tn = cvFloor( 2 * Pi * itheta );\r
 256 \r
 257     list = h_create_list__index( linesMax < 1000 ? linesMax : 1000 );\r
 258     vi.value = threshold;\r
 259     vi.rho = -1;\r
 260     h_add_head__index( list, &vi );\r
 261 \r
 262     /* Precalculating sin */\r
 263     CV_CALL( sinTable = (float*)cvAlloc( 5 * tn * stn * sizeof( float )));\r
 264 \r
 265     for( index = 0; index < 5 * tn * stn; index++ )\r
 266     {\r
 267         sinTable[index] = (float)cos( stheta * index * 0.2f );\r
 268     }\r
 269 \r
 270     CV_CALL( caccum = (uchar*)cvAlloc( rn * tn * sizeof( caccum[0] )));\r
 271     memset( caccum, 0, rn * tn * sizeof( caccum[0] ));\r
 272 \r
 273     /* Counting all feature pixels */\r
 274     for( row = 0; row < h; row++ )\r
 275         for( col = 0; col < w; col++ )\r
 276             fn += _POINT( row, col ) != 0;\r
 277 \r
 278     CV_CALL( x = (int*)cvAlloc( fn * sizeof(x[0])));\r
 279     CV_CALL( y = (int*)cvAlloc( fn * sizeof(y[0])));\r
 280 \r
 281     /* Full Hough Transform (it's accumulator update part) */\r
 282     fi = 0;\r
 283     for( row = 0; row < h; row++ )\r
 284     {\r
 285         for( col = 0; col < w; col++ )\r
 286         {\r
 287             if( _POINT( row, col ))\r
 288             {\r
 289                 int halftn;\r
 290                 float r0;\r
 291                 float scale_factor;\r
 292                 int iprev = -1;\r
 293                 float phi, phi1;\r
 294                 float theta_it;     /* Value of theta for iterating */\r
 295 \r
 296                 /* Remember the feature point */\r
 297                 x[fi] = col;\r
 298                 y[fi] = row;\r
 299                 fi++;\r
 300 \r
 301                 yc = (float) row + 0.5f;\r
 302                 xc = (float) col + 0.5f;\r
 303 \r
 304                 /* Update the accumulator */\r
 305                 t = (float) fabs( cvFastArctan( yc, xc ) * d2r );\r
 306                 r = (float) sqrt( (double)xc * xc + (double)yc * yc );\r
 307                 r0 = r * irho;\r
 308                 ti0 = cvFloor( (t + Pi / 2) * itheta );\r
 309 \r
 310                 caccum[ti0]++;\r
 311 \r
 312                 theta_it = rho / r;\r
 313                 theta_it = theta_it < theta ? theta_it : theta;\r
 314                 scale_factor = theta_it * itheta;\r
 315                 halftn = cvFloor( Pi / theta_it );\r
 316                 for( ti1 = 1, phi = theta_it - halfPi, phi1 = (theta_it + t) * itheta;\r
 317                      ti1 < halftn; ti1++, phi += theta_it, phi1 += scale_factor )\r
 318                 {\r
 319                     rv = r0 * _cos( phi );\r
 320                     i = cvFloor( rv ) * tn;\r
 321                     i += cvFloor( phi1 );\r
 322                     assert( i >= 0 );\r
 323                     assert( i < rn * tn );\r
 324                     caccum[i] = (uchar) (caccum[i] + ((i ^ iprev) != 0));\r
 325                     iprev = i;\r
 326                     if( cmax < caccum[i] )\r
 327                         cmax = caccum[i];\r
 328                 }\r
 329             }\r
 330         }\r
 331     }\r
 332 \r
 333     /* Starting additional analysis */\r
 334     count = 0;\r
 335     for( ri = 0; ri < rn; ri++ )\r
 336     {\r
 337         for( ti = 0; ti < tn; ti++ )\r
 338         {\r
 339             if( caccum[ri * tn + ti > threshold] )\r
 340             {\r
 341                 count++;\r
 342             }\r
 343         }\r
 344     }\r
 345 \r
 346     if( count * 100 > rn * tn )\r
 347     {\r
 348         icvHoughLinesStandard( img, rho, theta, threshold, lines, linesMax );\r
 349         EXIT;\r
 350     }\r
 351 \r
 352     CV_CALL( buffer = (uchar *) cvAlloc(srn * stn + 2));\r
 353     mcaccum = buffer + 1;\r
 354 \r
 355     count = 0;\r
 356     for( ri = 0; ri < rn; ri++ )\r
 357     {\r
 358         for( ti = 0; ti < tn; ti++ )\r
 359         {\r
 360             if( caccum[ri * tn + ti] > threshold )\r
 361             {\r
 362                 count++;\r
 363                 memset( mcaccum, 0, sfn * sizeof( uchar ));\r
 364 \r
 365                 for( index = 0; index < fn; index++ )\r
 366                 {\r
 367                     int ti2;\r
 368                     float r0;\r
 369 \r
 370                     yc = (float) y[index] + 0.5f;\r
 371                     xc = (float) x[index] + 0.5f;\r
 372 \r
 373                     /* Update the accumulator */\r
 374                     t = (float) fabs( cvFastArctan( yc, xc ) * d2r );\r
 375                     r = (float) sqrt( (double)xc * xc + (double)yc * yc ) * isrho;\r
 376                     ti0 = cvFloor( (t + Pi * 0.5f) * istheta );\r
 377                     ti2 = (ti * stn - ti0) * 5;\r
 378                     r0 = (float) ri *srn;\r
 379 \r
 380                     for( ti1 = 0 /*, phi = ti*theta - Pi/2 - t */ ; ti1 < stn; ti1++, ti2 += 5\r
 381                          /*phi += stheta */  )\r
 382                     {\r
 383                         /*rv = r*_cos(phi) - r0; */\r
 384                         rv = r * sinTable[(int) (abs( ti2 ))] - r0;\r
 385                         i = cvFloor( rv ) * stn + ti1;\r
 386 \r
 387                         i = CV_IMAX( i, -1 );\r
 388                         i = CV_IMIN( i, sfn );\r
 389                         mcaccum[i]++;\r
 390                         assert( i >= -1 );\r
 391                         assert( i <= sfn );\r
 392                     }\r
 393                 }\r
 394 \r
 395                 /* Find peaks in maccum... */\r
 396                 for( index = 0; index < sfn; index++ )\r
 397                 {\r
 398                     i = 0;\r
 399                     pos = h_get_tail_pos__index( list );\r
 400                     if( h_get_prev__index( &pos )->value < mcaccum[index] )\r
 401                     {\r
 402                         vi.value = mcaccum[index];\r
 403                         vi.rho = index / stn * srho + ri * rho;\r
 404                         vi.theta = index % stn * stheta + ti * theta - halfPi;\r
 405                         while( h_is_pos__index( pos ))\r
 406                         {\r
 407                             if( h_get__index( pos )->value > mcaccum[index] )\r
 408                             {\r
 409                                 h_insert_after__index( list, pos, &vi );\r
 410                                 if( h_get_count__index( list ) > linesMax )\r
 411                                 {\r
 412                                     h_remove_tail__index( list );\r
 413                                 }\r
 414                                 break;\r
 415                             }\r
 416                             h_get_prev__index( &pos );\r
 417                         }\r
 418                         if( !h_is_pos__index( pos ))\r
 419                         {\r
 420                             h_add_head__index( list, &vi );\r
 421                             if( h_get_count__index( list ) > linesMax )\r
 422                             {\r
 423                                 h_remove_tail__index( list );\r
 424                             }\r
 425                         }\r
 426                     }\r
 427                 }\r
 428             }\r
 429         }\r
 430     }\r
 431 \r
 432     pos = h_get_head_pos__index( list );\r
 433     if( h_get_count__index( list ) == 1 )\r
 434     {\r
 435         if( h_get__index( pos )->rho < 0 )\r
 436         {\r
 437             h_clear_list__index( list );\r
 438         }\r
 439     }\r
 440     else\r
 441     {\r
 442         while( h_is_pos__index( pos ))\r
 443         {\r
 444             CvLinePolar line;\r
 445             pindex = h_get__index( pos );\r
 446             if( pindex->rho < 0 )\r
 447             {\r
 448                 /* This should be the last element... */\r
 449                 h_get_next__index( &pos );\r
 450                 assert( !h_is_pos__index( pos ));\r
 451                 break;\r
 452             }\r
 453             line.rho = pindex->rho;\r
 454             line.angle = pindex->theta;\r
 455             cvSeqPush( lines, &line );\r
 456 \r
 457             if( lines->total >= linesMax )\r
 458                 EXIT;\r
 459             h_get_next__index( &pos );\r
 460         }\r
 461     }\r
 462 \r
 463     __END__;\r
 464 \r
 465     h_destroy_list__index( list );\r
 466     cvFree( &sinTable );\r
 467     cvFree( &x );\r
 468     cvFree( &y );\r
 469     cvFree( &caccum );\r
 470     cvFree( &buffer );\r
 471 }\r
 472 \r
 473 \r
 474 /****************************************************************************************\\r
 475 *                              Probabilistic Hough Transform                             *\r
 476 \****************************************************************************************/\r
 477 \r
 478 static void\r
 479 icvHoughLinesProbabalistic( CvMat* image,\r
 480                             float rho, float theta, int threshold,\r
 481                             int lineLength, int lineGap,\r
 482                             CvSeq *lines, int linesMax )\r
 483 {\r
 484     cv::Mat accum, mask;\r
 485     cv::vector<float> trigtab;\r
 486     cv::MemStorage storage(cvCreateMemStorage(0));\r
 487 \r
 488     CvSeq* seq;\r
 489     CvSeqWriter writer;\r
 490     int width, height;\r
 491     int numangle, numrho;\r
 492     float ang;\r
 493     int r, n, count;\r
 494     CvPoint pt;\r
 495     float irho = 1 / rho;\r
 496     CvRNG rng = cvRNG(-1);\r
 497     const float* ttab;\r
 498     uchar* mdata0;\r
 499 \r
 500     CV_Assert( CV_IS_MAT(image) && CV_MAT_TYPE(image->type) == CV_8UC1 );\r
 501 \r
 502     width = image->cols;\r
 503     height = image->rows;\r
 504 \r
 505     numangle = cvRound(CV_PI / theta);\r
 506     numrho = cvRound(((width + height) * 2 + 1) / rho);\r
 507 \r
 508     accum.create( numangle, numrho, CV_32SC1 );\r
 509     mask.create( height, width, CV_8UC1 );\r
 510     trigtab.resize(numangle*2);\r
 511     accum = cv::Scalar(0);\r
 512 \r
 513     for( ang = 0, n = 0; n < numangle; ang += theta, n++ )\r
 514     {\r
 515         trigtab[n*2] = (float)(cos(ang) * irho);\r
 516         trigtab[n*2+1] = (float)(sin(ang) * irho);\r
 517     }\r
 518     ttab = &trigtab[0];\r
 519     mdata0 = mask.data;\r
 520 \r
 521     cvStartWriteSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage, &writer );\r
 522 \r
 523     // stage 1. collect non-zero image points\r
 524     for( pt.y = 0, count = 0; pt.y < height; pt.y++ )\r
 525     {\r
 526         const uchar* data = image->data.ptr + pt.y*image->step;\r
 527         uchar* mdata = mdata0 + pt.y*width;\r
 528         for( pt.x = 0; pt.x < width; pt.x++ )\r
 529         {\r
 530             if( data[pt.x] )\r
 531             {\r
 532                 mdata[pt.x] = (uchar)1;\r
 533                 CV_WRITE_SEQ_ELEM( pt, writer );\r
 534             }\r
 535             else\r
 536                 mdata[pt.x] = 0;\r
 537         }\r
 538     }\r
 539 \r
 540     seq = cvEndWriteSeq( &writer );\r
 541     count = seq->total;\r
 542 \r
 543     // stage 2. process all the points in random order\r
 544     for( ; count > 0; count-- )\r
 545     {\r
 546         // choose random point out of the remaining ones\r
 547         int idx = cvRandInt(&rng) % count;\r
 548         int max_val = threshold-1, max_n = 0;\r
 549         CvPoint* pt = (CvPoint*)cvGetSeqElem( seq, idx );\r
 550         CvPoint line_end[2] = {{0,0}, {0,0}};\r
 551         float a, b;\r
 552         int* adata = (int*)accum.data;\r
 553         int i, j, k, x0, y0, dx0, dy0, xflag;\r
 554         int good_line;\r
 555         const int shift = 16;\r
 556 \r
 557         i = pt->y;\r
 558         j = pt->x;\r
 559 \r
 560         // "remove" it by overriding it with the last element\r
 561         *pt = *(CvPoint*)cvGetSeqElem( seq, count-1 );\r
 562 \r
 563         // check if it has been excluded already (i.e. belongs to some other line)\r
 564         if( !mdata0[i*width + j] )\r
 565             continue;\r
 566 \r
 567         // update accumulator, find the most probable line\r
 568         for( n = 0; n < numangle; n++, adata += numrho )\r
 569         {\r
 570             r = cvRound( j * ttab[n*2] + i * ttab[n*2+1] );\r
 571             r += (numrho - 1) / 2;\r
 572             int val = ++adata[r];\r
 573             if( max_val < val )\r
 574             {\r
 575                 max_val = val;\r
 576                 max_n = n;\r
 577             }\r
 578         }\r
 579 \r
 580         // if it is too "weak" candidate, continue with another point\r
 581         if( max_val < threshold )\r
 582             continue;\r
 583 \r
 584         // from the current point walk in each direction\r
 585         // along the found line and extract the line segment\r
 586         a = -ttab[max_n*2+1];\r
 587         b = ttab[max_n*2];\r
 588         x0 = j;\r
 589         y0 = i;\r
 590         if( fabs(a) > fabs(b) )\r
 591         {\r
 592             xflag = 1;\r
 593             dx0 = a > 0 ? 1 : -1;\r
 594             dy0 = cvRound( b*(1 << shift)/fabs(a) );\r
 595             y0 = (y0 << shift) + (1 << (shift-1));\r
 596         }\r
 597         else\r
 598         {\r
 599             xflag = 0;\r
 600             dy0 = b > 0 ? 1 : -1;\r
 601             dx0 = cvRound( a*(1 << shift)/fabs(b) );\r
 602             x0 = (x0 << shift) + (1 << (shift-1));\r
 603         }\r
 604 \r
 605         for( k = 0; k < 2; k++ )\r
 606         {\r
 607             int gap = 0, x = x0, y = y0, dx = dx0, dy = dy0;\r
 608 \r
 609             if( k > 0 )\r
 610                 dx = -dx, dy = -dy;\r
 611 \r
 612             // walk along the line using fixed-point arithmetics,\r
 613             // stop at the image border or in case of too big gap\r
 614             for( ;; x += dx, y += dy )\r
 615             {\r
 616                 uchar* mdata;\r
 617                 int i1, j1;\r
 618 \r
 619                 if( xflag )\r
 620                 {\r
 621                     j1 = x;\r
 622                     i1 = y >> shift;\r
 623                 }\r
 624                 else\r
 625                 {\r
 626                     j1 = x >> shift;\r
 627                     i1 = y;\r
 628                 }\r
 629 \r
 630                 if( j1 < 0 || j1 >= width || i1 < 0 || i1 >= height )\r
 631                     break;\r
 632 \r
 633                 mdata = mdata0 + i1*width + j1;\r
 634 \r
 635                 // for each non-zero point:\r
 636                 //    update line end,\r
 637                 //    clear the mask element\r
 638                 //    reset the gap\r
 639                 if( *mdata )\r
 640                 {\r
 641                     gap = 0;\r
 642                     line_end[k].y = i1;\r
 643                     line_end[k].x = j1;\r
 644                 }\r
 645                 else if( ++gap > lineGap )\r
 646                     break;\r
 647             }\r
 648         }\r
 649 \r
 650         good_line = abs(line_end[1].x - line_end[0].x) >= lineLength ||\r
 651                     abs(line_end[1].y - line_end[0].y) >= lineLength;\r
 652 \r
 653         for( k = 0; k < 2; k++ )\r
 654         {\r
 655             int x = x0, y = y0, dx = dx0, dy = dy0;\r
 656 \r
 657             if( k > 0 )\r
 658                 dx = -dx, dy = -dy;\r
 659 \r
 660             // walk along the line using fixed-point arithmetics,\r
 661             // stop at the image border or in case of too big gap\r
 662             for( ;; x += dx, y += dy )\r
 663             {\r
 664                 uchar* mdata;\r
 665                 int i1, j1;\r
 666 \r
 667                 if( xflag )\r
 668                 {\r
 669                     j1 = x;\r
 670                     i1 = y >> shift;\r
 671                 }\r
 672                 else\r
 673                 {\r
 674                     j1 = x >> shift;\r
 675                     i1 = y;\r
 676                 }\r
 677 \r
 678                 mdata = mdata0 + i1*width + j1;\r
 679 \r
 680                 // for each non-zero point:\r
 681                 //    update line end,\r
 682                 //    clear the mask element\r
 683                 //    reset the gap\r
 684                 if( *mdata )\r
 685                 {\r
 686                     if( good_line )\r
 687                     {\r
 688                         adata = (int*)accum.data;\r
 689                         for( n = 0; n < numangle; n++, adata += numrho )\r
 690                         {\r
 691                             r = cvRound( j1 * ttab[n*2] + i1 * ttab[n*2+1] );\r
 692                             r += (numrho - 1) / 2;\r
 693                             adata[r]--;\r
 694                         }\r
 695                     }\r
 696                     *mdata = 0;\r
 697                 }\r
 698 \r
 699                 if( i1 == line_end[k].y && j1 == line_end[k].x )\r
 700                     break;\r
 701             }\r
 702         }\r
 703 \r
 704         if( good_line )\r
 705         {\r
 706             CvRect lr = { line_end[0].x, line_end[0].y, line_end[1].x, line_end[1].y };\r
 707             cvSeqPush( lines, &lr );\r
 708             if( lines->total >= linesMax )\r
 709                 return;\r
 710         }\r
 711     }\r
 712 }\r
 713 \r
 714 /* Wrapper function for standard hough transform */\r
 715 CV_IMPL CvSeq*\r
 716 cvHoughLines2( CvArr* src_image, void* lineStorage, int method,\r
 717                double rho, double theta, int threshold,\r
 718                double param1, double param2 )\r
 719 {\r
 720     CvSeq* result = 0;\r
 721 \r
 722     CV_FUNCNAME( "cvHoughLines" );\r
 723 \r
 724     __BEGIN__;\r
 725 \r
 726     CvMat stub, *img = (CvMat*)src_image;\r
 727     CvMat* mat = 0;\r
 728     CvSeq* lines = 0;\r
 729     CvSeq lines_header;\r
 730     CvSeqBlock lines_block;\r
 731     int lineType, elemSize;\r
 732     int linesMax = INT_MAX;\r
 733     int iparam1, iparam2;\r
 734 \r
 735     CV_CALL( img = cvGetMat( img, &stub ));\r
 736 \r
 737     if( !CV_IS_MASK_ARR(img))\r
 738         CV_ERROR( CV_StsBadArg, "The source image must be 8-bit, single-channel" );\r
 739 \r
 740     if( !lineStorage )\r
 741         CV_ERROR( CV_StsNullPtr, "NULL destination" );\r
 742 \r
 743     if( rho <= 0 || theta <= 0 || threshold <= 0 )\r
 744         CV_ERROR( CV_StsOutOfRange, "rho, theta and threshold must be positive" );\r
 745 \r
 746     if( method != CV_HOUGH_PROBABILISTIC )\r
 747     {\r
 748         lineType = CV_32FC2;\r
 749         elemSize = sizeof(float)*2;\r
 750     }\r
 751     else\r
 752     {\r
 753         lineType = CV_32SC4;\r
 754         elemSize = sizeof(int)*4;\r
 755     }\r
 756 \r
 757     if( CV_IS_STORAGE( lineStorage ))\r
 758     {\r
 759         CV_CALL( lines = cvCreateSeq( lineType, sizeof(CvSeq), elemSize, (CvMemStorage*)lineStorage ));\r
 760     }\r
 761     else if( CV_IS_MAT( lineStorage ))\r
 762     {\r
 763         mat = (CvMat*)lineStorage;\r
 764 \r
 765         if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) )\r
 766             CV_ERROR( CV_StsBadArg,\r
 767             "The destination matrix should be continuous and have a single row or a single column" );\r
 768 \r
 769         if( CV_MAT_TYPE( mat->type ) != lineType )\r
 770             CV_ERROR( CV_StsBadArg,\r
 771             "The destination matrix data type is inappropriate, see the manual" );\r
 772 \r
 773         CV_CALL( lines = cvMakeSeqHeaderForArray( lineType, sizeof(CvSeq), elemSize, mat->data.ptr,\r
 774                                                   mat->rows + mat->cols - 1, &lines_header, &lines_block ));\r
 775         linesMax = lines->total;\r
 776         CV_CALL( cvClearSeq( lines ));\r
 777     }\r
 778     else\r
 779     {\r
 780         CV_ERROR( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );\r
 781     }\r
 782 \r
 783     iparam1 = cvRound(param1);\r
 784     iparam2 = cvRound(param2);\r
 785 \r
 786     switch( method )\r
 787     {\r
 788     case CV_HOUGH_STANDARD:\r
 789           CV_CALL( icvHoughLinesStandard( img, (float)rho,\r
 790                 (float)theta, threshold, lines, linesMax ));\r
 791           break;\r
 792     case CV_HOUGH_MULTI_SCALE:\r
 793           CV_CALL( icvHoughLinesSDiv( img, (float)rho, (float)theta,\r
 794                 threshold, iparam1, iparam2, lines, linesMax ));\r
 795           break;\r
 796     case CV_HOUGH_PROBABILISTIC:\r
 797           CV_CALL( icvHoughLinesProbabalistic( img, (float)rho, (float)theta,\r
 798                 threshold, iparam1, iparam2, lines, linesMax ));\r
 799           break;\r
 800     default:\r
 801         CV_ERROR( CV_StsBadArg, "Unrecognized method id" );\r
 802     }\r
 803 \r
 804     if( mat )\r
 805     {\r
 806         if( mat->cols > mat->rows )\r
 807             mat->cols = lines->total;\r
 808         else\r
 809             mat->rows = lines->total;\r
 810     }\r
 811     else\r
 812     {\r
 813         result = lines;\r
 814     }\r
 815 \r
 816     __END__;\r
 817 \r
 818     return result;\r
 819 }\r
 820 \r
 821 \r
 822 /****************************************************************************************\\r
 823 *                                     Circle Detection                                   *\r
 824 \****************************************************************************************/\r
 825 \r
 826 static void\r
 827 icvHoughCirclesGradient( CvMat* img, float dp, float min_dist,\r
 828                          int min_radius, int max_radius,\r
 829                          int canny_threshold, int acc_threshold,\r
 830                          CvSeq* circles, int circles_max )\r
 831 {\r
 832     const int SHIFT = 10, ONE = 1 << SHIFT, R_THRESH = 30;\r
 833     CvMat *dx = 0, *dy = 0;\r
 834     CvMat *edges = 0;\r
 835     CvMat *accum = 0;\r
 836     int* sort_buf = 0;\r
 837     CvMat* dist_buf = 0;\r
 838     CvMemStorage* storage = 0;\r
 839 \r
 840     CV_FUNCNAME( "icvHoughCirclesGradient" );\r
 841 \r
 842     __BEGIN__;\r
 843 \r
 844     int x, y, i, j, center_count, nz_count;\r
 845     int rows, cols, arows, acols;\r
 846     int astep, *adata;\r
 847     float* ddata;\r
 848     CvSeq *nz, *centers;\r
 849     float idp, dr;\r
 850     CvSeqReader reader;\r
 851 \r
 852     CV_CALL( edges = cvCreateMat( img->rows, img->cols, CV_8UC1 ));\r
 853     CV_CALL( cvCanny( img, edges, MAX(canny_threshold/2,1), canny_threshold, 3 ));\r
 854 \r
 855     CV_CALL( dx = cvCreateMat( img->rows, img->cols, CV_16SC1 ));\r
 856     CV_CALL( dy = cvCreateMat( img->rows, img->cols, CV_16SC1 ));\r
 857     CV_CALL( cvSobel( img, dx, 1, 0, 3 ));\r
 858     CV_CALL( cvSobel( img, dy, 0, 1, 3 ));\r
 859 \r
 860     if( dp < 1.f )\r
 861         dp = 1.f;\r
 862     idp = 1.f/dp;\r
 863     CV_CALL( accum = cvCreateMat( cvCeil(img->rows*idp)+2, cvCeil(img->cols*idp)+2, CV_32SC1 ));\r
 864     CV_CALL( cvZero(accum));\r
 865 \r
 866     CV_CALL( storage = cvCreateMemStorage() );\r
 867     CV_CALL( nz = cvCreateSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage ));\r
 868     CV_CALL( centers = cvCreateSeq( CV_32SC1, sizeof(CvSeq), sizeof(int), storage ));\r
 869 \r
 870     rows = img->rows;\r
 871     cols = img->cols;\r
 872     arows = accum->rows - 2;\r
 873     acols = accum->cols - 2;\r
 874     adata = accum->data.i;\r
 875     astep = accum->step/sizeof(adata[0]);\r
 876 \r
 877     for( y = 0; y < rows; y++ )\r
 878     {\r
 879         const uchar* edges_row = edges->data.ptr + y*edges->step;\r
 880         const short* dx_row = (const short*)(dx->data.ptr + y*dx->step);\r
 881         const short* dy_row = (const short*)(dy->data.ptr + y*dy->step);\r
 882 \r
 883         for( x = 0; x < cols; x++ )\r
 884         {\r
 885             float vx, vy;\r
 886             int sx, sy, x0, y0, x1, y1, r, k;\r
 887             CvPoint pt;\r
 888 \r
 889             vx = dx_row[x];\r
 890             vy = dy_row[x];\r
 891 \r
 892             if( !edges_row[x] || (vx == 0 && vy == 0) )\r
 893                 continue;\r
 894 \r
 895             float mag = sqrt(vx*vx+vy*vy);\r
 896             assert( mag >= 1 );\r
 897             sx = cvRound((vx*idp)*ONE/mag);\r
 898             sy = cvRound((vy*idp)*ONE/mag);\r
 899 \r
 900             x0 = cvRound((x*idp)*ONE);\r
 901             y0 = cvRound((y*idp)*ONE);\r
 902 \r
 903             for( k = 0; k < 2; k++ )\r
 904             {\r
 905                 x1 = x0 + min_radius * sx;\r
 906                 y1 = y0 + min_radius * sy;\r
 907 \r
 908                 for( r = min_radius; r <= max_radius; x1 += sx, y1 += sy, r++ )\r
 909                 {\r
 910                     int x2 = x1 >> SHIFT, y2 = y1 >> SHIFT;\r
 911                     if( (unsigned)x2 >= (unsigned)acols ||\r
 912                         (unsigned)y2 >= (unsigned)arows )\r
 913                         break;\r
 914                     adata[y2*astep + x2]++;\r
 915                 }\r
 916 \r
 917                 sx = -sx; sy = -sy;\r
 918             }\r
 919 \r
 920             pt.x = x; pt.y = y;\r
 921             cvSeqPush( nz, &pt );\r
 922         }\r
 923     }\r
 924 \r
 925     nz_count = nz->total;\r
 926     if( !nz_count )\r
 927         EXIT;\r
 928 \r
 929     for( y = 1; y < arows - 1; y++ )\r
 930     {\r
 931         for( x = 1; x < acols - 1; x++ )\r
 932         {\r
 933             int base = y*(acols+2) + x;\r
 934             if( adata[base] > acc_threshold &&\r
 935                 adata[base] > adata[base-1] && adata[base] > adata[base+1] &&\r
 936                 adata[base] > adata[base-acols-2] && adata[base] > adata[base+acols+2] )\r
 937                 cvSeqPush(centers, &base);\r
 938         }\r
 939     }\r
 940 \r
 941     center_count = centers->total;\r
 942     if( !center_count )\r
 943         EXIT;\r
 944 \r
 945     CV_CALL( sort_buf = (int*)cvAlloc( MAX(center_count,nz_count)*sizeof(sort_buf[0]) ));\r
 946     cvCvtSeqToArray( centers, sort_buf );\r
 947 \r
 948     icvHoughSortDescent32s( sort_buf, center_count, adata );\r
 949     cvClearSeq( centers );\r
 950     cvSeqPushMulti( centers, sort_buf, center_count );\r
 951 \r
 952     CV_CALL( dist_buf = cvCreateMat( 1, nz_count, CV_32FC1 ));\r
 953     ddata = dist_buf->data.fl;\r
 954 \r
 955     dr = dp;\r
 956     min_dist = MAX( min_dist, dp );\r
 957     min_dist *= min_dist;\r
 958 \r
 959     for( i = 0; i < centers->total; i++ )\r
 960     {\r
 961         int ofs = *(int*)cvGetSeqElem( centers, i );\r
 962         y = ofs/(acols+2) - 1;\r
 963         x = ofs - (y+1)*(acols+2) - 1;\r
 964         float cx = (float)(x*dp), cy = (float)(y*dp);\r
 965         int start_idx = nz_count - 1;\r
 966         float start_dist, dist_sum;\r
 967         float r_best = 0, c[3];\r
 968         int max_count = R_THRESH;\r
 969 \r
 970         for( j = 0; j < circles->total; j++ )\r
 971         {\r
 972             float* c = (float*)cvGetSeqElem( circles, j );\r
 973             if( (c[0] - cx)*(c[0] - cx) + (c[1] - cy)*(c[1] - cy) < min_dist )\r
 974                 break;\r
 975         }\r
 976 \r
 977         if( j < circles->total )\r
 978             continue;\r
 979 \r
 980         cvStartReadSeq( nz, &reader );\r
 981         for( j = 0; j < nz_count; j++ )\r
 982         {\r
 983             CvPoint pt;\r
 984             float _dx, _dy;\r
 985             CV_READ_SEQ_ELEM( pt, reader );\r
 986             _dx = cx - pt.x; _dy = cy - pt.y;\r
 987             ddata[j] = _dx*_dx + _dy*_dy;\r
 988             sort_buf[j] = j;\r
 989         }\r
 990 \r
 991         cvPow( dist_buf, dist_buf, 0.5 );\r
 992         icvHoughSortDescent32s( sort_buf, nz_count, (int*)ddata );\r
 993 \r
 994         dist_sum = start_dist = ddata[sort_buf[nz_count-1]];\r
 995         for( j = nz_count - 2; j >= 0; j-- )\r
 996         {\r
 997             float d = ddata[sort_buf[j]];\r
 998 \r
 999             if( d > max_radius )\r
1000                 break;\r
1001 \r
1002             if( d - start_dist > dr )\r
1003             {\r
1004                 float r_cur = ddata[sort_buf[(j + start_idx)/2]];\r
1005                 if( (start_idx - j)*r_best >= max_count*r_cur ||\r
1006                     (r_best < FLT_EPSILON && start_idx - j >= max_count) )\r
1007                 {\r
1008                     r_best = r_cur;\r
1009                     max_count = start_idx - j;\r
1010                 }\r
1011                 start_dist = d;\r
1012                 start_idx = j;\r
1013                 dist_sum = 0;\r
1014             }\r
1015             dist_sum += d;\r
1016         }\r
1017 \r
1018         if( max_count > R_THRESH )\r
1019         {\r
1020             c[0] = cx;\r
1021             c[1] = cy;\r
1022             c[2] = (float)r_best;\r
1023             cvSeqPush( circles, c );\r
1024             if( circles->total > circles_max )\r
1025                 EXIT;\r
1026         }\r
1027     }\r
1028 \r
1029     __END__;\r
1030 \r
1031     cvReleaseMat( &dist_buf );\r
1032     cvFree( &sort_buf );\r
1033     cvReleaseMemStorage( &storage );\r
1034     cvReleaseMat( &edges );\r
1035     cvReleaseMat( &dx );\r
1036     cvReleaseMat( &dy );\r
1037     cvReleaseMat( &accum );\r
1038 }\r
1039 \r
1040 CV_IMPL CvSeq*\r
1041 cvHoughCircles( CvArr* src_image, void* circle_storage,\r
1042                 int method, double dp, double min_dist,\r
1043                 double param1, double param2,\r
1044                 int min_radius, int max_radius )\r
1045 {\r
1046     CvSeq* result = 0;\r
1047 \r
1048     CV_FUNCNAME( "cvHoughCircles" );\r
1049 \r
1050     __BEGIN__;\r
1051 \r
1052     CvMat stub, *img = (CvMat*)src_image;\r
1053     CvMat* mat = 0;\r
1054     CvSeq* circles = 0;\r
1055     CvSeq circles_header;\r
1056     CvSeqBlock circles_block;\r
1057     int circles_max = INT_MAX;\r
1058     int canny_threshold = cvRound(param1);\r
1059     int acc_threshold = cvRound(param2);\r
1060 \r
1061     CV_CALL( img = cvGetMat( img, &stub ));\r
1062 \r
1063     if( !CV_IS_MASK_ARR(img))\r
1064         CV_ERROR( CV_StsBadArg, "The source image must be 8-bit, single-channel" );\r
1065 \r
1066     if( !circle_storage )\r
1067         CV_ERROR( CV_StsNullPtr, "NULL destination" );\r
1068 \r
1069     if( dp <= 0 || min_dist <= 0 || canny_threshold <= 0 || acc_threshold <= 0 )\r
1070         CV_ERROR( CV_StsOutOfRange, "dp, min_dist, canny_threshold and acc_threshold must be all positive numbers" );\r
1071 \r
1072     min_radius = MAX( min_radius, 0 );\r
1073     if( max_radius <= 0 )\r
1074         max_radius = MAX( img->rows, img->cols );\r
1075     else if( max_radius <= min_radius )\r
1076         max_radius = min_radius + 2;\r
1077 \r
1078     if( CV_IS_STORAGE( circle_storage ))\r
1079     {\r
1080         CV_CALL( circles = cvCreateSeq( CV_32FC3, sizeof(CvSeq),\r
1081             sizeof(float)*3, (CvMemStorage*)circle_storage ));\r
1082     }\r
1083     else if( CV_IS_MAT( circle_storage ))\r
1084     {\r
1085         mat = (CvMat*)circle_storage;\r
1086 \r
1087         if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) ||\r
1088             CV_MAT_TYPE(mat->type) != CV_32FC3 )\r
1089             CV_ERROR( CV_StsBadArg,\r
1090             "The destination matrix should be continuous and have a single row or a single column" );\r
1091 \r
1092         CV_CALL( circles = cvMakeSeqHeaderForArray( CV_32FC3, sizeof(CvSeq), sizeof(float)*3,\r
1093                 mat->data.ptr, mat->rows + mat->cols - 1, &circles_header, &circles_block ));\r
1094         circles_max = circles->total;\r
1095         CV_CALL( cvClearSeq( circles ));\r
1096     }\r
1097     else\r
1098     {\r
1099         CV_ERROR( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );\r
1100     }\r
1101 \r
1102     switch( method )\r
1103     {\r
1104     case CV_HOUGH_GRADIENT:\r
1105           CV_CALL( icvHoughCirclesGradient( img, (float)dp, (float)min_dist,\r
1106                                     min_radius, max_radius, canny_threshold,\r
1107                                     acc_threshold, circles, circles_max ));\r
1108           break;\r
1109     default:\r
1110         CV_ERROR( CV_StsBadArg, "Unrecognized method id" );\r
1111     }\r
1112 \r
1113     if( mat )\r
1114     {\r
1115         if( mat->cols > mat->rows )\r
1116             mat->cols = circles->total;\r
1117         else\r
1118             mat->rows = circles->total;\r
1119     }\r
1120     else\r
1121         result = circles;\r
1122 \r
1123     __END__;\r
1124 \r
1125     return result;\r
1126 }\r
1127 \r
1128 \r
1129 namespace cv\r
1130 {\r
1131 \r
1132 const int STORAGE_SIZE = 1 << 12;\r
1133 \r
1134 void HoughLines( const Mat& image, vector<Vec2f>& lines,\r
1135                  double rho, double theta, int threshold,\r
1136                  double srn, double stn )\r
1137 {\r
1138     CvMemStorage* storage = cvCreateMemStorage(STORAGE_SIZE);\r
1139     CvMat _image = image;\r
1140     CvSeq* seq = cvHoughLines2( &_image, storage, srn == 0 && stn == 0 ?\r
1141                     CV_HOUGH_STANDARD : CV_HOUGH_MULTI_SCALE,\r
1142                     rho, theta, threshold, srn, stn );\r
1143     Seq<Vec2f>(seq).copyTo(lines);\r
1144 }\r
1145 \r
1146 void HoughLinesP( Mat& image, vector<Vec4i>& lines,\r
1147                   double rho, double theta, int threshold,\r
1148                   double minLineLength, double maxGap )\r
1149 {\r
1150     CvMemStorage* storage = cvCreateMemStorage(STORAGE_SIZE);\r
1151     CvMat _image = image;\r
1152     CvSeq* seq = cvHoughLines2( &_image, storage, CV_HOUGH_PROBABILISTIC,\r
1153                     rho, theta, threshold, minLineLength, maxGap );\r
1154     Seq<Vec4i>(seq).copyTo(lines);\r
1155 }\r
1156 \r
1157 void HoughCircles( const Mat& image, vector<Vec3f>& circles,\r
1158                    int method, double dp, double min_dist,\r
1159                    double param1, double param2,\r
1160                    int minRadius, int maxRadius )\r
1161 {\r
1162     CvMemStorage* storage = cvCreateMemStorage(STORAGE_SIZE);\r
1163     CvMat _image = image;\r
1164     CvSeq* seq = cvHoughCircles( &_image, storage, method,\r
1165                     dp, min_dist, param1, param2, minRadius, maxRadius );\r
1166     Seq<Vec3f>(seq).copyTo(circles);\r
1167 }\r
1168 \r
1169 }\r
1170 \r
1171 /* End of file. */\r