VLFeat SIFT with OpenCV + Code

Sunday, February 26, 2012

So, you used the VLFeat  SIFT successfuly in Matlab but you need to use the library with C++ and you can't find the functions reference nor a tutorial? Well I have been there, done that and sharing the code for integrating VLFeat's SIFT with OpenCV.

 Actually I harvest the below code from the 'Toolbox' folder of VLFeat, from the source of the dot mex files. With minor modification the ' VLSIFT' function takes an OpenCV IplImage (grayscale) and returns pointers to the frames and descriptors found as well as the number of the frames. I removed the selection for exporting the descriptors in double; they always be exported as the vl_uint8 datatype.




void VLSIFT(IplImage* image, vl_uint8* DATAdescr, double* DATAframes, int* nframes){
    //Take IplImage -> convert to SINGLE (float):
    float* frame = (float*)malloc(image->height*image->width*sizeof(float));
    uchar* Ldata      = (uchar *)image->imageData;
    for(int i=0;i<image->height;i  )
        for(int j=0;j<image->width;j  )
            frame[j*image->height i*image->nChannels] = (float)Ldata[i*image->widthStep j*image->nChannels];

    // VL SIFT computation:
    vl_sift_pix const *data ;
    int                M, N ;
    data = (vl_sift_pix*)frame;
    M = image->height;
    N = image->width;
   
    // VL SIFT PARAMETERS
    int                verbose = 1 ; // change to 2 for more verbose..
    int                O     =   -1 ; //Octaves
    int                S     =   3 ; //Levels
    int                o_min =   0 ;

    double             edge_thresh = -1 ;  //-1 will use the default (as in matlab)
    double             peak_thresh =  -1 ;
    double             norm_thresh = -1 ;
    double             magnif      = -1 ;
    double             window_size = -1 ;

    double            *ikeys = 0 ; //?
    int                nikeys = -1 ; //?
    vl_bool            force_orientations = 0 ;
    vl_bool            floatDescriptors = 0 ;

    /* -----------------------------------------------------------------
    *                                                            Do job
    * -------------------------------------------------------------- */
    {
        VlSiftFilt        *filt ;
        vl_bool            first ;
        double            *frames = 0 ;
        vl_uint8              *descr  = 0 ;
        int                reserved = 0, i,j,q ;

        /* create a filter to process the image */
        filt = vl_sift_new (M, N, O, S, o_min) ;

        if (peak_thresh >= 0) vl_sift_set_peak_thresh (filt, peak_thresh) ;
        if (edge_thresh >= 0) vl_sift_set_edge_thresh (filt, edge_thresh) ;
        if (norm_thresh >= 0) vl_sift_set_norm_thresh (filt, norm_thresh) ;
        if (magnif      >= 0) vl_sift_set_magnif      (filt, magnif) ;
        if (window_size >= 0) vl_sift_set_window_size (filt, window_size) ;

        if (verbose) {
            printf("vl_sift: filter settings:\n") ;
            printf("vl_sift:   octaves      (O)      = %d\n",
                vl_sift_get_noctaves      (filt)) ;
            printf("vl_sift:   levels       (S)      = %d\n",
                vl_sift_get_nlevels       (filt)) ;
            printf("vl_sift:   first octave (o_min)  = %d\n",
                vl_sift_get_octave_first  (filt)) ;
            printf("vl_sift:   edge thresh           = %g\n",
                vl_sift_get_edge_thresh   (filt)) ;
            printf("vl_sift:   peak thresh           = %g\n",
                vl_sift_get_peak_thresh   (filt)) ;
            printf("vl_sift:   norm thresh           = %g\n",
                vl_sift_get_norm_thresh   (filt)) ;
            printf("vl_sift:   window size           = %g\n",
                vl_sift_get_window_size   (filt)) ;
            printf("vl_sift:   float descriptor      = %d\n",
                floatDescriptors) ;

            printf((nikeys >= 0) ?
                "vl_sift: will source frames? yes (%d read)\n" :
            "vl_sift: will source frames? no\n", nikeys) ;
            printf("vl_sift: will force orientations? %s\n",
                force_orientations ? "yes" : "no") ;
        }

        /* ...............................................................
        *                                             Process each octave
        * ............................................................ */
        i     = 0 ;
        first = 1 ;
        while (1) {
            int                   err ;
            VlSiftKeypoint const *keys  = 0 ;
            int                   nkeys = 0 ;

            if (verbose) {
                printf ("vl_sift: processing octave %d\n",
                    vl_sift_get_octave_index (filt)) ;
            }

            /* Calculate the GSS for the next octave .................... */
            if (first) {
                err   = vl_sift_process_first_octave (filt, data) ;
                first = 0 ;
            } else {
                err   = vl_sift_process_next_octave  (filt) ;
            }

            if (err) break ;

            if (verbose > 1) {
                printf("vl_sift: GSS octave %d computed\n",
                    vl_sift_get_octave_index (filt));
            }

            /* Run detector ............................................. */
            if (nikeys < 0) {
                vl_sift_detect (filt) ;

                keys  = vl_sift_get_keypoints  (filt) ;
                nkeys = vl_sift_get_nkeypoints (filt) ;
                i     = 0 ;

                if (verbose > 1) {
                    printf ("vl_sift: detected %d (unoriented) keypoints\n", nkeys) ;
                }
            } else {
                nkeys = nikeys ;
            }

            /* For each keypoint ........................................ */
            for (; i < nkeys ;   i) {
                double                angles [4] ;
                int                   nangles ;
                VlSiftKeypoint        ik ;
                VlSiftKeypoint const *k ;

                /* Obtain keypoint orientations ........................... */
                if (nikeys >= 0) {
                    vl_sift_keypoint_init (filt, &ik,
                        ikeys [4 * i   1] - 1,
                        ikeys [4 * i   0] - 1,
                        ikeys [4 * i   2]) ;

                    if (ik.o != vl_sift_get_octave_index (filt)) {
                        break ;
                    }

                    k = &ik ;

                    /* optionally compute orientations too */
                    if (force_orientations) {
                        nangles = vl_sift_calc_keypoint_orientations
                            (filt, angles, k) ;
                    } else {
                        angles [0] = VL_PI / 2 - ikeys [4 * i   3] ;
                        nangles    = 1 ;
                    }
                } else {
                    k = keys   i ;
                    nangles = vl_sift_calc_keypoint_orientations
                        (filt, angles, k) ;
                }

                /* For each orientation ................................... */
                for (q = 0 ; q < nangles ;   q) {
                    vl_sift_pix  buf [128] ;
                    vl_sift_pix rbuf [128] ;

                    /* compute descriptor (if necessary) */
                    vl_sift_calc_keypoint_descriptor (filt, buf, k, angles [q]) ;
                    transpose_descriptor (rbuf, buf) ;

                    /* make enough room for all these keypoints and more */
                    if (reserved < (*nframes)   1) {
                        reserved  = 2 * nkeys ;
                        frames = (double*)realloc (frames, 4 * sizeof(double) * reserved) ;
                        descr  = (vl_uint8*)realloc (descr,  128 * sizeof(vl_uint8) * reserved) ;
                    }

                    /* Save back with MATLAB conventions. Notice tha the input
                    * image was the transpose of the actual image. */
                    frames [4 * (*nframes)   0] = k -> y ;
                    frames [4 * (*nframes)   1] = k -> x ;
                    frames [4 * (*nframes)   2] = k -> sigma ;
                    frames [4 * (*nframes)   3] = VL_PI / 2 - angles [q] ;

                    for (j = 0 ; j < 128 ;   j) {
                        float x = 512.0F * rbuf [j] ;
                        x = (x < 255.0F) ? x : 255.0F ;
                        descr[128 * (*nframes)   j] = (vl_uint8)x ;
                    }

                       (*nframes) ;
                } /* next orientation */
            } /* next keypoint */
        } /* next octave */

        if (verbose) {
            printf ("vl_sift: found %d keypoints\n", (*nframes)) ;
        }
        // save variables:
        memcpy(DATAframes, frames, 4 * (*nframes ) * sizeof(double));
        memcpy(DATAdescr, descr, 128 * (*nframes ) * sizeof(vl_uint8));

        /* cleanup */
        vl_sift_delete (filt) ;
    } /* end: do job */


    return;
}



Concecuently the next function will compute matches between to pictures using the matlab's 'vl_ubcmatch' function transformed for C++. Again I harvest this from the source of the .mex file. I deleted the macros used to choose between datatypes, because I use the vl_uint8 type for my descriptors (I think is the vl's default). So if you want to use a different datatype, don'd forget to change this.




void VLMATCH(vl_uint8* L1_pt,vl_uint8* L2_pt, int K1, int K2, double thresh, int* nMatches, double* MATCHES ){
    //Match descriptors!

    int ND = 128;

    Pair* pairs_begin = (Pair*) malloc(sizeof(Pair) * (K1 K2)) ;
    Pair* pairs_iterator = pairs_begin ;

    int k1, k2 ;                                                      
    const int maxval = 0x7fffffff ;                        
    for(k1 = 0 ; k1 < K1 ;   k1, L1_pt  = ND ) {    //kalooo!                  

        int best = maxval ;                                    
        int second_best = maxval ;                            
        int bestk = -1 ;                                                

        /* For each point P2[k2] in the second image... */              
        for(k2 =  0 ; k2 < K2 ;   k2, L2_pt  = ND) {                    

            int bin ;                                                      
            int acc = 0 ;                                        
            for(bin = 0 ; bin < ND ;   bin) {                              
                int delta =                                        
                    ((int) L1_pt[bin]) -                            
                    ((int) L2_pt[bin]) ;                            
                acc  = delta*delta ;                                        
            }                                                              

            /* Filter the best and second best matching point. */          
            if(acc < best) {                                              
                second_best = best ;                                        
                best = acc ;                                                
                bestk = k2 ;                                                
            } else if(acc < second_best) {                                
                second_best = acc ;                                          
            }                                                              
        }                                                                

        L2_pt -= ND*K2 ;                                                

        /* Lowe's method: accept the match only if unique. */            
        if(thresh * (float) best < (float) second_best &&                
            bestk != -1) {                                                
                pairs_iterator->k1 = k1 ;                                      
                pairs_iterator->k2 = bestk ;                                  
                pairs_iterator->score = best ;                                
                pairs_iterator   ;
                (*nMatches)  ;
        }                                                                
    }                                                                  

    Pair* pairs_end = pairs_iterator ;
    //double* M_pt = (double*)calloc((pairs_end-pairs_begin)*2,sizeof(double));
    double* M_pt = (double*)calloc((*nMatches)*2,sizeof(double));
    //double* M_start = M_pt;

    for(pairs_iterator = pairs_begin ;
        pairs_iterator < pairs_end  ;
          pairs_iterator) {
            *M_pt   = pairs_iterator->k1 ;
            *M_pt   = pairs_iterator->k2 ;
    }
    M_pt -= (*nMatches)*2 ;
    memcpy(MATCHES,M_pt,(*nMatches) * 2 * sizeof(double));
    free(pairs_begin) ;
    free(M_pt);

    return;
}



Finally call the functions in the main. Optionally plot the image with the found frames overlapped using OpenCV.






typedef unsigned char uchar;
int main()
{
    // Load template image:
    IplImage* Timage = cvLoadImage("T.png",0);

    double*            TFrames = (double*)calloc ( 4 * 10000, sizeof(double) ) ;
    vl_uint8*          TDescr  = (vl_uint8*)calloc ( 128 * 10000, sizeof(vl_uint8) ) ;
    int                Tnframes = 0;
    VLSIFT(Timage, TDescr, TFrames, &Tnframes);
    TFrames = (double*)realloc (TFrames, 4 * sizeof(double) * Tnframes) ; // = Y X Scale Angle
    TDescr  = (vl_uint8*)realloc (TDescr,  128 * sizeof(vl_uint8) * Tnframes) ;
   
    /*
    for(int i=0;i<Tnframes;i  ){
        cvCircle(Timage,              
        cvPoint(TFrames[0 i*4], TFrames[1 i*4]), TFrames[2 i*4],  
        cvScalar(255, 0, 0, 0),
        1, 8, 0);
    }
    cvShowImage("FrameT", Timage);
    cvWaitKey(0);
    */
return(0);
}



Sorry for the extra spaces in the code. I'm still familiarizing with the syntaxhighlighter..


Steve

7 comments:

Meet Shah said...

Can you please first tell me how to set up VLfeat for use in OpenCV?

I got following error while trying to run this code error:
siftTest.cpp:1: error: variable or field ‘VLSIFT’ declared void
siftTest.cpp:1: error: ‘IplImage’ was not declared in this scope
siftTest.cpp:1: error: ‘image’ was not declared in this scope
siftTest.cpp:1: error: ‘vl_uint8’ was not declared in this scope
siftTest.cpp:1: error: ‘DATAdescr’ was not declared in this scope
siftTest.cpp:1: error: expected primary-expression before ‘double’
siftTest.cpp:1: error: expected primary-expression before ‘int’

steve said...

Hello and sorry for the delay. After setting up the OpenCV as described here (http://youtu.be/XeBhwbRoKvk) and creating your openCV project, you must add to project's additional include directories the 'vlfeat-version' folder (replace version with your vlfeat version ex C:\Libraries\vlfeat-0.9.14) containing the 'bin', 'doc', 'src', etc and under linker, add the library directory 'vlfeat-version/bin/arch' (ex C:\Libraries\vlfeat-0.9.14\bin\win64). Finally in Linker>Input, add vl.lib.
Use the code provided here:
http://dl.dropbox.com/u/218924/Blog/VLFeat.txt
it should compile just fine. It loads a template image and a query one and computes the matches given a threshold.
Don't forget to add in the same folder as your executable program all the necessary openCV .dlls and the appropriate vl.dll (found inside your vlfeat-/bin/ directory).

Unknown said...

Thank you very much for the great work. However, it seems that the VLMATCH function run much slower than vl_ubcmatch in Matlab (about five time). What might be the reason?

Unknown said...

Thanks, Steve! Your adaptation of SIFT form vlfeat toolbox to OpenCV works great. I was similarly trying to import Vlfeat covariace detector to OpenCV. I have added below its code snippet. But the number of features detected are zero. Please let me know, what should be corrected.


VlCovDet * covdet = vl_covdet_new(VL_COVDET_METHOD_HESSIAN_LAPLACE);

//set parameters
vl_covdet_set_first_octave(covdet,-1);
vl_covdet_set_octave_resolution(covdet,3);
vl_covdet_set_peak_threshold(covdet, 0.1) ;
vl_covdet_set_edge_threshold(covdet, 1) ;

// process the image and run the detector
vl_covdet_put_image(covdet, ImageData, Image->height, Image->width) ;
vl_covdet_detect(covdet) ;

// get feature frames back
vl_size numFeatures = vl_covdet_get_num_features(covdet) ;

VlCovDetFeature const * feature = (VlCovDetFeature*)vl_covdet_get_features(covdet) ;


for (int i=0;i<numFeatures;i++)
{

float x = feature[i].frame.x;
float y = feature[i].frame.y;
float sigma = 2;//
std::cout << x << "\t" << y << std::endl;
cvDrawCircle(Image, cvPoint(x,y), sigma/2, CV_RGB( 255,255,255));
}

JN said...

Hi steve and thanks for your code! But i don t understand the left curly bracket after the comment "do job", what it refers to?

JN said...

It doesn't work for me... it crashes inside VLSIFT function at
frame[j*image->height + i*image->nChannels] = (float)Ldata[i*image->widthStep + j*image->nChannels];

JN said...

Hi Steve, i have a dubt about your code... when you create a "frame" you put four information in this order [y,x,scale,angle] but when you use the cvPoint constructor which requires (x,y) you do cvPoint(TFrames[0 i*4], TFrames[1 i*4])

Post a Comment

 
Copyright © Synaptic Activity
Blogger Theme by BloggerThemes Sponsored by Busy Buzz Blogging