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