hough.cpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                           License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2000, Intel Corporation, all rights reserved.
00014 // Copyright (C) 2013, OpenCV Foundation, all rights reserved.
00015 // Copyright (C) 2014, Itseez, Inc, all rights reserved.
00016 // Third party copyrights are property of their respective owners.
00017 //
00018 // Redistribution and use in source and binary forms, with or without modification,
00019 // are permitted provided that the following conditions are met:
00020 //
00021 //   * Redistribution's of source code must retain the above copyright notice,
00022 //     this list of conditions and the following disclaimer.
00023 //
00024 //   * Redistribution's in binary form must reproduce the above copyright notice,
00025 //     this list of conditions and the following disclaimer in the documentation
00026 //     and/or other materials provided with the distribution.
00027 //
00028 //   * The name of the copyright holders may not be used to endorse or promote products
00029 //     derived from this software without specific prior written permission.
00030 //
00031 // This software is provided by the copyright holders and contributors "as is" and
00032 // any express or implied warranties, including, but not limited to, the implied
00033 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00034 // In no event shall the Intel Corporation or contributors be liable for any direct,
00035 // indirect, incidental, special, exemplary, or consequential damages
00036 // (including, but not limited to, procurement of substitute goods or services;
00037 // loss of use, data, or profits; or business interruption) however caused
00038 // and on any theory of liability, whether in contract, strict liability,
00039 // or tort (including negligence or otherwise) arising in any way out of
00040 // the use of this software, even if advised of the possibility of such damage.
00041 //
00042 //M*/
00043 
00044 #include "precomp.hpp"
00045 #include "opencl_kernels_imgproc.hpp"
00046 
00047 namespace cv
00048 {
00049 
00050 // Classical Hough Transform
00051 struct LinePolar
00052 {
00053     float rho;
00054     float angle;
00055 };
00056 
00057 
00058 struct hough_cmp_gt
00059 {
00060     hough_cmp_gt(const int* _aux) : aux(_aux) {}
00061     bool operator()(int l1, int l2) const
00062     {
00063         return aux[l1] > aux[l2] || (aux[l1] == aux[l2] && l1 < l2);
00064     }
00065     const int* aux;
00066 };
00067 
00068 
00069 /*
00070 Here image is an input raster;
00071 step is it's step; size characterizes it's ROI;
00072 rho and theta are discretization steps (in pixels and radians correspondingly).
00073 threshold is the minimum number of pixels in the feature for it
00074 to be a candidate for line. lines is the output
00075 array of (rho, theta) pairs. linesMax is the buffer size (number of pairs).
00076 Functions return the actual number of found lines.
00077 */
00078 static void
00079 HoughLinesStandard( const Mat& img, float rho, float theta,
00080                     int threshold, std::vector<Vec2f>& lines, int linesMax,
00081                     double min_theta, double max_theta )
00082 {
00083     int i, j;
00084     float irho = 1 / rho;
00085 
00086     CV_Assert( img.type() == CV_8UC1 );
00087 
00088     const uchar* image = img.ptr();
00089     int step = (int)img.step;
00090     int width = img.cols;
00091     int height = img.rows;
00092 
00093     if (max_theta < min_theta ) {
00094         CV_Error( CV_StsBadArg, "max_theta must be greater than min_theta" );
00095     }
00096     int numangle = cvRound((max_theta - min_theta) / theta);
00097     int numrho = cvRound(((width + height) * 2 + 1) / rho);
00098 
00099 #if defined HAVE_IPP && !defined(HAVE_IPP_ICV_ONLY) && IPP_VERSION_X100 >= 810 && IPP_DISABLE_BLOCK
00100     CV_IPP_CHECK()
00101     {
00102         IppiSize srcSize = { width, height };
00103         IppPointPolar delta = { rho, theta };
00104         IppPointPolar dstRoi[2] = {{(Ipp32f) -(width + height), (Ipp32f) min_theta},{(Ipp32f) (width + height), (Ipp32f) max_theta}};
00105         int bufferSize;
00106         int nz = countNonZero(img);
00107         int ipp_linesMax = std::min(linesMax, nz*numangle/threshold);
00108         int linesCount = 0;
00109         lines.resize(ipp_linesMax);
00110         IppStatus ok = ippiHoughLineGetSize_8u_C1R(srcSize, delta, ipp_linesMax, &bufferSize);
00111         Ipp8u* buffer = ippsMalloc_8u(bufferSize);
00112         if (ok >= 0) ok = ippiHoughLine_Region_8u32f_C1R(image, step, srcSize, (IppPointPolar*) &lines[0], dstRoi, ipp_linesMax, &linesCount, delta, threshold, buffer);
00113         ippsFree(buffer);
00114         if (ok >= 0)
00115         {
00116             lines.resize(linesCount);
00117             CV_IMPL_ADD(CV_IMPL_IPP);
00118             return;
00119         }
00120         lines.clear();
00121         setIppErrorStatus();
00122     }
00123 #endif
00124 
00125     AutoBuffer<int> _accum((numangle+2) * (numrho+2));
00126     std::vector<int> _sort_buf;
00127     AutoBuffer<float> _tabSin(numangle);
00128     AutoBuffer<float> _tabCos(numangle);
00129     int *accum = _accum;
00130     float *tabSin = _tabSin, *tabCos = _tabCos;
00131 
00132     memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );
00133 
00134     float ang = static_cast<float>(min_theta);
00135     for(int n = 0; n < numangle; ang += theta, n++ )
00136     {
00137         tabSin[n] = (float)(sin((double)ang) * irho);
00138         tabCos[n] = (float)(cos((double)ang) * irho);
00139     }
00140 
00141     // stage 1. fill accumulator
00142     for( i = 0; i < height; i++ )
00143         for( j = 0; j < width; j++ )
00144         {
00145             if( image[i * step + j] != 0 )
00146                 for(int n = 0; n < numangle; n++ )
00147                 {
00148                     int r = cvRound( j * tabCos[n] + i * tabSin[n] );
00149                     r += (numrho - 1) / 2;
00150                     accum[(n+1) * (numrho+2) + r+1]++;
00151                 }
00152         }
00153 
00154     // stage 2. find local maximums
00155     for(int r = 0; r < numrho; r++ )
00156         for(int n = 0; n < numangle; n++ )
00157         {
00158             int base = (n+1) * (numrho+2) + r+1;
00159             if( accum[base] > threshold &&
00160                 accum[base] > accum[base - 1] && accum[base] >= accum[base + 1] &&
00161                 accum[base] > accum[base - numrho - 2] && accum[base] >= accum[base + numrho + 2] )
00162                 _sort_buf.push_back(base);
00163         }
00164 
00165     // stage 3. sort the detected lines by accumulator value
00166     std::sort(_sort_buf.begin(), _sort_buf.end(), hough_cmp_gt(accum));
00167 
00168     // stage 4. store the first min(total,linesMax) lines to the output buffer
00169     linesMax = std::min(linesMax, (int)_sort_buf.size());
00170     double scale = 1./(numrho+2);
00171     for( i = 0; i < linesMax; i++ )
00172     {
00173         LinePolar line;
00174         int idx = _sort_buf[i];
00175         int n = cvFloor(idx*scale) - 1;
00176         int r = idx - (n+1)*(numrho+2) - 1;
00177         line.rho = (r - (numrho - 1)*0.5f) * rho;
00178         line.angle = static_cast<float>(min_theta) + n * theta;
00179         lines.push_back(Vec2f(line.rho, line.angle));
00180     }
00181 }
00182 
00183 
00184 // Multi-Scale variant of Classical Hough Transform
00185 
00186 struct hough_index
00187 {
00188     hough_index() : value(0), rho(0.f), theta(0.f) {}
00189     hough_index(int _val, float _rho, float _theta)
00190     : value(_val), rho(_rho), theta(_theta) {}
00191 
00192     int value;
00193     float rho, theta;
00194 };
00195 
00196 
00197 static void
00198 HoughLinesSDiv( const Mat& img,
00199                 float rho, float theta, int threshold,
00200                 int srn, int stn,
00201                 std::vector<Vec2f>& lines, int linesMax,
00202                 double min_theta, double max_theta )
00203 {
00204     #define _POINT(row, column)\
00205         (image_src[(row)*step+(column)])
00206 
00207     int index, i;
00208     int ri, ti, ti1, ti0;
00209     int row, col;
00210     float r, t;                 /* Current rho and theta */
00211     float rv;                   /* Some temporary rho value */
00212 
00213     int fn = 0;
00214     float xc, yc;
00215 
00216     const float d2r = (float)(CV_PI / 180);
00217     int sfn = srn * stn;
00218     int fi;
00219     int count;
00220     int cmax = 0;
00221 
00222     std::vector<hough_index> lst;
00223 
00224     CV_Assert( img.type() == CV_8UC1 );
00225     CV_Assert( linesMax > 0 );
00226 
00227     threshold = MIN( threshold, 255 );
00228 
00229     const uchar* image_src = img.ptr();
00230     int step = (int)img.step;
00231     int w = img.cols;
00232     int h = img.rows;
00233 
00234     float irho = 1 / rho;
00235     float itheta = 1 / theta;
00236     float srho = rho / srn;
00237     float stheta = theta / stn;
00238     float isrho = 1 / srho;
00239     float istheta = 1 / stheta;
00240 
00241     int rn = cvFloor( std::sqrt( (double)w * w + (double)h * h ) * irho );
00242     int tn = cvFloor( 2 * CV_PI * itheta );
00243 
00244     lst.push_back(hough_index(threshold, -1.f, 0.f));
00245 
00246     // Precalculate sin table
00247     std::vector<float> _sinTable( 5 * tn * stn );
00248     float* sinTable = &_sinTable[0];
00249 
00250     for( index = 0; index < 5 * tn * stn; index++ )
00251         sinTable[index] = (float)cos( stheta * index * 0.2f );
00252 
00253     std::vector<uchar> _caccum(rn * tn, (uchar)0);
00254     uchar* caccum = &_caccum[0];
00255 
00256     // Counting all feature pixels
00257     for( row = 0; row < h; row++ )
00258         for( col = 0; col < w; col++ )
00259             fn += _POINT( row, col ) != 0;
00260 
00261     std::vector<int> _x(fn), _y(fn);
00262     int* x = &_x[0], *y = &_y[0];
00263 
00264     // Full Hough Transform (it's accumulator update part)
00265     fi = 0;
00266     for( row = 0; row < h; row++ )
00267     {
00268         for( col = 0; col < w; col++ )
00269         {
00270             if( _POINT( row, col ))
00271             {
00272                 int halftn;
00273                 float r0;
00274                 float scale_factor;
00275                 int iprev = -1;
00276                 float phi, phi1;
00277                 float theta_it;     // Value of theta for iterating
00278 
00279                 // Remember the feature point
00280                 x[fi] = col;
00281                 y[fi] = row;
00282                 fi++;
00283 
00284                 yc = (float) row + 0.5f;
00285                 xc = (float) col + 0.5f;
00286 
00287                 /* Update the accumulator */
00288                 t = (float) fabs( cvFastArctan( yc, xc ) * d2r );
00289                 r = (float) std::sqrt( (double)xc * xc + (double)yc * yc );
00290                 r0 = r * irho;
00291                 ti0 = cvFloor( (t + CV_PI*0.5) * itheta );
00292 
00293                 caccum[ti0]++;
00294 
00295                 theta_it = rho / r;
00296                 theta_it = theta_it < theta ? theta_it : theta;
00297                 scale_factor = theta_it * itheta;
00298                 halftn = cvFloor( CV_PI / theta_it );
00299                 for( ti1 = 1, phi = theta_it - (float)(CV_PI*0.5), phi1 = (theta_it + t) * itheta;
00300                      ti1 < halftn; ti1++, phi += theta_it, phi1 += scale_factor )
00301                 {
00302                     rv = r0 * std::cos( phi );
00303                     i = (int)rv * tn;
00304                     i += cvFloor( phi1 );
00305                     assert( i >= 0 );
00306                     assert( i < rn * tn );
00307                     caccum[i] = (uchar) (caccum[i] + ((i ^ iprev) != 0));
00308                     iprev = i;
00309                     if( cmax < caccum[i] )
00310                         cmax = caccum[i];
00311                 }
00312             }
00313         }
00314     }
00315 
00316     // Starting additional analysis
00317     count = 0;
00318     for( ri = 0; ri < rn; ri++ )
00319     {
00320         for( ti = 0; ti < tn; ti++ )
00321         {
00322             if( caccum[ri * tn + ti] > threshold )
00323                 count++;
00324         }
00325     }
00326 
00327     if( count * 100 > rn * tn )
00328     {
00329         HoughLinesStandard( img, rho, theta, threshold, lines, linesMax, min_theta, max_theta );
00330         return;
00331     }
00332 
00333     std::vector<uchar> _buffer(srn * stn + 2);
00334     uchar* buffer = &_buffer[0];
00335     uchar* mcaccum = buffer + 1;
00336 
00337     count = 0;
00338     for( ri = 0; ri < rn; ri++ )
00339     {
00340         for( ti = 0; ti < tn; ti++ )
00341         {
00342             if( caccum[ri * tn + ti] > threshold )
00343             {
00344                 count++;
00345                 memset( mcaccum, 0, sfn * sizeof( uchar ));
00346 
00347                 for( index = 0; index < fn; index++ )
00348                 {
00349                     int ti2;
00350                     float r0;
00351 
00352                     yc = (float) y[index] + 0.5f;
00353                     xc = (float) x[index] + 0.5f;
00354 
00355                     // Update the accumulator
00356                     t = (float) fabs( cvFastArctan( yc, xc ) * d2r );
00357                     r = (float) std::sqrt( (double)xc * xc + (double)yc * yc ) * isrho;
00358                     ti0 = cvFloor( (t + CV_PI * 0.5) * istheta );
00359                     ti2 = (ti * stn - ti0) * 5;
00360                     r0 = (float) ri *srn;
00361 
00362                     for( ti1 = 0; ti1 < stn; ti1++, ti2 += 5 )
00363                     {
00364                         rv = r * sinTable[(int) (std::abs( ti2 ))] - r0;
00365                         i = cvFloor( rv ) * stn + ti1;
00366 
00367                         i = CV_IMAX( i, -1 );
00368                         i = CV_IMIN( i, sfn );
00369                         mcaccum[i]++;
00370                         assert( i >= -1 );
00371                         assert( i <= sfn );
00372                     }
00373                 }
00374 
00375                 // Find peaks in maccum...
00376                 for( index = 0; index < sfn; index++ )
00377                 {
00378                     i = 0;
00379                     int pos = (int)(lst.size() - 1);
00380                     if( pos < 0 || lst[pos].value < mcaccum[index] )
00381                     {
00382                         hough_index vi(mcaccum[index],
00383                                        index / stn * srho + ri * rho,
00384                                        index % stn * stheta + ti * theta - (float)(CV_PI*0.5));
00385                         lst.push_back(vi);
00386                         for( ; pos >= 0; pos-- )
00387                         {
00388                             if( lst[pos].value > vi.value )
00389                                 break;
00390                             lst[pos+1] = lst[pos];
00391                         }
00392                         lst[pos+1] = vi;
00393                         if( (int)lst.size() > linesMax )
00394                             lst.pop_back();
00395                     }
00396                 }
00397             }
00398         }
00399     }
00400 
00401     for( size_t idx = 0; idx < lst.size(); idx++ )
00402     {
00403         if( lst[idx].rho < 0 )
00404             continue;
00405         lines.push_back(Vec2f(lst[idx].rho, lst[idx].theta));
00406     }
00407 }
00408 
00409 
00410 /****************************************************************************************\
00411 *                              Probabilistic Hough Transform                             *
00412 \****************************************************************************************/
00413 
00414 static void
00415 HoughLinesProbabilistic( Mat& image,
00416                          float rho, float theta, int threshold,
00417                          int lineLength, int lineGap,
00418                          std::vector<Vec4i>& lines, int linesMax )
00419 {
00420     Point pt;
00421     float irho = 1 / rho;
00422     RNG rng((uint64)-1);
00423 
00424     CV_Assert( image.type() == CV_8UC1 );
00425 
00426     int width = image.cols;
00427     int height = image.rows;
00428 
00429     int numangle = cvRound(CV_PI / theta);
00430     int numrho = cvRound(((width + height) * 2 + 1) / rho);
00431 
00432 #if defined HAVE_IPP && !defined(HAVE_IPP_ICV_ONLY) && IPP_VERSION_X100 >= 810 && IPP_DISABLE_BLOCK
00433     CV_IPP_CHECK()
00434     {
00435         IppiSize srcSize = { width, height };
00436         IppPointPolar delta = { rho, theta };
00437         IppiHoughProbSpec* pSpec;
00438         int bufferSize, specSize;
00439         int ipp_linesMax = std::min(linesMax, numangle*numrho);
00440         int linesCount = 0;
00441         lines.resize(ipp_linesMax);
00442         IppStatus ok = ippiHoughProbLineGetSize_8u_C1R(srcSize, delta, &specSize, &bufferSize);
00443         Ipp8u* buffer = ippsMalloc_8u(bufferSize);
00444         pSpec = (IppiHoughProbSpec*) malloc(specSize);
00445         if (ok >= 0) ok = ippiHoughProbLineInit_8u32f_C1R(srcSize, delta, ippAlgHintNone, pSpec);
00446         if (ok >= 0) ok = ippiHoughProbLine_8u32f_C1R(image.data, image.step, srcSize, threshold, lineLength, lineGap, (IppiPoint*) &lines[0], ipp_linesMax, &linesCount, buffer, pSpec);
00447 
00448         free(pSpec);
00449         ippsFree(buffer);
00450         if (ok >= 0)
00451         {
00452             lines.resize(linesCount);
00453             CV_IMPL_ADD(CV_IMPL_IPP);
00454             return;
00455         }
00456         lines.clear();
00457         setIppErrorStatus();
00458     }
00459 #endif
00460 
00461     Mat accum = Mat::zeros( numangle, numrho, CV_32SC1 );
00462     Mat mask( height, width, CV_8UC1 );
00463     std::vector<float> trigtab(numangle*2);
00464 
00465     for( int n = 0; n < numangle; n++ )
00466     {
00467         trigtab[n*2] = (float)(cos((double)n*theta) * irho);
00468         trigtab[n*2+1] = (float)(sin((double)n*theta) * irho);
00469     }
00470     const float* ttab = &trigtab[0];
00471     uchar* mdata0 = mask.ptr();
00472     std::vector<Point> nzloc;
00473 
00474     // stage 1. collect non-zero image points
00475     for( pt.y = 0; pt.y < height; pt.y++ )
00476     {
00477         const uchar* data = image.ptr(pt.y);
00478         uchar* mdata = mask.ptr(pt.y);
00479         for( pt.x = 0; pt.x < width; pt.x++ )
00480         {
00481             if( data[pt.x] )
00482             {
00483                 mdata[pt.x] = (uchar)1;
00484                 nzloc.push_back(pt);
00485             }
00486             else
00487                 mdata[pt.x] = 0;
00488         }
00489     }
00490 
00491     int count = (int)nzloc.size();
00492 
00493     // stage 2. process all the points in random order
00494     for( ; count > 0; count-- )
00495     {
00496         // choose random point out of the remaining ones
00497         int idx = rng.uniform(0, count);
00498         int max_val = threshold-1, max_n = 0;
00499         Point point = nzloc[idx];
00500         Point line_end[2];
00501         float a, b;
00502         int* adata = accum.ptr<int>();
00503         int i = point.y, j = point.x, k, x0, y0, dx0, dy0, xflag;
00504         int good_line;
00505         const int shift = 16;
00506 
00507         // "remove" it by overriding it with the last element
00508         nzloc[idx] = nzloc[count-1];
00509 
00510         // check if it has been excluded already (i.e. belongs to some other line)
00511         if( !mdata0[i*width + j] )
00512             continue;
00513 
00514         // update accumulator, find the most probable line
00515         for( int n = 0; n < numangle; n++, adata += numrho )
00516         {
00517             int r = cvRound( j * ttab[n*2] + i * ttab[n*2+1] );
00518             r += (numrho - 1) / 2;
00519             int val = ++adata[r];
00520             if( max_val < val )
00521             {
00522                 max_val = val;
00523                 max_n = n;
00524             }
00525         }
00526 
00527         // if it is too "weak" candidate, continue with another point
00528         if( max_val < threshold )
00529             continue;
00530 
00531         // from the current point walk in each direction
00532         // along the found line and extract the line segment
00533         a = -ttab[max_n*2+1];
00534         b = ttab[max_n*2];
00535         x0 = j;
00536         y0 = i;
00537         if( fabs(a) > fabs(b) )
00538         {
00539             xflag = 1;
00540             dx0 = a > 0 ? 1 : -1;
00541             dy0 = cvRound( b*(1 << shift)/fabs(a) );
00542             y0 = (y0 << shift) + (1 << (shift-1));
00543         }
00544         else
00545         {
00546             xflag = 0;
00547             dy0 = b > 0 ? 1 : -1;
00548             dx0 = cvRound( a*(1 << shift)/fabs(b) );
00549             x0 = (x0 << shift) + (1 << (shift-1));
00550         }
00551 
00552         for( k = 0; k < 2; k++ )
00553         {
00554             int gap = 0, x = x0, y = y0, dx = dx0, dy = dy0;
00555 
00556             if( k > 0 )
00557                 dx = -dx, dy = -dy;
00558 
00559             // walk along the line using fixed-point arithmetics,
00560             // stop at the image border or in case of too big gap
00561             for( ;; x += dx, y += dy )
00562             {
00563                 uchar* mdata;
00564                 int i1, j1;
00565 
00566                 if( xflag )
00567                 {
00568                     j1 = x;
00569                     i1 = y >> shift;
00570                 }
00571                 else
00572                 {
00573                     j1 = x >> shift;
00574                     i1 = y;
00575                 }
00576 
00577                 if( j1 < 0 || j1 >= width || i1 < 0 || i1 >= height )
00578                     break;
00579 
00580                 mdata = mdata0 + i1*width + j1;
00581 
00582                 // for each non-zero point:
00583                 //    update line end,
00584                 //    clear the mask element
00585                 //    reset the gap
00586                 if( *mdata )
00587                 {
00588                     gap = 0;
00589                     line_end[k].y = i1;
00590                     line_end[k].x = j1;
00591                 }
00592                 else if( ++gap > lineGap )
00593                     break;
00594             }
00595         }
00596 
00597         good_line = std::abs(line_end[1].x - line_end[0].x) >= lineLength ||
00598                     std::abs(line_end[1].y - line_end[0].y) >= lineLength;
00599 
00600         for( k = 0; k < 2; k++ )
00601         {
00602             int x = x0, y = y0, dx = dx0, dy = dy0;
00603 
00604             if( k > 0 )
00605                 dx = -dx, dy = -dy;
00606 
00607             // walk along the line using fixed-point arithmetics,
00608             // stop at the image border or in case of too big gap
00609             for( ;; x += dx, y += dy )
00610             {
00611                 uchar* mdata;
00612                 int i1, j1;
00613 
00614                 if( xflag )
00615                 {
00616                     j1 = x;
00617                     i1 = y >> shift;
00618                 }
00619                 else
00620                 {
00621                     j1 = x >> shift;
00622                     i1 = y;
00623                 }
00624 
00625                 mdata = mdata0 + i1*width + j1;
00626 
00627                 // for each non-zero point:
00628                 //    update line end,
00629                 //    clear the mask element
00630                 //    reset the gap
00631                 if( *mdata )
00632                 {
00633                     if( good_line )
00634                     {
00635                         adata = accum.ptr<int>();
00636                         for( int n = 0; n < numangle; n++, adata += numrho )
00637                         {
00638                             int r = cvRound( j1 * ttab[n*2] + i1 * ttab[n*2+1] );
00639                             r += (numrho - 1) / 2;
00640                             adata[r]--;
00641                         }
00642                     }
00643                     *mdata = 0;
00644                 }
00645 
00646                 if( i1 == line_end[k].y && j1 == line_end[k].x )
00647                     break;
00648             }
00649         }
00650 
00651         if( good_line )
00652         {
00653             Vec4i lr(line_end[0].x, line_end[0].y, line_end[1].x, line_end[1].y);
00654             lines.push_back(lr);
00655             if( (int)lines.size() >= linesMax )
00656                 return;
00657         }
00658     }
00659 }
00660 
00661 #ifdef HAVE_OPENCL
00662 
00663 #define OCL_MAX_LINES 4096
00664 
00665 static bool ocl_makePointsList(InputArray _src, OutputArray _pointsList, InputOutputArray _counters)
00666 {
00667     UMat src = _src.getUMat();
00668     _pointsList.create(1, (int) src.total(), CV_32SC1);
00669     UMat pointsList = _pointsList.getUMat();
00670     UMat counters = _counters.getUMat();
00671     ocl::Device dev = ocl::Device::getDefault();
00672 
00673     const int pixPerWI = 16;
00674     int workgroup_size = min((int) dev.maxWorkGroupSize(), (src.cols + pixPerWI - 1)/pixPerWI);
00675     ocl::Kernel pointListKernel("make_point_list", ocl::imgproc::hough_lines_oclsrc,
00676                                 format("-D MAKE_POINTS_LIST -D GROUP_SIZE=%d -D LOCAL_SIZE=%d", workgroup_size, src.cols));
00677     if (pointListKernel.empty())
00678         return false;
00679 
00680     pointListKernel.args(ocl::KernelArg::ReadOnly(src), ocl::KernelArg::WriteOnlyNoSize(pointsList),
00681                          ocl::KernelArg::PtrWriteOnly(counters));
00682 
00683     size_t localThreads[2]  = { (size_t)workgroup_size, 1 };
00684     size_t globalThreads[2] = { (size_t)workgroup_size, (size_t)src.rows };
00685 
00686     return pointListKernel.run(2, globalThreads, localThreads, false);
00687 }
00688 
00689 static bool ocl_fillAccum(InputArray _pointsList, OutputArray _accum, int total_points, double rho, double theta, int numrho, int numangle)
00690 {
00691     UMat pointsList = _pointsList.getUMat();
00692     _accum.create(numangle + 2, numrho + 2, CV_32SC1);
00693     UMat accum = _accum.getUMat();
00694     ocl::Device dev = ocl::Device::getDefault();
00695 
00696     float irho = (float) (1 / rho);
00697     int workgroup_size = min((int) dev.maxWorkGroupSize(), total_points);
00698 
00699     ocl::Kernel fillAccumKernel;
00700     size_t localThreads[2];
00701     size_t globalThreads[2];
00702 
00703     size_t local_memory_needed = (numrho + 2)*sizeof(int);
00704     if (local_memory_needed > dev.localMemSize())
00705     {
00706         accum.setTo(Scalar::all(0));
00707         fillAccumKernel.create("fill_accum_global", ocl::imgproc::hough_lines_oclsrc,
00708                                 format("-D FILL_ACCUM_GLOBAL"));
00709         if (fillAccumKernel.empty())
00710             return false;
00711         globalThreads[0] = workgroup_size; globalThreads[1] = numangle;
00712         fillAccumKernel.args(ocl::KernelArg::ReadOnlyNoSize(pointsList), ocl::KernelArg::WriteOnlyNoSize(accum),
00713                         total_points, irho, (float) theta, numrho, numangle);
00714         return fillAccumKernel.run(2, globalThreads, NULL, false);
00715     }
00716     else
00717     {
00718         fillAccumKernel.create("fill_accum_local", ocl::imgproc::hough_lines_oclsrc,
00719                                 format("-D FILL_ACCUM_LOCAL -D LOCAL_SIZE=%d -D BUFFER_SIZE=%d", workgroup_size, numrho + 2));
00720         if (fillAccumKernel.empty())
00721             return false;
00722         localThreads[0] = workgroup_size; localThreads[1] = 1;
00723         globalThreads[0] = workgroup_size; globalThreads[1] = numangle+2;
00724         fillAccumKernel.args(ocl::KernelArg::ReadOnlyNoSize(pointsList), ocl::KernelArg::WriteOnlyNoSize(accum),
00725                         total_points, irho, (float) theta, numrho, numangle);
00726         return fillAccumKernel.run(2, globalThreads, localThreads, false);
00727     }
00728 }
00729 
00730 static bool ocl_HoughLines(InputArray _src, OutputArray _lines, double rho, double theta, int threshold,
00731                            double min_theta, double max_theta)
00732 {
00733     CV_Assert(_src.type() == CV_8UC1);
00734 
00735     if (max_theta < 0 || max_theta > CV_PI ) {
00736         CV_Error( CV_StsBadArg, "max_theta must fall between 0 and pi" );
00737     }
00738     if (min_theta < 0 || min_theta > max_theta ) {
00739         CV_Error( CV_StsBadArg, "min_theta must fall between 0 and max_theta" );
00740     }
00741     if (!(rho > 0 && theta > 0)) {
00742         CV_Error( CV_StsBadArg, "rho and theta must be greater 0" );
00743     }
00744 
00745     UMat src = _src.getUMat();
00746     int numangle = cvRound((max_theta - min_theta) / theta);
00747     int numrho = cvRound(((src.cols + src.rows) * 2 + 1) / rho);
00748 
00749     UMat pointsList;
00750     UMat counters(1, 2, CV_32SC1, Scalar::all(0));
00751 
00752     if (!ocl_makePointsList(src, pointsList, counters))
00753         return false;
00754 
00755     int total_points = counters.getMat(ACCESS_READ).at<int>(0, 0);
00756     if (total_points <= 0)
00757     {
00758         _lines.assign(UMat(0,0,CV_32FC2));
00759         return true;
00760     }
00761 
00762     UMat accum;
00763     if (!ocl_fillAccum(pointsList, accum, total_points, rho, theta, numrho, numangle))
00764         return false;
00765 
00766     const int pixPerWI = 8;
00767     ocl::Kernel getLinesKernel("get_lines", ocl::imgproc::hough_lines_oclsrc,
00768                                format("-D GET_LINES"));
00769     if (getLinesKernel.empty())
00770         return false;
00771 
00772     int linesMax = threshold > 0 ? min(total_points*numangle/threshold, OCL_MAX_LINES) : OCL_MAX_LINES;
00773     UMat lines(linesMax, 1, CV_32FC2);
00774 
00775     getLinesKernel.args(ocl::KernelArg::ReadOnly(accum), ocl::KernelArg::WriteOnlyNoSize(lines),
00776                         ocl::KernelArg::PtrWriteOnly(counters), linesMax, threshold, (float) rho, (float) theta);
00777 
00778     size_t globalThreads[2] = { ((size_t)numrho + pixPerWI - 1)/pixPerWI, (size_t)numangle };
00779     if (!getLinesKernel.run(2, globalThreads, NULL, false))
00780         return false;
00781 
00782     int total_lines = min(counters.getMat(ACCESS_READ).at<int>(0, 1), linesMax);
00783     if (total_lines > 0)
00784         _lines.assign(lines.rowRange(Range(0, total_lines)));
00785     else
00786         _lines.assign(UMat(0,0,CV_32FC2));
00787     return true;
00788 }
00789 
00790 static bool ocl_HoughLinesP(InputArray _src, OutputArray _lines, double rho, double theta, int threshold,
00791                            double minLineLength, double maxGap)
00792 {
00793     CV_Assert(_src.type() == CV_8UC1);
00794 
00795     if (!(rho > 0 && theta > 0)) {
00796         CV_Error( CV_StsBadArg, "rho and theta must be greater 0" );
00797     }
00798 
00799     UMat src = _src.getUMat();
00800     int numangle = cvRound(CV_PI / theta);
00801     int numrho = cvRound(((src.cols + src.rows) * 2 + 1) / rho);
00802 
00803     UMat pointsList;
00804     UMat counters(1, 2, CV_32SC1, Scalar::all(0));
00805 
00806     if (!ocl_makePointsList(src, pointsList, counters))
00807         return false;
00808 
00809     int total_points = counters.getMat(ACCESS_READ).at<int>(0, 0);
00810     if (total_points <= 0)
00811     {
00812         _lines.assign(UMat(0,0,CV_32SC4));
00813         return true;
00814     }
00815 
00816     UMat accum;
00817     if (!ocl_fillAccum(pointsList, accum, total_points, rho, theta, numrho, numangle))
00818         return false;
00819 
00820     ocl::Kernel getLinesKernel("get_lines", ocl::imgproc::hough_lines_oclsrc,
00821                                format("-D GET_LINES_PROBABOLISTIC"));
00822     if (getLinesKernel.empty())
00823         return false;
00824 
00825     int linesMax = threshold > 0 ? min(total_points*numangle/threshold, OCL_MAX_LINES) : OCL_MAX_LINES;
00826     UMat lines(linesMax, 1, CV_32SC4);
00827 
00828     getLinesKernel.args(ocl::KernelArg::ReadOnly(accum), ocl::KernelArg::ReadOnly(src),
00829                         ocl::KernelArg::WriteOnlyNoSize(lines), ocl::KernelArg::PtrWriteOnly(counters),
00830                         linesMax, threshold, (int) minLineLength, (int) maxGap, (float) rho, (float) theta);
00831 
00832     size_t globalThreads[2] = { (size_t)numrho, (size_t)numangle };
00833     if (!getLinesKernel.run(2, globalThreads, NULL, false))
00834         return false;
00835 
00836     int total_lines = min(counters.getMat(ACCESS_READ).at<int>(0, 1), linesMax);
00837     if (total_lines > 0)
00838         _lines.assign(lines.rowRange(Range(0, total_lines)));
00839     else
00840         _lines.assign(UMat(0,0,CV_32SC4));
00841 
00842     return true;
00843 }
00844 
00845 #endif /* HAVE_OPENCL */
00846 
00847 }
00848 
00849 void cv::HoughLines( InputArray _image, OutputArray _lines,
00850                     double rho, double theta, int threshold,
00851                     double srn, double stn, double min_theta, double max_theta )
00852 {
00853 #ifdef HAVE_OPENCL
00854     CV_OCL_RUN(srn == 0 && stn == 0 && _image.isUMat() && _lines.isUMat(),
00855                ocl_HoughLines(_image, _lines, rho, theta, threshold, min_theta, max_theta));
00856 #endif
00857 
00858     Mat image = _image.getMat();
00859     std::vector<Vec2f> lines;
00860 
00861     if( srn == 0 && stn == 0 )
00862         HoughLinesStandard(image, (float)rho, (float)theta, threshold, lines, INT_MAX, min_theta, max_theta );
00863     else
00864         HoughLinesSDiv(image, (float)rho, (float)theta, threshold, cvRound(srn), cvRound(stn), lines, INT_MAX, min_theta, max_theta);
00865 
00866     Mat(lines).copyTo(_lines);
00867 }
00868 
00869 
00870 void cv::HoughLinesP(InputArray _image, OutputArray _lines,
00871                      double rho, double theta, int threshold,
00872                      double minLineLength, double maxGap )
00873 {
00874 #ifdef HAVE_OPENCL
00875     CV_OCL_RUN(_image.isUMat() && _lines.isUMat(),
00876                ocl_HoughLinesP(_image, _lines, rho, theta, threshold, minLineLength, maxGap));
00877 #endif
00878 
00879     Mat image = _image.getMat();
00880     std::vector<Vec4i> lines;
00881     HoughLinesProbabilistic(image, (float)rho, (float)theta, threshold, cvRound(minLineLength), cvRound(maxGap), lines, INT_MAX);
00882     Mat(lines).copyTo(_lines);
00883 }
00884 
00885 
00886 
00887 /* Wrapper function for standard hough transform */
00888 CV_IMPL CvSeq*
00889 cvHoughLines2( CvArr* src_image, void* lineStorage, int method,
00890                double rho, double theta, int threshold,
00891                double param1, double param2,
00892                double min_theta, double max_theta )
00893 {
00894     cv::Mat image = cv::cvarrToMat(src_image);
00895     std::vector<cv::Vec2f> l2;
00896     std::vector<cv::Vec4i> l4;
00897     CvSeq* result = 0;
00898 
00899     CvMat* mat = 0;
00900     CvSeq* lines = 0;
00901     CvSeq lines_header;
00902     CvSeqBlock lines_block;
00903     int lineType, elemSize;
00904     int linesMax = INT_MAX;
00905     int iparam1, iparam2;
00906 
00907     if( !lineStorage )
00908         CV_Error( CV_StsNullPtr, "NULL destination" );
00909 
00910     if( rho <= 0 || theta <= 0 || threshold <= 0 )
00911         CV_Error( CV_StsOutOfRange, "rho, theta and threshold must be positive" );
00912 
00913     if( method != CV_HOUGH_PROBABILISTIC )
00914     {
00915         lineType = CV_32FC2;
00916         elemSize = sizeof(float)*2;
00917     }
00918     else
00919     {
00920         lineType = CV_32SC4;
00921         elemSize = sizeof(int)*4;
00922     }
00923 
00924     if( CV_IS_STORAGE( lineStorage ))
00925     {
00926         lines = cvCreateSeq( lineType, sizeof(CvSeq), elemSize, (CvMemStorage*)lineStorage );
00927     }
00928     else if( CV_IS_MAT( lineStorage ))
00929     {
00930         mat = (CvMat*)lineStorage;
00931 
00932         if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) )
00933             CV_Error( CV_StsBadArg,
00934             "The destination matrix should be continuous and have a single row or a single column" );
00935 
00936         if( CV_MAT_TYPE( mat->type ) != lineType )
00937             CV_Error( CV_StsBadArg,
00938             "The destination matrix data type is inappropriate, see the manual" );
00939 
00940         lines = cvMakeSeqHeaderForArray( lineType, sizeof(CvSeq), elemSize, mat->data.ptr,
00941                                          mat->rows + mat->cols - 1, &lines_header, &lines_block );
00942         linesMax = lines->total;
00943         cvClearSeq( lines );
00944     }
00945     else
00946         CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );
00947 
00948     iparam1 = cvRound(param1);
00949     iparam2 = cvRound(param2);
00950 
00951     switch( method )
00952     {
00953     case CV_HOUGH_STANDARD:
00954         HoughLinesStandard( image, (float)rho,
00955                 (float)theta, threshold, l2, linesMax, min_theta, max_theta );
00956         break;
00957     case CV_HOUGH_MULTI_SCALE:
00958         HoughLinesSDiv( image, (float)rho, (float)theta,
00959                 threshold, iparam1, iparam2, l2, linesMax, min_theta, max_theta );
00960         break;
00961     case CV_HOUGH_PROBABILISTIC:
00962         HoughLinesProbabilistic( image, (float)rho, (float)theta,
00963                 threshold, iparam1, iparam2, l4, linesMax );
00964         break;
00965     default:
00966         CV_Error( CV_StsBadArg, "Unrecognized method id" );
00967     }
00968 
00969     int nlines = (int)(l2.size() + l4.size());
00970 
00971     if( mat )
00972     {
00973         if( mat->cols > mat->rows )
00974             mat->cols = nlines;
00975         else
00976             mat->rows = nlines;
00977     }
00978 
00979     if( nlines )
00980     {
00981         cv::Mat lx = method == CV_HOUGH_STANDARD || method == CV_HOUGH_MULTI_SCALE ?
00982             cv::Mat(nlines, 1, CV_32FC2, &l2[0]) : cv::Mat(nlines, 1, CV_32SC4, &l4[0]);
00983 
00984         if( mat )
00985         {
00986             cv::Mat dst(nlines, 1, lx.type(), mat->data.ptr);
00987             lx.copyTo(dst);
00988         }
00989         else
00990         {
00991             cvSeqPushMulti(lines, lx.ptr(), nlines);
00992         }
00993     }
00994 
00995     if( !mat )
00996         result = lines;
00997     return result;
00998 }
00999 
01000 
01001 /****************************************************************************************\
01002 *                                     Circle Detection                                   *
01003 \****************************************************************************************/
01004 
01005 static void
01006 icvHoughCirclesGradient( CvMat* img, float dp, float min_dist,
01007                          int min_radius, int max_radius,
01008                          int canny_threshold, int acc_threshold,
01009                          CvSeq* circles, int circles_max )
01010 {
01011     const int SHIFT = 10, ONE = 1 << SHIFT;
01012     cv::Ptr<CvMat> dx, dy;
01013     cv::Ptr<CvMat> edges, accum, dist_buf;
01014     std::vector<int> sort_buf;
01015     cv::Ptr<CvMemStorage> storage;
01016 
01017     int x, y, i, j, k, center_count, nz_count;
01018     float min_radius2 = (float)min_radius*min_radius;
01019     float max_radius2 = (float)max_radius*max_radius;
01020     int rows, cols, arows, acols;
01021     int astep, *adata;
01022     float* ddata;
01023     CvSeq *nz, *centers;
01024     float idp, dr;
01025     CvSeqReader reader;
01026 
01027     edges.reset(cvCreateMat( img->rows, img->cols, CV_8UC1 ));
01028 
01029     // Use the Canny Edge Detector to detect all the edges in the image.
01030     cvCanny( img, edges, MAX(canny_threshold/2,1), canny_threshold, 3 );
01031 
01032     dx.reset(cvCreateMat( img->rows, img->cols, CV_16SC1 ));
01033     dy.reset(cvCreateMat( img->rows, img->cols, CV_16SC1 ));
01034 
01035     /*Use the Sobel Derivative to compute the local gradient of all the non-zero pixels in the edge image.*/
01036     cvSobel( img, dx, 1, 0, 3 );
01037     cvSobel( img, dy, 0, 1, 3 );
01038 
01039     if( dp < 1.f )
01040         dp = 1.f;
01041     idp = 1.f/dp;
01042     accum.reset(cvCreateMat( cvCeil(img->rows*idp)+2, cvCeil(img->cols*idp)+2, CV_32SC1 ));
01043     cvZero(accum);
01044 
01045     storage.reset(cvCreateMemStorage());
01046     /* Create sequences for the nonzero pixels in the edge image and the centers of circles
01047     which could be detected.*/
01048     nz = cvCreateSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage );
01049     centers = cvCreateSeq( CV_32SC1, sizeof(CvSeq), sizeof(int), storage );
01050 
01051     rows = img->rows;
01052     cols = img->cols;
01053     arows = accum->rows - 2;
01054     acols = accum->cols - 2;
01055     adata = accum->data.i;
01056     astep = accum->step/sizeof(adata[0]);
01057     // Accumulate circle evidence for each edge pixel
01058     for( y = 0; y < rows; y++ )
01059     {
01060         const uchar* edges_row = edges->data.ptr + y*edges->step;
01061         const short* dx_row = (const short*)(dx->data.ptr + y*dx->step);
01062         const short* dy_row = (const short*)(dy->data.ptr + y*dy->step);
01063 
01064         for( x = 0; x < cols; x++ )
01065         {
01066             float vx, vy;
01067             int sx, sy, x0, y0, x1, y1, r;
01068             CvPoint pt;
01069 
01070             vx = dx_row[x];
01071             vy = dy_row[x];
01072 
01073             if( !edges_row[x] || (vx == 0 && vy == 0) )
01074                 continue;
01075 
01076             float mag = std::sqrt(vx*vx+vy*vy);
01077             assert( mag >= 1 );
01078             sx = cvRound((vx*idp)*ONE/mag);
01079             sy = cvRound((vy*idp)*ONE/mag);
01080 
01081             x0 = cvRound((x*idp)*ONE);
01082             y0 = cvRound((y*idp)*ONE);
01083             // Step from min_radius to max_radius in both directions of the gradient
01084             for(int k1 = 0; k1 < 2; k1++ )
01085             {
01086                 x1 = x0 + min_radius * sx;
01087                 y1 = y0 + min_radius * sy;
01088 
01089                 for( r = min_radius; r <= max_radius; x1 += sx, y1 += sy, r++ )
01090                 {
01091                     int x2 = x1 >> SHIFT, y2 = y1 >> SHIFT;
01092                     if( (unsigned)x2 >= (unsigned)acols ||
01093                         (unsigned)y2 >= (unsigned)arows )
01094                         break;
01095                     adata[y2*astep + x2]++;
01096                 }
01097 
01098                 sx = -sx; sy = -sy;
01099             }
01100 
01101             pt.x = x; pt.y = y;
01102             cvSeqPush( nz, &pt );
01103         }
01104     }
01105 
01106     nz_count = nz->total;
01107     if( !nz_count )
01108         return;
01109     //Find possible circle centers
01110     for( y = 1; y < arows - 1; y++ )
01111     {
01112         for( x = 1; x < acols - 1; x++ )
01113         {
01114             int base = y*(acols+2) + x;
01115             if( adata[base] > acc_threshold &&
01116                 adata[base] > adata[base-1] && adata[base] > adata[base+1] &&
01117                 adata[base] > adata[base-acols-2] && adata[base] > adata[base+acols+2] )
01118                 cvSeqPush(centers, &base);
01119         }
01120     }
01121 
01122     center_count = centers->total;
01123     if( !center_count )
01124         return;
01125 
01126     sort_buf.resize( MAX(center_count,nz_count) );
01127     cvCvtSeqToArray( centers, &sort_buf[0] );
01128     /*Sort candidate centers in descending order of their accumulator values, so that the centers
01129     with the most supporting pixels appear first.*/
01130     std::sort(sort_buf.begin(), sort_buf.begin() + center_count, cv::hough_cmp_gt(adata));
01131     cvClearSeq( centers );
01132     cvSeqPushMulti( centers, &sort_buf[0], center_count );
01133 
01134     dist_buf.reset(cvCreateMat( 1, nz_count, CV_32FC1 ));
01135     ddata = dist_buf->data.fl;
01136 
01137     dr = dp;
01138     min_dist = MAX( min_dist, dp );
01139     min_dist *= min_dist;
01140     // For each found possible center
01141     // Estimate radius and check support
01142     for( i = 0; i < centers->total; i++ )
01143     {
01144         int ofs = *(int*)cvGetSeqElem( centers, i );
01145         y = ofs/(acols+2);
01146         x = ofs - (y)*(acols+2);
01147         //Calculate circle's center in pixels
01148         float cx = (float)((x + 0.5f)*dp), cy = (float)(( y + 0.5f )*dp);
01149         float start_dist, dist_sum;
01150         float r_best = 0;
01151         int max_count = 0;
01152         // Check distance with previously detected circles
01153         for( j = 0; j < circles->total; j++ )
01154         {
01155             float* c = (float*)cvGetSeqElem( circles, j );
01156             if( (c[0] - cx)*(c[0] - cx) + (c[1] - cy)*(c[1] - cy) < min_dist )
01157                 break;
01158         }
01159 
01160         if( j < circles->total )
01161             continue;
01162         // Estimate best radius
01163         cvStartReadSeq( nz, &reader );
01164         for( j = k = 0; j < nz_count; j++ )
01165         {
01166             CvPoint pt;
01167             float _dx, _dy, _r2;
01168             CV_READ_SEQ_ELEM( pt, reader );
01169             _dx = cx - pt.x; _dy = cy - pt.y;
01170             _r2 = _dx*_dx + _dy*_dy;
01171             if(min_radius2 <= _r2 && _r2 <= max_radius2 )
01172             {
01173                 ddata[k] = _r2;
01174                 sort_buf[k] = k;
01175                 k++;
01176             }
01177         }
01178 
01179         int nz_count1 = k, start_idx = nz_count1 - 1;
01180         if( nz_count1 == 0 )
01181             continue;
01182         dist_buf->cols = nz_count1;
01183         cvPow( dist_buf, dist_buf, 0.5 );
01184         // Sort non-zero pixels according to their distance from the center.
01185         std::sort(sort_buf.begin(), sort_buf.begin() + nz_count1, cv::hough_cmp_gt((int*)ddata));
01186 
01187         dist_sum = start_dist = ddata[sort_buf[nz_count1-1]];
01188         for( j = nz_count1 - 2; j >= 0; j-- )
01189         {
01190             float d = ddata[sort_buf[j]];
01191 
01192             if( d > max_radius )
01193                 break;
01194 
01195             if( d - start_dist > dr )
01196             {
01197                 float r_cur = ddata[sort_buf[(j + start_idx)/2]];
01198                 if( (start_idx - j)*r_best >= max_count*r_cur ||
01199                     (r_best < FLT_EPSILON && start_idx - j >= max_count) )
01200                 {
01201                     r_best = r_cur;
01202                     max_count = start_idx - j;
01203                 }
01204                 start_dist = d;
01205                 start_idx = j;
01206                 dist_sum = 0;
01207             }
01208             dist_sum += d;
01209         }
01210         // Check if the circle has enough support
01211         if( max_count > acc_threshold )
01212         {
01213             float c[3];
01214             c[0] = cx;
01215             c[1] = cy;
01216             c[2] = (float)r_best;
01217             cvSeqPush( circles, c );
01218             if( circles->total > circles_max )
01219                 return;
01220         }
01221     }
01222 }
01223 
01224 CV_IMPL CvSeq*
01225 cvHoughCircles( CvArr* src_image, void* circle_storage,
01226                 int method, double dp, double min_dist,
01227                 double param1, double param2,
01228                 int min_radius, int max_radius )
01229 {
01230     CvSeq* result = 0;
01231 
01232     CvMat stub, *img = (CvMat*)src_image;
01233     CvMat* mat = 0;
01234     CvSeq* circles = 0;
01235     CvSeq circles_header;
01236     CvSeqBlock circles_block;
01237     int circles_max = INT_MAX;
01238     int canny_threshold = cvRound(param1);
01239     int acc_threshold = cvRound(param2);
01240 
01241     img = cvGetMat( img, &stub );
01242 
01243     if( !CV_IS_MASK_ARR(img))
01244         CV_Error( CV_StsBadArg, "The source image must be 8-bit, single-channel" );
01245 
01246     if( !circle_storage )
01247         CV_Error( CV_StsNullPtr, "NULL destination" );
01248 
01249     if( dp <= 0 || min_dist <= 0 || canny_threshold <= 0 || acc_threshold <= 0 )
01250         CV_Error( CV_StsOutOfRange, "dp, min_dist, canny_threshold and acc_threshold must be all positive numbers" );
01251 
01252     min_radius = MAX( min_radius, 0 );
01253     if( max_radius <= 0 )
01254         max_radius = MAX( img->rows, img->cols );
01255     else if( max_radius <= min_radius )
01256         max_radius = min_radius + 2;
01257 
01258     if( CV_IS_STORAGE( circle_storage ))
01259     {
01260         circles = cvCreateSeq( CV_32FC3, sizeof(CvSeq),
01261             sizeof(float)*3, (CvMemStorage*)circle_storage );
01262     }
01263     else if( CV_IS_MAT( circle_storage ))
01264     {
01265         mat = (CvMat*)circle_storage;
01266 
01267         if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) ||
01268             CV_MAT_TYPE(mat->type) != CV_32FC3 )
01269             CV_Error( CV_StsBadArg,
01270             "The destination matrix should be continuous and have a single row or a single column" );
01271 
01272         circles = cvMakeSeqHeaderForArray( CV_32FC3, sizeof(CvSeq), sizeof(float)*3,
01273                 mat->data.ptr, mat->rows + mat->cols - 1, &circles_header, &circles_block );
01274         circles_max = circles->total;
01275         cvClearSeq( circles );
01276     }
01277     else
01278         CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );
01279 
01280     switch( method )
01281     {
01282     case CV_HOUGH_GRADIENT:
01283         icvHoughCirclesGradient( img, (float)dp, (float)min_dist,
01284                                 min_radius, max_radius, canny_threshold,
01285                                 acc_threshold, circles, circles_max );
01286           break;
01287     default:
01288         CV_Error( CV_StsBadArg, "Unrecognized method id" );
01289     }
01290 
01291     if( mat )
01292     {
01293         if( mat->cols > mat->rows )
01294             mat->cols = circles->total;
01295         else
01296             mat->rows = circles->total;
01297     }
01298     else
01299         result = circles;
01300 
01301     return result;
01302 }
01303 
01304 
01305 namespace cv
01306 {
01307 
01308 const int STORAGE_SIZE = 1 << 12;
01309 
01310 static void seqToMat(const CvSeq* seq, OutputArray _arr)
01311 {
01312     if( seq && seq->total > 0 )
01313     {
01314         _arr.create(1, seq->total, seq->flags, -1, true);
01315         Mat arr = _arr.getMat();
01316         cvCvtSeqToArray(seq, arr.ptr());
01317     }
01318     else
01319         _arr.release();
01320 }
01321 
01322 }
01323 
01324 void cv::HoughCircles( InputArray _image, OutputArray _circles,
01325                        int method, double dp, double min_dist,
01326                        double param1, double param2,
01327                        int minRadius, int maxRadius )
01328 {
01329     Ptr<CvMemStorage> storage(cvCreateMemStorage(STORAGE_SIZE));
01330     Mat image = _image.getMat();
01331     CvMat c_image = image;
01332     CvSeq* seq = cvHoughCircles( &c_image, storage, method,
01333                     dp, min_dist, param1, param2, minRadius, maxRadius );
01334     seqToMat(seq, _circles);
01335 }
01336 
01337 /* End of file. */
01338