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;

};

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;

};

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;

}

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;

}

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*