Thursday, July 5, 2012

SURF Souce Code (Part-1)

When you install openCV, openCV samples folder also gets installed. The samples will be present in your InstallationDirectory/samples/C. For me it's present in OpenCV-2.1.0/samples/c.This folder contains the sample codes of many good openCV programmes that can be used for a wide variety of purposes. One more thing is that, it also contains the compiled object files along with the source code for each programme. The programme , we will be looking is find_obj.cpp and the compiled code will be with the name find_obj. This programme uses SURF to do an object detection. The original code, written by Liu, is modified by me to give the below code. Note that the comments made by me, start with Dileep:
 * A Demo to OpenCV Implementation of SURF
 * Further Information Refer to "SURF: Speed-Up Robust Feature"
 * Author: Liu Liu
 * Modified by :
 * Modifying Author : Dileep Kumar Kotha

#include <cv.h>
#include <highgui.h>
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>

#include <iostream>
#include <vector>

using namespace std;

static double dis=0;//For calculating the average descriptor distance
// define whether to use approximate nearest-neighbor search
//#define USE_FLANN

IplImage *image = 0;

compareSURFDescriptors( const float* d1, const float* d2, double best, int length )
    double total_cost = 0;
    assert( length % 4 == 0 );
    for( int i = 0; i < length; i += 4 )
        double t0 = d1[i] - d2[i];
        double t1 = d1[i+1] - d2[i+1];
        double t2 = d1[i+2] - d2[i+2];
        double t3 = d1[i+3] - d2[i+3];
        total_cost += t0*t0 + t1*t1 + t2*t2 + t3*t3;
We are sending a total cost, that's slightly greater or smaller than the best 

      if( total_cost > best )
    return total_cost;

naiveNearestNeighbor( const float* vec, int laplacian,
                      const CvSeq* model_keypoints,
                      const CvSeq* model_descriptors )
    int length = (int)(model_descriptors->elem_size/sizeof(float));
    int i, neighbor = -1;
    double d, dist1 = 1e6, dist2 = 1e6;
    CvSeqReader reader, kreader;
    cvStartReadSeq( model_keypoints, &kreader, 0 );
    cvStartReadSeq( model_descriptors, &reader, 0 );

    for( i = 0; i < model_descriptors->total; i++ )
        const CvSURFPoint* kp = (const CvSURFPoint*)kreader.ptr;
        const float* mvec = (const float*)reader.ptr;
     CV_NEXT_SEQ_ELEM( kreader.seq->elem_size, kreader );
        CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
        if( laplacian != kp->laplacian )
        d = compareSURFDescriptors( vec, mvec, dist2, length );

        if( d < dist1 )
            dist2 = dist1;
            dist1 = d;
            neighbor = i;
        else if ( d < dist2 )
            dist2 = d;

    if ( dist1 < 0.6*dist2 )
        return neighbor;
    return -1;

findPairs( const CvSeq* objectKeypoints, const CvSeq* objectDescriptors,
           const CvSeq* imageKeypoints, const CvSeq* imageDescriptors, vector<int>& ptpairs )
    int i;
    CvSeqReader reader, kreader;
    cvStartReadSeq( objectKeypoints, &kreader );
    cvStartReadSeq( objectDescriptors, &reader );

    for( i = 0; i < objectDescriptors->total; i++ )
        const CvSURFPoint* kp = (const CvSURFPoint*)kreader.ptr;
        const float* descriptor = (const float*)reader.ptr;
        CV_NEXT_SEQ_ELEM( kreader.seq->elem_size, kreader );
        CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
        int nearest_neighbor = naiveNearestNeighbor( descriptor, kp->laplacian, imageKeypoints, imageDescriptors);
        if( nearest_neighbor >= 0 )

printf("\n%lf\n",(dis/objectDescriptors->total));////Dileep:Here's where I am outputting the distance between the images
/*Dileep: If you are using this for recognition, write this distance to a file along with the name of the image you are matching against. After doing this for several images, you can then sort them in ascending order to find the best possible match - the one with the smallest distance. Here, I am outputting the distance to stdout
////Dileep:I am commenting the below one since it's of flann method and what ever the method we choose to find distance, the result is the same
flannFindPairs( const CvSeq*, const CvSeq* objectDescriptors,
           const CvSeq*, const CvSeq* imageDescriptors, vector<int>& ptpairs )
 int length = (int)(objectDescriptors->elem_size/sizeof(float));

    cv::Mat m_object(objectDescriptors->total, length, CV_32F);
 cv::Mat m_image(imageDescriptors->total, length, CV_32F);

 // copy descriptors
    CvSeqReader obj_reader;
 float* obj_ptr = m_object.ptr<float>(0);
    cvStartReadSeq( objectDescriptors, &obj_reader );
    for(int i = 0; i < objectDescriptors->total; i++ )
        const float* descriptor = (const float*)obj_reader.ptr;
        CV_NEXT_SEQ_ELEM( obj_reader.seq->elem_size, obj_reader );
        memcpy(obj_ptr, descriptor, length*sizeof(float));
        obj_ptr += length;
    CvSeqReader img_reader;
 float* img_ptr = m_image.ptr<float>(0);
    cvStartReadSeq( imageDescriptors, &img_reader );
    for(int i = 0; i < imageDescriptors->total; i++ )
        const float* descriptor = (const float*)img_reader.ptr;
        CV_NEXT_SEQ_ELEM( img_reader.seq->elem_size, img_reader );
        memcpy(img_ptr, descriptor, length*sizeof(float));
        img_ptr += length;

    // find nearest neighbors using FLANN
    cv::Mat m_indices(objectDescriptors->total, 2, CV_32S);
    cv::Mat m_dists(objectDescriptors->total, 2, CV_32F);
    cv::flann::Index flann_index(m_image, cv::flann::KDTreeIndexParams(4));  // using 4 randomized kdtrees
    flann_index.knnSearch(m_object, m_indices, m_dists, 2, cv::flann::SearchParams(64) ); // maximum number of leafs checked

    int* indices_ptr = m_indices.ptr<int>(0);
    float* dists_ptr = m_dists.ptr<float>(0);
    for (int i=0;i<m_indices.rows;++i) {
     if (dists_ptr[2*i]<0.6*dists_ptr[2*i+1]) {

// a rough implementation for object location 
locatePlanarObject( const CvSeq* objectKeypoints, const CvSeq* objectDescriptors,
                    const CvSeq* imageKeypoints, const CvSeq* imageDescriptors,
                    const CvPoint src_corners[4], CvPoint dst_corners[4] )
    double h[9];
    CvMat _h = cvMat(3, 3, CV_64F, h);
    vector<int> ptpairs;
    vector<CvPoint2D32f> pt1, pt2;
    CvMat _pt1, _pt2;
    int i, n;

#ifdef USE_FLANN
    flannFindPairs( objectKeypoints, objectDescriptors, imageKeypoints, imageDescriptors, ptpairs );
    findPairs( objectKeypoints, objectDescriptors, imageKeypoints, imageDescriptors, ptpairs );

    n = ptpairs.size()/2;
    if( n < 4 )
        return 0;

    for( i = 0; i < n; i++ )
        pt1[i] = ((CvSURFPoint*)cvGetSeqElem(objectKeypoints,ptpairs[i*2]))->pt;
        pt2[i] = ((CvSURFPoint*)cvGetSeqElem(imageKeypoints,ptpairs[i*2+1]))->pt;

    _pt1 = cvMat(1, n, CV_32FC2, &pt1[0] );
    _pt2 = cvMat(1, n, CV_32FC2, &pt2[0] );
    if( !cvFindHomography( &_pt1, &_pt2, &_h, CV_RANSAC, 5 ))
        return 0;

    for( i = 0; i < 4; i++ )
        double x = src_corners[i].x, y = src_corners[i].y;
        double Z = 1./(h[6]*x + h[7]*y + h[8]);
        double X = (h[0]*x + h[1]*y + h[2])*Z;
        double Y = (h[3]*x + h[4]*y + h[5])*Z;
        dst_corners[i] = cvPoint(cvRound(X), cvRound(Y));

    return 1;

int main(int argc, char** argv)
    const char* object_filename = argc == 3 ? argv[1] : "box.png";
    const char* scene_filename = argc == 3 ? argv[2] : "box_in_scene.png";

    CvMemStorage* storage = cvCreateMemStorage(0);

    cvNamedWindow("Object", 1);
    cvNamedWindow("Object Correspond", 1);

    static CvScalar colors[] = 

    IplImage* object = cvLoadImage( object_filename, CV_LOAD_IMAGE_GRAYSCALE );
    IplImage* image = cvLoadImage( scene_filename, CV_LOAD_IMAGE_GRAYSCALE );
    if( !object || !image )
        fprintf( stderr, "Can not load %s and/or %s\n"
            "Usage: find_obj [<object_filename> <scene_filename>]\n",
            object_filename, scene_filename );
    IplImage* object_color = cvCreateImage(cvGetSize(object), 8, 3);
    cvCvtColor( object, object_color, CV_GRAY2BGR );
    CvSeq *objectKeypoints = 0, *objectDescriptors = 0;
    CvSeq *imageKeypoints = 0, *imageDescriptors = 0;
    int i;
    CvSURFParams params = cvSURFParams(500, 1);

    double tt = (double)cvGetTickCount();
    cvExtractSURF( object, 0, &objectKeypoints, &objectDescriptors, storage, params );
    printf("Object Descriptors: %d\n", objectDescriptors->total);
    cvExtractSURF( image, 0, &imageKeypoints, &imageDescriptors, storage, params );
    printf("Image Descriptors: %d\n", imageDescriptors->total);
    tt = (double)cvGetTickCount() - tt;
    printf( "Extraction time = %gms\n", tt/(cvGetTickFrequency()*1000.));
    CvPoint src_corners[4] = {{0,0}, {object->width,0}, {object->width, object->height}, {0, object->height}};
    CvPoint dst_corners[4];
    IplImage* correspond = cvCreateImage( cvSize(image->width, object->height+image->height), 8, 1 );
    cvSetImageROI( correspond, cvRect( 0, 0, object->width, object->height ) );
    cvCopy( object, correspond );
    cvSetImageROI( correspond, cvRect( 0, object->height, correspond->width, correspond->height ) );
    cvCopy( image, correspond );
    cvResetImageROI( correspond );

#ifdef USE_FLANN
    printf("Using approximate nearest neighbor search\n");
    if( locatePlanarObject( objectKeypoints, objectDescriptors, imageKeypoints,
        imageDescriptors, src_corners, dst_corners ))
        for( i = 0; i < 4; i++ )
            CvPoint r1 = dst_corners[i%4];
            CvPoint r2 = dst_corners[(i+1)%4];
            cvLine( correspond, cvPoint(r1.x, r1.y+object->height ),
                cvPoint(r2.x, r2.y+object->height ), colors[8] );

    vector<int> ptpairs;
#ifdef USE_FLANN
    flannFindPairs( objectKeypoints, objectDescriptors, imageKeypoints, imageDescriptors, ptpairs );
    findPairs( objectKeypoints, objectDescriptors, imageKeypoints, imageDescriptors, ptpairs );

////Dileep:It is used to draw the the lines and rectangles in the windows that display images. If you want just the distance, you can comment the whole 
////below code
    for( i = 0; i < (int)ptpairs.size(); i += 2 )
        CvSURFPoint* r1 = (CvSURFPoint*)cvGetSeqElem( objectKeypoints, ptpairs[i] );
        CvSURFPoint* r2 = (CvSURFPoint*)cvGetSeqElem( imageKeypoints, ptpairs[i+1] );
        cvLine( correspond, cvPointFrom32f(r1->pt),
            cvPoint(cvRound(r2->pt.x), cvRound(r2->pt.y+object->height)), colors[8] );

    cvShowImage( "Object Correspond", correspond );
    for( i = 0; i < objectKeypoints->total; i++ )
        CvSURFPoint* r = (CvSURFPoint*)cvGetSeqElem( objectKeypoints, i );
        CvPoint center;
        int radius;
        center.x = cvRound(r->pt.x);
        center.y = cvRound(r->pt.y);
        radius = cvRound(r->size*1.2/9.*2);
        cvCircle( object_color, center, radius, colors[0], 1, 8, 0 );
    cvShowImage( "Object", object_color );


    cvDestroyWindow("Object SURF");
    cvDestroyWindow("Object Correspond");

    return 0;

There are enough comments in it, explaining what each part does. compile it. If you compile it to the object file find_obj, the command below works.


Otherwise, replace it with a.out. Also, copy the two images in the samples folder named box.png and box_in_scene.png into the folder where you run the above command.

Now for ./find_obj Since we are not supplying arguements, it will take the default box.png and box_in_scene.png and tries to find the first object inside the second object. The below figure appears.

SURF example

What's actually happening here is object detection. But I included a "dis" variable and also printing it to standard ouput, so that you will also know the average distance of all the descriptors. If you want, you can keep two photographs of the same person with different facial expressions in the same folder and supply them as arg1 and arg2 (in which case you will get distance as zero, ofcourse :P)

./findobj arg1 arg2 

If we keep on changing the second arguement(the image we are supplying as arg2), we can actually use it for recognition. Then the one with the least distance is the best possible match for the first argument (the image we are supplying as arg1). In the next post, I will strip all the unnecessary components of this code and make this usable for recognition. You can try to do it yourself.


محمد حسام الدين دويدري said...

Thanks a lot for this useful post,
Actually I used compiled the same code using Opencv 2.1, I always get an exception error: at the function cvExtractSURF !! even though I load the same images mentioned in the opencv sample??
could you please help to solve this problem

Dileep Kumar said...

That's an inbuilt function of openCV and it should not cause any problem. May be the problem is with openCV installation in your system. Can you try in one of your friend's system and revert back. Or try in a linux m/c. Also there is one post in this blog, that explains about proper installation of openCV. Try looking into that, if the problem is not resolved

محمد حسام الدين دويدري said...

Thanks a lot
I've tested the same code on another pc which has the same version of opencv, and the same windows and VS, it works well, but unfortunatly it does not work on mine, even I developped alot of project on mine without any problem !
what is your suggestion ?

Dileep Kumar said...

Can you try reinstalling openCV with the help of your friend, in whom the project worked.

محمد حسام الدين دويدري said...

I did but in vain !!
I confirm that I already acheived a lot of projects using the current configurations without any problem. I wonder when I get an error only at cvExtractSURF !
hope to find a solution

siya said...

hi i had a doubt with opencv little coding!
could i email u my doubt!
my email id is

Dileep Kumar said...

Yeah, u can mail me. Check my profile to get the mail ID.

But it's better to comment, since it will be useful to the ones, who got a similar problem.

srushti said...

Hi dileep. There is a problem in cvCopy.. if the sizes are different, it results in exceptions. can you plz fix it

stark said...

hello dileep sir. i m engineering student and my final year project is on object recognition system so i need some help so plz help me..mail id

Dileep Kumar said...

Hi Srushti,
cvCopy is a library of openCV. I cannot fix them.

Hi stark,
I am totally busy with my office work these days. Please mention specific points where you don't understand. This blog is small enough and concise to read and start working with openCV

Nikhil Nayanar said...

Hi the vector bearing location of each feature is ptpairs ?

Post a Comment