SURF 중요 코드

File List

o Main.cpp
o integral.cpp
o fasthessian.cpp
o surf.cpp
o surflib.h
o utils.cpp
o ipoint.h
o fasthessian.h
o utils.h
o surf.h

Main.cpp

#include "surflib.h"

//-------------------------------------------------------
// Define PROCEDURE as:
//  - 1 and supply image path to run on static image
//  - 2 to capture from a webcam
#define PROCEDURE 1

//-------------------------------------------------------

int mainImage(const char *filename);
int mainVideo(void);

int main() 
{
  if (PROCEDURE == 1) return mainImage("IMG_4207.png");
  if (PROCEDURE == 2) return mainVideo();
}

//-------------------------------------------------------

int mainImage(const char *filename)
{
  // Declare Ipoints and other stuff
  IpVec ipts;
  IplImage *img=cvLoadImage(filename);

  // Detect and describe interest points in the image
  surfDetDes(img, ipts, false, 3, 4, 1, 0.0004f);

  // Draw the detected points
  drawIpoints(img, ipts);

  // Display the result
  showImage(img);

  return 0;
}

//-------------------------------------------------------

int mainVideo(void)
{
  // Initialise capture device
  CvCapture* capture = cvCaptureFromCAM( CV_CAP_ANY );
  if(!capture) error("No Capture");

  // Create a window 
  cvNamedWindow("OpenSURF", CV_WINDOW_AUTOSIZE );

  // Declare Ipoints and other stuff
  IpVec ipts;
  IplImage *img=NULL;

  // Main capture loop
  while( 1 ) 
  {
    // Grab frame from the capture source
    img = cvQueryFrame(capture);

    // Detect and describe interest points in the image
    surfDetDes(img, ipts, false, 3, 4, 2, 0.002f);

    // Draw the detected points
    drawIpoints(img, ipts);

    // Draw the FPS figure
    drawFPS(img);

    // Display the result
    cvShowImage("OpenSURF", img);

    // If ESC key pressed exit loop
    if( (cvWaitKey(10) & 255) == 27 ) break;
  }

  // Release the capture device
  cvReleaseCapture( &capture );
  cvDestroyWindow( "OpenSURF" );
  return 0;
}

//-------------------------------------------------------

integral.cpp

#include "cv.h"
#include "highgui.h"
#include "integral.h"
#include "utils.h"


//! Computes the integral image of image img.  Assumes source image to be a 
//! 32-bit floating point.  Returns IplImage of 32-bit float form.
IplImage *Integral(IplImage *source)
{
  // convert the image to single channel 32f
  IplImage *img = getGray(source);
  IplImage *int_img = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, 1);

  // set up variables for data access
  int height = img->height;
  int width = img->width;
  int step = img->widthStep/sizeof(float);
  float *data   = (float *) img->imageData;  
  float *i_data = (float *) int_img->imageData;  

  // calculate integral image values
  for(int i=0; i<height; i++) 
  {
    float rs = 0.0f;
    for(int j=0; j<width; j++) 
    {
      // cumulative row sum
      rs = rs+data[i*step+j]; 
                  
      // integral image sum
      if (i*step+j >= step) 
      {
        i_data[i*step+j] = i_data[(i-1)*step+j]+rs; 
      }
      else 
      {
        // these values computed for first row only
        i_data[i*step+j] = rs; 
      }
    }
  }

  // release the gray image
  cvReleaseImage(&img);

  // return the integral image
  return int_img;
}

//-------------------------------------------------------

//! Computes the sum of pixels within the rectangle specified by the top-left start
//! co-ordinate (x,y) with width w and height h.
float Area(IplImage *img, int x, int y, int w, int h)
{
  float *data = (float *) img->imageData;
  int step = img->widthStep/sizeof(float);
  int x1=x-1, y1=y-1, x2=x1+w, y2=y1+h;

  // if all values are within bounds return BR + TL - (BL + TR)
  if (x1>=0 && y1 >=0 && x2<img->width && y2<img->height) 
  {
    return data[y2*step+x2] + data[y1*step+x1] 
          -(data[y2*step+x1] + data[y1*step+x2]);
  }
  else 
    return 0;
}

fasthessian.cpp

#include "cv.h"
#include "integral.h"
#include "fasthessian.h"
#include "ipoint.h"
#include "utils.h"
#include <vector>
#include <iostream>
using namespace std;



//! Destructor
FastHessian::~FastHessian() 
{
  // free the det array
  if (m_det) delete [] m_det;
}

//-------------------------------------------------------

//! Constructor without image
FastHessian::FastHessian(std::vector<Ipoint> &ipts, 
                         const int octaves, const int intervals, const int init_sample, 
                         const float thres) 
                         : ipts(ipts), m_det(NULL), i_width(0), i_height(0)
{
  // Initialise variables and construct Det array
  this->octaves = octaves;
  this->intervals = intervals;
  this->init_sample = init_sample;
  this->thres = thres;
}

//-------------------------------------------------------

//! Constructor with image
FastHessian::FastHessian(IplImage *img, std::vector<Ipoint> &ipts, 
                         const int octaves, const int intervals, const int init_sample, 
                         const float thres) 
                         : ipts(ipts), m_det(NULL), i_width(0), i_height(0)
{
  // Initialise variables and construct Det array
  this->octaves = octaves;
  this->intervals = intervals;
  this->init_sample = init_sample;
  this->thres = thres;

  // Set the current image
  setIntImage(img);
}

//-------------------------------------------------------

//! Set or re-set the integral image source
void FastHessian::setIntImage(IplImage *img)
{
  // Change the source image
  this->img = img;

  // Redefine width, height and det map only if image has changed size
  if (img->width != i_width || img->height != i_height) 
  {
    i_width = img->width;
    i_height = img->height;

    // Allocate space for determinant of hessian pyramid 
    if (m_det) delete [] m_det;
    m_det = new float [octaves*intervals*i_width*i_height];
    memset(m_det,0,sizeof(m_det));
  }
}

//-------------------------------------------------------

//! Find the image features and write into vector of features
void FastHessian::getIpoints()
{
  // Clear the vector of exisiting ipts
  ipts.clear();

  // Calculate approximated determinant of hessian values
  buildDet();

  for(int o=0; o < octaves; o++) 
  {
    // for each octave double the sampling step of the previous
    int step = init_sample * cvRound(pow(2.0f,o));

    // determine border width for the largest filter for each ocave
    int border = (3 * cvRound(pow(2.0,o+1)*(intervals)+1) + 1)/2;

    // check for maxima across the scale space
    for(int i = 1; i < intervals-1; ++i) 
      for(int r = (border); r < i_height - (border); r += step)
        for(int c = (border); c < i_width - (border); c += step) 
          if(isExtremum(o, i, c, r)) 
            interp_extremum(o,i,r,c);
  } 
}

//-------------------------------------------------------

//! Calculate determinant of hessian responses
void FastHessian::buildDet()
{
  int lobe, border, step;
  float Dxx, Dyy, Dxy, scale;

  for(int o=0; o<octaves; o++) 
  {
    // calculate filter border for this octave
    border = (3 * fRound(pow(2.0f,o+1)*(intervals)+1) + 1)/2;
    step = init_sample * fRound(pow(2.0f,o));

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

      // calculate lobe length (filter side length/3)
      lobe = fRound(pow(2.0f,o+1)*(i+1)+1);
      scale = 1.0f/pow((float)(3*lobe),2);

      for(int r=border; r<i_height-border; r+=step) 
      {
        for(int c=border; c<i_width-border; c+=step) 
        {
          Dyy = Area(img, c-(lobe-1), r-((3*lobe-1)/2), (2*lobe-1), lobe)
             -2*Area(img, c-(lobe-1), r-((lobe-1)/2), (2*lobe-1), lobe)
               +Area(img, c-(lobe-1), r+((lobe+1)/2), (2*lobe-1), lobe);

          Dxx = Area(img, c-((3*lobe-1)/2), r-(lobe-1), lobe, (2*lobe-1))
             -2*Area(img, c-((lobe-1)/2), r-(lobe-1), lobe, (2*lobe-1))
               +Area(img, c+((lobe+1)/2), r-(lobe-1), lobe, (2*lobe-1));

          Dxy =  Area(img, c-lobe, r-lobe, lobe, lobe)
                +Area(img, c+1, r+1, lobe, lobe)
                -Area(img, c-lobe, r+1, lobe, lobe)
                -Area(img, c+1, r-lobe, lobe, lobe);

          // Normalise the filter responses with respect to their size
          Dxx *= scale;
          Dyy *= scale;
          Dxy *= scale;

          // Get the sign of the laplacian
          int lap_sign = (Dxx+Dyy >= 0 ? 1 : -1);

          // Get the determinant of hessian response
          float res = (Dxx*Dyy - pow(0.9f*Dxy,2));
          res = (res < thres ? 0 : lap_sign * res);

          // calculate approximated determinant of hessian value
          m_det[(o*intervals+i)*(i_width*i_height) + (r*i_width+c)] = res;
        }
      }
    }
  }
}   

//-------------------------------------------------------

//! Non Maximal Suppression function
int FastHessian::isExtremum(int octave, int interval, int c, int r)
{
  float val = getVal(octave, interval, c, r);
  int step = init_sample * cvRound(pow(2.0f,octave));

  // reject points with low response to the det of hessian function
  if(val < thres) return 0;

  // check for maximum 
  for(int i = -1; i <= 1; i++ )
    for(int j = -step; j <= step; j+=step )
      for(int k = -step; k <= step; k+=step ) 
        if (i != 0 || j != 0 || k != 0)
          if(getVal(octave, interval+i, c+j, r+k) > val)
            return 0;

  return 1;
}

//-------------------------------------------------------

//! Return the value of the approximated determinant of hessian
inline float FastHessian::getVal(int o, int i, int c, int r)
{
  return fabs(m_det[(o*intervals+i)*(i_width*i_height) + (r*i_width+c)]);
}

//-------------------------------------------------------

//! Return the sign of the laplacian (trace of the hessian)
inline int FastHessian::getLaplacian(int o, int i, int c, int r)
{
  float res = (m_det[(o*intervals+i)*(i_width*i_height) + (r*i_width+c)]);

  return (res >= 0 ? 1 : -1);
}

//-------------------------------------------------------

//! Round an int to a float
inline int FastHessian::fRound(float flt)
{
  return (int) floor(flt+0.5f);
}

//-------------------------------------------------------

//! Return the value of the approximated determinant of hessian
inline float FastHessian::getValLowe(int o, int i, int r, int c)
{
  return fabs(m_det[(o*intervals+i)*(i_width*i_height) + (r*i_width+c)]);
}

//-------------------------------------------------------

//Interpolates a scale-space extremum's location and scale to subpixel
//accuracy to form an image feature.   
void FastHessian::interp_extremum(int octv, int intvl, int r, int c)
{
	double xi = 0, xr = 0, xc = 0;
	int i = 0;
  int step = init_sample * cvRound(pow(2.0f,octv));

  // Get the offsets to the actual location of the extremum
  interp_step( octv, intvl, r, c, &xi, &xr, &xc );
  
  // If point is sufficiently close to the actual extremum
  if( fabs( xi ) <= 0.5  &&  fabs( xr ) <= 0.5  &&  fabs( xc ) <= 0.5 )
  {
    // Create Ipoint and push onto Ipoints vector
    Ipoint ipt;
    ipt.x = (float) (c + step*xc);
    ipt.y = (float) (r + step*xr);
    ipt.scale = (float) ((1.2f/9.0f) * (3*(pow(2.0f, octv+1) * (intvl+xi+1)+1)));
    ipt.laplacian = getLaplacian(octv, intvl, c, r);
    ipts.push_back(ipt);
  }
}

//-------------------------------------------------------

//Performs one step of extremum interpolation. 
void FastHessian::interp_step( int octv, int intvl, int r, int c, double* xi, double* xr, double* xc )
{
	CvMat* dD, * H, * H_inv, X;
	double x[3] = { 0 };

	dD = deriv_3D( octv, intvl, r, c );
	H = hessian_3D( octv, intvl, r, c );
	H_inv = cvCreateMat( 3, 3, CV_64FC1 );
	cvInvert( H, H_inv, CV_SVD );
	cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );
	cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );

	cvReleaseMat( &dD );
	cvReleaseMat( &H );
	cvReleaseMat( &H_inv );

	*xi = x[2];
	*xr = x[1];
	*xc = x[0];
}

//-------------------------------------------------------

//Computes the partial derivatives in x, y, and scale of a pixel in the DoG
//scale space pyramid.
CvMat* FastHessian::deriv_3D( int octv, int intvl, int r, int c )
{
	CvMat* dI;
	double dx, dy, ds;
  int step = init_sample * cvRound(pow(2.0f,octv));

  dx = ( getValLowe(octv,intvl, r, c+step ) -
		getValLowe( octv,intvl, r, c-step ) ) / 2.0;
	dy = ( getValLowe( octv,intvl, r+step, c ) -
		getValLowe( octv,intvl, r-step, c ) ) / 2.0;
	ds = ( getValLowe( octv,intvl+1, r, c ) -
		getValLowe( octv,intvl-1, r, c ) ) / 2.0;

	dI = cvCreateMat( 3, 1, CV_64FC1 );
	cvmSet( dI, 0, 0, dx );
	cvmSet( dI, 1, 0, dy );
	cvmSet( dI, 2, 0, ds );

	return dI;
}

//-------------------------------------------------------

//Computes the 3D Hessian matrix for a pixel in the DoG scale space pyramid.
CvMat* FastHessian::hessian_3D(int octv, int intvl, int r, int c )
{
	CvMat* H;
	double v, dxx, dyy, dss, dxy, dxs, dys;
  int step = init_sample * cvRound(pow(2.0f,octv));

	v = getValLowe( octv,intvl, r, c );
	dxx = ( getValLowe( octv,intvl, r, c+step ) + 
			getValLowe( octv,intvl, r, c-step ) - 2 * v );
	dyy = ( getValLowe( octv,intvl, r+step, c ) +
			getValLowe( octv,intvl, r-step, c ) - 2 * v );
	dss = ( getValLowe( octv,intvl+1, r, c ) +
			getValLowe( octv,intvl-1, r, c ) - 2 * v );
	dxy = ( getValLowe( octv,intvl, r+step, c+step ) -
			getValLowe( octv,intvl, r+step, c-step ) -
			getValLowe( octv,intvl, r-step, c+step ) +
			getValLowe( octv,intvl, r-step, c-step ) ) / 4.0;
	dxs = ( getValLowe( octv,intvl+1, r, c+step ) -
			getValLowe( octv,intvl+1, r, c-step ) -
			getValLowe( octv,intvl-1, r, c+step ) +
			getValLowe( octv,intvl-1, r, c-step ) ) / 4.0;
	dys = ( getValLowe( octv,intvl+1, r+step, c ) -
			getValLowe( octv,intvl+1, r-step, c ) -
			getValLowe( octv,intvl-1, r+step, c ) +
			getValLowe( octv,intvl-1, r-step, c ) ) / 4.0;

	H = cvCreateMat( 3, 3, CV_64FC1 );
	cvmSet( H, 0, 0, dxx );
	cvmSet( H, 0, 1, dxy );
	cvmSet( H, 0, 2, dxs );
	cvmSet( H, 1, 0, dxy );
	cvmSet( H, 1, 1, dyy );
	cvmSet( H, 1, 2, dys );
	cvmSet( H, 2, 0, dxs );
	cvmSet( H, 2, 1, dys );
	cvmSet( H, 2, 2, dss );

	return H;
}

//-------------------------------------------------------

surf.cpp

#include "cv.h"
#include "highgui.h"
#include "integral.h"
#include "surf.h"
#include "utils.h"
#include <fstream>
#include <iostream>
using namespace std;

const float pi = float(CV_PI);

//! Destructor
Surf::~Surf()
{
}

//-------------------------------------------------------

//! Constructor
Surf::Surf(IplImage *img, IpVec &ipts)
: ipts(ipts)
{
  this->img = img;
}

//-------------------------------------------------------

//! Describe all features in the supplied vector
void Surf::getDescriptors(bool upright)
{
  // Check there are Ipoints to be described
  if (!ipts.size()) return;

  // Get the size of the vector for fixed loop bounds
  int ipts_size = ipts.size();

  if (upright)
  {
    // U-SURF loop just gets descriptors
    for (unsigned int i = 0; i < ipts_size; ++i)
    {
      // Set the Ipoint to be described
      index = i;

      // Extract upright (i.e. not rotation invariant) descriptors
      getUprightDescriptor();
    }
  }
  else
  {
    // Main SURF-64 loop assigns orientations and gets descriptors
    for (unsigned int i = 0; i < ipts_size; ++i)
    {
      // Set the Ipoint to be described
      index = i;

      // Assign Orientations and extract rotation invariant descriptors
      getOrientation();
      getDescriptor();
    }
  }
}

//-------------------------------------------------------

//! Assign the supplied Ipoint an orientation
void Surf::getOrientation()
{
  Ipoint *ipt = &ipts.at(index);
  float gauss, scale = ipt->scale;
  int s = cvRound(scale), r = cvRound(ipt->y), c = cvRound(ipt->x);
  std::vector<float> resX, resY, Ang;

  // calculate haar responses for points within radius of 6*scale
  for(int i = -6*s; i <= 6*s; i+=s) 
  {
    for(int j = -6*s; j <= 6*s; j+=s) 
    {
      if( i*i + j*j  < 36*s*s ) 
      {
        gauss = gaussian(i, j, 2.5f*s);
        resX.push_back( gauss * haarX(r+j, c+i, 4*s) );
        resY.push_back( gauss * haarY(r+j, c+i, 4*s) );
        Ang.push_back(getAngle(resX.back(), resY.back()));
      }
    }
  }

  // calculate the dominant direction 
  float sumX, sumY;
  float max=0, old_max = 0, orientation = 0, old_orientation = 0;
  float ang1, ang2, ang;

  // loop slides pi/3 window around feature point
  for(ang1 = 0; ang1 < 2*pi;  ang1+=0.2f) {
    ang2 = ( ang1+pi/3.0f > 2*pi ? ang1-5.0f*pi/3.0f : ang1+pi/3.0f);
    sumX = sumY = 0; 
    for(unsigned int k = 0; k < Ang.size(); k++) 
    {
      // get angle from the x-axis of the sample point
      ang = Ang.at(k);

      // determine whether the point is within the window
      if (ang1 < ang2 && ang1 < ang && ang < ang2) 
      {
        sumX+=resX.at(k);  
        sumY+=resY.at(k);
      } 
      else if (ang2 < ang1 && 
        ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2*pi) )) 
      {
        sumX+=resX.at(k);  
        sumY+=resY.at(k);
      }
    }

    // if the vector produced from this window is longer than all 
    // previous vectors then this forms the new dominant direction
    if (sumX*sumX + sumY*sumY > max) 
    {
      // store second largest orientation
      old_max = max;
      old_orientation = orientation;

      // store largest orientation
      max = sumX*sumX + sumY*sumY;
      orientation = getAngle(sumX, sumY);
    }
  }

  // check whether there are two dominant orientations based on 0.8 threshold
  if (old_max >= 0.8*max)
  {
    // assign second largest orientation and push copy onto vector
    ipt->orientation = old_orientation;
    ipts.push_back(*ipt);
    
    // Reset ipt to point to correct Ipoint in the vector
    ipt = &ipts.at(index);
  }

  // assign orientation of the dominant response vector
  ipt->orientation = orientation;
}

//-------------------------------------------------------

//! Get the descriptor vector of the provided Ipoint
void Surf::getDescriptor()
{
  int y, x, count=0;
  float scale, *desc, dx, dy, mdx, mdy, co, si;
  float gauss, rx, ry, rrx, rry, sample_x, sample_y, len=0;
  
  Ipoint *ipt = &ipts.at(index);
  scale = ipt->scale;
  x = cvRound(ipt->x);
  y = cvRound(ipt->y);  
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;

  // Calculate descriptor for this interest point
  for (int i = -10; i < 10; i+=5)
  {
    for (int j = -10; j < 10; j+=5) 
    {
      dx=dy=mdx=mdy=0;
      for (int k = i; k < i + 5; ++k) 
      {
        for (int l = j; l < j + 5; ++l) 
        {
          // Get coords of sample point on the rotated axis
          sample_x = cvRound(x + (-l*scale*si + k*scale*co));
          sample_y = cvRound(y + ( l*scale*co + k*scale*si));

          // Get the gaussian weighted x and y responses
          gauss = gaussian(k*scale, l*scale, 3.3f*scale);  
          rx = gauss * haarX(sample_y, sample_x, 2*scale);
          ry = gauss * haarY(sample_y, sample_x, 2*scale);

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);
  for( i = 0; i < 64; i++)
    desc[i] /= len;

}

//-------------------------------------------------------

//! Get the upright descriptor vector of the provided Ipoint
void Surf::getUprightDescriptor()
{
  int y, x, count=0;
  float scale, *desc, dx, dy, mdx, mdy;
  float gauss, rx, ry, len = 0.0f;

  Ipoint *ipt = &ipts.at(index);
  scale = ipt->scale;
  y = cvRound(ipt->y);  
  x = cvRound(ipt->x);
  desc = ipt->descriptor;

  // Calculate descriptor for this interest point
  for (int i = -10; i < 10; i+=5)
  {
    for (int j = -10; j < 10; j+=5) 
    {
      dx=dy=mdx=mdy=0;
      for (int k = i; k < i + 5; ++k) 
      {
        for (int l = j; l < j + 5; ++l) 
        {
          // get Gaussian weighted x and y responses
          gauss = gaussian(k*scale, l*scale, 3.3f*scale);  
          rx = gauss * haarX(k*scale+y, l*scale+x, 2*scale);
          ry = gauss * haarY(k*scale+y, l*scale+x, 2*scale);

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);
  for(i = 0; i < 64; i++)
    desc[i] /= len;

}

//-------------------------------------------------------

//! Calculate the value of the 2d gaussian at x,y
inline float Surf::gaussian(int x, int y, float sig)
{
  return 1.0f/(2.0f*pi*sig*sig) * exp( -(x*x+y*y)/(2.0f*sig*sig));
}

//-------------------------------------------------------

//! Calculate the value of the 2d gaussian at x,y
inline float Surf::gaussian(float x, float y, float sig)
{
  return 1.0f/(2.0f*pi*sig*sig) * exp( -(x*x+y*y)/(2.0f*sig*sig));
}

//-------------------------------------------------------

//! Calculate Haar wavelet responses in x direction
inline float Surf::haarX(int row, int column, int s)
{
  return Area(img, column, row-s/2, s/2, s) 
    -1 * Area(img, column-s/2, row-s/2, s/2, s);
}

//-------------------------------------------------------

//! Calculate Haar wavelet responses in y direction
inline float Surf::haarY(int row, int column, int s)
{
  return Area(img, column-s/2, row, s, s/2) 
    -1 * Area(img, column-s/2, row-s/2, s, s/2);
}

//-------------------------------------------------------

//! Get the angle from the +ve x-axis of the vector given by (X Y)
float Surf::getAngle(float X, float Y)
{
  if(X >= 0 && Y >= 0)
    return atan(Y/X);

  if(X < 0 && Y >= 0)
    return pi - atan(-Y/X);

  if(X < 0 && Y < 0)
    return pi + atan(Y/X);

  if(X >= 0 && Y < 0)
    return 2*pi - atan(-Y/X);

  return 0;
}

surflib.h

#ifndef SURFLIB_H
#define SURFLIB_H

#include "cv.h"
#include "highgui.h"
#include "integral.h"
#include "fasthessian.h"
#include "surf.h"
#include "ipoint.h"
#include "utils.h"

// typedef away STL mess!
typedef std::vector<Ipoint> IpVec;


//! Library function builds vector of described interest points
inline void surfDetDes(IplImage *img,  /* image to find Ipoints in */
                       std::vector<Ipoint> &ipts, /* reference to vector of Ipoints */
                       bool upright = false, /* run in rotation invariant mode? */
                       int octaves = OCTAVES, /* number of octaves to calculate */
                       int intervals = INTERVALS, /* number of intervals per octave */
                       int init_sample = INIT_SAMPLE, /* initial sampling step */
                       float thres = THRES /* blob response threshold */)
{
  // Create integral-image representation of the image
  IplImage *int_img = Integral(img);
  
  // Create Fast Hessian Object
  FastHessian fh(int_img, ipts, octaves, intervals, init_sample, thres);
 
  // Extract interest points and store in vector ipts
  fh.getIpoints();
  
  // Create Surf Descriptor Object
  Surf des(int_img, ipts);

  // Extract the descriptors for the ipts
  des.getDescriptors(upright);

  // Deallocate the integral image
  cvReleaseImage(&int_img);
}


//! Library function builds vector of interest points
inline void surfDet(IplImage *img,  /* image to find Ipoints in */
                    std::vector<Ipoint> &ipts, /* reference to vector of Ipoints */
                    int octaves = OCTAVES, /* number of octaves to calculate */
                    int intervals = INTERVALS, /* number of intervals per octave */
                    int init_sample = INIT_SAMPLE, /* initial sampling step */
                    float thres = THRES /* blob response threshold */)
{
  // Create integral image representation of the image
  IplImage *int_img = Integral(img);

  // Create Fast Hessian Object
  FastHessian fh(int_img, ipts, octaves, intervals, init_sample, thres);

  // Extract interest points and store in vector ipts
  fh.getIpoints();

  // Deallocate the integral image
  cvReleaseImage(&int_img);
}




//! Library function describes interest points in vector
inline void surfDes(IplImage *img,  /* image to find Ipoints in */
                    std::vector<Ipoint> &ipts, /* reference to vector of Ipoints */
                    bool upright = false) /* run in rotation invariant mode? */
{ 
  // Create integral image representation of the image
  IplImage *int_img = Integral(img);

  // Create Surf Descriptor Object
  Surf des(int_img, ipts);

  // Extract the descriptors for the ipts
  des.getDescriptors(upright);
  
  // Deallocate the integral image
  cvReleaseImage(&int_img);
}


#endif

utils.cpp

#include "cv.h"
#include "highgui.h"
#include "utils.h"
#include <iostream>
#include <fstream>
#include <time.h>
using namespace std;

//-------------------------------------------------------

//! Display error message and terminate program
void error(char *msg) 
{
  cout << "\nError: " << msg;
  getchar();
  exit(0);
}

//-------------------------------------------------------

//! Show the provided image and wait for keypress
void showImage(IplImage *img)
{
  cvNamedWindow("Surf", CV_WINDOW_AUTOSIZE); 
  cvShowImage("Surf", img);  
  cvWaitKey(0);
}

//-------------------------------------------------------

//! Show the provided image in titled window and wait for keypress
void showImage(char *title, IplImage *img)
{
  cvNamedWindow(title, CV_WINDOW_AUTOSIZE); 
  cvShowImage(title, img);  
  cvWaitKey(0);
}

//-------------------------------------------------------

// Convert image to single channel 32F
IplImage* getGray(IplImage *img)
{
  IplImage* gray8, * gray32;

  gray8 = cvCreateImage( cvGetSize(img), IPL_DEPTH_8U, 1 );
  gray32 = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

  if( img->nChannels == 1 )
    gray8 = (IplImage *) cvClone( img );
  else
    cvCvtColor( img, gray8, CV_BGR2GRAY );

  cvConvertScale( gray8, gray32, 1.0 / 255.0, 0 );

  cvReleaseImage( &gray8 );
  return gray32;
}

//-------------------------------------------------------

//! Draw all the Ipoints in the provided vector
void drawIpoints(IplImage *img, vector<Ipoint> &ipts, int tailSize)
{
  Ipoint *ipt;
  float s, o;
  int r1, c1, r2, c2, lap;

  for(unsigned int i = 0; i < ipts.size(); i++) 
  {
    ipt = &ipts.at(i);
    s = ((9.0f/1.2f) * ipt->scale) / 3.0f;
    o = ipt->orientation;
    lap = ipt->laplacian;
    r1 = cvRound(ipt->y);
    c1 = cvRound(ipt->x);
    c2 = cvRound(s * cos(o)) + c1;
    r2 = cvRound(s * sin(o)) + r1;

    // Green line indicates orientation
    cvLine(img, cvPoint(c1, r1), cvPoint(c2, r2), cvScalar(0, 255, 0));

    if (lap >= 0)
    { // Blue circles indicate light blobs on dark backgrounds
      cvCircle(img, cvPoint(c1,r1), cvRound(s), cvScalar(255, 0, 0),1);
    }
    else
    { // Red circles indicate light blobs on dark backgrounds
      cvCircle(img, cvPoint(c1,r1), cvRound(s), cvScalar(0, 0, 255),1);
    }

    // Draw motion from ipoint dx and dy
    if (tailSize)
    {
      cvLine(img, cvPoint(c1,r1),
        cvPoint(int(c1+ipt->dx*tailSize), int(r1+ipt->dy*tailSize)),
        cvScalar(255,255,255), 1);
    }
  }
}

//-------------------------------------------------------

//! Draw a single feature on the image
void drawIpoint(IplImage *img, Ipoint &ipt)
{
  float s, o;
  int r1, c1, r2, c2, lap;

  s = ((9.0f/1.2f) * ipt.scale) / 3.0f;
  s*=3;
  o = ipt.orientation;
  lap = ipt.laplacian;
  r1 = cvRound(ipt.y);
  c1 = cvRound(ipt.x);
  c2 = cvRound(s * cos(o)) + c1;
  r2 = cvRound(s * sin(o)) + r1;

  // Green line indicates orientation
  cvLine(img, cvPoint(c1, r1), cvPoint(c2, r2), cvScalar(0, 255, 0));

  if (lap >= 0)
  { // Blue circles indicate light blobs on dark backgrounds
    cvCircle(img, cvPoint(c1,r1), cvRound(s), cvScalar(255, 0, 0),1);
  }
  else
  { // Red circles indicate light blobs on dark backgrounds
    cvCircle(img, cvPoint(c1,r1), cvRound(s), cvScalar(0, 0, 255),1);
  }
}

//-------------------------------------------------------

//! Draw a single feature on the image
void drawPoint(IplImage *img, Ipoint &ipt)
{
  float s, o;
  int r1, c1;

  s = 4;
  o = ipt.orientation;
  r1 = cvRound(ipt.y);
  c1 = cvRound(ipt.x);

  cvCircle(img, cvPoint(c1,r1), cvRound(s), cvScalar(0, 255, 0), -1);
  cvCircle(img, cvPoint(c1,r1), cvRound(s+1), cvScalar(255, 0,0), 2);

}

//-------------------------------------------------------

//! Draw a single feature on the image
void drawPoints(IplImage *img, vector<Ipoint> &ipts)
{
  float s, o;
  int r1, c1;

  for(unsigned int i = 0; i < ipts.size(); i++) 
  {
    s = 10;
    o = ipts[i].orientation;
    r1 = cvRound(ipts[i].y);
    c1 = cvRound(ipts[i].x);

    cvCircle(img, cvPoint(c1,r1), cvRound(s), cvScalar(0, 255, 0), -1);
    cvCircle(img, cvPoint(c1,r1), cvRound(s+1), cvScalar(255, 0,0), 2);
  }
}

//-------------------------------------------------------

//! Draw descriptor windows around Ipoints in the provided vector
void drawWindows(IplImage *img, vector<Ipoint> &ipts)
{
  Ipoint *ipt;
  float s, o, cd, sd;
  int x, y;
  CvPoint2D32f src[4];

  for(int i = 0; i < cvRound(ipts.size()); i++) 
  {
    ipt = &ipts.at(i);
    s = (10 * ipt->scale);
    o = ipt->orientation;
    y = cvRound(ipt->y);
    x = cvRound(ipt->x);
    cd = cos(o);
    sd = sin(o);

    src[0].x=sd*s+cd*s+x;   src[0].y=-cd*s+sd*s+y;
    src[1].x=sd*s+cd*-s+x;  src[1].y=-cd*s+sd*-s+y;
    src[2].x=sd*-s+cd*-s+x; src[2].y=-cd*-s+sd*-s+y;
    src[3].x=sd*-s+cd*s+x;  src[3].y=-cd*-s+sd*s+y;

    // Draw orientation line
    cvLine(img, cvPoint(cvRound(x), cvRound(y)), 
      cvPoint(cvRound(s*cd + x), cvRound(s*sd + y)), cvScalar(0, 255, 0),1);

    // Draw box window around the point
    cvLine(img, cvPoint(cvRound(src[0].x), cvRound(src[0].y)), 
      cvPoint(cvRound(src[1].x), cvRound(src[1].y)), cvScalar(255, 0, 0),2);
    cvLine(img, cvPoint(cvRound(src[1].x), cvRound(src[1].y)), 
      cvPoint(cvRound(src[2].x), cvRound(src[2].y)), cvScalar(255, 0, 0),2);
    cvLine(img, cvPoint(cvRound(src[2].x), cvRound(src[2].y)), 
      cvPoint(cvRound(src[3].x), cvRound(src[3].y)), cvScalar(255, 0, 0),2);
    cvLine(img, cvPoint(cvRound(src[3].x), cvRound(src[3].y)), 
      cvPoint(cvRound(src[0].x), cvRound(src[0].y)), cvScalar(255, 0, 0),2);

  }
}

//-------------------------------------------------------

// Draw the FPS figure on the image (requires at least 2 calls)
void drawFPS(IplImage *img)
{
  static clock_t t;
  char fps_text[20];
  CvFont font;
  cvInitFont(&font,CV_FONT_HERSHEY_SIMPLEX|CV_FONT_ITALIC, 1.0,1.0,0,2);

  // Add fps figure
  float fps = (1.0f/(clock()-t) * CLOCKS_PER_SEC);

  // Update variables
  t=clock(); 

  // Get the figure as a string
  sprintf(fps_text,"FPS: %.2f",fps);

  // Draw the string on the image
  cvPutText (img,fps_text,cvPoint(10,25), &font, cvScalar(255,255,0));
}

//-------------------------------------------------------

//! Save the SURF features to file
void saveSurf(char *filename, vector<Ipoint> &ipts)
{
  ofstream outfile(filename);

  // output descriptor length
  outfile << "64\n";
  outfile << ipts.size() << "\n";

  // create output line as:  scale  x  y  des
  for(unsigned int i=0; i < ipts.size(); i++) 
  {
    outfile << ipts.at(i).scale << "  ";
    outfile << ipts.at(i).x << " ";
    outfile << ipts.at(i).y << " ";
    for(int j=0; j<64; j++)
      outfile << ipts.at(i).descriptor[j] << " ";

    outfile << "\n";
  }

  outfile.close();
}

//-------------------------------------------------------

//! Load the SURF features from file
void loadSurf(char *filename, vector<Ipoint> &ipts)
{
  int count;
  float temp;
  ifstream infile(filename);

  // output descriptor length
  infile >> count;
  infile >> count;

  for (int i = 0; i < count; i++) 
  {
    Ipoint ipt;

    infile >> ipt.x;
    infile >> ipt.y;
    infile >> temp;
    ipt.scale = temp; 
    ipt.laplacian = 1;

    for (int j = 0; j < 3; j++)
      infile >> temp;

    for (j = 0; j < 64; j++)
      infile >> ipt.descriptor[j];

    ipts.push_back(ipt);

  }
}

ipoint.h

#ifndef IPOINT_H
#define IPOINT_H

class Ipoint {

public:

  //! Destructor
  ~Ipoint() {};

  //! Constructor
  Ipoint() : orientation(0) {};

  //! Coordinates of the detected interest point
  float x, y;

  //! Detected scale
  float scale;

  //! Orientation measured anti-clockwise from +ve x-axis
  float orientation;

  //! Sign of laplacian for fast matching purposes
  int laplacian;

  //! Vector of descriptor components
  float descriptor[64];

  //! Placeholds for point motion (can be used for frame to frame motion analysis)
  float dx, dy;

  //! Used to store cluster index
  int clusterIndex;
};



#endif

fasthessian.h

#ifndef FASTHESSIAN_H
#define FASTHESSIAN_H

#include "cv.h"
#include "integral.h"
#include "ipoint.h"
#include <vector>

static const int OCTAVES = 3;
static const int INTERVALS = 4;
static const float THRES = 0.002f;
static const int INIT_SAMPLE = 2;


class FastHessian {
  
  public:
    
    //! Destructor
    ~FastHessian();

    //! Constructor without image
    FastHessian(std::vector<Ipoint> &ipts, 
                const int octaves = OCTAVES, 
                const int intervals = INTERVALS, 
                const int init_sample = INIT_SAMPLE, 
                const float thres = THRES);

    //! Constructor with image
    FastHessian(IplImage *img, 
                std::vector<Ipoint> &ipts, int octaves = OCTAVES, 
                const int intervals = INTERVALS, 
                const int init_sample = INIT_SAMPLE, 
                const float thres = THRES);

    //! Set or re-set the integral image source
    void setIntImage(IplImage *img);

    //! Find the image features and write into vector of features
    void getIpoints();
    
  private:

    //---------------- Private Functions -----------------//

    //! Calculate determinant of hessian responses
    void buildDet();
   
    //! Non Maximal Suppression function
    int isExtremum(int octave, int interval, int column, int row);    
    
    //! Return the value of the approximated determinant of hessian
    inline float getVal(int octave, int interval, int column, int row);

    //! Return the sign of the laplacian (trace of the hessian)
    inline int getLaplacian(int o, int i, int c, int r);

    //! Round an int to a float
    inline int fRound(float flt);

    //! Interpolation functions - adapted from Lowe's SIFT interpolation
    void interp_extremum(int octv, int intvl, int r, int c);
    void interp_step( int octv, int intvl, int r, int c, double* xi, double* xr, double* xc );
    CvMat* deriv_3D( int octv, int intvl, int r, int c );
    CvMat* hessian_3D(int octv, int intvl, int r, int c );
    inline float getValLowe(int octave, int interval, int column, int row);

    //---------------- Private Variables -----------------//

    //! Pointer to the integral Image, and its attributes 
    IplImage *img;
    int i_width, i_height;

    //! Reference to vector of features passed from outside 
    std::vector<Ipoint> &ipts;

    //! Number of Octaves
    int octaves;

    //! Number of Intervals per octave
    int intervals;

    //! Initial sampling step for Ipoint detection
    int init_sample;

    //! Threshold value for blob resonses
    float thres;

    //! Array stack of determinant of hessian values
    float *m_det;
};



#endif

utils.h

#ifndef UTILS_H
#define UTILS_H

#include "cv.h"
#include "ipoint.h"
#include <vector>


//! Display error message and terminate program
void error(char *msg);

//! Show the provided image and wait for keypress
void showImage(IplImage *img);

//! Show the provided image in titled window and wait for keypress
void showImage(char *title, IplImage *img);

// Convert image to single channel 32F
IplImage* getGray(IplImage *img);

//! Draw a single feature on the image
void drawIpoint(IplImage *img, Ipoint &ipt);

//! Draw all the Ipoints in the provided vector
void drawIpoints(IplImage *img, std::vector<Ipoint> &ipts, int tailSize = 0);

//! Draw descriptor windows around Ipoints in the provided vector
void drawWindows(IplImage *img, std::vector<Ipoint> &ipts);

// Draw the FPS figure on the image (requires at least 2 calls)
void drawFPS(IplImage *img);

//! Draw a Point at feature location
void drawPoint(IplImage *img, Ipoint &ipt);

//! Draw a Point at all features
void drawPoints(IplImage *img, std::vector<Ipoint> &ipts);

//! Save the SURF features to file
void saveSurf(char *filename, std::vector<Ipoint> &ipts);

//! Load the SURF features from file
void loadSurf(char *filename, std::vector<Ipoint> &ipts);

#endif

surf.h

#ifndef SURF_H
#define SURF_H

#include "cv.h"
#include "ipoint.h"
#include "integral.h"
#include <vector>

typedef std::vector<Ipoint> IpVec;

class Surf {
  
  public:
    
    //! Destructor
    ~Surf();

    //! Standard Constructor (img is an integral image)
    Surf(IplImage *img, std::vector<Ipoint> &ipts);

    //! Describe all features in the supplied vector
    void getDescriptors(bool upright);
  
  private:
    
    //---------------- Private Functions -----------------//

    //! Assign the current Ipoint an orientation
    void getOrientation();
    
    //! Get the descriptor vector of the current Ipoint
	  void getDescriptor();

    //! Get the upright descriptor vector of the current Ipoint
    void getUprightDescriptor();

    //! Calculate the value of the 2d gaussian at x,y
    inline float gaussian(int x, int y, float sig);
    inline float gaussian(float x, float y, float sig);

    //! Calculate Haar wavelet responses in x and y directions
    inline float haarX(int row, int column, int size);
    inline float haarY(int row, int column, int size);

    //! Get the angle from the +ve x-axis of the vector given by [X Y]
    float getAngle(float X, float Y);


    //---------------- Private Variables -----------------//

    //! Integral image where Ipoints have been detected
    IplImage *img;

    //! Ipoints vector
    IpVec &ipts;

    //! Index of current Ipoint in the vector
    int index;
};


#endif