돛단배

void matchXY ( vector<CvPoint2D32f> &p, vector<CvPoint2D32f> &CRimage, int numIx, int numIy, vector<CvPoint3D32f> &world, vector<CvPoint2D32f> &CRworld, int numPx, int numPy, IplImage* iplMatch, CvPoint2D32f scale )
{   
    for( int i = 0; i < numIx; i++ ) // points in x-direction on the input image
    {
        // check if x-component of the point is valid
        if( -1 == CRimage[i*numIy+0].x )
        {
            cout << endl << "could not make matching in x-direction" << endl;
            continue;
        }
       
        // CvScalar generateRandomColor(unsigned char thR, unsigned char thG, unsigned char thB) defined in "matching.h"
        CvScalar colorMatch = generateRandomColor(50,50,50); 
       
        for( int j = 0; j < numIy; j++ ) // points in y-direction on the input image
        {       
            // check if y-component of the point is valid
            if( -1 == CRimage[i*numIy+j].y )
            {
                   cout << endl << "could not make matching in y-direction" << endl;
                continue;
            }
           
            // to find the x-index of the corresponding point
            int indexPx = 0;                   
            float errX_min = fabs( CRimage[i*numIy+j].x- CRworld[indexPx*numPy+0].x );
            // search points in x-direction on the real pattern
            for( int wx = 0; wx < numPx; wx++ )
            {
                float errX = CRimage[i*numIy+j].x - CRworld[wx*numPy+0].x;
                if ( fabs(errX) < errX_min )
                {
                    errX_min = fabs(errX);
                    indexPx = wx;
                }
            }
           
            // to find the y-index of the corresponding point
            int indexPy = 0;
            float errY_min = fabs( CRimage[i*numIy+j].y - CRworld[0*numPy+indexPy].y );           
            // search points in y-direction on the real pattern    
            for( int wy = 0; wy < numPy; wy++ )
            {
                float errY = CRimage[i*numIy+j].y - CRworld[0*numPy+wy].y;
                if ( fabs(errY) < errY_min )
                {
                    errY_min = fabs(errY);
                    indexPy = wy;
                }
            }
           
//            cout << endl << i << ", " << j << " point in the input frame is matched with "
//            << indexPy << "-th point in the real pattern" << endl;
           
            // draw the line to connect "world" point and "image" point
            CvPoint pointImage = cvPoint(cvRound(p[i*numIy+j].x), cvRound(IMG_HEIGHT + p[i*numIy+j].y));
            CvPoint pointPattern = cvPoint(cvRound(world[indexPx*numPy+indexPy].x*scale.x), cvRound(world[indexPx*numPy+indexPy].y*scale.y));
           
            cvLine( iplMatch, pointImage, pointPattern, colorMatch, 1 );
            cvCircle( iplMatch, pointImage, 3, colorMatch, 2 );  
            cvCircle( iplMatch, pointPattern, 3, colorMatch, 2 );       
           
        } // end for j
    } // end for i
}

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

#define BLOCK 16

char *fgetline(FILE *fp) {
   
    char *s, *tmp, buf[BLOCK];
    size_t count;
   
    for(count = 0, s = NULL; fgets(buf, sizeof(buf), fp); count++)
    {       
        if((tmp = (char *)realloc(s, (count + 1) * BLOCK)) == NULL) {
            free(s);
            return NULL;
        }       
        s = tmp;       
        if(count == 0) *s = '\0';       
        strcat(s,buf);       
        if((tmp = strrchr(s,'\n')) != NULL) {           
            *tmp = '\0';
            break;
        }
    }
    return s;
}


#define IMG_WIDTH 400
#define IMG_HEIGHT 300

int main (int argc, char * const argv[]) {
    // input file
    //    FILE* fp = fopen(argv[1], "r");
    FILE* fp1 = fopen("bestPatternHorizonal.txt", "r");
    FILE* fp2 = fopen("bestPatternVertical.txt", "r");   

    vector<double> crList1, crList2;
    char *mystring1, *mystring2;
   
    double width = 0.0, widthSum = 0.0;
    double height = 0.0, heightSum = 0.0;
    double scaleX = 0.0, scaleY = 0.0;
   
    IplImage* img = cvCreateImage(cvSize(IMG_WIDTH, IMG_HEIGHT), 8, 1); // image to draw pattern on
    cvNamedWindow("pattern");
   
    while( mystring1 = fgetline(fp1) )
    {
        width = atof(mystring1);
        crList1.push_back(width);
        widthSum += width;
    }
    scaleX = IMG_WIDTH / widthSum;
    while( mystring2 = fgetline(fp2) )
    {
        height = atof(mystring2);
        crList2.push_back(height);
        heightSum += height;
    }   
    scaleY = IMG_HEIGHT / heightSum;
   
/*    for(int i = 0; i < crList1.size(); i++) {
        printf("%f\n", crList1[i]);
    }
*/
   
    cvZero(img);
    CvPoint2D32f p;
    CvPoint2D32f l;
   
    p.x = 0.0, p.y = 0.0;
    l.x = 0.0, l.y = 0.0;
   
    cout << "totla # = " << (crList1.size()-1) * (crList2.size()-1) << endl;
   
    for( int i = 0; i < crList1.size(); i++ )       
    {
        l.x = crList1[i];       
        p.y = 0.0;

        for( int j = 0; j < crList2.size(); j++ )
        {
            int number =  i*(crList2.size()-1) + j;
            cout << "# " << number <<  endl;
           
            l.y = crList2[j];
           
            if( 0 == number % 2) {       
                CvRect grid = cvRect( cvRound(p.x), cvRound(p.y), cvRound(l.x), cvRound(l.y) ); // cvRect(<#int x#>, <#int y#>, <#int width#>, <#int height#>)
                CvRect pattern = cvRect( cvRound(p.x*scaleX), cvRound(p.y*scaleY), cvRound(l.x*scaleX), cvRound(l.y*scaleY) );
               
                cvSetImageROI(img, pattern);
                cvSet(img, cvScalarAll(255));
                cvResetImageROI(img);
//                cvRectangle(<#CvArr * img#>, <#CvPoint pt1#>, <#CvPoint pt2#>, <#CvScalar color#>, <#int thickness#>, <#int line_type#>, <#int shift#>)(<#int x#>, <#int y#>, <#int width#>, <#int height#>);
                p.y = p.y + crList2[j];   
            }
           
            else {
                p.y = p.y + crList2[j];   
                continue;
            }
        } // end for j
        p.x = p.x + crList1[i];       
    } // end for i
   
    cvShowImage("pattern", img);   
    cvSaveImage("pattern.bmp", img);

/*  Draw grid pattern
 using kyu's pattern generator with optimal cross ratios
 2010, lym
 
 cf.
 cvRectangle(<#CvArr * img#>, <#CvPoint pt1#>, <#CvPoint pt2#>, <#CvScalar color#>, <#int thickness#>, <#int line_type#>, <#int shift#>)(<#int x#>, <#int y#>, <#int width#>, <#int height#>);
 cvRect(<#int x#>, <#int y#>, <#int width#>, <#int height#>)
 */

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

#define BLOCK 16

char *fgetline(FILE *fp) {
   
    char *s, *tmp, buf[BLOCK];
    size_t count;
   
    for(count = 0, s = NULL; fgets(buf, sizeof(buf), fp); count++)
    {       
        if((tmp = (char *)realloc(s, (count + 1) * BLOCK)) == NULL) {
            free(s);
            return NULL;
        }       
        s = tmp;       
        if(count == 0) *s = '\0';       
        strcat(s,buf);       
        if((tmp = strrchr(s,'\n')) != NULL) {           
            *tmp = '\0';
            break;
        }
    }
    return s;
}


#define IMG_WIDTH 400
#define IMG_HEIGHT 300

int main (int argc, char * const argv[]) {
    // input file
    //    FILE* fp = fopen(argv[1], "r");
    FILE* fp1 = fopen("bestPatternHorizonal.txt", "r");
    FILE* fp2 = fopen("bestPatternVertical.txt", "r");   
   
    vector<double> crList1, crList2;
    char *mystring1, *mystring2;
   
    double width = 0.0, widthSum = 0.0;
    double height = 0.0, heightSum = 0.0;
    double scaleX = 0.0, scaleY = 0.0;
   
    IplImage* img = cvCreateImage(cvSize(IMG_WIDTH, IMG_HEIGHT), 8, 1); // image to draw pattern on
    cvNamedWindow("pattern");
   
    while( mystring1 = fgetline(fp1) )
    {
        width = atof(mystring1);
        crList1.push_back(width);
        widthSum += width;
    }
    scaleX = IMG_WIDTH / widthSum;
    while( mystring2 = fgetline(fp2) )
    {
        height = atof(mystring2);
        crList2.push_back(height);
        heightSum += height;
    }   
    scaleY = IMG_HEIGHT / heightSum;
   
    cvZero(img);
    CvPoint2D32f p;
    CvPoint2D32f l;
   
    p.x = 0.0, p.y = 0.0;
    l.x = 0.0, l.y = 0.0;
   
//    cout << "totla # = " << (crList1.size()-1) * (crList2.size()-1) << endl;
   
    // draw the vertical bars in white
    for( int i = 0; i < crList1.size(); i++ )       
    {
        l.x = crList1[i];                   
        if ( 0 == i%2 ) {
            CvRect patternX = cvRect( cvRound(p.x*scaleX), 0, cvRound(l.x*scaleX), IMG_HEIGHT );
            cvSetImageROI(img, patternX);
            cvSet(img, cvScalarAll(255));
            cvResetImageROI(img);
        }
        p.x = p.x + crList1[i];               
    } // end for i
   
    // draw the horizontal bars in white
     for( int j = 0; j < crList2.size(); j++ )
     {
         l.y = crList2[j];         
         if ( 0 == j%2 ) {
             CvRect patternY = cvRect( 0, cvRound(p.y*scaleY), IMG_WIDTH, cvRound(l.y*scaleY) );
             cvSetImageROI(img, patternY);
             cvSet(img, cvScalarAll(255));
             cvResetImageROI(img);
         }
         p.y = p.y + crList2[j];       
     } // end for j
        
    // fill the overlapped parts in black
    p.x = 0.0;
    for( int i = 0; i < crList1.size(); i++ )
    {
        l.x = crList1[i];
        p.y = 0.0;
        for( int j = 0; j < crList2.size(); j++ )
        {
            l.y = crList2[j];
            if ( 0 == i%2 && 0 == j%2 )
            {
                CvRect overlap = cvRect( cvRound(p.x*scaleX), cvRound(p.y*scaleY), cvRound(l.x*scaleX), cvRound(l.y*scaleY) );
                cvSetImageROI(img, overlap);
                cvSet(img, cvScalarAll(0));
                cvResetImageROI(img);
            }
            p.y = p.y + crList2[j];
        } // end for j
        p.x = p.x + crList1[i];   
    } // end for i

    cvShowImage("pattern", img);   
    cvSaveImage("pattern.bmp", img);
    cvWaitKey(0);
   
    return 0;
}

/* Test: feature points identification in implementing a virtual studio
 1) grid pattern design with cross ratios
 2) lines detection by Hough transform with Non Maximum Suppression,
 modifying 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
 4) http://leeway.tistory.com/813
 camera: Logitech QuickCam Pro 4000
 2010, lym
 */


#include <OpenCV/OpenCV.h> // framework on Mac
//#include <cv.h>
//#include <highgui.h>
//#include <cxmisc.h>
#include <iostream>
#include <vector>
using namespace std;

//#include "nms.h" // Non Maximum Suppression to extract vertical and horizontal edges separately
//#include "nmshough.h" // Hough transform with Non Maximum Suppression to detect lines


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

struct line_param // structure to contain parameters to define a line
{
    // eqn of a line: a*x + b*y = rho, when a = cos(theta) & b = sin(theta)
    float rho, theta;
    float a, b;
    float x0, y0; // x-intercept, y-intercepts
};


// rearrange vertical lines from left to right
void indexLinesY (vector<line_param>& line, CvSeq* lines, IplImage* image )
{   
    int numLines = lines->total; // total number of detected lines
    line.resize(numLines); // define the size of "line" vector as "numLines"
    char txt[100]; // text to represent the index number of an ordered line
   
    // get "rho" and "theta" values of lines detected by Hough transform and NMS
    for( int n = 0; n < numLines; n++ )
    {
        float* newline = (float*)cvGetSeqElem(lines,n);
        line[n].rho = newline[0];
        line[n].theta = newline[1];
    }
   
    // rearrange "line" array in geometrical order, that is, by values of x-intercept in the image frame coordinate
    float temp_rho, temp_theta;
    for( int n = 0; n < numLines-1; n++ )
    {
        for ( int k = n+1; k < numLines; k++ )
        {
            float x0_here = line[n].rho / cos(line[n].theta);
            float x0_next = line[k].rho / cos(line[k].theta);
           
            if( x0_here > x0_next ) {
                temp_rho = line[n].rho;        temp_theta = line[n].theta;
                line[n].rho = line[k].rho;        line[n].theta = line[k].theta;
                line[k].rho = temp_rho;        line[k].theta = temp_theta;
            }
        }
    }
   
    // calculate the other parameters of the rearranged lines
    for( int n = 0; n < numLines; n++ )
    {
        line[n].a = cos(line[n].theta);
        line[n].b = sin(line[n].theta);
        line[n].x0 = line[n].rho / line[n].a;
        line[n].y0 = line[n].rho / line[n].b;
       
        // display the ordered vertical lines
        cout << "x[" << n << "] = " << line[n].x0 << "    y[" << n << "] = " << line[n].y0 ;
        cout << "    rho[" << n << "] = " << line[n].rho << "    theta[" << n << "] = " << line[n].theta << endl;
       
        sprintf(txt, "%d", n);
        cvPutText(image, txt, cvPoint(line[n].x0, 10+n*10), &cvFont(0.8), CV_RGB(255,50,50));
    }
    return;
}


// rearrange horizontal lines from up to bottom
void indexLinesX (vector<line_param>& line, CvSeq* lines, IplImage* image )
{     
    int numLines = lines->total; // total number of detected lines
    line.resize(numLines); // define the size of "line" vector as "numLines"
    char txt[100]; // text to represent the index number of an ordered line
   
    // get "rho" and "theta" values of lines detected by Hough transform and NMS
    for( int n = 0; n < numLines; n++ )
    {
        float* newline = (float*)cvGetSeqElem(lines,n);
        line[n].rho = newline[0];
        line[n].theta = newline[1];
    }
   
    // rearrange "line" array in geometrical order, that is, by values of y-intercept in the image frame coordinate
    float temp_rho, temp_theta;
    for( int n = 0; n < numLines-1; n++ )
    {
        for ( int k = n+1; k < numLines; k++ )
        {
            float y0_here = line[n].rho / sin(line[n].theta);
            float y0_next = line[k].rho / sin(line[k].theta);
           
            if( y0_here > y0_next ) {
                temp_rho = line[n].rho;        temp_theta = line[n].theta;
                line[n].rho = line[k].rho;        line[n].theta = line[k].theta;
                line[k].rho = temp_rho;        line[k].theta = temp_theta;
            }
        }
    }
   
    // calculate the other parameters of the rearranged lines
    for( int n = 0; n < numLines; n++ )
    {
        line[n].a = cos(line[n].theta);
        line[n].b = sin(line[n].theta);
        line[n].x0 = line[n].rho / line[n].a;
        line[n].y0 = line[n].rho / line[n].b;
       
        // display the ordered horizontal lines
        cout << "y[" << n << "] = " << line[n].y0 << "    x[" << n << "] = " << line[n].x0 ;
        cout << "    rho[" << n << "] = " << line[n].rho << "    theta[" << n << "] = " << line[n].theta << endl;
   
        sprintf(txt, "%d", n);
        cvPutText(image, txt, cvPoint(10+n*10, line[n].y0), &cvFont(0.8), CV_RGB(50,50,250));
    }
    return;
}


// index intersection points of lines in X and Y
void indexIntersections (vector<CvPoint>& p, vector<line_param>& lineX, vector<line_param>& lineY, IplImage* image )
{ // find "p", intersections of lines, "linesX" & "linesY", and draw them in "image"   
    int numLinesX = lineX.size(), numLinesY = lineY.size(); // total number of detected lines
    int numPoints = numLinesX * numLinesY; // total number of intersection of the lines
    p.resize(numPoints); // define the size of "p" vector as "numPoints"
    char txt[100]; // text to represent the index number of an intersection point
   
    // calculate intersection points of vertical and horizontal lines
    for( int i = 0; i < numLinesX; i++ )
    {
        for( int j = 0; j < numLinesY; j++ )
        {
            int indexP = numLinesY * i + j; // index number of intersection points in geometrical order
            float Px = ( lineX[i].rho*lineY[j].b - lineY[j].rho*lineX[i].b ) / ( lineX[i].a*lineY[j].b - lineX[i].b*lineY[j].a ) ;
            float Py = ( lineX[i].rho - lineX[i].a*Px ) / lineX[i].b ;
            p[indexP].x = cvRound(Px);
            p[indexP].y = cvRound(Py);
           
            // display the points in an image
            cvCircle( image, p[indexP], 3, CV_RGB(0,255,50) /* , <#int line_type#>, <#int shift#> */ );   
            sprintf(txt, "%d-%d", i, j);    cvPutText(image, txt, p[indexP], &cvFont(0.7), CV_RGB(50,255,250));           
        }
    }       
    return;
}


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 by filtering in each direction, to be used as input in line-fitting
   
    double minValx, maxValx; // minimum & maximum of pixel intensity values
    double minValy, maxValy;
    double minValt, maxValt; 
   
    int width, height; // window size of input frame 
    int kernel = 1; float edgethres; // parameters of NMS function
   
    double rho = 0.8; // distance resolution in pixel-related units
    double theta = 0.8; // angle resolution measured in radians
    // "A line is returned by the function if the corresponding accumulator value is greater than threshold."
    //    int threshold = 24, rN = 5, tN = 5; // for grid pattern of 11x7 squares
    int threshold = 20, rN = 5, tN = 5; // for grid pattern of lines with cross ratios    
   
    double h[] = { -1, -7, -15, 0, 15, 7, 1 }; // 1-D kernel of DoG filter   
    CvMat DoGx = cvMat( 1, 7, CV_64FC1, h ); // DoG filter in x-direction
    CvMat* DoGy = cvCreateMat( 7, 1, CV_64FC1 ); // DoG filter in y-direction
    cvTranspose( &DoGx, DoGy ); // transpose(&DoGx) -> DoGy
   
    // output information of lines found by Hough transform with NMS
    CvMemStorage* storageX = cvCreateMemStorage(0), *storageY = cvCreateMemStorage(0);
    CvSeq* linesX = 0, *linesY = 0;
   
    // arranged detected lines according the image frame coordinate
    vector<line_param> linesXorder; // ordered horizontal lines
    vector<line_param> linesYorder; // ordered vertical lines
   
    vector<CvPoint> p; // ordered intersection points on the "linesXorder" & "linesYorder"
   
   
    // 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);
   
   
    // initialize capture from a camera
    CvCapture* capture = cvCaptureFromCAM(0); // capture from video device #0
    int frame = 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           
            cvSaveImage("original.bmp", iplInput);
            if(iplInput) {
                if(0 == frame) {   
                    // 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);   
                    width = iplInput->width, height = iplInput->height;
                    cvMoveWindow( "temp", 100+width+10, 100 );
                    cvMoveWindow( title_fx, 100, 100+height+30 );
                    cvMoveWindow( title_fy, 100+width+10, 100+height+30 );
                    cvMoveWindow( title_ey, 100, 100+(height+30)*2 );
                    cvMoveWindow( title_ex, 100+width+10, 100+(height+30)*2 );
                }
                // convert the input color image to gray one   
                cvCvtColor(iplInput, iplGray, CV_BGR2GRAY); // convert an image from one color space to another
                cvSaveImage("gray.bmp", iplGray);       
                // convert one array to another with optional linear transformation
                cvConvert(iplGray, iplTemp);
                // increase the frame number
                frame++;
            }       
           
            // #1. DoG filtering
           
            // convolve an image with the DoG 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         

            // convert negative values to positive to filter the image in reverse direction
            cvAbs(iplDoGx, iplDoGx);            cvAbs(iplDoGy, iplDoGy);
           
            // normalize the pixel values
            cvMinMaxLoc( iplDoGx, &minValx, &maxValx ); // find global minimum and maximum in image array
            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("temp.bmp", iplTemp);           
            cvSaveImage("DoGx.bmp", iplDoGx);            cvSaveImage("DoGy.bmp", iplDoGy);
           
           
            // #2. separate selected edges into vertical and horizontal
           
            cvCopy(iplDoGy, iplDoGyClone), cvCopy(iplDoGx, iplDoGxClone);
            // non-Maximum suppression (NMS)
            edgethres = 0.5;
//            edgethres = 200;
            nonMaximumSuppression2 ( iplDoGx, iplDoGyClone, kernel, edgethres );           
            nonMaximumSuppression2 ( iplDoGy, iplDoGxClone, kernel, edgethres );

           
            // 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( iplTemp, &minValt, &maxValt );
            cvScale( iplTemp, iplTemp, 1.0 / maxValt );
            cvMinMaxLoc( iplDoGx, &minValx, &maxValx );        cvMinMaxLoc( iplDoGy, &minValy, &maxValy );
            cvScale( iplDoGx, iplDoGx, 1.0 / maxValx );        cvScale( iplDoGy, iplDoGy, 1.0 / maxValy );
            // display windows 
            cvShowImage( title_ex, iplDoGx );            cvShowImage( title_ey, iplDoGy );  
            // save images to files
            cvSaveImage("edgeX.bmp", iplDoGx);            cvSaveImage("edgeY.bmp", iplDoGy);
           
           
            // #3. line detection

//            cvConvert(iplDoGx, iplEdgeY);    cvConvert(iplDoGy, iplEdgeX);           
            cvMinMaxLoc( iplDoGx, &minValx, &maxValx );            cvMinMaxLoc( iplDoGy, &minValy, &maxValy );       
            cvConvertScale(iplDoGx, iplEdgeY, 255.0 / maxValx);            cvConvertScale(iplDoGy, iplEdgeX, 255.0 / maxValy);
//            cvConvertScale(iplDoGx, iplEdgeY, 255.0);            cvConvertScale(iplDoGy, iplEdgeX, 255.0)

            // Hough transform with Non Maximal Suppression
            linesX = cvHoughLinesLYM(iplEdgeX, storageX, CV_HOUGH_STANDARD, 1.0*rho, CV_PI/180*theta, threshold, rN, tN);
            linesY = cvHoughLinesLYM(iplEdgeY, storageY, CV_HOUGH_STANDARD, 1.0*rho, CV_PI/180*theta, threshold, rN, tN);           
           
            cout << endl << "frame # " << frame << "  ---------------------------" << endl;
            cout << "# of found lines = " << linesY->total << " vertical, " << linesX->total << " horizontal " << endl;
           
            // draw found lines
            drawLines( linesX, iplInput );            drawLines( linesY, iplInput );           
       
           
            // #4. intersection points calculation
           
            // arrange vertical lines from left to right           
            cout << "vertical" << endl;           
            indexLinesY(linesYorder, linesY, iplInput );
            // arrange horizontal lines    from up to bottom           
            cout << "horizontal" << endl;       
            indexLinesX(linesXorder, linesX, iplInput );           
           
            // calculate and index intersection points
            indexIntersections(p, linesXorder, linesYorder, iplInput);       
           
            cvShowImage( "input", iplInput );  
            cvSaveImage( "input.bmp", iplInput );
           
            if( cvWaitKey(10) >= 0 )
                break;
        } // end of else
    } // end of while
   
    cvReleaseCapture( &capture ); // release the capture source
    cvDestroyWindow( "input" );
    cvDestroyWindow( "temp" );
    cvDestroyWindow(title_fx);
    cvDestroyWindow(title_fy);
    cvDestroyWindow(title_ex);
    cvDestroyWindow(title_ey);
   
    return 0;
}

/* Test: feature point extraction in implementing virtual studio
 : using Gaussin gradient filter, first-order Derivative of Gaussian (DoG) filter
 with a kernel h = [-1, -7, -15, 0, 15, 7, 1]
 ref. swPark_2000rti.pdf:440p
 2010, lym
 camera: Logitech QuickCam Pro 4000
 */
// + non maximum suppression


#include <OpenCV/OpenCV.h>
#include <iostream>s
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++ )
        {
           
            if( x == image->width - 1 ) {
                x = 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;
            }
        }    
    }
}


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; // filtered image by DoG in x-direction
    IplImage* iplDoGy = 0; // filtered image by DoG in y-direction
    
    double minValx, maxValx;
    double minValy, maxValy;
    double minValt, maxValt;    
    
    double threshold = 100.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 titlef[200], titlex[200], titley[200];
    sprintf(titlef, "filtered image by DoG");   
    sprintf(titlex, "filtered image by DoGx");
    sprintf(titley, "filtered image by DoGy");    
    
    cvNamedWindow("input");
    cvNamedWindow( "temp" );
    cvNamedWindow(titlex);
    cvNamedWindow(titley);
//    cvNamedWindow(titlef);
    

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

            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);     
                }
                // 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);
           
            // non-Maximum suppression (NMS)
            // void NonMaximumSuppression ( IplImage* image, int kernel, int threshold )
            NonMaximumSuppression ( iplDoGx, kernel, threshold );
            NonMaximumSuppression ( iplDoGy, kernel, threshold );
           
            // 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( titlex, iplDoGx );
            cvShowImage( titley, iplDoGy );            
//            cvShowImage( titlef, iplDoG );
           
            // save images to files
            cvSaveImage("input.bmp", iplInput);
            cvSaveImage("gray.bmp", iplGray);
            cvSaveImage("temp.bmp", iplTemp);
            cvSaveImage("DoGx.bmp", iplDoGx);
            cvSaveImage("DoGy.bmp", iplDoGy);    
//            cvSaveImage("edgesDoG.bmp", iplDoG);    
           
            if( cvWaitKey(10) >= 0 )
                break;
        }
    }  
    
    cvReleaseCapture( &capture ); // release the capture source
    cvDestroyWindow( "input" );
    cvDestroyWindow( "temp" );
    cvDestroyWindow(titlex);
    cvDestroyWindow(titley);
//    cvDestroyWindow(titlef);
    
    return 0;  
}

/* Test: feature point extraction in implementing virtual studio
 : using Gaussin gradient filter, first-order Derivative of Gaussian (DoG) filter
 with a kernel h = [-1, -7, -15, 0, 15, 7, 1]
 ref. swPark_2000rti.pdf:440p
 2010, lym
 camera: Logitech QuickCam Pro 4000
 */

#include <OpenCV/OpenCV.h>
#include <iostream>s
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++ )
        {
           
            if( x == image->width - 1 ) {
                x = 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;
            }
        }    
    }
}


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; // filtered image by DoG in x-direction
    IplImage* iplDoGy = 0; // filtered image by DoG in y-direction
    
    double minValx, maxValx;
    double minValy, maxValy;
    double minValt, maxValt;    
    
    double threshold = 100.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 titlef[200], titlex[200], titley[200];
    sprintf(titlef, "filtered image by DoG");   
    sprintf(titlex, "filtered image by DoGx");
    sprintf(titley, "filtered image by DoGy");    
    
    cvNamedWindow("input");
    cvNamedWindow( "temp" );
    cvNamedWindow(titlex);
    cvNamedWindow(titley);
//    cvNamedWindow(titlef);
    

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

            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);     
                }
                // 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);
           
            // non-Maximum suppression (NMS)
            // void NonMaximumSuppression ( IplImage* image, int kernel, int threshold )
            NonMaximumSuppression ( iplDoGx, kernel, threshold );
            NonMaximumSuppression ( iplDoGy, kernel, threshold );
           
            // 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( titlex, iplDoGx );
            cvShowImage( titley, iplDoGy );            
//            cvShowImage( titlef, iplDoG );
           
            // save images to files
            cvSaveImage("input.bmp", iplInput);
            cvSaveImage("gray.bmp", iplGray);
            cvSaveImage("temp.bmp", iplTemp);
            cvSaveImage("DoGx.bmp", iplDoGx);
            cvSaveImage("DoGy.bmp", iplDoGy);    
//            cvSaveImage("edgesDoG.bmp", iplDoG);    
           
            if( cvWaitKey(10) >= 0 )
                break;
        }
    }  
    
    cvReleaseCapture( &capture ); // release the capture source
    cvDestroyWindow( "input" );
    cvDestroyWindow( "temp" );
    cvDestroyWindow(titlex);
    cvDestroyWindow(titley);
//    cvDestroyWindow(titlef);
    
    return 0;  
}

/* Test: feature point extraction in implementing virtual studio
 : using Gaussin gradient filter, first-order Derivative of Gaussian (DoG) filter
 with a kernel h = [-1, -7, -15, 0, 15, 7, 1]
 ref. swPark_2000rti.pdf:440p
 camera: Logitech QuickCam Pro 4000
 2010, lym & kyu
 */

// + non maximum suppression

#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;
                }
            }
        }
    }
}


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, "edges in x direction");
    sprintf(title_ey, "edges in y direction"); 
   
    cvNamedWindow(title_ex);
    cvNamedWindow(title_ey);
    //    cvNamedWindow(title_e);
   
   
    // 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);
                   
                }
                // 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( "temp", iplTemp ); 
            cvShowImage( title_ex, iplDoGx );
            cvShowImage( title_ey, iplDoGy );         
           
            // save images to files
            cvSaveImage("input.bmp", iplInput);
            cvSaveImage("gray.bmp", iplGray);
            cvSaveImage("temp.bmp", iplTemp);
            cvSaveImage("edgeX.bmp", iplDoGx);
            cvSaveImage("edgeY.bmp", iplDoGy); 
            
            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;
}

// cross ratio test 
// to design a pattern to be used in real-time virtual studio
// lym 2010-02-25

#include <iostream>
#include <cstring>
using namespace std;


int main (int argc, char * const argv[]) {
   
    string pi ("314159265358979323846264338327950288");
//    cout << "The size of str is " << pi.size() << endl;
    int number = pi.size();
   
    int p[number];
    for( int n = 0; n < number; n++ )
    {
        p[n] = (pi[n]) - 48 + 1;
        cout << p[n] << " ";
    }
    cout << endl;   

    int a = 0;
    int b = 0;
    int c = 0;
    int d = 0;
    double cross_ratio = 0.0;
   
    for( int n = 0; n < number-3; n++ )
    {
//        cout << "n = " << n << endl;
       
        a = p[n];
        b = a + p[n+1];
        c = b + p[n+2];
        d = c + p[n+3];

//        cout << "a=" << a << " b=" << b << " c=" << c << " d=" << d << endl;

        cross_ratio = double( (c-a)*(d-b) ) / double( (d-a)*(c-b) );
       
//        cout << "cross ratio = " << cross_ratio << endl;
        printf("cross ratio = %f\n", cross_ratio);
    }
       
    return 0;
}

/* Test: pattern design in implementing a virtual studio
 to design grid lines using cross-ratio as represented in "swPark_2000rti": (16)
 C(x1,x2,x3,x4) = (x2-x1)(x4-x3) / (x4-x2)(x3-x1)
 based on two plots in figure 7, 439p
 2010, lym
 */
 
#include <iostream>
using namespace std;

int main (int argc, char * const argv[]) {

    // ref. swPark_2000rti: 43p: figure 7
    // 37 consecutive cross-ratios for 40 vertical lines
    double crossratio_vertical[] = { 0.47, 0.11, 0.32, 0.17, 0.44, 0.08, 0.42, 0.25, 0.24, 0.13,
        0.46, 0.18, 0.19, 0.29, 0.21, 0.37, 0.16, 0.38, 0.23, 0.09, 0.37, 0.26, 0.31, 0.18,
    0.30, 0.15, 0.39, 0.16, 0.32, 0.27, 0.20, 0.28, 0.39, 0.12, 0.23, 0.28, 0.35 };   
    // 17 consecutive cross-ratios for 20 horizontal lines
    double crossratio_horizontal[] = { 0.42, 0.13, 0.32, 0.16, 0.49, 0.08, 0.40, 0.20, 0.29, 0.19,
    0.37, 0.13, 0.26, 0.38, 0.21, 0.16, 0.42 };
   
    cout << "# of cross-ratios in vertical lines = " << sizeof(crossratio_vertical) / sizeof(double) << endl;
    cout << "# of cross-ratios in horizontal lines = " << sizeof(crossratio_horizontal) / sizeof(double) << endl;

    // initialize grid lines of the pattern
    double x[40] = { 0.0 }; // 40 vertical lines
    double y[20] = { 0.0 }; // 20 horizontal lines
    double r = 0.0; // cross ratio of 4 consecutive lines
    double a,b,c,d; // temporary variables to be used in calculating a cross ratio 

    // set the positions of the first three vertical lines
    x[0] = 1.0;
    x[1] = x[0] + 1.0;
    x[2] = x[1] + 2.0;
   
    // set the positions of the first three horizontal lines
    y[0] = 1.0;
    y[1] = y[0] + 2.0;
    y[2] = y[1] + 1.0;
   
    // for vertical lines
    int number = 40;
    cout << endl << "x[0]=" << x[0] << "  x[1]=" << x[1] << "  x[2]=" << x[2] << endl;
    for( int n = 0; n < number-3; n++ )
    {
//        cout << "n = " << n << endl;
       
        r = crossratio_vertical[n];
 /*       a = x[n];
        b = x[n+1];
        c = x[n+2];
       
//        cout << "a=" << a << " b=" << b << " c=" << c << " ratio =" << r;
       
        d = ( (b-a)*c - (c-a)*b*r ) / ( (b-a) - (c-a)*r ) ;
       
//        cout << "  ==> d= " << d << endl;
   
        x[n+3] = d;
*/   
        double son = (x[n+1]-x[n])*x[n+2] - r*(x[n+2]-x[n])*x[n+1];
        double mom = x[n+1] - x[n] - r*(x[n+2]-x[n]);
        x[n+3] = son / mom;
       
        cout << "x[" << n+3 << "]=" << x[n+3] << "  ";
    }
    return 0;
}