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