Wednesday, April 25, 2012

SURF correspondence and Matching

Today I'm going through using OpenCV (2.1) or Blepo to create a nice and easy to use class for feature correspondence and correspondence matching. SURF stances for Speeded up Robust Features and is an improvement to SIFT (Scale Invariant Feature Transform). It isn't too bad to implement (and a good exercise) if you have the time but OpenCV has a nice version that does multiple octaves and octave layers.
One of the drawbacks of the OpenCV implementation is that they do not have a good correspondence matching algorithm, which is a shame, since a decent one can be done extremely easily. The code I have was written courtesy of Brian Peasley, which I have edited heavily (much to his probable dismay).

Click here for full code (blepo people look in the blepo folder)
Make sure your OpenCV lib and header files are set up properly

I build two classes for this. The first is a class for the feature, descriptor, and match information and is shown below:

class SURFDescriptor{
public:
    bool            valid;
    float             x, y, radius;
    //128 or 64 floats depending on whether extended is chosen or not
    float             *descriptor;   
    SURFDescriptor        *match;
};


The second class is the main SURF class for the feature extraction and putative matching. Note that this is not a necessary class and can be done with a series of functions instead. However, since I use a very large library, it is often more convenient for me to group these within classes. Feel free to take code out of context (though be careful about memory leaks).

class SURF{
public:
    typedef SURFDescriptor* Iterator;
    typedef SURFDescriptor* ConstIterator;

    SURF(double hessianThreshold = 500, bool extended = true) //good values are (500, true)
        : m_threshold(hessianThreshold), m_extended(extended), m_CvSeqKeypoints(NULL), m_CvSeqDescriptors(NULL), m_descriptors(NULL), m_size(0) {}
    ~SURF() {}

    //must release before done using, not included in the destructor for shallow copying purposes
    void Release() {
        if(m_storage != NULL) {
            cvClearMemStorage(m_storage);
            cvReleaseMemStorage(&m_storage);
            m_storage = NULL;
        }
        if(m_descriptors != NULL) {
            free(m_descriptors);
            m_descriptors = NULL;
        };
        m_size = 0;
    }
    inline Iterator Begin() { return m_descriptors; }
    inline Iterator End() { return m_descriptors + m_size; }
    inline ConstIterator Begin() const { return m_descriptors; }
    inline ConstIterator End() const { return m_descriptors + m_size; }
    int ExtractFeatures(IplImage *img);
    void OverlayFeatures(IplImage *img);
    IplImage* DisplayMatches(IplImage *img, IplImage *img2);
    int size() const {return m_size;}
    SURFDescriptor &operator()(int index){ return m_descriptors[index]; }
    SURFDescriptor operator()(int index) const { return m_descriptors[index]; }
    SURFDescriptor *operator[](int index) const { return &(m_descriptors[index]); }
    bool extended() const { return m_extended; }
    float PutativeMatch(const SURF &surf2, int prune_dist, float percent = 0.1f);

private:
    double         m_threshold;
    bool        m_extended;
    SURFDescriptor    *m_descriptors;
    CvSeq         *m_CvSeqKeypoints, *m_CvSeqDescriptors;
    CvMemStorage* m_storage;
    int        m_size;
};


There are a couple of nice feature display and matching display functions in the class as well for debugging. If you want those, just download the original file at the top or bottom.
For now, I will just show the extract features and matching algorithms:

#define SQR(X) ((X) * (X))

inline float SAD_Descriptors(const SURFDescriptor &desc1, const SURFDescriptor &desc2, int length){
    float SAD = 0;
    float *p1 = desc1.descriptor, *p2 = desc2.descriptor;
    for(int i = 0; i < length; i++)
        SAD += fabs( *p1++ - *p2++ );
   
    return SAD;
}

int SURF::ExtractFeatures(IplImage *img){
    m_storage = cvCreateMemStorage(0);
    CvMat *gray_img = cvCreateMat(img->height, img->width, CV_8UC1); //create matrix which is what CV SURF operates on
    cvCvtColor(img, gray_img, CV_BGR2GRAY); //convert color IplImage to a gray scale matrix
   
    //Find surf features and calculate their descriptors
    CvSURFParams params = cvSURFParams(m_threshold, m_extended ? 1 : 0);       
    cvExtractSURF(gray_img, 0, &m_CvSeqKeypoints, &m_CvSeqDescriptors, m_storage, params );

    m_descriptors = (SURFDescriptor *)malloc(m_CvSeqDescriptors->total*sizeof(SURFDescriptor));
           
    CvSURFPoint *r;
    SURFDescriptor *desc = m_descriptors;
    for(int i = 0; i < m_CvSeqKeypoints->total; i++){
        r = (CvSURFPoint*)cvGetSeqElem( m_CvSeqKeypoints, i );
       
        desc->x = r->pt.x;
        desc->y = r->pt.y;
        desc->radius = (float)r->size;

        desc->descriptor =    (float *)cvGetSeqElem( m_CvSeqDescriptors, i);
        desc++;       
    }
    m_size = m_CvSeqDescriptors->total;
    cvReleaseMat(&gray_img);
    return m_CvSeqDescriptors->total;
}


float SURF::PutativeMatch(const SURF &surf2, int prune_dist, float percent) {
    vector<int>     index1(m_size);
    vector<float>    minval1(m_size);
   
    vector<int>     index2(surf2.size());
    vector<float>    minval2(surf2.size());

    float val;
    //for every feature in surf1, compare descriptor to every feature in surf2 and save the index
    //of the minimum SAD.
    vector<float>::iterator min1 = minval1.begin(), min2 = minval2.begin();
    vector<int>::iterator in1 = index1.begin(), in2 = index2.begin();
    SURFDescriptor *desc1 = m_descriptors;
    for(unsigned int i = 0; i < minval2.size(); i++)
        *min2++ = 999.0;

    int tempcnt,ext=0;
    while(1)
    {
        tempcnt=0;
        min1 = minval1.begin();
        in1 = index1.begin();
        for(int i = 0; i < m_size; i++){
            *min1 = 999.0;
            SURF::Iterator desc2 = surf2.Begin();
            min2 = minval2.begin();
            in2 = index2.begin();
            for(int j = 0; j < surf2.size(); j++) {
                if(sqrt(SQR(desc1->x-desc2->x) + SQR(desc1->y-desc2->y)) <= (prune_dist+ext))
                    tempcnt++;
                val = SAD_Descriptors( *desc1, *desc2, m_extended ? 128 : 64 );   
           
                //comparing ith point in surf1 to every point in surf2
                if(val < *min1){
                    *min1 = val;
                    *in1 = j;
                }

                //comparing jth point in surf2 to ith point in surf1
                if(val < *min2){
                    *min2 = val;
                    *in2 = i;
                }
                desc2++;
                in2++;
                min2++;
            }
            desc1++;
            in1++;
            min1++;
        }
        if(tempcnt==0)
            ext++;
        else break;
    }

    //sort minvals to get top percent
    int minIndex;
    float temp;
    float avg = 0;
    min1 = minval1.begin();
    //find the average of minval1
    while(min1 != minval1.end())
        avg += *min1++;
    avg /= minval1.size();
    assert(percent > 0 && percent < 1);
    int endSort = (int)(percent*minval1.size());
    min1 = minval1.begin();
    for(unsigned int i = 0; i < endSort; i++) {
        minIndex = i;
        float currMin = *min1;
        vector<float>::const_iterator tmpMin = (minval1.begin()+i+1);
        for(unsigned int j = i+1; j < minval1.size(); j++) {
            if(*tmpMin < currMin) {
                minIndex = j;
                currMin =*tmpMin;
            }
            tmpMin++;
        }
        if(minIndex != i) {
            temp = *min1;
            *min1 = currMin;
            minval1[minIndex] = temp;
        }
        min1++;
    }

    float T = minval1[endSort];
    //see if minimum comparison is putative between feature point sets
    //dirty code, a little bit better now
    desc1 = m_descriptors;
    in1 = index1.begin();
    for(unsigned int i = 0; i < index1.size(); i++)
    {
        if(index2[*in1] == i && minval2[*in1] < T)       
            desc1->match = surf2[*in1];
        else
            desc1->match = NULL;
        in1++;
        desc1++;
    }
    return avg;
}


 Consider making the matching algorithm better or more efficient, possibly with an implementation involving RANSAC.




If you want the full source, including an example of how it works, download the source below:

Click here for full code (blepo people look in the blepo folder)
Make sure your OpenCV lib and header files are set up properly

Consider donating to further my tinkering.


Places you can find me

No comments:

Post a Comment