OpenCV: cvHoughLines2() 연습 코드

OpenCV 라이브러리의 Hough transform에 의한 직선 찾기 함수

CvSeq* cvHoughLines2(CvArr* image, void* storage, int method, double rho, double theta, int threshold, double param1=0, double param2=0)¶

Finds lines in a binary image using a Hough transform.

Parameters:

image – The 8-bit, single-channel, binary source image. In the case of a probabilistic method, the image is modified by the function storage – The storage for the lines that are detected. It can be a memory storage (in this case a sequence of lines is created in the storage and returned by the function) or single row/single column matrix (CvMat*) of a particular type (see below) to which the lines’ parameters are written. The matrix header is modified by the function so its cols or rows will contain the number of lines detected. If storage is a matrix and the actual number of lines exceeds the matrix size, the maximum possible number of lines is returned (in the case of standard hough transform the lines are sorted by the accumulator value) method – The Hough transform variant, one of the following: CV_HOUGH_STANDARD - classical or standard Hough transform. Every line is represented by two floating-point numbers , where is a distance between (0,0) point and the line, and is the angle between x-axis and the normal to the line. Thus, the matrix must be (the created sequence will be) of CV_32FC2 type CV_HOUGH_PROBABILISTIC - probabilistic Hough transform (more efficient in case if picture contains a few long linear segments). It returns line segments rather than the whole line. Each segment is represented by starting and ending points, and the matrix must be (the created sequence will be) of CV_32SC4 type CV_HOUGH_MULTI_SCALE - multi-scale variant of the classical Hough transform. The lines are encoded the same way as CV_HOUGH_STANDARD rho – Distance resolution in pixel-related units theta – Angle resolution measured in radians threshold – Threshold parameter. A line is returned by the function if the corresponding accumulator value is greater than threshold param1 – The first method-dependent parameter: For the classical Hough transform it is not used (0). For the probabilistic Hough transform it is the minimum line length. For the multi-scale Hough transform it is the divisor for the distance resolution . (The coarse distance resolution will be and the accurate resolution will be ). param2 – The second method-dependent parameter: For the classical Hough transform it is not used (0). For the probabilistic Hough transform it is the maximum gap between line segments lying on the same line to treat them as a single line segment (i.e. to join them). For the multi-scale Hough transform it is the divisor for the angle resolution . (The coarse angle resolution will be and the accurate resolution will be ).

Memory storage is a low-level structure used to store dynamicly growing data structures such as sequences, contours, graphs, subdivisions, etc.

입력 이미지가 8비트 단일 채널이어야 하므로,
다음과 같이 "IPL_DEPTH_32F"로 생성했던 입력 이미지 (iplDoGx)를 바꾸어 "8" 비트 depth짜리 새로운 이미지 (iplEdgeY)에 저장한다.

            cvConvert(iplDoGx, iplEdgeY);

두번째 인자 " void* storage" 는 탐지된 직선을 저장할 메모리. 이 함수의 아웃풋에 해당한다.

CvMemStorage¶

Growing memory storage.

typedef struct CvMemStorage

{

    struct CvMemBlock* bottom;/* first allocated block */

    struct CvMemBlock* top; /* the current memory block - top of the stack */

    struct CvMemStorage* parent; /* borrows new blocks from */

    int block\_size; /* block size */

    int free\_space; /* free space in the \texttt{top} block (in bytes) */

} CvMemStorage;

CvMemStorage* cvCreateMemStorage(int blockSize=0)¶

Creates memory storage.

Parameter:	blockSize – Size of the storage blocks in bytes. If it is 0, the block size is set to a default value - currently it is about 64K.

 그 아웃풋을 다음의 CvSeq 형태의 자료 구조체 안에 저장한다.

CvSeq¶

Growable sequence of elements.

#define CV_SEQUENCE\_FIELDS() \

    int flags; /* micsellaneous flags */ \

    int header_size; /* size of sequence header */ \

    struct CvSeq* h_prev; /* previous sequence */ \

    struct CvSeq* h_next; /* next sequence */ \

    struct CvSeq* v_prev; /* 2nd previous sequence */ \

    struct CvSeq* v_next; /* 2nd next sequence */ \

    int total; /* total number of elements */ \

    int elem_size;/* size of sequence element in bytes */ \

    char* block_max;/* maximal bound of the last block */ \

    char* ptr; /* current write pointer */ \

    int delta_elems; /* how many elements allocated when the sequence grows

                  (sequence granularity) */ \

    CvMemStorage* storage; /* where the seq is stored */ \

    CvSeqBlock* free_blocks; /* free blocks list */ \

    CvSeqBlock* first; /* pointer to the first sequence block */



typedef struct CvSeq

{

    CV_SEQUENCE_FIELDS()

} CvSeq;

The structure CvSeq is a base for all of OpenCV dynamic data structures.

그 저장된 값을 읽는 함수

char* cvGetSeqElem(const CvSeq* seq, int index)¶

Returns a pointer to a sequence element according to its index.

#define CV_GET_SEQ_ELEM( TYPE, seq, index )  (TYPE*)cvGetSeqElem( (CvSeq*)(seq), (index) )

Parameters:	seq – Sequence index – Index of element

accumulator value 란?

/* Test: line fitting in implementing a virtual studio
 using cvHoughLines2() function in OpenCV library
 ref.
 1) swPark_2000rti.pdf
 2) 박승우_1999대한전자공학회지 제36권 S편 제7호
 3) http://opencv.willowgarage.com/documentation/feature_detection.html?highlight=cvhoughlines#cvHoughLines2
 camera: Logitech QuickCam Pro 4000
 2010, lym
 */

#include <OpenCV/OpenCV.h>
#include <iostream>
using namespace std;

// non-maximum suppression (NMS)
void nonMaximumSuppression ( IplImage* image, int kernel, int threshold )
{
    for ( int y = 0; y < image->height; y++ )
    {
        //        cout << "y = " << y << endl;
        for ( int x = 0; x < image->width; x++ )
        {
            float intensity = CV_IMAGE_ELEM( image, float, y, x );
            if ( intensity > threshold ) {
                float neighbor;
                int flag = 0;
               
                for ( int ky = -kernel; ky <= kernel; ky++ ) // in y-direction
                {
                    if ( y+ky < 0 || y+ky >= image->height ) { // border check
                        continue;
                    }
                    for ( int kx = -kernel; kx <= kernel; kx++ ) // in x-direction
                    {
                        if ( x+kx < 0 || x+kx >= image->width ) {  // border check
                            continue;
                        }
                        neighbor = CV_IMAGE_ELEM( image, float, y+ky, x+kx );
                        if ( intensity < neighbor ) {
                            CV_IMAGE_ELEM( image, float, y, x ) = 0.0;
                            flag = 1;
                            break;
                        }
                    }
                    if ( 1 == flag ) {
                        break;
                    }
                }
            }
           
            else {
                CV_IMAGE_ELEM( image, float, y, x ) = 0.0;
            }
        } 
    }
}

// non-maximum suppression (NMS)
void nonMaximumSuppression2 ( IplImage* image, IplImage* image2, int kernel)
{
    float neighbor, neighbor2;
    for ( int y = 0; y < image->height; y++ )
    {
        //        cout << "y = " << y << endl;
        for ( int x = 0; x < image->width; x++ )
        {
            float intensity = CV_IMAGE_ELEM( image, float, y, x );
            //            if ( intensity > threshold ) {
            if (intensity > 0) {
                int flag = 0;
               
                for ( int ky = -kernel; ky <= kernel; ky++ ) // in y-direction
                {
                    if ( y+ky < 0 || y+ky >= image->height ) { // border check
                        continue;
                    }
                    for ( int kx = -kernel; kx <= kernel; kx++ ) // in x-direction
                    {
                        if ( x+kx < 0 || x+kx >= image->width ) {  // border check
                            continue;
                        }
                        neighbor = CV_IMAGE_ELEM( image, float, y+ky, x+kx );
                        neighbor2 = CV_IMAGE_ELEM( image2, float, y+ky, x+kx );
                        //                        if ( intensity < neighbor ) {
                        if ( intensity < neighbor || intensity < neighbor2) {
                            CV_IMAGE_ELEM( image, float, y, x ) = 0.0;
                            flag = 1;
                            break;
                        }
                    }
                    if ( 1 == flag ) {
                        break;
                    }
                }
            }
           
            else {
                CV_IMAGE_ELEM( image, float, y, x ) = 0.0;
            }
        } 
    }
}


// in the gradient direction
void selectEdges( IplImage* image1, IplImage* image2 )
{
    for ( int y = 0; y < image1->height; y++ )
    {
        //        cout << "y = " << y << endl;
        for ( int x = 0; x < image1->width; x++ )
        {
           
            if( x == image1->width - 1 ) {
                x = x;
            }
           
            float intensity1 = CV_IMAGE_ELEM( image1, float, y, x );
            if ( intensity1 > 0.0 ) { // if the pixel is a edge point surviving NMS
                float intensity2 = CV_IMAGE_ELEM( image2, float, y, x );
                // compare it with the gradient value in the other direction
                if ( intensity1 < intensity2 ) {
                    CV_IMAGE_ELEM( image1, float, y, x ) = 0.0;
                }
            }
        }
    }
}


// draw found lines
void drawLines ( CvSeq* lines, IplImage* image )
{
    for( int i = 0; i < MIN(lines->total,100); i++ )
    {
        float* line = (float*)cvGetSeqElem(lines,i);
        float rho = line[0];
        float theta = line[1];
        CvPoint pt1, pt2;
        double a = cos(theta), b = sin(theta);
        double x0 = a*rho, y0 = b*rho;
        pt1.x = cvRound(x0 + 1000*(-b));
        pt1.y = cvRound(y0 + 1000*(a));
        pt2.x = cvRound(x0 - 1000*(-b));
        pt2.y = cvRound(y0 - 1000*(a));
        // cvLine(<#CvArr * img#>, <#CvPoint pt1#>, <#CvPoint pt2#>, <#CvScalar color#>, <#int thickness#>, <#int line_type#>, <#int shift#>)
        cvLine( image, pt1, pt2, CV_RGB(255,0,0), 1, 8 );
    }
}


int main()
{
    IplImage* iplInput = 0; // input image
    IplImage* iplGray = 0; // grey image converted from input image
    IplImage *iplTemp = 0; // converted image from input image with a change of bit depth 
//    IplImage* iplDoG = 0; // filtered image by DoG
    IplImage* iplDoGx = 0, *iplDoGxClone; // filtered image by DoG in x-direction
    IplImage* iplDoGy = 0, *iplDoGyClone; // filtered image by DoG in y-direction
   
    double minValx, maxValx;
    double minValy, maxValy;
    double minValt, maxValt; 
   
    //    double threshold = 100.0;
    double threshold = 0.0;
    int kernel = 1;
   
    double h[] = { -1, -7, -15, 0, 15, 7, 1 };
   
    CvMat DoGx = cvMat( 1, 7, CV_64FC1, h );
    CvMat* DoGy = cvCreateMat( 7, 1, CV_64FC1 );
    cvTranspose( &DoGx, DoGy ); // transpose(&DoGx) -> DoGy
   
    char title_fx[200], title_fy[200];
    sprintf(title_fx, "filtered image by DoGx");
    sprintf(title_fy, "filtered image by DoGy"); 
   
    cvNamedWindow("input");
    cvNamedWindow( "temp" );
    cvNamedWindow(title_fx);
    cvNamedWindow(title_fy);
    //    cvNamedWindow(title_f);
   
    char title_ex[200], title_ey[200];
    sprintf(title_ex, "detected edges in x direction");
    sprintf(title_ey, "detected edges in y direction"); 
   
    cvNamedWindow(title_ex);
    cvNamedWindow(title_ey);
    //    cvNamedWindow(title_e);
   
    // Hough transform
    IplImage* iplEdgeX = 0, *iplEdgeY = 0; // edge-detected image in each direction
    CvMemStorage* storageX = cvCreateMemStorage(0), *storageY = cvCreateMemStorage(0);
    CvSeq* linesX = 0, *linesY = 0;
   
   
   
    // initialize capture from a camera
    CvCapture* capture = cvCaptureFromCAM(0); // capture from video device #0
    int count = 0; // number of grabbed frames
   
    while(1) {
        // get video frames from the camera
        //        if (0) {
        if ( !cvGrabFrame(capture) ) {
            printf("Could not grab a frame\n\7");
            exit(0);
        }
        else {
            cvGrabFrame( capture ); // capture a frame
            iplInput = cvRetrieveFrame(capture); // retrieve the caputred frame
            //            iplInput = cvLoadImage("P:/input.bmp"); // retrieve the caputred frame
           
            if(iplInput) {
                if(0 == count) {   
                    // create an image header and allocate the image data
                    /*                   iplGray = cvCreateImage(cvGetSize(iplInput), 8, 1);
                     iplDoGx = cvCreateImage(cvGetSize(iplInput), 8, 1);
                     iplDoGy = cvCreateImage(cvGetSize(iplInput), 8, 1);              
                     */                 
                    iplGray = cvCreateImage(cvGetSize(iplInput), 8, 1);
                    iplTemp = cvCreateImage(cvGetSize(iplInput), IPL_DEPTH_32F, 1);                  
                    iplDoGx = cvCreateImage(cvGetSize(iplInput), IPL_DEPTH_32F, 1); 
                    iplDoGy = cvCreateImage(cvGetSize(iplInput), IPL_DEPTH_32F, 1); 
                    iplDoGyClone = cvCloneImage(iplDoGy), iplDoGxClone = cvCloneImage(iplDoGx);
                    iplEdgeX = cvCreateImage(cvGetSize(iplInput), 8, 1);
                    iplEdgeY = cvCreateImage(cvGetSize(iplInput), 8, 1);                   
                }
                // convert the input color image to gray one
                cvCvtColor(iplInput, iplGray, CV_BGR2GRAY); // convert an image from one color space to another
                // convert one array to another with optional linear transformation
                cvConvert(iplGray, iplTemp);
                // increase the frame number
                count++;
            }       
            //            cvShowImage( "input", iplInput ); 
           
            // convolve an image with the kernel
            // void cvFilter2D(const CvArr* src, CvArr* dst, const CvMat* kernel, CvPoint  anchor=cvPoint(-1, -1)
            cvFilter2D( iplTemp, iplDoGx, &DoGx ); // convolve an image with the DoG kernel in x-direction
            cvFilter2D( iplTemp, iplDoGy, DoGy ); // convolve an image with the DoG kernel in y-direction         
            // ref. http://opencv.willowgarage.com/documentation/operations_on_arrays.html?highlight=cvabs#cvAbsDiffS
            cvAbs(iplDoGx, iplDoGx);            cvAbs(iplDoGy, iplDoGy);
           
            // normalize the pixel values
            // http://opencv.willowgarage.com/documentation/operations_on_arrays.html?highlight=cvminmax#cvMinMaxLoc
            // void cvMinMaxLoc(const CvArr* arr, double* minVal, double* maxVal, CvPoint* minLoc=NULL, CvPoint* maxLoc=NULL, const CvArr* mask=NULL)¶
            cvMinMaxLoc( iplDoGx, &minValx, &maxValx );
            cvMinMaxLoc( iplDoGy, &minValy, &maxValy );
            cvMinMaxLoc( iplTemp, &minValt, &maxValt );
            cvScale( iplDoGx, iplDoGx, 1.0 / maxValx );
            cvScale( iplDoGy, iplDoGy, 1.0 / maxValy );
            cvScale( iplTemp, iplTemp, 1.0 / maxValt );
           
            // display windows
            cvShowImage( "temp", iplTemp ); 
            cvShowImage( title_fx, iplDoGx );
            cvShowImage( title_fy, iplDoGy );         
           
            // save images to files
            cvSaveImage("input.bmp", iplInput);
            cvSaveImage("gray.bmp", iplGray);
            cvSaveImage("temp.bmp", iplTemp);
            cvSaveImage("DoGx.bmp", iplDoGx);
            cvSaveImage("DoGy.bmp", iplDoGy);
           
            // non-Maximum suppression (NMS)
            // void NonMaximumSuppression ( IplImage* image, int kernel, int threshold )
            //            nonMaximumSuppression ( iplDoGx, kernel, threshold ); nonMaximumSuppression ( iplDoGy, kernel, threshold );
            //           cvCopyImage(iplDoGy, iplDoGyClone), cvCopyImage(iplDoGx, iplDoGxClone);
            cvCopy(iplDoGy, iplDoGyClone), cvCopy(iplDoGx, iplDoGxClone);
           
            nonMaximumSuppression2 ( iplDoGx, iplDoGyClone, kernel );
            nonMaximumSuppression2 ( iplDoGy, iplDoGxClone, kernel );
           
            // normalize the pixel values
            // http://opencv.willowgarage.com/documentation/operations_on_arrays.html?highlight=cvminmax#cvMinMaxLoc
            // void cvMinMaxLoc(const CvArr* arr, double* minVal, double* maxVal, CvPoint* minLoc=NULL, CvPoint* maxLoc=NULL, const CvArr* mask=NULL)¶
            cvMinMaxLoc( iplDoGx, &minValx, &maxValx );
            cvMinMaxLoc( iplDoGy, &minValy, &maxValy );
            cvMinMaxLoc( iplTemp, &minValt, &maxValt );
            cvScale( iplDoGx, iplDoGx, 1.0 / maxValx );
            cvScale( iplDoGy, iplDoGy, 1.0 / maxValy );
            cvScale( iplTemp, iplTemp, 1.0 / maxValt );
           
            // display windows
            cvShowImage( title_ex, iplDoGx );
            cvShowImage( title_ey, iplDoGy );         
           
            // save images to files
            cvSaveImage("edgeX.bmp", iplDoGx);
            cvSaveImage("edgeY.bmp", iplDoGy); 
           
            cvConvert(iplDoGx, iplEdgeY);
            cvConvert(iplDoGy, iplEdgeX);
           
            // line fitting
            // ref. http://opencv.willowgarage.com/documentation/feature_detection.html?highlight=cvhoughlines#cvHoughLines2
            // CvSeq* cvHoughLines2(CvArr*  image, void* storage, int method, double rho, double theta, int threshold, double param1=0, double param2=0)
            double rho = 1.0; // distance resolution in pixel-related units
            double theta = 1.0; // angle resolution measured in radians
            int threshold = 40; // ("A line is returned by the function if the corresponding accumulator value is greater than threshold.")
            linesX = cvHoughLines2(iplEdgeX, storageX, CV_HOUGH_STANDARD, 1.0*rho, CV_PI/180*theta, 20, 0, 0);
//            linesX = cvHoughLines2(iplLineX, storage, CV_HOUGH_PROBABILISTIC, 2.0, CV_PI/180*1.0, 10, 1.0, 1.0);
            linesY = cvHoughLines2(iplEdgeY, storageY, CV_HOUGH_STANDARD, 1.0*rho, CV_PI/180*theta, 20, 0, 0);
           
            cout << "# of found lines = " << linesX->total << " horizontal, " << linesY->total << " vertical " << endl;
           
            // draw found lines
            drawLines ( linesX, iplInput );
            drawLines ( linesY, iplInput );           
           
            cvShowImage( "input", iplInput ); 
           
            if( cvWaitKey(10) >= 0 )
                break;
        }
    }
   
    cvReleaseCapture( &capture ); // release the capture source
    cvDestroyWindow( "input" );
    cvDestroyWindow( "temp" );
    cvDestroyWindow(title_fx);
    cvDestroyWindow(title_fy);
    cvDestroyWindow(title_ex);
    cvDestroyWindow(title_ey);
   
    return 0;
}

"detected edges" 이미지에 대해 Hough transform에 의한 line fitting 한 결과를 "input" 이미지에 그리고 있음

opencv/opencv/src/cv/cvhough.cpp 를 열면, 다음의 네 부분으로 나뉘어 정의되어 있다.

Classical Hough Transform
Multi-Scale variant of Classical Hough Transform 
Probabilistic Hough Transform      
Circle Detection

이 중 "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 );
}