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