cvhough.cpp

opencv库在TI DM6437上的移植,目前包括两个库cv.lib和cxcore.lib的工程

        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. */