在vs环境中跑动sift特征提取（代码部分）

因为在前两天的学习中发现。在opencv环境中跑动sift特征点提取还是比较困难的。

所以在此，进行记述。

遇到的问题分别有，csdn不愿意花费积分、配置gtk困难、教程海量然而能跑者鲜。描述不详尽等。

【然后我却是发现这个borwhess实在是不知道叫先生何名为好。】

话归正题。

以下跑动具体过程：

首先去：
http://blog.csdn.net/masibuaa/article/details/9246493
发现main.cpp

也就是：检测sift的部分。

这个回头慢慢凿。先跑起来：

sift主体：

img1_Feat = cvCloneImage(img1);//复制图1，深拷贝，用来画特征点
img2_Feat = cvCloneImage(img2);//复制图2，深拷贝，用来画特征点

//默认提取的是LOWE格式的SIFT特征点
//提取并显示第1幅图片上的特征点
n1 = sift_features( img1, &feat1 );//检测图1中的SIFT特征点,n1是图1的特征点个数
export_features("feature1.txt",feat1,n1);//将特征向量数据写入到文件
draw_features( img1_Feat, feat1, n1 );//画出特征点
cvNamedWindow(IMG1_FEAT);//创建窗口
cvShowImage(IMG1_FEAT,img1_Feat);//显示

//提取并显示第2幅图片上的特征点
n2 = sift_features( img2, &feat2 );//检测图2中的SIFT特征点，n2是图2的特征点个数
export_features("feature2.txt",feat2,n2);//将特征向量数据写入到文件
draw_features( img2_Feat, feat2, n2 );//画出特征点
cvNamedWindow(IMG2_FEAT);//创建窗口
cvShowImage(IMG2_FEAT,img2_Feat);//显示

这个基本算是函数主体，也就是要运行的主程。

由于调试过程中，调动多次：次序依然蒙蔽。

主体来自：
https://github.com/robwhess/opensift
然而，现在的win7装载gtk并不容易，我在读了2/5之后没有打算继续安装的想法。

遂西凑东拼，

勉强实现这样的结果：



算是sift，检测出来了。

调用函数主体：
http://www.cnblogs.com/freedomshe/archive/2012/04/28/2475057.html
一共七个文件：imgfeatures.h sift.h utils.h imgfeatures.c main.cpp sift.c utils.c 以下源代码。

sift.h:

/**@file
Functions for detecting SIFT image features.

For more information, refer to:

Lowe, D.  Distinctive image features from scale-invariant keypoints.
<EM>International Journal of Computer Vision, 60</EM>, 2 (2004),
pp.91--110.

Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>
Note: The SIFT algorithm is patented in the United States and cannot be
used in commercial products without a license from the University of
British Columbia.  For more information, refer to the file LICENSE.ubc
that accompanied this distribution.

@version 1.1.2-20100521
*/

#ifndef SIFT_H
#define SIFT_H

#include "cxcore.h"

/******************************** Structures *********************************/

/** holds feature data relevant to detection */
struct detection_data
{
int r;
int c;
int octv;
int intvl;
double subintvl;
double scl_octv;
};

struct feature;

/******************************* Defs and macros *****************************/

/** default number of sampled intervals per octave */
#define SIFT_INTVLS 3

/** default sigma for initial gaussian smoothing */
#define SIFT_SIGMA 1.6

/** default threshold on keypoint contrast |D(x)| */
#define SIFT_CONTR_THR 0.04

/** default threshold on keypoint ratio of principle curvatures */
#define SIFT_CURV_THR 10

/** double image size before pyramid construction? */
#define SIFT_IMG_DBL 1

/** default width of descriptor histogram array */
#define SIFT_DESCR_WIDTH 4

/** default number of bins per histogram in descriptor array */
#define SIFT_DESCR_HIST_BINS 8

/* assumed gaussian blur for input image */
#define SIFT_INIT_SIGMA 0.5

/* width of border in which to ignore keypoints */
#define SIFT_IMG_BORDER 5

/* maximum steps of keypoint interpolation before failure */
#define SIFT_MAX_INTERP_STEPS 5

/* default number of bins in histogram for orientation assignment */
#define SIFT_ORI_HIST_BINS 36

/* determines gaussian sigma for orientation assignment */
#define SIFT_ORI_SIG_FCTR 1.5

/* determines the radius of the region used in orientation assignment */
#define SIFT_ORI_RADIUS 3.0 * SIFT_ORI_SIG_FCTR

/* number of passes of orientation histogram smoothing */
#define SIFT_ORI_SMOOTH_PASSES 2

/* orientation magnitude relative to max that results in new feature */
#define SIFT_ORI_PEAK_RATIO 0.8

/* determines the size of a single descriptor orientation histogram */
#define SIFT_DESCR_SCL_FCTR 3.0

/* threshold on magnitude of elements of descriptor vector */
#define SIFT_DESCR_MAG_THR 0.2

/* factor used to convert floating-point descriptor to unsigned char */
#define SIFT_INT_DESCR_FCTR 512.0

/* returns a feature's detection data */
#define feat_detection_data(f) ( (struct detection_data*)(f->feature_data) )

/*************************** Function Prototypes *****************************/

#ifdef __cplusplus
extern "C" {
#endif

/**
Finds SIFT features in an image using default parameter values.  All
detected features are stored in the array pointed to by \a feat.

@param img the image in which to detect features
@param feat a pointer to an array in which to store detected features; memory
for this array is allocated by this function and must be freed by the caller
using free(*feat)

@return Returns the number of features stored in \a feat or -1 on failure
@see _sift_features()
*/
extern int sift_features(IplImage* img, struct feature** feat);

/**
Finda SIFT features in an image using user-specified parameter values.  All
detected features are stored in the array pointed to by \a feat.

@param img the image in which to detect features
@param feat a pointer to an array in which to store detected features; memory
for this array is allocated by this function and must be freed by the caller
using free(*feat)
@param intvls the number of intervals sampled per octave of scale space
@param sigma the amount of Gaussian smoothing applied to each image level
before building the scale space representation for an octave
@param contr_thr a threshold on the value of the scale space function
\f$\left|D(\hat{x})\right|\f$, where \f$\hat{x}\f$ is a vector specifying
feature location and scale, used to reject unstable features;  assumes
pixel values in the range [0, 1]
@param curv_thr threshold on a feature's ratio of principle curvatures
used to reject features that are too edge-like
@param img_dbl should be 1 if image doubling prior to scale space
construction is desired or 0 if not
@param descr_width the width, \f$n\f$, of the \f$n \times n\f$ array of
orientation histograms used to compute a feature's descriptor
@param descr_hist_bins the number of orientations in each of the
histograms in the array used to compute a feature's descriptor

@return Returns the number of keypoints stored in \a feat or -1 on failure
@see sift_features()
*/
extern int _sift_features(IplImage* img, struct feature** feat, int intvls,
double sigma, double contr_thr, int curv_thr,
int img_dbl, int descr_width, int descr_hist_bins);

#ifdef __cplusplus
}
#endif

#endif

imgfeatures.h

/**@file
Functions and structures for dealing with image features

Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>

@version 1.1.2-20100521
*/

#ifndef IMGFEATURES_H
#define IMGFEATURES_H

#include "cxcore.h"

/** FEATURE_OXFD <BR> FEATURE_LOWE */
enum feature_type
{
FEATURE_OXFD,
FEATURE_LOWE,
};

/** FEATURE_FWD_MATCH <BR> FEATURE_BCK_MATCH <BR> FEATURE_MDL_MATCH */
enum feature_match_type
{
FEATURE_FWD_MATCH,
FEATURE_BCK_MATCH,
FEATURE_MDL_MATCH,
};

/* colors in which to display different feature types */
#define FEATURE_OXFD_COLOR CV_RGB(255,255,0)
#define FEATURE_LOWE_COLOR CV_RGB(255,0,255)

/** max feature descriptor length */
#define FEATURE_MAX_D 128

/**
Structure to represent an affine invariant image feature.  The fields
x, y, a, b, c represent the affine region around the feature:

a(x-u)(x-u) + 2b(x-u)(y-v) + c(y-v)(y-v) = 1
*/
struct feature
{
double x;                      /**< x coord */
double y;                      /**< y coord */
double a;                      /**< Oxford-type affine region parameter */
double b;                      /**< Oxford-type affine region parameter */
double c;                      /**< Oxford-type affine region parameter */
double scl;                    /**< scale of a Lowe-style feature */
double ori;                    /**< orientation of a Lowe-style feature */
int d;                         /**< descriptor length */
double descr[FEATURE_MAX_D];   /**< descriptor */
int type;                      /**< feature type, OXFD or LOWE */
int category;                  /**< all-purpose feature category */
struct feature* fwd_match;     /**< matching feature from forward image */
struct feature* bck_match;     /**< matching feature from backmward image */
struct feature* mdl_match;     /**< matching feature from model */
CvPoint2D64f img_pt;           /**< location in image */
CvPoint2D64f mdl_pt;           /**< location in model */
void* feature_data;            /**< user-definable data */
};

#ifdef __cplusplus
extern "C" {
#endif

/**
Reads image features from file.  The file should be formatted as from
the code provided by the Visual Geometry Group at Oxford or from the
code provided by David Lowe.

@param filename location of a file containing image features
@param type determines how features are input.  If \a type is FEATURE_OXFD,
the input file is treated as if it is from the code provided by the VGG
at Oxford: http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html <BR><BR>
If \a type is FEATURE_LOWE, the input file is treated as if it is from
David Lowe's SIFT code: http://www.cs.ubc.ca/~lowe/keypoints @param feat pointer to an array in which to store imported features; memory for
this array is allocated by this function and must be freed by the caller using
free(*feat)

@return Returns the number of features imported from filename or -1 on error
*/
extern int import_features(char* filename, int type, struct feature** feat);

/**
Exports a feature set to a file formatted depending on the type of
features, as specified in the feature struct's type field.

@param filename name of file to which to export features
@param feat feature array
@param n number of features

@return Returns 0 on success or 1 on error
*/
extern int export_features(char* filename, struct feature* feat, int n);

/**
Displays a set of features on an image

@param img image on which to display features
@param feat array of Oxford-type features
@param n number of features
*/
extern void draw_features(IplImage* img, struct feature* feat, int n);

/**
Calculates the squared Euclidian distance between two feature descriptors.

@param f1 first feature
@param f2 second feature

@return Returns the squared Euclidian distance between the descriptors of
\a f1 and \a f2.
*/
extern double descr_dist_sq(struct feature* f1, struct feature* f2);

#ifdef __cplusplus
}
#endif

#endif

utils.h

/**@file
Miscellaneous utility functions.

Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>

@version 1.1.2-20100521
*/

#ifndef UTILS_H
#define UTILS_H

#include "cxcore.h"

#include <stdio.h>

/* absolute value */
#ifndef ABS
#define ABS(x) ( ( (x) < 0 )? -(x) : (x) )
#endif

/***************************** Inline Functions ******************************/

/**
A function to get a pixel value from an 8-bit unsigned image.

@param img an image
@param r row
@param c column
@return Returns the value of the pixel at (\a r, \a c) in \a img
*/
static __inline int pixval8(IplImage* img, int r, int c)
{
return (int)(((uchar*)(img->imageData + img->widthStep*r))[c]);
}

/**
A function to set a pixel value in an 8-bit unsigned image.

@param img an image
@param r row
@param c column
@param val pixel value
*/
static __inline void setpix8(IplImage* img, int r, int c, uchar val)
{
((uchar*)(img->imageData + img->widthStep*r))[c] = val;
}

/**
A function to get a pixel value from a 32-bit floating-point image.

@param img an image
@param r row
@param c column
@return Returns the value of the pixel at (\a r, \a c) in \a img
*/
static __inline float pixval32f(IplImage* img, int r, int c)
{
return ((float*)(img->imageData + img->widthStep*r))[c];
}

/**
A function to set a pixel value in a 32-bit floating-point image.

@param img an image
@param r row
@param c column
@param val pixel value
*/
static __inline void setpix32f(IplImage* img, int r, int c, float val)
{
((float*)(img->imageData + img->widthStep*r))[c] = val;
}

/**
A function to get a pixel value from a 64-bit floating-point image.

@param img an image
@param r row
@param c column
@return Returns the value of the pixel at (\a r, \a c) in \a img
*/
static __inline double pixval64f(IplImage* img, int r, int c)
{
return (double)(((double*)(img->imageData + img->widthStep*r))[c]);
}

/**
A function to set a pixel value in a 64-bit floating-point image.

@param img an image
@param r row
@param c column
@param val pixel value
*/
static __inline void setpix64f(IplImage* img, int r, int c, double val)
{
((double*)(img->imageData + img->widthStep*r))[c] = val;
}

/**************************** Function Prototypes ****************************/

/**
Prints an error message and aborts the program.  The error message is
of the form "Error: ...", where the ... is specified by the \a format
argument

@param format an error message format string (as with \c printf(3)).
*/
extern void fatal_error(char* format, ...);

/**
Replaces a file's extension, which is assumed to be everything after the
last dot ('.') character.

@param file the name of a file

@param extn a new extension for \a file; should not include a dot (i.e.
\c "jpg", not \c ".jpg") unless the new file extension should contain
two dots.

@return Returns a new string formed as described above.  If \a file does
not have an extension, this function simply adds one.
*/
extern char* replace_extension(const char* file, const char* extn);

/**
A function that removes the path from a filename.  Similar to the Unix
basename command.

@param pathname a (full) path name

@return Returns the basename of \a pathname.
*/
extern char* basename(const char* pathname);

/**
Displays progress in the console with a spinning pinwheel.  Every time this
function is called, the state of the pinwheel is incremented.  The pinwheel
has four states that loop indefinitely: '|', '/', '-', '\'.

@param done if 0, this function simply increments the state of the pinwheel;
otherwise it prints "done"
*/
extern void progress(int done);

/**
Erases a specified number of characters from a stream.

@param stream the stream from which to erase characters
@param n the number of characters to erase
*/
extern void erase_from_stream(FILE* stream, int n);

/**
Doubles the size of an array with error checking

@param array pointer to an array whose size is to be doubled
@param n number of elements allocated for \a array
@param size size in bytes of elements in \a array

@return Returns the new number of elements allocated for \a array.  If no
memory is available, returns 0 and frees array.
*/
extern int array_double(void** array, int n, int size);

/**
Calculates the squared distance between two points.

@param p1 a point
@param p2 another point
*/
extern double dist_sq_2D(CvPoint2D64f p1, CvPoint2D64f p2);

/**
Draws an x on an image.

@param img an image
@param pt the center point of the x
@param r the x's radius
@param w the x's line weight
@param color the color of the x
*/
extern void draw_x(IplImage* img, CvPoint pt, int r, int w, CvScalar color);

/**
Combines two images by scacking one on top of the other

@param img1 top image
@param img2 bottom image

@return Returns the image resulting from stacking \a img1 on top if \a img2
*/
extern IplImage* stack_imgs(IplImage* img1, IplImage* img2);

/**
Allows user to view an array of images as a video.  Keyboard controls
are as follows:

<ul>
<li>Space - start and pause playback</li>
<li>Page Up - skip forward 10 frames</li>
<li>Page Down - jump back 10 frames</li>
<li>Right Arrow - skip forward 1 frame</li>
<li>Left Arrow - jump back 1 frame</li>
<li>Backspace - jump back to beginning</li>
<li>Esc - exit playback</li>
<li>Closing the window also exits playback</li>
</ul>

@param imgs an array of images
@param n number of images in \a imgs
@param win_name name of window in which images are displayed
*/
extern void vid_view(IplImage** imgs, int n, char* win_name);

/**
Checks if a HighGUI window is still open or not

@param name the name of the window we're checking

@return Returns 1 if the window named \a name has been closed or 0 otherwise
*/
extern int win_closed(char* name);

#endif

imgfeatures.c

/*
Functions and structures for dealing with image features

Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>

@version 1.1.2-20100521
*/

#include "utils.h"
#include "imgfeatures.h"

#include <cxcore.h>

#include <math.h>

static int import_oxfd_features(char*, struct feature**);
static int export_oxfd_features(char*, struct feature*, int);
static void draw_oxfd_features(IplImage*, struct feature*, int);
static void draw_oxfd_feature(IplImage*, struct feature*, CvScalar);

static int import_lowe_features(char*, struct feature**);
static int export_lowe_features(char*, struct feature*, int);
static void draw_lowe_features(IplImage*, struct feature*, int);
static void draw_lowe_feature(IplImage*, struct feature*, CvScalar);

/*
Reads image features from file.  The file should be formatted as from
the code provided by the Visual Geometry Group at Oxford:

@param filename location of a file containing image features
@param type determines how features are input.  If \a type is FEATURE_OXFD,
the input file is treated as if it is from the code provided by the VGG
at Oxford:
 http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html 
If \a type is FEATURE_LOWE, the input file is treated as if it is from
David Lowe's SIFT code:
 http://www.cs.ubc.ca/~lowe/keypoints @param features pointer to an array in which to store features

@return Returns the number of features imported from filename or -1 on error
*/
int import_features(char* filename, int type, struct feature** feat)
{
int n;

switch (type)
{
case FEATURE_OXFD:
n = import_oxfd_features(filename, feat);
break;
case FEATURE_LOWE:
n = import_lowe_features(filename, feat);
break;
default:
fprintf(stderr, "Warning: import_features(): unrecognized feature" \
"type, %s, line %d\n", __FILE__, __LINE__);
return -1;
}

if (n == -1)
fprintf(stderr, "Warning: unable to import features from %s,"    \
" %s, line %d\n", filename, __FILE__, __LINE__);
return n;
}

/*
Exports a feature set to a file formatted depending on the type of
features, as specified in the feature struct's type field.

@param filename name of file to which to export features
@param feat feature array
@param n number of features

@return Returns 0 on success or 1 on error
*/
int export_features(char* filename, struct feature* feat, int n)
{
int r, type;

if (n <= 0 || !feat)
{
fprintf(stderr, "Warning: no features to export, %s line %d\n",
__FILE__, __LINE__);
return 1;
}
type = feat[0].type;
switch (type)
{
case FEATURE_OXFD:
r = export_oxfd_features(filename, feat, n);
break;
case FEATURE_LOWE:
r = export_lowe_features(filename, feat, n);
break;
default:
fprintf(stderr, "Warning: export_features(): unrecognized feature" \
"type, %s, line %d\n", __FILE__, __LINE__);
return -1;
}

if (r)
fprintf(stderr, "Warning: unable to export features to %s,"    \
" %s, line %d\n", filename, __FILE__, __LINE__);
return r;
}

/*
Draws a set of features on an image

@param img image on which to draw features
@param feat array of Oxford-type features
@param n number of features
*/
void draw_features(IplImage* img, struct feature* feat, int n)
{
int type;

if (n <= 0 || !feat)
{
fprintf(stderr, "Warning: no features to draw, %s line %d\n",
__FILE__, __LINE__);
return;
}
type = feat[0].type;
switch (type)
{
case FEATURE_OXFD:
draw_oxfd_features(img, feat, n);
break;
case FEATURE_LOWE:
draw_lowe_features(img, feat, n);
break;
default:
fprintf(stderr, "Warning: draw_features(): unrecognized feature" \
" type, %s, line %d\n", __FILE__, __LINE__);
break;
}
}

/*
Calculates the squared Euclidian distance between two feature descriptors.

@param f1 first feature
@param f2 second feature

@return Returns the squared Euclidian distance between the descriptors of
f1 and f2.
*/
double descr_dist_sq(struct feature* f1, struct feature* f2)
{
double diff, dsq = 0;
double* descr1, *descr2;
int i, d;

d = f1->d;
if (f2->d != d)
return DBL_MAX;
descr1 = f1->descr;
descr2 = f2->descr;

for (i = 0; i < d; i++)
{
diff = descr1[i] - descr2[i];
dsq += diff*diff;
}
return dsq;
}

/***************************** Local Functions *******************************/

/*
Reads image features from file.  The file should be formatted as from
the code provided by the Visual Geometry Group at Oxford:
 http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html 
@param filename location of a file containing image features
@param features pointer to an array in which to store features

@return Returns the number of features imported from filename or -1 on error
*/
static int import_oxfd_features(char* filename, struct feature** features)
{
struct feature* f;
int i, j, n, d;
double x, y, a, b, c, dv;
FILE* file;

if (!features)
fatal_error("NULL pointer error, %s, line %d", __FILE__, __LINE__);

if (!(file = fopen(filename, "r")))
{
fprintf(stderr, "Warning: error opening %s, %s, line %d\n",
filename, __FILE__, __LINE__);
return -1;
}

/* read dimension and number of features */
if (fscanf(file, " %d %d ", &d, &n) != 2)
{
fprintf(stderr, "Warning: file read error, %s, line %d\n",
__FILE__, __LINE__);
return -1;
}
if (d > FEATURE_MAX_D)
{
fprintf(stderr, "Warning: descriptor too long, %s, line %d\n",
__FILE__, __LINE__);
return -1;
}

f = calloc(n, sizeof(struct feature));
for (i = 0; i < n; i++)
{
/* read affine region parameters */
if (fscanf(file, " %lf %lf %lf %lf %lf ", &x, &y, &a, &b, &c) != 5)
{
fprintf(stderr, "Warning: error reading feature #%d, %s, line %d\n",
i + 1, __FILE__, __LINE__);
free(f);
return -1;
}
f[i].img_pt.x = f[i].x = x;
f[i].img_pt.y = f[i].y = y;
f[i].a = a;
f[i].b = b;
f[i].c = c;
f[i].d = d;
f[i].type = FEATURE_OXFD;

/* read descriptor */
for (j = 0; j < d; j++)
{
if (!fscanf(file, " %lf ", &dv))
{
fprintf(stderr, "Warning: error reading feature descriptor" \
" #%d, %s, line %d\n", i + 1, __FILE__, __LINE__);
free(f);
return -1;
}
f[i].descr[j] = dv;
}

f[i].scl = f[i].ori = 0;
f[i].category = 0;
f[i].fwd_match = f[i].bck_match = f[i].mdl_match = NULL;
f[i].mdl_pt.x = f[i].mdl_pt.y = -1;
f[i].feature_data = NULL;
}

if (fclose(file))
{
fprintf(stderr, "Warning: file close error, %s, line %d\n",
__FILE__, __LINE__);
free(f);
return -1;
}

*features = f;
return n;
}

/*
Exports a feature set to a file formatted as one from the code provided
by the Visual Geometry Group at Oxford:
 http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html 
@param filename name of file to which to export features
@param feat feature array
@param n number of features

@return Returns 0 on success or 1 on error
*/
static int export_oxfd_features(char* filename, struct feature* feat, int n)
{
FILE* file;
int i, j, d;

if (n <= 0)
{
fprintf(stderr, "Warning: feature count %d, %s, line %s\n",
n, __FILE__, __LINE__);
return 1;
}
if (!(file = fopen(filename, "w")))
{
fprintf(stderr, "Warning: error opening %s, %s, line %d\n",
filename, __FILE__, __LINE__);
return 1;
}

d = feat[0].d;
fprintf(file, "%d\n%d\n", d, n);
for (i = 0; i < n; i++)
{
fprintf(file, "%f %f %f %f %f", feat[i].x, feat[i].y, feat[i].a,
feat[i].b, feat[i].c);
for (j = 0; j < d; j++)
fprintf(file, " %f", feat[i].descr[j]);
fprintf(file, "\n");
}

if (fclose(file))
{
fprintf(stderr, "Warning: file close error, %s, line %d\n",
__FILE__, __LINE__);
return 1;
}

return 0;
}

/*
Draws Oxford-type affine features

@param img image on which to draw features
@param feat array of Oxford-type features
@param n number of features
*/
static void draw_oxfd_features(IplImage* img, struct feature* feat, int n)
{
CvScalar color = CV_RGB(255, 255, 255);
int i;

if (img->nChannels > 1)
color = FEATURE_OXFD_COLOR;
for (i = 0; i < n; i++)
draw_oxfd_feature(img, feat + i, color);
}

/*
Draws a single Oxford-type feature

@param img image on which to draw
@param feat feature to be drawn
@param color color in which to draw
*/
static void draw_oxfd_feature(IplImage* img, struct feature* feat, CvScalar color)
{
double m[4] = { feat->a, feat->b, feat->b, feat->c };
double v[4] = { 0 };
double e[2] = { 0 };
CvMat M, V, E;
double alpha, l1, l2;

/* compute axes and orientation of ellipse surrounding affine region */
cvInitMatHeader(&M, 2, 2, CV_64FC1, m, CV_AUTOSTEP);
cvInitMatHeader(&V, 2, 2, CV_64FC1, v, CV_AUTOSTEP);
cvInitMatHeader(&E, 2, 1, CV_64FC1, e, CV_AUTOSTEP);
cvEigenVV(&M, &V, &E, DBL_EPSILON, 0, 0);
l1 = 1 / sqrt(e[1]);
l2 = 1 / sqrt(e[0]);
alpha = -atan2(v[1], v[0]);
alpha *= 180 / CV_PI;

cvEllipse(img, cvPoint(feat->x, feat->y), cvSize(l2, l1), alpha,
0, 360, CV_RGB(0, 0, 0), 3, 8, 0);
cvEllipse(img, cvPoint(feat->x, feat->y), cvSize(l2, l1), alpha,
0, 360, color, 1, 8, 0);
cvLine(img, cvPoint(feat->x + 2, feat->y), cvPoint(feat->x - 2, feat->y),
color, 1, 8, 0);
cvLine(img, cvPoint(feat->x, feat->y + 2), cvPoint(feat->x, feat->y - 2),
color, 1, 8, 0);
}

/*
Reads image features from file.  The file should be formatted as from
the code provided by David Lowe:
 http://www.cs.ubc.ca/~lowe/keypoints/ 
@param filename location of a file containing image features
@param features pointer to an array in which to store features

@return Returns the number of features imported from filename or -1 on error
*/
static int import_lowe_features(char* filename, struct feature** features)
{
struct feature* f;
int i, j, n, d;
double x, y, s, o, dv;
FILE* file;

if (!features)
fatal_error("NULL pointer error, %s, line %d", __FILE__, __LINE__);

if (!(file = fopen(filename, "r")))
{
fprintf(stderr, "Warning: error opening %s, %s, line %d\n",
filename, __FILE__, __LINE__);
return -1;
}

/* read number of features and dimension */
if (fscanf(file, " %d %d ", &n, &d) != 2)
{
fprintf(stderr, "Warning: file read error, %s, line %d\n",
__FILE__, __LINE__);
return -1;
}
if (d > FEATURE_MAX_D)
{
fprintf(stderr, "Warning: descriptor too long, %s, line %d\n",
__FILE__, __LINE__);
return -1;
}

f = calloc(n, sizeof(struct feature));
for (i = 0; i < n; i++)
{
/* read affine region parameters */
if (fscanf(file, " %lf %lf %lf %lf ", &y, &x, &s, &o) != 4)
{
fprintf(stderr, "Warning: error reading feature #%d, %s, line %d\n",
i + 1, __FILE__, __LINE__);
free(f);
return -1;
}
f[i].img_pt.x = f[i].x = x;
f[i].img_pt.y = f[i].y = y;
f[i].scl = s;
f[i].ori = o;
f[i].d = d;
f[i].type = FEATURE_LOWE;

/* read descriptor */
for (j = 0; j < d; j++)
{
if (!fscanf(file, " %lf ", &dv))
{
fprintf(stderr, "Warning: error reading feature descriptor" \
" #%d, %s, line %d\n", i + 1, __FILE__, __LINE__);
free(f);
return -1;
}
f[i].descr[j] = dv;
}

f[i].a = f[i].b = f[i].c = 0;
f[i].category = 0;
f[i].fwd_match = f[i].bck_match = f[i].mdl_match = NULL;
f[i].mdl_pt.x = f[i].mdl_pt.y = -1;
}

if (fclose(file))
{
fprintf(stderr, "Warning: file close error, %s, line %d\n",
__FILE__, __LINE__);
free(f);
return -1;
}

*features = f;
return n;
}

/*
Exports a feature set to a file formatted as one from the code provided
by David Lowe:
 http://www.cs.ubc.ca/~lowe/keypoints/ 
@param filename name of file to which to export features
@param feat feature array
@param n number of features

@return Returns 0 on success or 1 on error
*/
static int export_lowe_features(char* filename, struct feature* feat, int n)
{
FILE* file;
int i, j, d;

if (n <= 0)
{
fprintf(stderr, "Warning: feature count %d, %s, line %s\n",
n, __FILE__, __LINE__);
return 1;
}
if (!(file = fopen(filename, "w")))
{
fprintf(stderr, "Warning: error opening %s, %s, line %d\n",
filename, __FILE__, __LINE__);
return 1;
}

d = feat[0].d;
fprintf(file, "%d %d\n", n, d);
for (i = 0; i < n; i++)
{
fprintf(file, "%f %f %f %f", feat[i].y, feat[i].x,
feat[i].scl, feat[i].ori);
for (j = 0; j < d; j++)
{
/* write 20 descriptor values per line */
if (j % 20 == 0)
fprintf(file, "\n");
fprintf(file, " %d", (int)(feat[i].descr[j]));
}
fprintf(file, "\n");
}

if (fclose(file))
{
fprintf(stderr, "Warning: file close error, %s, line %d\n",
__FILE__, __LINE__);
return 1;
}

return 0;
}

/*
Draws Lowe-type features

@param img image on which to draw features
@param feat array of Oxford-type features
@param n number of features
*/
static void draw_lowe_features(IplImage* img, struct feature* feat, int n)
{
CvScalar color = CV_RGB(255, 255, 255);
int i;

if (img->nChannels > 1)
color = FEATURE_LOWE_COLOR;
for (i = 0; i < n; i++)
draw_lowe_feature(img, feat + i, color);
}

/*
Draws a single Lowe-type feature

@param img image on which to draw
@param feat feature to be drawn
@param color color in which to draw
*/
static void draw_lowe_feature(IplImage* img, struct feature* feat, CvScalar color)
{
int len, hlen, blen, start_x, start_y, end_x, end_y, h1_x, h1_y, h2_x, h2_y;
double scl, ori;
double scale = 5.0;
double hscale = 0.75;
CvPoint start, end, h1, h2;

/* compute points for an arrow scaled and rotated by feat's scl and ori */
start_x = cvRound(feat->x);
start_y = cvRound(feat->y);
scl = feat->scl;
ori = feat->ori;
len = cvRound(scl * scale);
hlen = cvRound(scl * hscale);
blen = len - hlen;
end_x = cvRound(len *  cos(ori)) + start_x;
end_y = cvRound(len * -sin(ori)) + start_y;
h1_x = cvRound(blen *  cos(ori + CV_PI / 18.0)) + start_x;
h1_y = cvRound(blen * -sin(ori + CV_PI / 18.0)) + start_y;
h2_x = cvRound(blen *  cos(ori - CV_PI / 18.0)) + start_x;
h2_y = cvRound(blen * -sin(ori - CV_PI / 18.0)) + start_y;
start = cvPoint(start_x, start_y);
end = cvPoint(end_x, end_y);
h1 = cvPoint(h1_x, h1_y);
h2 = cvPoint(h2_x, h2_y);

cvLine(img, start, end, color, 1, 8, 0);
cvLine(img, end, h1, color, 1, 8, 0);
cvLine(img, end, h2, color, 1, 8, 0);
}

main.cpp

#include "opencv2/highgui/highgui.hpp"//这个是基本类库，高级图形用户接口，必须要引入
#include "sift.h"//解决了没有sift_features函数的问题
#include "imgfeatures.h"//解决了找不到export_features 和draw_features 的问题
#include "utils.h"

int main(){

IplImage * img1 = cvLoadImage("a.png");
IplImage * img2 = cvLoadImage("b.png");

IplImage *img1_Feat = cvCloneImage(img1);//复制图1，深拷贝，用来画特征点
IplImage *img2_Feat = cvCloneImage(img2);//复制图2，深拷贝，用来画特征点

img1_Feat = cvCloneImage(img1);//复制图1，深拷贝，用来画特征点
img2_Feat = cvCloneImage(img2);//复制图2，深拷贝，用来画特征点

//默认提取的是LOWE格式的SIFT特征点
//提取并显示第1幅图片上的特征点

feature *feat1, *feat2;

int n1 = sift_features(img1, &feat1);//检测图1中的SIFT特征点,n1是图1的特征点个数
export_features("feature1.txt", feat1, n1);//将特征向量数据写入到文件
draw_features(img1_Feat, feat1, n1);//画出特征点
cvNamedWindow("IMG1_FEAT");//创建窗口
cvShowImage("IMG1_FEAT", img1_Feat);//显示

//提取并显示第2幅图片上的特征点
int n2 = sift_features(img2, &feat2);//检测图2中的SIFT特征点，n2是图2的特征点个数
export_features("feature2.txt", feat2, n2);//将特征向量数据写入到文件
draw_features(img2_Feat, feat2, n2);//画出特征点
cvNamedWindow("IMG2_FEAT");//创建窗口
cvShowImage("IMG2_FEAT", img2_Feat);//显示

cvWaitKey(0);

}

sift.c

/*
Functions for detecting SIFT image features.

For more information, refer to:

Lowe, D.  Distinctive image features from scale-invariant keypoints.
<EM>International Journal of Computer Vision, 60</EM>, 2 (2004),
pp.91--110.

Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>

Note: The SIFT algorithm is patented in the United States and cannot be
used in commercial products without a license from the University of
British Columbia.  For more information, refer to the file LICENSE.ubc
that accompanied this distribution.

@version 1.1.2-20100521
*/

#include "sift.h"
#include "imgfeatures.h"
#include "utils.h"

#include <cxcore.h>
#include <cv.h>

/************************* Local Function Prototypes *************************/

static IplImage* create_init_img(IplImage*, int, double);
static IplImage* convert_to_gray32(IplImage*);
static IplImage*** build_gauss_pyr(IplImage*, int, int, double);
static IplImage* downsample(IplImage*);
static IplImage*** build_dog_pyr(IplImage***, int, int);
static CvSeq* scale_space_extrema(IplImage***, int, int, double, int, CvMemStorage*);
static int is_extremum(IplImage***, int, int, int, int);
static struct feature* interp_extremum(IplImage***, int, int, int, int, int, double);
static void interp_step(IplImage***, int, int, int, int, double*, double*, double*);
static CvMat* deriv_3D(IplImage***, int, int, int, int);
static CvMat* hessian_3D(IplImage***, int, int, int, int);
static double interp_contr(IplImage***, int, int, int, int, double, double, double);
static struct feature* new_feature(void);
static int is_too_edge_like(IplImage*, int, int, int);
static void calc_feature_scales(CvSeq*, double, int);
static void adjust_for_img_dbl(CvSeq*);
static void calc_feature_oris(CvSeq*, IplImage***);
static double* ori_hist(IplImage*, int, int, int, int, double);
static int calc_grad_mag_ori(IplImage*, int, int, double*, double*);
static void smooth_ori_hist(double*, int);
static double dominant_ori(double*, int);
static void add_good_ori_features(CvSeq*, double*, int, double, struct feature*);
static struct feature* clone_feature(struct feature*);
static void compute_descriptors(CvSeq*, IplImage***, int, int);
static double*** descr_hist(IplImage*, int, int, double, double, int, int);
static void interp_hist_entry(double***, double, double, double, double, int, int);
static void hist_to_descr(double***, int, int, struct feature*);
static void normalize_descr(struct feature*);
static int feature_cmp(void*, void*, void*);
static void release_descr_hist(double****, int);
static void release_pyr(IplImage****, int, int);

/*********************** Functions prototyped in sift.h **********************/

/**
Finds SIFT features in an image using default parameter values.  All
detected features are stored in the array pointed to by \a feat.

@param img the image in which to detect features
@param feat a pointer to an array in which to store detected features

@return Returns the number of features stored in \a feat or -1 on failure
@see _sift_features()
*/
int sift_features(IplImage* img, struct feature** feat)
{
return _sift_features(img, feat, SIFT_INTVLS, SIFT_SIGMA, SIFT_CONTR_THR,
SIFT_CURV_THR, SIFT_IMG_DBL, SIFT_DESCR_WIDTH,
SIFT_DESCR_HIST_BINS);
}

/**
Finds SIFT features in an image using user-specified parameter values.  All
detected features are stored in the array pointed to by \a feat.

@param img the image in which to detect features
@param fea a pointer to an array in which to store detected features
@param intvls the number of intervals sampled per octave of scale space
@param sigma the amount of Gaussian smoothing applied to each image level
before building the scale space representation for an octave
@param cont_thr a threshold on the value of the scale space function
\f$\left|D(\hat{x})\right|\f$, where \f$\hat{x}\f$ is a vector specifying
feature location and scale, used to reject unstable features;  assumes
pixel values in the range [0, 1]
@param curv_thr threshold on a feature's ratio of principle curvatures
used to reject features that are too edge-like
@param img_dbl should be 1 if image doubling prior to scale space
construction is desired or 0 if not
@param descr_width the width, \f$n\f$, of the \f$n \times n\f$ array of
orientation histograms used to compute a feature's descriptor
@param descr_hist_bins the number of orientations in each of the
histograms in the array used to compute a feature's descriptor

@return Returns the number of keypoints stored in \a feat or -1 on failure
@see sift_keypoints()
*/
int _sift_features(IplImage* img, struct feature** feat, int intvls,
double sigma, double contr_thr, int curv_thr,
int img_dbl, int descr_width, int descr_hist_bins)
{
IplImage* init_img;
IplImage*** gauss_pyr, *** dog_pyr;
CvMemStorage* storage;
CvSeq* features;
int octvs, i, n = 0;

/* check arguments */
if (!img)
fatal_error("NULL pointer error, %s, line %d", __FILE__, __LINE__);

if (!feat)
fatal_error("NULL pointer error, %s, line %d", __FILE__, __LINE__);

/* build scale space pyramid; smallest dimension of top level is ~4 pixels */
init_img = create_init_img(img, img_dbl, sigma);
octvs = log(MIN(init_img->width, init_img->height)) / log(2) - 2;
gauss_pyr = build_gauss_pyr(init_img, octvs, intvls, sigma);
dog_pyr = build_dog_pyr(gauss_pyr, octvs, intvls);

storage = cvCreateMemStorage(0);
features = scale_space_extrema(dog_pyr, octvs, intvls, contr_thr,
curv_thr, storage);
calc_feature_scales(features, sigma, intvls);
if (img_dbl)
adjust_for_img_dbl(features);
calc_feature_oris(features, gauss_pyr);
compute_descriptors(features, gauss_pyr, descr_width, descr_hist_bins);

/* sort features by decreasing scale and move from CvSeq to array */
cvSeqSort(features, (CvCmpFunc)feature_cmp, NULL);
n = features->total;
*feat = calloc(n, sizeof(struct feature));
*feat = cvCvtSeqToArray(features, *feat, CV_WHOLE_SEQ);
for (i = 0; i < n; i++)
{
free((*feat)[i].feature_data);
(*feat)[i].feature_data = NULL;
}

cvReleaseMemStorage(&storage);
cvReleaseImage(&init_img);
release_pyr(&gauss_pyr, octvs, intvls + 3);
release_pyr(&dog_pyr, octvs, intvls + 2);
return n;
}

/************************ Functions prototyped here **************************/

/*
Converts an image to 8-bit grayscale and Gaussian-smooths it.  The image is
optionally doubled in size prior to smoothing.

@param img input image
@param img_dbl if true, image is doubled in size prior to smoothing
@param sigma total std of Gaussian smoothing
*/
static IplImage* create_init_img(IplImage* img, int img_dbl, double sigma)
{
IplImage* gray, *dbl;
float sig_diff;

gray = convert_to_gray32(img);
if (img_dbl)
{
sig_diff = sqrt(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4);
dbl = cvCreateImage(cvSize(img->width * 2, img->height * 2),
IPL_DEPTH_32F, 1);
cvResize(gray, dbl, CV_INTER_CUBIC);
cvSmooth(dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff);
cvReleaseImage(&gray);
return dbl;
}
else
{
sig_diff = sqrt(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA);
cvSmooth(gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff);
return gray;
}
}

/*
Converts an image to 32-bit grayscale

@param img a 3-channel 8-bit color (BGR) or 8-bit gray image

@return Returns a 32-bit grayscale image
*/
static IplImage* convert_to_gray32(IplImage* img)
{
IplImage* gray8, *gray32;

gray32 = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, 1);
if (img->nChannels == 1)
gray8 = cvClone(img);
else
{
gray8 = cvCreateImage(cvGetSize(img), IPL_DEPTH_8U, 1);
cvCvtColor(img, gray8, CV_BGR2GRAY);
}
cvConvertScale(gray8, gray32, 1.0 / 255.0, 0);

cvReleaseImage(&gray8);
return gray32;
}

/*
Builds Gaussian scale space pyramid from an image

@param base base image of the pyramid
@param octvs number of octaves of scale space
@param intvls number of intervals per octave
@param sigma amount of Gaussian smoothing per octave

@return Returns a Gaussian scale space pyramid as an octvs x (intvls + 3) array
*/
static IplImage*** build_gauss_pyr(IplImage* base, int octvs,
int intvls, double sigma)
{
IplImage*** gauss_pyr;
double* sig = calloc(intvls + 3, sizeof(double));
double sig_total, sig_prev, k;
int i, o;

gauss_pyr = calloc(octvs, sizeof(IplImage**));
for (i = 0; i < octvs; i++)
gauss_pyr[i] = calloc(intvls + 3, sizeof(IplImage*));

/*
precompute Gaussian sigmas using the following formula:

\sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2
*/
sig[0] = sigma;
k = pow(2.0, 1.0 / intvls);
for (i = 1; i < intvls + 3; i++)
{
sig_prev = pow(k, i - 1) * sigma;
sig_total = sig_prev * k;
sig[i] = sqrt(sig_total * sig_total - sig_prev * sig_prev);
}

for (o = 0; o < octvs; o++)
for (i = 0; i < intvls + 3; i++)
{
if (o == 0 && i == 0)
gauss_pyr[o][i] = cvCloneImage(base);

/* base of new octvave is halved image from end of previous octave */
else if (i == 0)
gauss_pyr[o][i] = downsample(gauss_pyr[o - 1][intvls]);

/* blur the current octave's last image to create the next one */
else
{
gauss_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i - 1]),
IPL_DEPTH_32F, 1);
cvSmooth(gauss_pyr[o][i - 1], gauss_pyr[o][i],
CV_GAUSSIAN, 0, 0, sig[i], sig[i]);
}
}

free(sig);
return gauss_pyr;
}

/*
Downsamples an image to a quarter of its size (half in each dimension)
using nearest-neighbor interpolation

@param img an image

@return Returns an image whose dimensions are half those of img
*/
static IplImage* downsample(IplImage* img)
{
IplImage* smaller = cvCreateImage(cvSize(img->width / 2, img->height / 2),
img->depth, img->nChannels);
cvResize(img, smaller, CV_INTER_NN);

return smaller;
}

/*
Builds a difference of Gaussians scale space pyramid by subtracting adjacent
intervals of a Gaussian pyramid

@param gauss_pyr Gaussian scale-space pyramid
@param octvs number of octaves of scale space
@param intvls number of intervals per octave

@return Returns a difference of Gaussians scale space pyramid as an
octvs x (intvls + 2) array
*/
static IplImage*** build_dog_pyr(IplImage*** gauss_pyr, int octvs, int intvls)
{
IplImage*** dog_pyr;
int i, o;

dog_pyr = calloc(octvs, sizeof(IplImage**));
for (i = 0; i < octvs; i++)
dog_pyr[i] = calloc(intvls + 2, sizeof(IplImage*));

for (o = 0; o < octvs; o++)
for (i = 0; i < intvls + 2; i++)
{
dog_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i]),
IPL_DEPTH_32F, 1);
cvSub(gauss_pyr[o][i + 1], gauss_pyr[o][i], dog_pyr[o][i], NULL);
}

return dog_pyr;
}

/*
Detects features at extrema in DoG scale space.  Bad features are discarded
based on contrast and ratio of principal curvatures.

@param dog_pyr DoG scale space pyramid
@param octvs octaves of scale space represented by dog_pyr
@param intvls intervals per octave
@param contr_thr low threshold on feature contrast
@param curv_thr high threshold on feature ratio of principal curvatures
@param storage memory storage in which to store detected features

@return Returns an array of detected features whose scales, orientations,
and descriptors are yet to be determined.
*/
static CvSeq* scale_space_extrema(IplImage*** dog_pyr, int octvs, int intvls,
double contr_thr, int curv_thr,
CvMemStorage* storage)
{
CvSeq* features;
double prelim_contr_thr = 0.5 * contr_thr / intvls;
struct feature* feat;
struct detection_data* ddata;
int o, i, r, c;

features = cvCreateSeq(0, sizeof(CvSeq), sizeof(struct feature), storage);
for (o = 0; o < octvs; o++)
for (i = 1; i <= intvls; i++)
for (r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height - SIFT_IMG_BORDER; r++)
for (c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width - SIFT_IMG_BORDER; c++)
/* perform preliminary check on contrast */
if (ABS(pixval32f(dog_pyr[o][i], r, c)) > prelim_contr_thr)
if (is_extremum(dog_pyr, o, i, r, c))
{
feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr);
if (feat)
{
ddata = feat_detection_data(feat);
if (!is_too_edge_like(dog_pyr[ddata->octv][ddata->intvl],
ddata->r, ddata->c, curv_thr))
{
cvSeqPush(features, feat);
}
else
free(ddata);
free(feat);
}
}

return features;
}

/*
Determines whether a pixel is a scale-space extremum by comparing it to it's
3x3x3 pixel neighborhood.

@param dog_pyr DoG scale space pyramid
@param octv pixel's scale space octave
@param intvl pixel's within-octave interval
@param r pixel's image row
@param c pixel's image col

@return Returns 1 if the specified pixel is an extremum (max or min) among
it's 3x3x3 pixel neighborhood.
*/
static int is_extremum(IplImage*** dog_pyr, int octv, int intvl, int r, int c)
{
float val = pixval32f(dog_pyr[octv][intvl], r, c);
int i, j, k;

/* check for maximum */
if (val > 0)
{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val < pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))
return 0;
}

/* check for minimum */
else
{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val > pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))
return 0;
}

return 1;
}

/*
Interpolates a scale-space extremum's location and scale to subpixel
accuracy to form an image feature.  Rejects features with low contrast.
Based on Section 4 of Lowe's paper.

@param dog_pyr DoG scale space pyramid
@param octv feature's octave of scale space
@param intvl feature's within-octave interval
@param r feature's image row
@param c feature's image column
@param intvls total intervals per octave
@param contr_thr threshold on feature contrast

@return Returns the feature resulting from interpolation of the given
parameters or NULL if the given location could not be interpolated or
if contrast at the interpolated loation was too low.  If a feature is
returned, its scale, orientation, and descriptor are yet to be determined.
*/
static struct feature* interp_extremum(IplImage*** dog_pyr, int octv, int intvl,
int r, int c, int intvls, double contr_thr)
{
struct feature* feat;
struct detection_data* ddata;
double xi, xr, xc, contr;
int i = 0;

while (i < SIFT_MAX_INTERP_STEPS)
{
interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);
if (ABS(xi) < 0.5  &&  ABS(xr) < 0.5  &&  ABS(xc) < 0.5)
break;

c += cvRound(xc);
r += cvRound(xr);
intvl += cvRound(xi);

if (intvl < 1 ||
intvl > intvls ||
c < SIFT_IMG_BORDER ||
r < SIFT_IMG_BORDER ||
c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||
r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER)
{
return NULL;
}

i++;
}

/* ensure convergence of interpolation */
if (i >= SIFT_MAX_INTERP_STEPS)
return NULL;

contr = interp_contr(dog_pyr, octv, intvl, r, c, xi, xr, xc);
if (ABS(contr) < contr_thr / intvls)
return NULL;

feat = new_feature();
ddata = feat_detection_data(feat);
feat->img_pt.x = feat->x = (c + xc) * pow(2.0, octv);
feat->img_pt.y = feat->y = (r + xr) * pow(2.0, octv);
ddata->r = r;
ddata->c = c;
ddata->octv = octv;
ddata->intvl = intvl;
ddata->subintvl = xi;

return feat;
}

/*
Performs one step of extremum interpolation.  Based on Eqn. (3) in Lowe's
paper.

@param dog_pyr difference of Gaussians scale space pyramid
@param octv octave of scale space
@param intvl interval being interpolated
@param r row being interpolated
@param c column being interpolated
@param xi output as interpolated subpixel increment to interval
@param xr output as interpolated subpixel increment to row
@param xc output as interpolated subpixel increment to col
*/

static void interp_step(IplImage*** dog_pyr, 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(dog_pyr, octv, intvl, r, c);
H = hessian_3D(dog_pyr, 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.

@param dog_pyr DoG scale space pyramid
@param octv pixel's octave in dog_pyr
@param intvl pixel's interval in octv
@param r pixel's image row
@param c pixel's image col

@return Returns the vector of partial derivatives for pixel I
{ dI/dx, dI/dy, dI/ds }^T as a CvMat*
*/
static CvMat* deriv_3D(IplImage*** dog_pyr, int octv, int intvl, int r, int c)
{
CvMat* dI;
double dx, dy, ds;

dx = (pixval32f(dog_pyr[octv][intvl], r, c + 1) -
pixval32f(dog_pyr[octv][intvl], r, c - 1)) / 2.0;
dy = (pixval32f(dog_pyr[octv][intvl], r + 1, c) -
pixval32f(dog_pyr[octv][intvl], r - 1, c)) / 2.0;
ds = (pixval32f(dog_pyr[octv][intvl + 1], r, c) -
pixval32f(dog_pyr[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.

@param dog_pyr DoG scale space pyramid
@param octv pixel's octave in dog_pyr
@param intvl pixel's interval in octv
@param r pixel's image row
@param c pixel's image col

@return Returns the Hessian matrix (below) for pixel I as a CvMat*

/ Ixx  Ixy  Ixs \ <BR>
| Ixy  Iyy  Iys | <BR>
\ Ixs  Iys  Iss /
*/
static CvMat* hessian_3D(IplImage*** dog_pyr, int octv, int intvl, int r, int c)
{
CvMat* H;
double v, dxx, dyy, dss, dxy, dxs, dys;

v = pixval32f(dog_pyr[octv][intvl], r, c);
dxx = (pixval32f(dog_pyr[octv][intvl], r, c + 1) +
pixval32f(dog_pyr[octv][intvl], r, c - 1) - 2 * v);
dyy = (pixval32f(dog_pyr[octv][intvl], r + 1, c) +
pixval32f(dog_pyr[octv][intvl], r - 1, c) - 2 * v);
dss = (pixval32f(dog_pyr[octv][intvl + 1], r, c) +
pixval32f(dog_pyr[octv][intvl - 1], r, c) - 2 * v);
dxy = (pixval32f(dog_pyr[octv][intvl], r + 1, c + 1) -
pixval32f(dog_pyr[octv][intvl], r + 1, c - 1) -
pixval32f(dog_pyr[octv][intvl], r - 1, c + 1) +
pixval32f(dog_pyr[octv][intvl], r - 1, c - 1)) / 4.0;
dxs = (pixval32f(dog_pyr[octv][intvl + 1], r, c + 1) -
pixval32f(dog_pyr[octv][intvl + 1], r, c - 1) -
pixval32f(dog_pyr[octv][intvl - 1], r, c + 1) +
pixval32f(dog_pyr[octv][intvl - 1], r, c - 1)) / 4.0;
dys = (pixval32f(dog_pyr[octv][intvl + 1], r + 1, c) -
pixval32f(dog_pyr[octv][intvl + 1], r - 1, c) -
pixval32f(dog_pyr[octv][intvl - 1], r + 1, c) +
pixval32f(dog_pyr[octv][intvl - 1], r - 1, 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;
}

/*
Calculates interpolated pixel contrast.  Based on Eqn. (3) in Lowe's paper.

@param dog_pyr difference of Gaussians scale space pyramid
@param octv octave of scale space
@param intvl within-octave interval
@param r pixel row
@param c pixel column
@param xi interpolated subpixel increment to interval
@param xr interpolated subpixel increment to row
@param xc interpolated subpixel increment to col

@param Returns interpolated contrast.
*/
static double interp_contr(IplImage*** dog_pyr, int octv, int intvl, int r,
int c, double xi, double xr, double xc)
{
CvMat* dD, X, T;
double t[1], x[3] = { xc, xr, xi };

cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP);
cvInitMatHeader(&T, 1, 1, CV_64FC1, t, CV_AUTOSTEP);
dD = deriv_3D(dog_pyr, octv, intvl, r, c);
cvGEMM(dD, &X, 1, NULL, 0, &T, CV_GEMM_A_T);
cvReleaseMat(&dD);

return pixval32f(dog_pyr[octv][intvl], r, c) + t[0] * 0.5;
}

/*
Allocates and initializes a new feature

@return Returns a pointer to the new feature
*/
static struct feature* new_feature(void)
{
struct feature* feat;
struct detection_data* ddata;

feat = malloc(sizeof(struct feature));
memset(feat, 0, sizeof(struct feature));
ddata = malloc(sizeof(struct detection_data));
memset(ddata, 0, sizeof(struct detection_data));
feat->feature_data = ddata;
feat->type = FEATURE_LOWE;

return feat;
}

/*
Determines whether a feature is too edge like to be stable by computing the
ratio of principal curvatures at that feature.  Based on Section 4.1 of
Lowe's paper.

@param dog_img image from the DoG pyramid in which feature was detected
@param r feature row
@param c feature col
@param curv_thr high threshold on ratio of principal curvatures

@return Returns 0 if the feature at (r,c) in dog_img is sufficiently
corner-like or 1 otherwise.
*/
static int is_too_edge_like(IplImage* dog_img, int r, int c, int curv_thr)
{
double d, dxx, dyy, dxy, tr, det;

/* principal curvatures are computed using the trace and det of Hessian */
d = pixval32f(dog_img, r, c);
dxx = pixval32f(dog_img, r, c + 1) + pixval32f(dog_img, r, c - 1) - 2 * d;
dyy = pixval32f(dog_img, r + 1, c) + pixval32f(dog_img, r - 1, c) - 2 * d;
dxy = (pixval32f(dog_img, r + 1, c + 1) - pixval32f(dog_img, r + 1, c - 1) -
pixval32f(dog_img, r - 1, c + 1) + pixval32f(dog_img, r - 1, c - 1)) / 4.0;
tr = dxx + dyy;
det = dxx * dyy - dxy * dxy;

/* negative determinant -> curvatures have different signs; reject feature */
if (det <= 0)
return 1;

if (tr * tr / det < (curv_thr + 1.0)*(curv_thr + 1.0) / curv_thr)
return 0;
return 1;
}

/*
Calculates characteristic scale for each feature in an array.

@param features array of features
@param sigma amount of Gaussian smoothing per octave of scale space
@param intvls intervals per octave of scale space
*/
static void calc_feature_scales(CvSeq* features, double sigma, int intvls)
{
struct feature* feat;
struct detection_data* ddata;
double intvl;
int i, n;

n = features->total;
for (i = 0; i < n; i++)
{
feat = CV_GET_SEQ_ELEM(struct feature, features, i);
ddata = feat_detection_data(feat);
intvl = ddata->intvl + ddata->subintvl;
feat->scl = sigma * pow(2.0, ddata->octv + intvl / intvls);
ddata->scl_octv = sigma * pow(2.0, intvl / intvls);
}
}

/*
Halves feature coordinates and scale in case the input image was doubled
prior to scale space construction.

@param features array of features
*/
static void adjust_for_img_dbl(CvSeq* features)
{
struct feature* feat;
int i, n;

n = features->total;
for (i = 0; i < n; i++)
{
feat = CV_GET_SEQ_ELEM(struct feature, features, i);
feat->x /= 2.0;
feat->y /= 2.0;
feat->scl /= 2.0;
feat->img_pt.x /= 2.0;
feat->img_pt.y /= 2.0;
}
}

/*
Computes a canonical orientation for each image feature in an array.  Based
on Section 5 of Lowe's paper.  This function adds features to the array when
there is more than one dominant orientation at a given feature location.

@param features an array of image features
@param gauss_pyr Gaussian scale space pyramid
*/
static void calc_feature_oris(CvSeq* features, IplImage*** gauss_pyr)
{
struct feature* feat;
struct detection_data* ddata;
double* hist;
double omax;
int i, j, n = features->total;

for (i = 0; i < n; i++)
{
feat = malloc(sizeof(struct feature));
cvSeqPopFront(features, feat);
ddata = feat_detection_data(feat);
hist = ori_hist(gauss_pyr[ddata->octv][ddata->intvl],
ddata->r, ddata->c, SIFT_ORI_HIST_BINS,
cvRound(SIFT_ORI_RADIUS * ddata->scl_octv),
SIFT_ORI_SIG_FCTR * ddata->scl_octv);
for (j = 0; j < SIFT_ORI_SMOOTH_PASSES; j++)
smooth_ori_hist(hist, SIFT_ORI_HIST_BINS);
omax = dominant_ori(hist, SIFT_ORI_HIST_BINS);
add_good_ori_features(features, hist, SIFT_ORI_HIST_BINS,
omax * SIFT_ORI_PEAK_RATIO, feat);
free(ddata);
free(feat);
free(hist);
}
}

/*
Computes a gradient orientation histogram at a specified pixel.

@param img image
@param r pixel row
@param c pixel col
@param n number of histogram bins
@param rad radius of region over which histogram is computed
@param sigma std for Gaussian weighting of histogram entries

@return Returns an n-element array containing an orientation histogram
representing orientations between 0 and 2 PI.
*/
static double* ori_hist(IplImage* img, int r, int c, int n, int rad, double sigma)
{
double* hist;
double mag, ori, w, exp_denom, PI2 = CV_PI * 2.0;
int bin, i, j;

hist = calloc(n, sizeof(double));
exp_denom = 2.0 * sigma * sigma;
for (i = -rad; i <= rad; i++)
for (j = -rad; j <= rad; j++)
if (calc_grad_mag_ori(img, r + i, c + j, &mag, &ori))
{
w = exp(-(i*i + j*j) / exp_denom);
bin = cvRound(n * (ori + CV_PI) / PI2);
bin = (bin < n) ? bin : 0;
hist[bin] += w * mag;
}

return hist;
}

/*
Calculates the gradient magnitude and orientation at a given pixel.

@param img image
@param r pixel row
@param c pixel col
@param mag output as gradient magnitude at pixel (r,c)
@param ori output as gradient orientation at pixel (r,c)

@return Returns 1 if the specified pixel is a valid one and sets mag and
ori accordingly; otherwise returns 0
*/
static int calc_grad_mag_ori(IplImage* img, int r, int c, double* mag, double* ori)
{
double dx, dy;

if (r > 0 && r < img->height - 1 && c > 0 && c < img->width - 1)
{
dx = pixval32f(img, r, c + 1) - pixval32f(img, r, c - 1);
dy = pixval32f(img, r - 1, c) - pixval32f(img, r + 1, c);
*mag = sqrt(dx*dx + dy*dy);
*ori = atan2(dy, dx);
return 1;
}

else
return 0;
}

/*
Gaussian smooths an orientation histogram.

@param hist an orientation histogram
@param n number of bins
*/
static void smooth_ori_hist(double* hist, int n)
{
double prev, tmp, h0 = hist[0];
int i;

prev = hist[n - 1];
for (i = 0; i < n; i++)
{
tmp = hist[i];
hist[i] = 0.25 * prev + 0.5 * hist[i] +
0.25 * ((i + 1 == n) ? h0 : hist[i + 1]);
prev = tmp;
}
}

/*
Finds the magnitude of the dominant orientation in a histogram

@param hist an orientation histogram
@param n number of bins

@return Returns the value of the largest bin in hist
*/
static double dominant_ori(double* hist, int n)
{
double omax;
int maxbin, i;

omax = hist[0];
maxbin = 0;
for (i = 1; i < n; i++)
if (hist[i] > omax)
{
omax = hist[i];
maxbin = i;
}
return omax;
}

/*
Interpolates a histogram peak from left, center, and right values
*/
#define interp_hist_peak( l, c, r ) ( 0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)) )

/*
Adds features to an array for every orientation in a histogram greater than
a specified threshold.

@param features new features are added to the end of this array
@param hist orientation histogram
@param n number of bins in hist
@param mag_thr new features are added for entries in hist greater than this
@param feat new features are clones of this with different orientations
*/
static void add_good_ori_features(CvSeq* features, double* hist, int n,
double mag_thr, struct feature* feat)
{
struct feature* new_feat;
double bin, PI2 = CV_PI * 2.0;
int l, r, i;

for (i = 0; i < n; i++)
{
l = (i == 0) ? n - 1 : i - 1;
r = (i + 1) % n;

if (hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr)
{
bin = i + interp_hist_peak(hist[l], hist[i], hist[r]);
bin = (bin < 0) ? n + bin : (bin >= n) ? bin - n : bin;
new_feat = clone_feature(feat);
new_feat->ori = ((PI2 * bin) / n) - CV_PI;
cvSeqPush(features, new_feat);
free(new_feat);
}
}
}

/*
Makes a deep copy of a feature

@param feat feature to be cloned

@return Returns a deep copy of feat
*/
static struct feature* clone_feature(struct feature* feat)
{
struct feature* new_feat;
struct detection_data* ddata;

new_feat = new_feature();
ddata = feat_detection_data(new_feat);
memcpy(new_feat, feat, sizeof(struct feature));
memcpy(ddata, feat_detection_data(feat), sizeof(struct detection_data));
new_feat->feature_data = ddata;

return new_feat;
}

/*
Computes feature descriptors for features in an array.  Based on Section 6
of Lowe's paper.

@param features array of features
@param gauss_pyr Gaussian scale space pyramid
@param d width of 2D array of orientation histograms
@param n number of bins per orientation histogram
*/
static void compute_descriptors(CvSeq* features, IplImage*** gauss_pyr, int d, int n)
{
struct feature* feat;
struct detection_data* ddata;
double*** hist;
int i, k = features->total;

for (i = 0; i < k; i++)
{
feat = CV_GET_SEQ_ELEM(struct feature, features, i);
ddata = feat_detection_data(feat);
hist = descr_hist(gauss_pyr[ddata->octv][ddata->intvl], ddata->r,
ddata->c, feat->ori, ddata->scl_octv, d, n);
hist_to_descr(hist, d, n, feat);
release_descr_hist(&hist, d);
}
}

/*
Computes the 2D array of orientation histograms that form the feature
descriptor.  Based on Section 6.1 of Lowe's paper.

@param img image used in descriptor computation
@param r row coord of center of orientation histogram array
@param c column coord of center of orientation histogram array
@param ori canonical orientation of feature whose descr is being computed
@param scl scale relative to img of feature whose descr is being computed
@param d width of 2d array of orientation histograms
@param n bins per orientation histogram

@return Returns a d x d array of n-bin orientation histograms.
*/
static double*** descr_hist(IplImage* img, int r, int c, double ori,
double scl, int d, int n)
{
double*** hist;
double cos_t, sin_t, hist_width, exp_denom, r_rot, c_rot, grad_mag,
grad_ori, w, rbin, cbin, obin, bins_per_rad, PI2 = 2.0 * CV_PI;
int radius, i, j;

hist = calloc(d, sizeof(double**));
for (i = 0; i < d; i++)
{
hist[i] = calloc(d, sizeof(double*));
for (j = 0; j < d; j++)
hist[i][j] = calloc(n, sizeof(double));
}

cos_t = cos(ori);
sin_t = sin(ori);
bins_per_rad = n / PI2;
exp_denom = d * d * 0.5;
hist_width = SIFT_DESCR_SCL_FCTR * scl;
radius = hist_width * sqrt(2) * (d + 1.0) * 0.5 + 0.5;
for (i = -radius; i <= radius; i++)
for (j = -radius; j <= radius; j++)
{
/*
Calculate sample's histogram array coords rotated relative to ori.
Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.
r_rot = 1.5) have full weight placed in row 1 after interpolation.
*/
c_rot = (j * cos_t - i * sin_t) / hist_width;
r_rot = (j * sin_t + i * cos_t) / hist_width;
rbin = r_rot + d / 2 - 0.5;
cbin = c_rot + d / 2 - 0.5;

if (rbin > -1.0  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d)
if (calc_grad_mag_ori(img, r + i, c + j, &grad_mag, &grad_ori))
{
grad_ori -= ori;
while (grad_ori < 0.0)
grad_ori += PI2;
while (grad_ori >= PI2)
grad_ori -= PI2;

obin = grad_ori * bins_per_rad;
w = exp(-(c_rot * c_rot + r_rot * r_rot) / exp_denom);
interp_hist_entry(hist, rbin, cbin, obin, grad_mag * w, d, n);
}
}

return hist;
}

/*
Interpolates an entry into the array of orientation histograms that form
the feature descriptor.

@param hist 2D array of orientation histograms
@param rbin sub-bin row coordinate of entry
@param cbin sub-bin column coordinate of entry
@param obin sub-bin orientation coordinate of entry
@param mag size of entry
@param d width of 2D array of orientation histograms
@param n number of bins per orientation histogram
*/
static void interp_hist_entry(double*** hist, double rbin, double cbin,
double obin, double mag, int d, int n)
{
double d_r, d_c, d_o, v_r, v_c, v_o;
double** row, *h;
int r0, c0, o0, rb, cb, ob, r, c, o;

r0 = cvFloor(rbin);
c0 = cvFloor(cbin);
o0 = cvFloor(obin);
d_r = rbin - r0;
d_c = cbin - c0;
d_o = obin - o0;

/*
The entry is distributed into up to 8 bins.  Each entry into a bin
is multiplied by a weight of 1 - d for each dimension, where d is the
distance from the center value of the bin measured in bin units.
*/
for (r = 0; r <= 1; r++)
{
rb = r0 + r;
if (rb >= 0 && rb < d)
{
v_r = mag * ((r == 0) ? 1.0 - d_r : d_r);
row = hist[rb];
for (c = 0; c <= 1; c++)
{
cb = c0 + c;
if (cb >= 0 && cb < d)
{
v_c = v_r * ((c == 0) ? 1.0 - d_c : d_c);
h = row[cb];
for (o = 0; o <= 1; o++)
{
ob = (o0 + o) % n;
v_o = v_c * ((o == 0) ? 1.0 - d_o : d_o);
h[ob] += v_o;
}
}
}
}
}
}

/*
Converts the 2D array of orientation histograms into a feature's descriptor
vector.

@param hist 2D array of orientation histograms
@param d width of hist
@param n bins per histogram
@param feat feature into which to store descriptor
*/
static void hist_to_descr(double*** hist, int d, int n, struct feature* feat)
{
int int_val, i, r, c, o, k = 0;

for (r = 0; r < d; r++)
for (c = 0; c < d; c++)
for (o = 0; o < n; o++)
feat->descr[k++] = hist[r][c][o];

feat->d = k;
normalize_descr(feat);
for (i = 0; i < k; i++)
if (feat->descr[i] > SIFT_DESCR_MAG_THR)
feat->descr[i] = SIFT_DESCR_MAG_THR;
normalize_descr(feat);

/* convert floating-point descriptor to integer valued descriptor */
for (i = 0; i < k; i++)
{
int_val = SIFT_INT_DESCR_FCTR * feat->descr[i];
feat->descr[i] = MIN(255, int_val);
}
}

/*
Normalizes a feature's descriptor vector to unitl length

@param feat feature
*/
static void normalize_descr(struct feature* feat)
{
double cur, len_inv, len_sq = 0.0;
int i, d = feat->d;

for (i = 0; i < d; i++)
{
cur = feat->descr[i];
len_sq += cur*cur;
}
len_inv = 1.0 / sqrt(len_sq);
for (i = 0; i < d; i++)
feat->descr[i] *= len_inv;
}

/*
Compares features for a decreasing-scale ordering.  Intended for use with
CvSeqSort

@param feat1 first feature
@param feat2 second feature
@param param unused

@return Returns 1 if feat1's scale is greater than feat2's, -1 if vice versa,
and 0 if their scales are equal
*/
static int feature_cmp(void* feat1, void* feat2, void* param)
{
struct feature* f1 = (struct feature*) feat1;
struct feature* f2 = (struct feature*) feat2;

if (f1->scl < f2->scl)
return 1;
if (f1->scl > f2->scl)
return -1;
return 0;
}

/*
De-allocates memory held by a descriptor histogram

@param hist pointer to a 2D array of orientation histograms
@param d width of hist
*/
static void release_descr_hist(double**** hist, int d)
{
int i, j;

for (i = 0; i < d; i++)
{
for (j = 0; j < d; j++)
free((*hist)[i][j]);
free((*hist)[i]);
}
free(*hist);
*hist = NULL;
}

/*
De-allocates memory held by a scale space pyramid

@param pyr scale space pyramid
@param octvs number of octaves of scale space
@param n number of images per octave
*/
static void release_pyr(IplImage**** pyr, int octvs, int n)
{
int i, j;
for (i = 0; i < octvs; i++)
{
for (j = 0; j < n; j++)
cvReleaseImage(&(*pyr)[i][j]);
free((*pyr)[i]);
}
free(*pyr);
*pyr = NULL;
}

utils.c

/*
Miscellaneous utility functions.

Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>

@version 1.1.2-20100521
*/

#include "utils.h"
#include <stdarg.h>

#include <cv.h>
#include <cxcore.h>
#include <highgui.h>

#include <errno.h>
#include <string.h>
#include <stdlib.h>

/*************************** Function Definitions ****************************/

/*
Prints an error message and aborts the program.  The error message is
of the form "Error: ...", where the ... is specified by the \a format
argument

@param format an error message format string (as with \c printf(3)).
*/
void fatal_error(char* format, ...)
{
va_list ap;

fprintf(stderr, "Error: ");

va_start(ap, format);
vfprintf(stderr, format, ap);
va_end(ap);
fprintf(stderr, "\n");
abort();
}

/*
Replaces a file's extension, which is assumed to be everything after the
last dot ('.') character.

@param file the name of a file

@param extn a new extension for \a file; should not include a dot (i.e.
\c "jpg", not \c ".jpg") unless the new file extension should contain
two dots.

@return Returns a new string formed as described above.  If \a file does
not have an extension, this function simply adds one.
*/
char* replace_extension(const char* file, const char* extn)
{
char* new_file, *lastdot;

new_file = calloc(strlen(file) + strlen(extn) + 2, sizeof(char));
strcpy(new_file, file);
lastdot = strrchr(new_file, '.');
if (lastdot)
*(lastdot + 1) = '\0';
else
strcat(new_file, ".");
strcat(new_file, extn);

return new_file;
}

/*
A function that removes the path from a filename.  Similar to the Unix
basename command.

@param pathname a (full) path name

@return Returns the basename of \a pathname.
*/
char* basename(const char* pathname)
{
char* base, *last_slash;

last_slash = strrchr(pathname, '/');
if (!last_slash)
{
base = calloc(strlen(pathname) + 1, sizeof(char));
strcpy(base, pathname);
}
else
{
base = calloc(strlen(last_slash++), sizeof(char));
strcpy(base, last_slash);
}

return base;
}

/*
Displays progress in the console with a spinning pinwheel.  Every time this
function is called, the state of the pinwheel is incremented.  The pinwheel
has four states that loop indefinitely: '|', '/', '-', '\'.

@param done if 0, this function simply increments the state of the pinwheel;
otherwise it prints "done"
*/
void progress(int done)
{
char state[4] = { '|', '/', '-', '\\' };
static int cur = -1;

if (cur == -1)
fprintf(stderr, "  ");

if (done)
{
fprintf(stderr, "\b\bdone\n");
cur = -1;
}
else
{
cur = (cur + 1) % 4;
fprintf(stdout, "\b\b%c ", state[cur]);
fflush(stderr);
}
}

/*
Erases a specified number of characters from a stream.

@param stream the stream from which to erase characters
@param n the number of characters to erase
*/
void erase_from_stream(FILE* stream, int n)
{
int j;
for (j = 0; j < n; j++)
fprintf(stream, "\b");
for (j = 0; j < n; j++)
fprintf(stream, " ");
for (j = 0; j < n; j++)
fprintf(stream, "\b");
}

/*
Doubles the size of an array with error checking

@param array pointer to an array whose size is to be doubled
@param n number of elements allocated for \a array
@param size size in bytes of elements in \a array

@return Returns the new number of elements allocated for \a array.  If no
memory is available, returns 0 and frees array.
*/
int array_double(void** array, int n, int size)
{
void* tmp;

tmp = realloc(*array, 2 * n * size);
if (!tmp)
{
fprintf(stderr, "Warning: unable to allocate memory in array_double(),"
" %s line %d\n", __FILE__, __LINE__);
if (*array)
free(*array);
*array = NULL;
return 0;
}
*array = tmp;
return n * 2;
}

/*
Calculates the squared distance between two points.

@param p1 a point
@param p2 another point
*/
double dist_sq_2D(CvPoint2D64f p1, CvPoint2D64f p2)
{
double x_diff = p1.x - p2.x;
double y_diff = p1.y - p2.y;

return x_diff * x_diff + y_diff * y_diff;
}

/*
Draws an x on an image.

@param img an image
@param pt the center point of the x
@param r the x's radius
@param w the x's line weight
@param color the color of the x
*/
void draw_x(IplImage* img, CvPoint pt, int r, int w, CvScalar color)
{
cvLine(img, pt, cvPoint(pt.x + r, pt.y + r), color, w, 8, 0);
cvLine(img, pt, cvPoint(pt.x - r, pt.y + r), color, w, 8, 0);
cvLine(img, pt, cvPoint(pt.x + r, pt.y - r), color, w, 8, 0);
cvLine(img, pt, cvPoint(pt.x - r, pt.y - r), color, w, 8, 0);
}

/*
Combines two images by scacking one on top of the other

@param img1 top image
@param img2 bottom image

@return Returns the image resulting from stacking \a img1 on top if \a img2
*/
extern IplImage* stack_imgs(IplImage* img1, IplImage* img2)
{
IplImage* stacked = cvCreateImage(cvSize(MAX(img1->width, img2->width),
img1->height + img2->height),
IPL_DEPTH_8U, 3);

cvZero(stacked);
cvSetImageROI(stacked, cvRect(0, 0, img1->width, img1->height));
cvAdd(img1, stacked, stacked, NULL);
cvSetImageROI(stacked, cvRect(0, img1->height, img2->width, img2->height));
cvAdd(img2, stacked, stacked, NULL);
cvResetImageROI(stacked);

return stacked;
}

/*
Allows user to view an array of images as a video.  Keyboard controls
are as follows:

<ul>
<li>Space - start and pause playback</li>
<li>Page Down - skip forward 10 frames</li>
<li>Page Up - jump back 10 frames</li>
<li>Right Arrow - skip forward 1 frame</li>
<li>Left Arrow - jump back 1 frame</li>
<li>Backspace - jump back to beginning</li>
<li>Esc - exit playback</li>
<li>Closing the window also exits playback</li>
</ul>

@param imgs an array of images
@param n number of images in \a imgs
@param win_name name of window in which images are displayed
*/
void vid_view(IplImage** imgs, int n, char* win_name)
{
int k, i = 0, playing = 0;

cvNamedWindow(win_name, 1);
cvShowImage(win_name, imgs[i]);
while (!win_closed(win_name))
{
/* if already playing, advance frame and check for pause */
if (playing)
{
i = MIN(i + 1, n - 1);
cvNamedWindow(win_name, 1);
cvShowImage(win_name, imgs[i]);
k = cvWaitKey(33);
if (k == ' ' || i == n - 1)
playing = 0;
}

else
{
k = cvWaitKey(0);
switch (k)
{
/* space */
case ' ':
playing = 1;
break;

/* esc */
case 27:
case 1048603:
cvDestroyWindow(win_name);
break;

/* backspace */
case '\b':
i = 0;
cvNamedWindow(win_name, 1);
cvShowImage(win_name, imgs[i]);
break;

/* left arrow */
case 65288:
case 1113937:
i = MAX(i - 1, 0);
cvNamedWindow(win_name, 1);
cvShowImage(win_name, imgs[i]);
break;

/* right arrow */
case 65363:
case 1113939:
i = MIN(i + 1, n - 1);
cvNamedWindow(win_name, 1);
cvShowImage(win_name, imgs[i]);
break;

/* page up */
case 65365:
case 1113941:
i = MAX(i - 10, 0);
cvNamedWindow(win_name, 1);
cvShowImage(win_name, imgs[i]);
break;

/* page down */
case 65366:
case 1113942:
i = MIN(i + 10, n - 1);
cvNamedWindow(win_name, 1);
cvShowImage(win_name, imgs[i]);
break;
}
}
}
}

/*
Checks if a HighGUI window is still open or not

@param name the name of the window we're checking

@return Returns 1 if the window named \a name has been closed or 0 otherwise
*/
int win_closed(char* win_name)
{
if (!cvGetWindowHandle(win_name))
return 1;
return 0;
}