virtual studio 구현: line fitting test

overview:
렌즈의 왜곡 현상 때문에 이미지 상에 검출된 edge points들은 (직선에서 왜곡된) 2차 곡선을 그리게 된다. 이 곡선의 방정식을 먼저 구한 후에, 렌즈 왜곡 변수를 "0"으로 두고 나오는 직선들로부터 비로소 cross-ratio 값이 보존된다.

ref.  2010/02/10 - [Visual Information Processing Lab] - Seong-Woo Park & Yongduek Seo & Ki-Sang Hong

swPark_2000rti 440쪽: "The cross-ratio is not preserved for the (image) frame coordinate, positions of the feature points in an image, or for the distorted image coordinate. Cross-ratio is invariant only for the undistorted coordinate." (swPark_20

박승우_1999전자공학회지 96쪽: "이렇게 곡선으로 나타난 가로선과 세로선을 직선으로 피팅할 경우 cross-ratio는 왜곡 현상 때문에 이 선들에 대해서는 보존되지 않게 된다. 따라서 정확한 피팅을 위해서는 아래와 같이 렌즈의 왜곡변수(k1)를 고려한 이차곡선으로의 피팅이 필요하다.

Y = a*X + b/(1+k1*R^2) = a*X + b/(1+k1*(X^2+Y^2)) <--- 이 식은 영어 논문 (19)식과 한글 논문 (15)식을 조합, 수정한 식임. 확인 필요

이 식을 피팅해서 계수 a, b, k1를 구하고, 여기서 k1=0을 두면 왜곡이 보상된 점에 대한 직선식을 구할 수 있다. 이렇게 구해진 직선들을 패턴의 가로선들과 세로선들의 cross-ratio와 비교함으로써 영상에서 찾아진 선들을 인식할 수 있다. 또한 영상에서의 특징점은 이 식에 의해 피팅된 가로선들과 세로선들의 교점으로 정확하게 구할 수 있다."

그런데, 
현재 시험용 패턴과 코드로부터 촬영, 검출된 이미지 상의 점들은 거의 직선에 가깝다. 우선 OpenCV 라이브러리의 cvHoughLines2() 함수에 의한 직선 찾기를 해 보자.

2010/04/07 - [Visual Information Processing Lab] - OpenCV: cvHoughLines2() 연습 코드


1) 교점 구하기 테스트
line fitting을 통해 찾은 직선들로부터 패턴 격자의 corner points를 구하는 것을 시험해 본다.

/* 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 );
    }
}

// calculate intersection points of lines in X and Y to find corner points
void findIntersections ( CvSeq* linesX, CvSeq* linesY, IplImage* image )
// find intersections of lines, "linesX" & "linesY", and draw them in "image"
{
    for( int i = 0; i < MIN(linesX->total,100); i++ )
    {
        float* lineX = (float*)cvGetSeqElem(linesX,i);
        for( int j = 0; j < MIN(linesY->total,100); j++ )
        {
            float* lineY = (float*)cvGetSeqElem(linesY,j);               
            float rhoX = lineX[0];
            float rhoY = lineY[0];
            float thetaX = lineX[1], thetaY = lineY[1];

            double aX = cos(thetaX), bX = sin(thetaX);
            double aY = cos(thetaY), bY = sin(thetaY);

            CvPoint c; // the intersection point of lineX[i] and lineY[j]
            double Cx = ( rhoX*bY - rhoY*bX ) / ( aX*bY - bX*aY ) ;
            double Cy = ( rhoX - aX*Cx ) / bX ;
            c.x = cvRound(Cx);
            c.y = cvRound(Cy);
           
            // display the points in an image
            cvCircle( image, c, 3, CV_RGB(0,255,100), 2 /* , <#int line_type#>, <#int shift#> */ );
        }
    }       
}

int width = 320;
int height = 240;

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* iplDoGx = 0, *iplDoGxClone; // filtered image by DoG in x-direction
    IplImage* iplDoGy = 0, *iplDoGyClone; // filtered image by DoG in y-direction
    IplImage* iplEdgeX = 0, *iplEdgeY = 0; // edge-detected image in each direction, to be used as input in line-fitting
   
    double minValx, maxValx;
    double minValy, maxValy;
    double minValt, maxValt; 
   
    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
   
    // create windows
    cvNamedWindow("input");
    cvNamedWindow( "temp" );
   
    char title_fx[200], title_fy[200];
    sprintf(title_fx, "filtered image by DoGx");
    sprintf(title_fy, "filtered image by DoGy"); 
   
    cvNamedWindow(title_fx);
    cvNamedWindow(title_fy);
   
    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);
   
    cvMoveWindow( "temp", 110+width, 100 );
    cvMoveWindow( title_fx, 100, 100+height+30 );
    cvMoveWindow( title_fy, 110+width, 100+height+30 );
    cvMoveWindow( title_ex, 100, 100+(height+30)*2 );
    cvMoveWindow( title_ey, 110+width, 100+(height+30)*2 );
   
   
    // output information of lines found by Hough transform
    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 ( !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);
                    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("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); 
           
            // line fitting
            cvConvert(iplDoGx, iplEdgeY);
            cvConvert(iplDoGy, iplEdgeX);
           
            // 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 = 20; // ("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, threshold, 0, 0);
            linesY = cvHoughLines2(iplEdgeY, storageY, CV_HOUGH_STANDARD, 1.0*rho, CV_PI/180*theta, threshold, 0, 0);
           
//            linesX = cvHoughLines2(iplEdgeX, storageX, CV_HOUGH_PROBABILISTIC, 1.0*rho, CV_PI/180*1.0, 10, 1.0, 1.0);
//            linesY = cvHoughLines2(iplEdgeY, storageY, CV_HOUGH_PROBABILISTIC, 1.0*rho, CV_PI/180*theta, 10, 1.0, 1.0);
           
            cout << "# of found lines = " << linesY->total << " vertical " << linesX->total << " horizontal, " << endl;
           
            // draw found lines
            drawLines ( linesX, iplInput );
            drawLines ( linesY, iplInput );           
           
            // find intersections of lines
            // void findIntersections ( CvSeq* linesX, CvSeq* linesY, IplImage* image )
            findIntersections( linesX, linesY, iplInput );
           
            cvShowImage( "input", iplInput ); 
            cvSaveImage( "input.bmp", 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;
}

실시간으로 산출하는 데 무리가 없음이 확인되었다.

2)
그러나, line fitting의 결과가 깔끔하지 않은 문제를 우선 해결해야 한다. (rho, theta, threshold 등의 함수 매개변수 값을 조정하는 것을 포함하여 사용 중인 웹캠에 적합한 데이터 처리가 필요하다.)

현재의 코드로부터 나오는 결과를 정리해 두면 아래와 같다.

NMS와 동시에 수평선 또는 수직선 위의 점들을 따로 추출한 결과 이미지 ("iplEdgeX"와  "iplEdgeY")를 cvHoughLines2() 함수의 입력으로 하고,

double rho = 1.0; // distance resolution in pixel-related units
double theta = 1.0; // angle resolution measured in radians
int threshold = 20; // (A line is returned by the function if the corresponding accumulator value is greater than threshold)

위와 같이 매개변수 값을 주면 검출된 직선들과 그로부터 계산한 교점들은 다음과 같이 나타난다.

수직선 상의 edges만 검출한 영상

수평선 상의 edges만 검출한 영상

Hough transform에 의한 line fitting 한 결과

(Non Maximal suppression (NMS)을 하기 전에) 1차 DoG 필터를 이미지 프레임의 x 방향, y 방향으로 적용한 결과 이미지 ("iplDoGx"와 "iplDoGy")를 cvHoughLines2() 함수의 입력으로 하고,

double rho = 1.0; // distance resolution in pixel-related units
double theta = 1.0; // angle resolution measured in radians
int threshold = 20; // (A line is returned by the function if the corresponding accumulator value is greater than threshold)

위와 같이 매개변수 값들을 주면 검출된 직선들과 그로부터 계산한 교점들은 다음과 같이 나타난다.

x방향으로 DoG 필터를 적용한 이미지

y방향으로 DoG 필터를 적용한 이미지

Hough transform에 의한 line fitting한 결과

그러니까... 실제로 한 직선 상의 점들로부터 여러 개의 직선을 찾게 되는 것은 edge points로 detection된 (흰색으로 보이는) 픽셀 부분의 세기값이 약하거나 일정하지 않기 때문인 것 같다. 입력 이미지를 binary로 바꾸고 cvHoughLines2()의 입력으로 accumulator value에 기준값을 주는 파라미터 threshold를 증가시키면 될 것 같다.



Try #1. 입력 이미지 이진화

NMS와 동시에 수평선 또는 수직선 위의 점들을 따로 추출한 결과 이미지 ("iplEdgeX"와  "iplEdgeY")를 이진화하고, 

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.")

위와 같이 매개변수 값들을 주면 검출된 직선들과 그로부터 계산한 교점들은 다음과 같이 나타난다.

수직선 상의 edges만 검출하여 이진화한 영상

수평선 상의 edges만 검출하여 이진화한 영상

Hough transform에 의한 line fitting한 결과

실제로 한 직선에 여러 개의 직선이 검출되는 빈도는 현저히 줄지만 대신 실제 직선 중에 검출되지 않는 것이 생긴다.


Try #2. line fitting의 입력 이미지 처리 & 매개변수 조정



Try #3. 실제로 하나인데 여러 개로 겹쳐서 나오는 직선들의 평균을 취해 하나로 합침

다음과 같은 입력 영상에 대해 탐지된 직선들의 방정식을 정의하는 매개변수 (rho와 theta) 값을 출력해 보면 아래와 같이 나온다.

# of found lines = 8 vertical   22 horizontal
vertical
rho = 172.6    theta = 0
rho = 133    theta = 0.139626
rho = -240.2    theta = 2.84489
rho = -209    theta = 2.98451
rho = 91.8    theta = 0.279253
rho = 173.8    theta = 0
rho = 52.6    theta = 0.401426
rho = 53.8    theta = 0.418879
horizontal
rho = 81    theta = 1.55334
rho = 53.4    theta = 1.55334
rho = 155    theta = 1.55334
rho = 114.6    theta = 1.55334
rho = 50.6    theta = 1.5708
rho = 29.8    theta = 1.55334
rho = 76.6    theta = 1.5708
rho = 112.6    theta = 1.5708
rho = 9.8    theta = 1.55334
rho = 152.6    theta = 1.5708
rho = 153.8    theta = 1.5708
rho = 150.6    theta = 1.5708
rho = 6.6    theta = 1.5708
rho = 78.6    theta = 1.5708
rho = 205.4    theta = 1.55334
rho = 27.8    theta = 1.5708
rho = 8.6    theta = 1.5708
rho = 201.8    theta = 1.5708
rho = 110.6    theta = 1.5708
rho = 49.8    theta = 1.5708
rho = 48.6    theta = 1.5708
rho = 111.8    theta = 1.5708

잠시 현재 상태 기록: cross ratios를 이용한 격자 무늬 패턴과 line detection 시험 + feature points matching을 위한 교점 찾기와 순번 부여 시험

to do next:
1) line detection의 error 교정
2) (rhoX, thetaX, rhoY, thetaY로 정의되는) 교점 indexing