modules/imgproc/src/hough.cpp

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