--- /dev/null
+/*M///////////////////////////////////////////////////////////////////////////////////////\r
+//\r
+// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
+//\r
+// By downloading, copying, installing or using the software you agree to this license.\r
+// If you do not agree to this license, do not download, install,\r
+// copy or use the software.\r
+//\r
+//\r
+// Intel License Agreement\r
+// For Open Source Computer Vision Library\r
+//\r
+// Copyright (C) 2000, Intel Corporation, all rights reserved.\r
+// Third party copyrights are property of their respective owners.\r
+//\r
+// Redistribution and use in source and binary forms, with or without modification,\r
+// are permitted provided that the following conditions are met:\r
+//\r
+// * Redistribution's of source code must retain the above copyright notice,\r
+// this list of conditions and the following disclaimer.\r
+//\r
+// * Redistribution's in binary form must reproduce the above copyright notice,\r
+// this list of conditions and the following disclaimer in the documentation\r
+// and/or other materials provided with the distribution.\r
+//\r
+// * The name of Intel Corporation may not be used to endorse or promote products\r
+// derived from this software without specific prior written permission.\r
+//\r
+// This software is provided by the copyright holders and contributors "as is" and\r
+// any express or implied warranties, including, but not limited to, the implied\r
+// warranties of merchantability and fitness for a particular purpose are disclaimed.\r
+// In no event shall the Intel Corporation or contributors be liable for any direct,\r
+// indirect, incidental, special, exemplary, or consequential damages\r
+// (including, but not limited to, procurement of substitute goods or services;\r
+// loss of use, data, or profits; or business interruption) however caused\r
+// and on any theory of liability, whether in contract, strict liability,\r
+// or tort (including negligence or otherwise) arising in any way out of\r
+// the use of this software, even if advised of the possibility of such damage.\r
+//\r
+//M*/\r
+\r
+#include "_cv.h"\r
+#include "_cvlist.h"\r
+\r
+#define halfPi ((float)(CV_PI*0.5))\r
+#define Pi ((float)CV_PI)\r
+#define a0 0 /*-4.172325e-7f*/ /*(-(float)0x7)/((float)0x1000000); */\r
+#define a1 1.000025f /*((float)0x1922253)/((float)0x1000000)*2/Pi; */\r
+#define a2 -2.652905e-4f /*(-(float)0x2ae6)/((float)0x1000000)*4/(Pi*Pi); */\r
+#define a3 -0.165624f /*(-(float)0xa45511)/((float)0x1000000)*8/(Pi*Pi*Pi); */\r
+#define a4 -1.964532e-3f /*(-(float)0x30fd3)/((float)0x1000000)*16/(Pi*Pi*Pi*Pi); */\r
+#define a5 1.02575e-2f /*((float)0x191cac)/((float)0x1000000)*32/(Pi*Pi*Pi*Pi*Pi); */\r
+#define a6 -9.580378e-4f /*(-(float)0x3af27)/((float)0x1000000)*64/(Pi*Pi*Pi*Pi*Pi*Pi); */\r
+\r
+#define _sin(x) ((((((a6*(x) + a5)*(x) + a4)*(x) + a3)*(x) + a2)*(x) + a1)*(x) + a0)\r
+#define _cos(x) _sin(halfPi - (x))\r
+\r
+/****************************************************************************************\\r
+* Classical Hough Transform *\r
+\****************************************************************************************/\r
+\r
+typedef struct CvLinePolar\r
+{\r
+ float rho;\r
+ float angle;\r
+}\r
+CvLinePolar;\r
+\r
+/*=====================================================================================*/\r
+\r
+#define hough_cmp_gt(l1,l2) (aux[l1] > aux[l2])\r
+\r
+static CV_IMPLEMENT_QSORT_EX( icvHoughSortDescent32s, int, hough_cmp_gt, const int* )\r
+\r
+/*\r
+Here image is an input raster;\r
+step is it's step; size characterizes it's ROI;\r
+rho and theta are discretization steps (in pixels and radians correspondingly).\r
+threshold is the minimum number of pixels in the feature for it\r
+to be a candidate for line. lines is the output\r
+array of (rho, theta) pairs. linesMax is the buffer size (number of pairs).\r
+Functions return the actual number of found lines.\r
+*/\r
+static void\r
+icvHoughLinesStandard( const CvMat* img, float rho, float theta,\r
+ int threshold, CvSeq *lines, int linesMax )\r
+{\r
+ int *accum = 0;\r
+ int *sort_buf=0;\r
+ float *tabSin = 0;\r
+ float *tabCos = 0;\r
+\r
+ CV_FUNCNAME( "icvHoughLinesStandard" );\r
+\r
+ __BEGIN__;\r
+\r
+ const uchar* image;\r
+ int step, width, height;\r
+ int numangle, numrho;\r
+ int total = 0;\r
+ float ang;\r
+ int r, n;\r
+ int i, j;\r
+ float irho = 1 / rho;\r
+ double scale;\r
+\r
+ CV_ASSERT( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1 );\r
+\r
+ image = img->data.ptr;\r
+ step = img->step;\r
+ width = img->cols;\r
+ height = img->rows;\r
+\r
+ numangle = cvRound(CV_PI / theta);\r
+ numrho = cvRound(((width + height) * 2 + 1) / rho);\r
+\r
+ CV_CALL( accum = (int*)cvAlloc( sizeof(accum[0]) * (numangle+2) * (numrho+2) ));\r
+ CV_CALL( sort_buf = (int*)cvAlloc( sizeof(accum[0]) * numangle * numrho ));\r
+ CV_CALL( tabSin = (float*)cvAlloc( sizeof(tabSin[0]) * numangle ));\r
+ CV_CALL( tabCos = (float*)cvAlloc( sizeof(tabCos[0]) * numangle ));\r
+ memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );\r
+\r
+ for( ang = 0, n = 0; n < numangle; ang += theta, n++ )\r
+ {\r
+ tabSin[n] = (float)(sin(ang) * irho);\r
+ tabCos[n] = (float)(cos(ang) * irho);\r
+ }\r
+\r
+ // stage 1. fill accumulator\r
+ for( i = 0; i < height; i++ )\r
+ for( j = 0; j < width; j++ )\r
+ {\r
+ if( image[i * step + j] != 0 )\r
+ for( n = 0; n < numangle; n++ )\r
+ {\r
+ r = cvRound( j * tabCos[n] + i * tabSin[n] );\r
+ r += (numrho - 1) / 2;\r
+ accum[(n+1) * (numrho+2) + r+1]++;\r
+ }\r
+ }\r
+\r
+ // stage 2. find local maximums\r
+ for( r = 0; r < numrho; r++ )\r
+ for( n = 0; n < numangle; n++ )\r
+ {\r
+ int base = (n+1) * (numrho+2) + r+1;\r
+ if( accum[base] > threshold &&\r
+ accum[base] > accum[base - 1] && accum[base] >= accum[base + 1] &&\r
+ accum[base] > accum[base - numrho - 2] && accum[base] >= accum[base + numrho + 2] )\r
+ sort_buf[total++] = base;\r
+ }\r
+\r
+ // stage 3. sort the detected lines by accumulator value\r
+ icvHoughSortDescent32s( sort_buf, total, accum );\r
+\r
+ // stage 4. store the first min(total,linesMax) lines to the output buffer\r
+ linesMax = MIN(linesMax, total);\r
+ scale = 1./(numrho+2);\r
+ for( i = 0; i < linesMax; i++ )\r
+ {\r
+ CvLinePolar line;\r
+ int idx = sort_buf[i];\r
+ int n = cvFloor(idx*scale) - 1;\r
+ int r = idx - (n+1)*(numrho+2) - 1;\r
+ line.rho = (r - (numrho - 1)*0.5f) * rho;\r
+ line.angle = n * theta;\r
+ cvSeqPush( lines, &line );\r
+ }\r
+\r
+ __END__;\r
+\r
+ cvFree( &sort_buf );\r
+ cvFree( &tabSin );\r
+ cvFree( &tabCos );\r
+ cvFree( &accum );\r
+}\r
+\r
+\r
+/****************************************************************************************\\r
+* Multi-Scale variant of Classical Hough Transform *\r
+\****************************************************************************************/\r
+\r
+#if defined _MSC_VER && _MSC_VER >= 1200\r
+#pragma warning( disable: 4714 )\r
+#endif\r
+\r
+//DECLARE_AND_IMPLEMENT_LIST( _index, h_ );\r
+IMPLEMENT_LIST( _index, h_ )\r
+\r
+static void\r
+icvHoughLinesSDiv( const CvMat* img,\r
+ float rho, float theta, int threshold,\r
+ int srn, int stn,\r
+ CvSeq* lines, int linesMax )\r
+{\r
+ uchar *caccum = 0;\r
+ uchar *buffer = 0;\r
+ float *sinTable = 0;\r
+ int *x = 0;\r
+ int *y = 0;\r
+ _CVLIST *list = 0;\r
+\r
+ CV_FUNCNAME( "icvHoughLinesSDiv" );\r
+\r
+ __BEGIN__;\r
+\r
+#define _POINT(row, column)\\r
+ (image_src[(row)*step+(column)])\r
+\r
+ uchar *mcaccum = 0;\r
+ int rn, tn; /* number of rho and theta discrete values */\r
+ int index, i;\r
+ int ri, ti, ti1, ti0;\r
+ int row, col;\r
+ float r, t; /* Current rho and theta */\r
+ float rv; /* Some temporary rho value */\r
+ float irho;\r
+ float itheta;\r
+ float srho, stheta;\r
+ float isrho, istheta;\r
+\r
+ const uchar* image_src;\r
+ int w, h, step;\r
+ int fn = 0;\r
+ float xc, yc;\r
+\r
+ const float d2r = (float)(Pi / 180);\r
+ int sfn = srn * stn;\r
+ int fi;\r
+ int count;\r
+ int cmax = 0;\r
+\r
+ CVPOS pos;\r
+ _index *pindex;\r
+ _index vi;\r
+\r
+ CV_ASSERT( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1 );\r
+ CV_ASSERT( linesMax > 0 && rho > 0 && theta > 0 );\r
+\r
+ threshold = MIN( threshold, 255 );\r
+\r
+ image_src = img->data.ptr;\r
+ step = img->step;\r
+ w = img->cols;\r
+ h = img->rows;\r
+\r
+ irho = 1 / rho;\r
+ itheta = 1 / theta;\r
+ srho = rho / srn;\r
+ stheta = theta / stn;\r
+ isrho = 1 / srho;\r
+ istheta = 1 / stheta;\r
+\r
+ rn = cvFloor( sqrt( (double)w * w + (double)h * h ) * irho );\r
+ tn = cvFloor( 2 * Pi * itheta );\r
+\r
+ list = h_create_list__index( linesMax < 1000 ? linesMax : 1000 );\r
+ vi.value = threshold;\r
+ vi.rho = -1;\r
+ h_add_head__index( list, &vi );\r
+\r
+ /* Precalculating sin */\r
+ CV_CALL( sinTable = (float*)cvAlloc( 5 * tn * stn * sizeof( float )));\r
+\r
+ for( index = 0; index < 5 * tn * stn; index++ )\r
+ {\r
+ sinTable[index] = (float)cos( stheta * index * 0.2f );\r
+ }\r
+\r
+ CV_CALL( caccum = (uchar*)cvAlloc( rn * tn * sizeof( caccum[0] )));\r
+ memset( caccum, 0, rn * tn * sizeof( caccum[0] ));\r
+\r
+ /* Counting all feature pixels */\r
+ for( row = 0; row < h; row++ )\r
+ for( col = 0; col < w; col++ )\r
+ fn += _POINT( row, col ) != 0;\r
+\r
+ CV_CALL( x = (int*)cvAlloc( fn * sizeof(x[0])));\r
+ CV_CALL( y = (int*)cvAlloc( fn * sizeof(y[0])));\r
+\r
+ /* Full Hough Transform (it's accumulator update part) */\r
+ fi = 0;\r
+ for( row = 0; row < h; row++ )\r
+ {\r
+ for( col = 0; col < w; col++ )\r
+ {\r
+ if( _POINT( row, col ))\r
+ {\r
+ int halftn;\r
+ float r0;\r
+ float scale_factor;\r
+ int iprev = -1;\r
+ float phi, phi1;\r
+ float theta_it; /* Value of theta for iterating */\r
+\r
+ /* Remember the feature point */\r
+ x[fi] = col;\r
+ y[fi] = row;\r
+ fi++;\r
+\r
+ yc = (float) row + 0.5f;\r
+ xc = (float) col + 0.5f;\r
+\r
+ /* Update the accumulator */\r
+ t = (float) fabs( cvFastArctan( yc, xc ) * d2r );\r
+ r = (float) sqrt( (double)xc * xc + (double)yc * yc );\r
+ r0 = r * irho;\r
+ ti0 = cvFloor( (t + Pi / 2) * itheta );\r
+\r
+ caccum[ti0]++;\r
+\r
+ theta_it = rho / r;\r
+ theta_it = theta_it < theta ? theta_it : theta;\r
+ scale_factor = theta_it * itheta;\r
+ halftn = cvFloor( Pi / theta_it );\r
+ for( ti1 = 1, phi = theta_it - halfPi, phi1 = (theta_it + t) * itheta;\r
+ ti1 < halftn; ti1++, phi += theta_it, phi1 += scale_factor )\r
+ {\r
+ rv = r0 * _cos( phi );\r
+ i = cvFloor( rv ) * tn;\r
+ i += cvFloor( phi1 );\r
+ assert( i >= 0 );\r
+ assert( i < rn * tn );\r
+ caccum[i] = (uchar) (caccum[i] + ((i ^ iprev) != 0));\r
+ iprev = i;\r
+ if( cmax < caccum[i] )\r
+ cmax = caccum[i];\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ /* Starting additional analysis */\r
+ count = 0;\r
+ for( ri = 0; ri < rn; ri++ )\r
+ {\r
+ for( ti = 0; ti < tn; ti++ )\r
+ {\r
+ if( caccum[ri * tn + ti > threshold] )\r
+ {\r
+ count++;\r
+ }\r
+ }\r
+ }\r
+\r
+ if( count * 100 > rn * tn )\r
+ {\r
+ icvHoughLinesStandard( img, rho, theta, threshold, lines, linesMax );\r
+ EXIT;\r
+ }\r
+\r
+ CV_CALL( buffer = (uchar *) cvAlloc(srn * stn + 2));\r
+ mcaccum = buffer + 1;\r
+\r
+ count = 0;\r
+ for( ri = 0; ri < rn; ri++ )\r
+ {\r
+ for( ti = 0; ti < tn; ti++ )\r
+ {\r
+ if( caccum[ri * tn + ti] > threshold )\r
+ {\r
+ count++;\r
+ memset( mcaccum, 0, sfn * sizeof( uchar ));\r
+\r
+ for( index = 0; index < fn; index++ )\r
+ {\r
+ int ti2;\r
+ float r0;\r
+\r
+ yc = (float) y[index] + 0.5f;\r
+ xc = (float) x[index] + 0.5f;\r
+\r
+ /* Update the accumulator */\r
+ t = (float) fabs( cvFastArctan( yc, xc ) * d2r );\r
+ r = (float) sqrt( (double)xc * xc + (double)yc * yc ) * isrho;\r
+ ti0 = cvFloor( (t + Pi * 0.5f) * istheta );\r
+ ti2 = (ti * stn - ti0) * 5;\r
+ r0 = (float) ri *srn;\r
+\r
+ for( ti1 = 0 /*, phi = ti*theta - Pi/2 - t */ ; ti1 < stn; ti1++, ti2 += 5\r
+ /*phi += stheta */ )\r
+ {\r
+ /*rv = r*_cos(phi) - r0; */\r
+ rv = r * sinTable[(int) (abs( ti2 ))] - r0;\r
+ i = cvFloor( rv ) * stn + ti1;\r
+\r
+ i = CV_IMAX( i, -1 );\r
+ i = CV_IMIN( i, sfn );\r
+ mcaccum[i]++;\r
+ assert( i >= -1 );\r
+ assert( i <= sfn );\r
+ }\r
+ }\r
+\r
+ /* Find peaks in maccum... */\r
+ for( index = 0; index < sfn; index++ )\r
+ {\r
+ i = 0;\r
+ pos = h_get_tail_pos__index( list );\r
+ if( h_get_prev__index( &pos )->value < mcaccum[index] )\r
+ {\r
+ vi.value = mcaccum[index];\r
+ vi.rho = index / stn * srho + ri * rho;\r
+ vi.theta = index % stn * stheta + ti * theta - halfPi;\r
+ while( h_is_pos__index( pos ))\r
+ {\r
+ if( h_get__index( pos )->value > mcaccum[index] )\r
+ {\r
+ h_insert_after__index( list, pos, &vi );\r
+ if( h_get_count__index( list ) > linesMax )\r
+ {\r
+ h_remove_tail__index( list );\r
+ }\r
+ break;\r
+ }\r
+ h_get_prev__index( &pos );\r
+ }\r
+ if( !h_is_pos__index( pos ))\r
+ {\r
+ h_add_head__index( list, &vi );\r
+ if( h_get_count__index( list ) > linesMax )\r
+ {\r
+ h_remove_tail__index( list );\r
+ }\r
+ }\r
+ }\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ pos = h_get_head_pos__index( list );\r
+ if( h_get_count__index( list ) == 1 )\r
+ {\r
+ if( h_get__index( pos )->rho < 0 )\r
+ {\r
+ h_clear_list__index( list );\r
+ }\r
+ }\r
+ else\r
+ {\r
+ while( h_is_pos__index( pos ))\r
+ {\r
+ CvLinePolar line;\r
+ pindex = h_get__index( pos );\r
+ if( pindex->rho < 0 )\r
+ {\r
+ /* This should be the last element... */\r
+ h_get_next__index( &pos );\r
+ assert( !h_is_pos__index( pos ));\r
+ break;\r
+ }\r
+ line.rho = pindex->rho;\r
+ line.angle = pindex->theta;\r
+ cvSeqPush( lines, &line );\r
+\r
+ if( lines->total >= linesMax )\r
+ EXIT;\r
+ h_get_next__index( &pos );\r
+ }\r
+ }\r
+\r
+ __END__;\r
+\r
+ h_destroy_list__index( list );\r
+ cvFree( &sinTable );\r
+ cvFree( &x );\r
+ cvFree( &y );\r
+ cvFree( &caccum );\r
+ cvFree( &buffer );\r
+}\r
+\r
+\r
+/****************************************************************************************\\r
+* Probabilistic Hough Transform *\r
+\****************************************************************************************/\r
+\r
+static void\r
+icvHoughLinesProbabalistic( CvMat* image,\r
+ float rho, float theta, int threshold,\r
+ int lineLength, int lineGap,\r
+ CvSeq *lines, int linesMax )\r
+{\r
+ cv::Mat accum, mask;\r
+ cv::vector<float> trigtab;\r
+ cv::MemStorage storage(cvCreateMemStorage(0));\r
+\r
+ CvSeq* seq;\r
+ CvSeqWriter writer;\r
+ int width, height;\r
+ int numangle, numrho;\r
+ float ang;\r
+ int r, n, count;\r
+ CvPoint pt;\r
+ float irho = 1 / rho;\r
+ CvRNG rng = cvRNG(-1);\r
+ const float* ttab;\r
+ uchar* mdata0;\r
+\r
+ CV_Assert( CV_IS_MAT(image) && CV_MAT_TYPE(image->type) == CV_8UC1 );\r
+\r
+ width = image->cols;\r
+ height = image->rows;\r
+\r
+ numangle = cvRound(CV_PI / theta);\r
+ numrho = cvRound(((width + height) * 2 + 1) / rho);\r
+\r
+ accum.create( numangle, numrho, CV_32SC1 );\r
+ mask.create( height, width, CV_8UC1 );\r
+ trigtab.resize(numangle*2);\r
+ accum = cv::Scalar(0);\r
+\r
+ for( ang = 0, n = 0; n < numangle; ang += theta, n++ )\r
+ {\r
+ trigtab[n*2] = (float)(cos(ang) * irho);\r
+ trigtab[n*2+1] = (float)(sin(ang) * irho);\r
+ }\r
+ ttab = &trigtab[0];\r
+ mdata0 = mask.data;\r
+\r
+ cvStartWriteSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage, &writer );\r
+\r
+ // stage 1. collect non-zero image points\r
+ for( pt.y = 0, count = 0; pt.y < height; pt.y++ )\r
+ {\r
+ const uchar* data = image->data.ptr + pt.y*image->step;\r
+ uchar* mdata = mdata0 + pt.y*width;\r
+ for( pt.x = 0; pt.x < width; pt.x++ )\r
+ {\r
+ if( data[pt.x] )\r
+ {\r
+ mdata[pt.x] = (uchar)1;\r
+ CV_WRITE_SEQ_ELEM( pt, writer );\r
+ }\r
+ else\r
+ mdata[pt.x] = 0;\r
+ }\r
+ }\r
+\r
+ seq = cvEndWriteSeq( &writer );\r
+ count = seq->total;\r
+\r
+ // stage 2. process all the points in random order\r
+ for( ; count > 0; count-- )\r
+ {\r
+ // choose random point out of the remaining ones\r
+ int idx = cvRandInt(&rng) % count;\r
+ int max_val = threshold-1, max_n = 0;\r
+ CvPoint* pt = (CvPoint*)cvGetSeqElem( seq, idx );\r
+ CvPoint line_end[2] = {{0,0}, {0,0}};\r
+ float a, b;\r
+ int* adata = (int*)accum.data;\r
+ int i, j, k, x0, y0, dx0, dy0, xflag;\r
+ int good_line;\r
+ const int shift = 16;\r
+\r
+ i = pt->y;\r
+ j = pt->x;\r
+\r
+ // "remove" it by overriding it with the last element\r
+ *pt = *(CvPoint*)cvGetSeqElem( seq, count-1 );\r
+\r
+ // check if it has been excluded already (i.e. belongs to some other line)\r
+ if( !mdata0[i*width + j] )\r
+ continue;\r
+\r
+ // update accumulator, find the most probable line\r
+ for( n = 0; n < numangle; n++, adata += numrho )\r
+ {\r
+ r = cvRound( j * ttab[n*2] + i * ttab[n*2+1] );\r
+ r += (numrho - 1) / 2;\r
+ int val = ++adata[r];\r
+ if( max_val < val )\r
+ {\r
+ max_val = val;\r
+ max_n = n;\r
+ }\r
+ }\r
+\r
+ // if it is too "weak" candidate, continue with another point\r
+ if( max_val < threshold )\r
+ continue;\r
+\r
+ // from the current point walk in each direction\r
+ // along the found line and extract the line segment\r
+ a = -ttab[max_n*2+1];\r
+ b = ttab[max_n*2];\r
+ x0 = j;\r
+ y0 = i;\r
+ if( fabs(a) > fabs(b) )\r
+ {\r
+ xflag = 1;\r
+ dx0 = a > 0 ? 1 : -1;\r
+ dy0 = cvRound( b*(1 << shift)/fabs(a) );\r
+ y0 = (y0 << shift) + (1 << (shift-1));\r
+ }\r
+ else\r
+ {\r
+ xflag = 0;\r
+ dy0 = b > 0 ? 1 : -1;\r
+ dx0 = cvRound( a*(1 << shift)/fabs(b) );\r
+ x0 = (x0 << shift) + (1 << (shift-1));\r
+ }\r
+\r
+ for( k = 0; k < 2; k++ )\r
+ {\r
+ int gap = 0, x = x0, y = y0, dx = dx0, dy = dy0;\r
+\r
+ if( k > 0 )\r
+ dx = -dx, dy = -dy;\r
+\r
+ // walk along the line using fixed-point arithmetics,\r
+ // stop at the image border or in case of too big gap\r
+ for( ;; x += dx, y += dy )\r
+ {\r
+ uchar* mdata;\r
+ int i1, j1;\r
+\r
+ if( xflag )\r
+ {\r
+ j1 = x;\r
+ i1 = y >> shift;\r
+ }\r
+ else\r
+ {\r
+ j1 = x >> shift;\r
+ i1 = y;\r
+ }\r
+\r
+ if( j1 < 0 || j1 >= width || i1 < 0 || i1 >= height )\r
+ break;\r
+\r
+ mdata = mdata0 + i1*width + j1;\r
+\r
+ // for each non-zero point:\r
+ // update line end,\r
+ // clear the mask element\r
+ // reset the gap\r
+ if( *mdata )\r
+ {\r
+ gap = 0;\r
+ line_end[k].y = i1;\r
+ line_end[k].x = j1;\r
+ }\r
+ else if( ++gap > lineGap )\r
+ break;\r
+ }\r
+ }\r
+\r
+ good_line = abs(line_end[1].x - line_end[0].x) >= lineLength ||\r
+ abs(line_end[1].y - line_end[0].y) >= lineLength;\r
+\r
+ for( k = 0; k < 2; k++ )\r
+ {\r
+ int x = x0, y = y0, dx = dx0, dy = dy0;\r
+\r
+ if( k > 0 )\r
+ dx = -dx, dy = -dy;\r
+\r
+ // walk along the line using fixed-point arithmetics,\r
+ // stop at the image border or in case of too big gap\r
+ for( ;; x += dx, y += dy )\r
+ {\r
+ uchar* mdata;\r
+ int i1, j1;\r
+\r
+ if( xflag )\r
+ {\r
+ j1 = x;\r
+ i1 = y >> shift;\r
+ }\r
+ else\r
+ {\r
+ j1 = x >> shift;\r
+ i1 = y;\r
+ }\r
+\r
+ mdata = mdata0 + i1*width + j1;\r
+\r
+ // for each non-zero point:\r
+ // update line end,\r
+ // clear the mask element\r
+ // reset the gap\r
+ if( *mdata )\r
+ {\r
+ if( good_line )\r
+ {\r
+ adata = (int*)accum.data;\r
+ for( n = 0; n < numangle; n++, adata += numrho )\r
+ {\r
+ r = cvRound( j1 * ttab[n*2] + i1 * ttab[n*2+1] );\r
+ r += (numrho - 1) / 2;\r
+ adata[r]--;\r
+ }\r
+ }\r
+ *mdata = 0;\r
+ }\r
+\r
+ if( i1 == line_end[k].y && j1 == line_end[k].x )\r
+ break;\r
+ }\r
+ }\r
+\r
+ if( good_line )\r
+ {\r
+ CvRect lr = { line_end[0].x, line_end[0].y, line_end[1].x, line_end[1].y };\r
+ cvSeqPush( lines, &lr );\r
+ if( lines->total >= linesMax )\r
+ return;\r
+ }\r
+ }\r
+}\r
+\r
+/* Wrapper function for standard hough transform */\r
+CV_IMPL CvSeq*\r
+cvHoughLines2( CvArr* src_image, void* lineStorage, int method,\r
+ double rho, double theta, int threshold,\r
+ double param1, double param2 )\r
+{\r
+ CvSeq* result = 0;\r
+\r
+ CV_FUNCNAME( "cvHoughLines" );\r
+\r
+ __BEGIN__;\r
+\r
+ CvMat stub, *img = (CvMat*)src_image;\r
+ CvMat* mat = 0;\r
+ CvSeq* lines = 0;\r
+ CvSeq lines_header;\r
+ CvSeqBlock lines_block;\r
+ int lineType, elemSize;\r
+ int linesMax = INT_MAX;\r
+ int iparam1, iparam2;\r
+\r
+ CV_CALL( img = cvGetMat( img, &stub ));\r
+\r
+ if( !CV_IS_MASK_ARR(img))\r
+ CV_ERROR( CV_StsBadArg, "The source image must be 8-bit, single-channel" );\r
+\r
+ if( !lineStorage )\r
+ CV_ERROR( CV_StsNullPtr, "NULL destination" );\r
+\r
+ if( rho <= 0 || theta <= 0 || threshold <= 0 )\r
+ CV_ERROR( CV_StsOutOfRange, "rho, theta and threshold must be positive" );\r
+\r
+ if( method != CV_HOUGH_PROBABILISTIC )\r
+ {\r
+ lineType = CV_32FC2;\r
+ elemSize = sizeof(float)*2;\r
+ }\r
+ else\r
+ {\r
+ lineType = CV_32SC4;\r
+ elemSize = sizeof(int)*4;\r
+ }\r
+\r
+ if( CV_IS_STORAGE( lineStorage ))\r
+ {\r
+ CV_CALL( lines = cvCreateSeq( lineType, sizeof(CvSeq), elemSize, (CvMemStorage*)lineStorage ));\r
+ }\r
+ else if( CV_IS_MAT( lineStorage ))\r
+ {\r
+ mat = (CvMat*)lineStorage;\r
+\r
+ if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) )\r
+ CV_ERROR( CV_StsBadArg,\r
+ "The destination matrix should be continuous and have a single row or a single column" );\r
+\r
+ if( CV_MAT_TYPE( mat->type ) != lineType )\r
+ CV_ERROR( CV_StsBadArg,\r
+ "The destination matrix data type is inappropriate, see the manual" );\r
+\r
+ CV_CALL( lines = cvMakeSeqHeaderForArray( lineType, sizeof(CvSeq), elemSize, mat->data.ptr,\r
+ mat->rows + mat->cols - 1, &lines_header, &lines_block ));\r
+ linesMax = lines->total;\r
+ CV_CALL( cvClearSeq( lines ));\r
+ }\r
+ else\r
+ {\r
+ CV_ERROR( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );\r
+ }\r
+\r
+ iparam1 = cvRound(param1);\r
+ iparam2 = cvRound(param2);\r
+\r
+ switch( method )\r
+ {\r
+ case CV_HOUGH_STANDARD:\r
+ CV_CALL( icvHoughLinesStandard( img, (float)rho,\r
+ (float)theta, threshold, lines, linesMax ));\r
+ break;\r
+ case CV_HOUGH_MULTI_SCALE:\r
+ CV_CALL( icvHoughLinesSDiv( img, (float)rho, (float)theta,\r
+ threshold, iparam1, iparam2, lines, linesMax ));\r
+ break;\r
+ case CV_HOUGH_PROBABILISTIC:\r
+ CV_CALL( icvHoughLinesProbabalistic( img, (float)rho, (float)theta,\r
+ threshold, iparam1, iparam2, lines, linesMax ));\r
+ break;\r
+ default:\r
+ CV_ERROR( CV_StsBadArg, "Unrecognized method id" );\r
+ }\r
+\r
+ if( mat )\r
+ {\r
+ if( mat->cols > mat->rows )\r
+ mat->cols = lines->total;\r
+ else\r
+ mat->rows = lines->total;\r
+ }\r
+ else\r
+ {\r
+ result = lines;\r
+ }\r
+\r
+ __END__;\r
+\r
+ return result;\r
+}\r
+\r
+\r
+/****************************************************************************************\\r
+* Circle Detection *\r
+\****************************************************************************************/\r
+\r
+static void\r
+icvHoughCirclesGradient( CvMat* img, float dp, float min_dist,\r
+ int min_radius, int max_radius,\r
+ int canny_threshold, int acc_threshold,\r
+ CvSeq* circles, int circles_max )\r
+{\r
+ const int SHIFT = 10, ONE = 1 << SHIFT, R_THRESH = 30;\r
+ CvMat *dx = 0, *dy = 0;\r
+ CvMat *edges = 0;\r
+ CvMat *accum = 0;\r
+ int* sort_buf = 0;\r
+ CvMat* dist_buf = 0;\r
+ CvMemStorage* storage = 0;\r
+\r
+ CV_FUNCNAME( "icvHoughCirclesGradient" );\r
+\r
+ __BEGIN__;\r
+\r
+ int x, y, i, j, center_count, nz_count;\r
+ int rows, cols, arows, acols;\r
+ int astep, *adata;\r
+ float* ddata;\r
+ CvSeq *nz, *centers;\r
+ float idp, dr;\r
+ CvSeqReader reader;\r
+\r
+ CV_CALL( edges = cvCreateMat( img->rows, img->cols, CV_8UC1 ));\r
+ CV_CALL( cvCanny( img, edges, MAX(canny_threshold/2,1), canny_threshold, 3 ));\r
+\r
+ CV_CALL( dx = cvCreateMat( img->rows, img->cols, CV_16SC1 ));\r
+ CV_CALL( dy = cvCreateMat( img->rows, img->cols, CV_16SC1 ));\r
+ CV_CALL( cvSobel( img, dx, 1, 0, 3 ));\r
+ CV_CALL( cvSobel( img, dy, 0, 1, 3 ));\r
+\r
+ if( dp < 1.f )\r
+ dp = 1.f;\r
+ idp = 1.f/dp;\r
+ CV_CALL( accum = cvCreateMat( cvCeil(img->rows*idp)+2, cvCeil(img->cols*idp)+2, CV_32SC1 ));\r
+ CV_CALL( cvZero(accum));\r
+\r
+ CV_CALL( storage = cvCreateMemStorage() );\r
+ CV_CALL( nz = cvCreateSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage ));\r
+ CV_CALL( centers = cvCreateSeq( CV_32SC1, sizeof(CvSeq), sizeof(int), storage ));\r
+\r
+ rows = img->rows;\r
+ cols = img->cols;\r
+ arows = accum->rows - 2;\r
+ acols = accum->cols - 2;\r
+ adata = accum->data.i;\r
+ astep = accum->step/sizeof(adata[0]);\r
+\r
+ for( y = 0; y < rows; y++ )\r
+ {\r
+ const uchar* edges_row = edges->data.ptr + y*edges->step;\r
+ const short* dx_row = (const short*)(dx->data.ptr + y*dx->step);\r
+ const short* dy_row = (const short*)(dy->data.ptr + y*dy->step);\r
+\r
+ for( x = 0; x < cols; x++ )\r
+ {\r
+ float vx, vy;\r
+ int sx, sy, x0, y0, x1, y1, r, k;\r
+ CvPoint pt;\r
+\r
+ vx = dx_row[x];\r
+ vy = dy_row[x];\r
+\r
+ if( !edges_row[x] || (vx == 0 && vy == 0) )\r
+ continue;\r
+\r
+ float mag = sqrt(vx*vx+vy*vy);\r
+ assert( mag >= 1 );\r
+ sx = cvRound((vx*idp)*ONE/mag);\r
+ sy = cvRound((vy*idp)*ONE/mag);\r
+\r
+ x0 = cvRound((x*idp)*ONE);\r
+ y0 = cvRound((y*idp)*ONE);\r
+\r
+ for( k = 0; k < 2; k++ )\r
+ {\r
+ x1 = x0 + min_radius * sx;\r
+ y1 = y0 + min_radius * sy;\r
+\r
+ for( r = min_radius; r <= max_radius; x1 += sx, y1 += sy, r++ )\r
+ {\r
+ int x2 = x1 >> SHIFT, y2 = y1 >> SHIFT;\r
+ if( (unsigned)x2 >= (unsigned)acols ||\r
+ (unsigned)y2 >= (unsigned)arows )\r
+ break;\r
+ adata[y2*astep + x2]++;\r
+ }\r
+\r
+ sx = -sx; sy = -sy;\r
+ }\r
+\r
+ pt.x = x; pt.y = y;\r
+ cvSeqPush( nz, &pt );\r
+ }\r
+ }\r
+\r
+ nz_count = nz->total;\r
+ if( !nz_count )\r
+ EXIT;\r
+\r
+ for( y = 1; y < arows - 1; y++ )\r
+ {\r
+ for( x = 1; x < acols - 1; x++ )\r
+ {\r
+ int base = y*(acols+2) + x;\r
+ if( adata[base] > acc_threshold &&\r
+ adata[base] > adata[base-1] && adata[base] > adata[base+1] &&\r
+ adata[base] > adata[base-acols-2] && adata[base] > adata[base+acols+2] )\r
+ cvSeqPush(centers, &base);\r
+ }\r
+ }\r
+\r
+ center_count = centers->total;\r
+ if( !center_count )\r
+ EXIT;\r
+\r
+ CV_CALL( sort_buf = (int*)cvAlloc( MAX(center_count,nz_count)*sizeof(sort_buf[0]) ));\r
+ cvCvtSeqToArray( centers, sort_buf );\r
+\r
+ icvHoughSortDescent32s( sort_buf, center_count, adata );\r
+ cvClearSeq( centers );\r
+ cvSeqPushMulti( centers, sort_buf, center_count );\r
+\r
+ CV_CALL( dist_buf = cvCreateMat( 1, nz_count, CV_32FC1 ));\r
+ ddata = dist_buf->data.fl;\r
+\r
+ dr = dp;\r
+ min_dist = MAX( min_dist, dp );\r
+ min_dist *= min_dist;\r
+\r
+ for( i = 0; i < centers->total; i++ )\r
+ {\r
+ int ofs = *(int*)cvGetSeqElem( centers, i );\r
+ y = ofs/(acols+2) - 1;\r
+ x = ofs - (y+1)*(acols+2) - 1;\r
+ float cx = (float)(x*dp), cy = (float)(y*dp);\r
+ int start_idx = nz_count - 1;\r
+ float start_dist, dist_sum;\r
+ float r_best = 0, c[3];\r
+ int max_count = R_THRESH;\r
+\r
+ for( j = 0; j < circles->total; j++ )\r
+ {\r
+ float* c = (float*)cvGetSeqElem( circles, j );\r
+ if( (c[0] - cx)*(c[0] - cx) + (c[1] - cy)*(c[1] - cy) < min_dist )\r
+ break;\r
+ }\r
+\r
+ if( j < circles->total )\r
+ continue;\r
+\r
+ cvStartReadSeq( nz, &reader );\r
+ for( j = 0; j < nz_count; j++ )\r
+ {\r
+ CvPoint pt;\r
+ float _dx, _dy;\r
+ CV_READ_SEQ_ELEM( pt, reader );\r
+ _dx = cx - pt.x; _dy = cy - pt.y;\r
+ ddata[j] = _dx*_dx + _dy*_dy;\r
+ sort_buf[j] = j;\r
+ }\r
+\r
+ cvPow( dist_buf, dist_buf, 0.5 );\r
+ icvHoughSortDescent32s( sort_buf, nz_count, (int*)ddata );\r
+\r
+ dist_sum = start_dist = ddata[sort_buf[nz_count-1]];\r
+ for( j = nz_count - 2; j >= 0; j-- )\r
+ {\r
+ float d = ddata[sort_buf[j]];\r
+\r
+ if( d > max_radius )\r
+ break;\r
+\r
+ if( d - start_dist > dr )\r
+ {\r
+ float r_cur = ddata[sort_buf[(j + start_idx)/2]];\r
+ if( (start_idx - j)*r_best >= max_count*r_cur ||\r
+ (r_best < FLT_EPSILON && start_idx - j >= max_count) )\r
+ {\r
+ r_best = r_cur;\r
+ max_count = start_idx - j;\r
+ }\r
+ start_dist = d;\r
+ start_idx = j;\r
+ dist_sum = 0;\r
+ }\r
+ dist_sum += d;\r
+ }\r
+\r
+ if( max_count > R_THRESH )\r
+ {\r
+ c[0] = cx;\r
+ c[1] = cy;\r
+ c[2] = (float)r_best;\r
+ cvSeqPush( circles, c );\r
+ if( circles->total > circles_max )\r
+ EXIT;\r
+ }\r
+ }\r
+\r
+ __END__;\r
+\r
+ cvReleaseMat( &dist_buf );\r
+ cvFree( &sort_buf );\r
+ cvReleaseMemStorage( &storage );\r
+ cvReleaseMat( &edges );\r
+ cvReleaseMat( &dx );\r
+ cvReleaseMat( &dy );\r
+ cvReleaseMat( &accum );\r
+}\r
+\r
+CV_IMPL CvSeq*\r
+cvHoughCircles( CvArr* src_image, void* circle_storage,\r
+ int method, double dp, double min_dist,\r
+ double param1, double param2,\r
+ int min_radius, int max_radius )\r
+{\r
+ CvSeq* result = 0;\r
+\r
+ CV_FUNCNAME( "cvHoughCircles" );\r
+\r
+ __BEGIN__;\r
+\r
+ CvMat stub, *img = (CvMat*)src_image;\r
+ CvMat* mat = 0;\r
+ CvSeq* circles = 0;\r
+ CvSeq circles_header;\r
+ CvSeqBlock circles_block;\r
+ int circles_max = INT_MAX;\r
+ int canny_threshold = cvRound(param1);\r
+ int acc_threshold = cvRound(param2);\r
+\r
+ CV_CALL( img = cvGetMat( img, &stub ));\r
+\r
+ if( !CV_IS_MASK_ARR(img))\r
+ CV_ERROR( CV_StsBadArg, "The source image must be 8-bit, single-channel" );\r
+\r
+ if( !circle_storage )\r
+ CV_ERROR( CV_StsNullPtr, "NULL destination" );\r
+\r
+ if( dp <= 0 || min_dist <= 0 || canny_threshold <= 0 || acc_threshold <= 0 )\r
+ CV_ERROR( CV_StsOutOfRange, "dp, min_dist, canny_threshold and acc_threshold must be all positive numbers" );\r
+\r
+ min_radius = MAX( min_radius, 0 );\r
+ if( max_radius <= 0 )\r
+ max_radius = MAX( img->rows, img->cols );\r
+ else if( max_radius <= min_radius )\r
+ max_radius = min_radius + 2;\r
+\r
+ if( CV_IS_STORAGE( circle_storage ))\r
+ {\r
+ CV_CALL( circles = cvCreateSeq( CV_32FC3, sizeof(CvSeq),\r
+ sizeof(float)*3, (CvMemStorage*)circle_storage ));\r
+ }\r
+ else if( CV_IS_MAT( circle_storage ))\r
+ {\r
+ mat = (CvMat*)circle_storage;\r
+\r
+ if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) ||\r
+ CV_MAT_TYPE(mat->type) != CV_32FC3 )\r
+ CV_ERROR( CV_StsBadArg,\r
+ "The destination matrix should be continuous and have a single row or a single column" );\r
+\r
+ CV_CALL( circles = cvMakeSeqHeaderForArray( CV_32FC3, sizeof(CvSeq), sizeof(float)*3,\r
+ mat->data.ptr, mat->rows + mat->cols - 1, &circles_header, &circles_block ));\r
+ circles_max = circles->total;\r
+ CV_CALL( cvClearSeq( circles ));\r
+ }\r
+ else\r
+ {\r
+ CV_ERROR( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );\r
+ }\r
+\r
+ switch( method )\r
+ {\r
+ case CV_HOUGH_GRADIENT:\r
+ CV_CALL( icvHoughCirclesGradient( img, (float)dp, (float)min_dist,\r
+ min_radius, max_radius, canny_threshold,\r
+ acc_threshold, circles, circles_max ));\r
+ break;\r
+ default:\r
+ CV_ERROR( CV_StsBadArg, "Unrecognized method id" );\r
+ }\r
+\r
+ if( mat )\r
+ {\r
+ if( mat->cols > mat->rows )\r
+ mat->cols = circles->total;\r
+ else\r
+ mat->rows = circles->total;\r
+ }\r
+ else\r
+ result = circles;\r
+\r
+ __END__;\r
+\r
+ return result;\r
+}\r
+\r
+\r
+namespace cv\r
+{\r
+\r
+const int STORAGE_SIZE = 1 << 12;\r
+\r
+void HoughLines( const Mat& image, vector<Vec2f>& lines,\r
+ double rho, double theta, int threshold,\r
+ double srn, double stn )\r
+{\r
+ CvMemStorage* storage = cvCreateMemStorage(STORAGE_SIZE);\r
+ CvMat _image = image;\r
+ CvSeq* seq = cvHoughLines2( &_image, storage, srn == 0 && stn == 0 ?\r
+ CV_HOUGH_STANDARD : CV_HOUGH_MULTI_SCALE,\r
+ rho, theta, threshold, srn, stn );\r
+ Seq<Vec2f>(seq).copyTo(lines);\r
+}\r
+\r
+void HoughLinesP( Mat& image, vector<Vec4i>& lines,\r
+ double rho, double theta, int threshold,\r
+ double minLineLength, double maxGap )\r
+{\r
+ CvMemStorage* storage = cvCreateMemStorage(STORAGE_SIZE);\r
+ CvMat _image = image;\r
+ CvSeq* seq = cvHoughLines2( &_image, storage, CV_HOUGH_PROBABILISTIC,\r
+ rho, theta, threshold, minLineLength, maxGap );\r
+ Seq<Vec4i>(seq).copyTo(lines);\r
+}\r
+\r
+void HoughCircles( const Mat& image, vector<Vec3f>& circles,\r
+ int method, double dp, double min_dist,\r
+ double param1, double param2,\r
+ int minRadius, int maxRadius )\r
+{\r
+ CvMemStorage* storage = cvCreateMemStorage(STORAGE_SIZE);\r
+ CvMat _image = image;\r
+ CvSeq* seq = cvHoughCircles( &_image, storage, method,\r
+ dp, min_dist, param1, param2, minRadius, maxRadius );\r
+ Seq<Vec3f>(seq).copyTo(circles);\r
+}\r
+\r
+}\r
+\r
+/* End of file. */\r
projects / opencv / blobdiff