sift是图像匹配的非常经典的算法，但是很复杂，要想自己拿C或C++实现很麻烦，如果只是使用的话，有国外某高人维护的sift库，前期只要自己能够调用即可，关键是要熟悉大致的流程，对sift库有个了解，具体的工作只要调用其中的函数即可。匹配效果：
sift是图像匹配的非常经典的算法，但是很复杂，要想自己拿C或C++实现很麻烦，如果只是使用的话，有国外某高人维护的sift库，前期只要自己能够调用即可，关键是要熟悉大致的流程，对sift库有个了解，具体的工作只要调用其中的函数即可。
一、sift简介
1、sift算法应用典型场合：
物体识别
机器人定位与导航
图像拼接
三维建模
手势识别
视频跟踪
笔记鉴定
指纹与人脸识别
犯罪现场特征提取
2、sift算法解决的问题
目标的旋转、缩放、平移（RST）
图像仿射/投影变换（视点viewpoint）
光照影响（illumination）
目标遮挡（occlusion）
杂物场景（clutter）
噪声
二、sift基本概念
1、SIFT匹配的定义
SIFT匹配（Scale-invariant feature transform，尺度不变特征转换)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由 David Lowe 在1999年所发表，2004年完善总结。其应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对。
局部影像特征的描述与侦测可以帮助辨识物体，SIFT 特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、些微视角改变的容忍度也相当高。基于这些特性，它们是高度显著而且相对容易撷取，在母数庞大的特征数据库中，很容易辨识物体而且鲜有误认。使用 SIFT特征描述对于部分物体遮蔽的侦测率也相当高，甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下，辨识速度可接近即时运算。SIFT特征的信息量大，适合在海量数据库中快速准确匹配。
2、SIFT特征的主要特点
从理论上说，SIFT是一种相似不变量，即对图像尺度变化和旋转是不变量。然而，由于构造SIFT特征时，在很多细节上进行了特殊处理，使得SIFT对图像的复杂变形和光照变化具有了较强的适应性，同时运算速度比较快，定位精度比较高。如：
在多尺度空间采用DOG算子检测关键点，相比传统的基于LOG算子的检测方法，运算速度大大加快；
关键点的精确定位不仅提高了精度，而且大大提高了关键点的稳定性；
在构造描述子时，以子区域的统计特性，而不是以单个像素作为研究对象，提高了对图像局部变形的适应能力；
对于16*16的关键点邻域和4*4的子区域，在处理梯度幅度时都进行了类似于高斯函数的加权处理，强化了中心区域，淡化了边缘区域的影响，从而提高了算法对几何变形的适应性；
该方法不仅对通用的线性光照模型具有不变性，而且对复杂的光照变化亦具有一定的适应性。关于这部分内容的细节，可参看文献“Distinctive Image Features from Scale-Invariant Keypoints”
SIFT算法的特点：
a.SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程 度的稳定性；
b.独特性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配；
c.多量性，即使少数的几 个物体也可以产生大量的SIFT特征向量；
d.高速性，经优化的SIFT匹配算法甚至可以达到实时的要求；
e.可扩展性，可以很方便的与其他形式 的特征向量进行联合。
三、sift算法的应用
1、sift的main函数
#include "sift.h"
#include "imgfeatures.h"
#include "kdtree.h"
#include "utils.h"
#include "xform.h"
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
#include <stdio.h>
/* the maximum number of keypoint NN candidates to check during BBF search */
#define KDTREE_BBF_MAX_NN_CHKS 200
/* threshold on squared ratio of distances between NN and 2nd NN */
#define NN_SQ_DIST_RATIO_THR 0.49
/******************************** Globals ************************************/
char img1_file[] = "beaver.png";
char img2_file[] = "beaver_xform.png";
/********************************** Main *************************************/
int main( int argc, char** argv )
{
IplImage* img1, * img2, * stacked;
struct feature* feat1, * feat2, * feat;
struct feature** nbrs;
struct kd_node* kd_root;
CvPoint pt1, pt2;
double d0, d1;
int n1, n2, k, i, m = 0;
img1 = cvLoadImage( img1_file, 1 );
if( ! img1 )
fatal_error( "unable to load image from %s", img1_file );
img2 = cvLoadImage( img2_file, 1 );
if( ! img2 )
fatal_error( "unable to load image from %s", img2_file );
stacked = stack_imgs( img1, img2 );
fprintf( stderr, "Finding features in %s...\n", img1_file );
n1 = sift_features( img1, &feat1 );
fprintf( stderr, "Finding features in %s...\n", img2_file );
n2 = sift_features( img2, &feat2 );
kd_root = kdtree_build( feat2, n2 );
for( i = 0; i < n1; i++ )
{
feat = feat1 + i;
k = kdtree_bbf_knn( kd_root, feat, 2, &nbrs, KDTREE_BBF_MAX_NN_CHKS );
if( k == 2 )
{
d0 = descr_dist_sq( feat, nbrs[0] );
d1 = descr_dist_sq( feat, nbrs[1] );
if( d0 < d1 * NN_SQ_DIST_RATIO_THR )
{
pt1 = cvPoint( cvRound( feat->x ), cvRound( feat->y ) );
pt2 = cvPoint( cvRound( nbrs[0]->x ), cvRound( nbrs[0]->y ) );
pt2.y += img1->height;
cvLine( stacked, pt1, pt2, CV_RGB(255,0,255), 1, 8, 0 );
m++;
feat1[i].fwd_match = nbrs[0];
}
}
free( nbrs );
}
fprintf( stderr, "Found %d total matches\n", m );
cvNamedWindow( "Matches", 1 );
cvShowImage( "Matches", stacked );
cvSaveImage("1.bmp",stacked,0);
cvWaitKey( 0 );
cvReleaseImage( &stacked );
cvReleaseImage( &img1 );
cvReleaseImage( &img2 );
kdtree_release( kd_root );
free( feat1 );
free( feat2 );
return 0;
}
Sift（Scale InvariantFeature Transform）是一个很好的图像匹配算法，同时能处理亮度、平移、旋转、尺度的变化，利用特征点来提取特征描述符，最后在特征描述符之间寻找匹配。
sift库分析：
总共包含6 个文件再加上一个main.c文件。
untils.c
minpq.c
imgfeature.c
sift.c
xform.c
kdtree.c
以及其对应的.h文件
最核心的文件是sift.c 和 kdtree.c
1 sift.c 寻找图片中的特征点
主要步骤：
a、构建尺度空间，检测极值点，获得尺度不变性；
b、特征点过滤并进行精确定位，剔除不稳定的特征点；
c、在特征点处提取特征描述符，为特征点分配方向值；
d、生成特征描述子，利用特征描述符寻找匹配点；
以特征点为中心取16*16的邻域作为采样窗口，
将采样点与特征点的相对方向通过高斯加权后归入包含8个bin的方向直方图，
最后获得4*4*8的128维特征描述子。
2 kdtree.c 对两幅图像进行匹配
主要步骤：
当两幅图像的Sift特征向量生成以后，下一步就可以采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。
取图1的某个关键点，通过遍历找到图像2中的距离最近的两个关键点。在这两个关键点中，如果次近距离除以最近距离小于某个阙值，则判定为一对匹配点。
3 imgfeature.c在以上图片中画上标记。连接对应的匹配的点。
四、具体文件分析
1、utils.h
utils是sift算法的一个比较基础的头文件，主要是对图像的基本操纵：
获取某位置的像素点
设置某位置的像素点（8位，32位和64位），
计算两点之间的距离的平方
在图片某一点画一个“X”
将两张图片合成为一个，高是二者之和，宽是二者的较大者。img1 在左上角，img2在右下角。
函数的解释如下：
static int pixval8 (IplImage *img, int r,int c)
从灰度图像中获取某点像素
static void setpix8 (IplImage *img, int r,int c, uchar val)
设置灰度图像中某点的像素
static float pixval32f (IplImage *img, intr, int c)
static void setpix32f (IplImage *img, intr, int c, float val)
static double pixval64f (IplImage *img, intr, int c)
static void setpix64f (IplImage *img, intr, int c, double val)
在32位和64位图像中获取和设置某点的像素值
void fatal_error (char *format,...)
错误处理，用VS2010+opencv2.3报错，直接把内部实现注释掉即可，无影响。
char * replace_extension (const char *file,const char *extn)
获取一个文件全名，将文件名和扩展名连接到一起如 traffic + jpg => traffic.jpg
char * prepend_path (const char *path,const char *file)
获取一个文件的全路径，将路径名添加到文件名之前 C:\\ + traffic.jpg => c:\\traffic.jpg
char * basename (const char *pathname)
文件名中去掉路径c:\\traffic.jpg => traffic.jpg
void progress (int done)
显示程序进展，用|/-\表示
int array_double(void **array, int n, int size)
数组长度加倍
double dist_sq_2D(CvPoint2D64f p1, CvPoint2D64f p2)
计算两点的对角线距离
void draw_x(IplImage *img, CvPoint pt, int r, int w, CvScalar color)
在点pt画个叉，本质就是在那个点画四条线
IplImage * stack_imgs (IplImage *img1,IplImage *img2)
两张图片生成一张图片，高是二者之和，宽是二者的较大者
int win_closed(char *name)
查看某个窗口是否已经关闭
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
###utils.c
#include "utils.h"
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
#include <errno.h>
#include <string.h>
#include <stdlib.h>
void fatal_error(char* format, ...)
{
}
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;
}
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;
}
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);
}
}
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" );
}
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;
}
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;
}
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 );
}
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;
}
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;
}
2、imfeature.h
int import_features (char *filename, inttype, struct feature **feat)
从一个文件中（txt）读入几个数据，作为特征的初始化参数
type指的是特征的类型FEATURE_OXFD或者FEATURE_LOWE
int export_features (char *filename, structfeature *feat, int n)
将几个参数存到一个文件中去，参数指的是feature 这个结构，feature**指的是这个结构的数组的指针 ，传递指针可以避免数据的复制
void draw_features (IplImage *img, structfeature *feat, int n)
将特征绘制到一个图片上
double descr_dist_sq(struct feature *f1, struct feature *f2)
计算两个特征描绘子之间的欧式距离
feature数据结构分析：
struct feature
{
double x; /**< x coord */
double y; /**< y coord */
double a; /**< Oxford-typeaffine region parameter */
double b; /**< Oxford-typeaffine region parameter */
double c; /**< Oxford-type affineregion parameter */
double scl; /**< scale of aLowe-style feature */
double ori; /**< orientation of aLowe-style feature */
int d; /**< descriptorlength */
double descr[FEATURE_MAX_D]; /**< descriptor */
int type; /**< feature type,OXFD or LOWE */
int category; /**< all-purpose featurecategory */
struct feature* fwd_match; /**< matching feature from forwardimage */
struct feature* bck_match; /**< matching feature from backmwardimage */
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 */
};
imfeature.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 */
};
/**
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 );
#endif
imfeature.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 );
}
3、xforms.h
CvMat * ransac_xform (struct feature*features, int n, int mtype, ransac_xform_fn xform_fn, int m, doublep_badxform, ransac_err_fn err_fn, double err_tol, struct feature ***inliers,int *n_in)
运用RANSAC方法进行特征匹配，计算出匹配最佳的图像转化矩阵
CvMat * dlt_homog (CvPoint2D64f *pts,CvPoint2D64f *mpts, int n)
运用线性变换，进行点匹配计算平面单应性
CvMat * lsq_homog (CvPoint2D64f *pts,CvPoint2D64f *mpts, int n)
进行点匹配，计算对角线最短的平面单应性
double homog_xfer_err(CvPoint2D64f pt, CvPoint2D64f mpt, CvMat *H)
对于给定的单应性，计算某点及其匹配的转换误差
CvPoint2D64f persp_xform_pt (CvPoint2D64fpt, CvMat *T)
透视变换
xforms.h
/**@file
Functions for computing transforms from image feature correspondences
Copyright (C) 2006-2010 Rob Hess <hess@eecs.oregonstate.edu>
@version 1.1.2-20100521
*/
#ifndef XFORM_H
#define XFORM_H
#include <cxcore.h>
/********************************** Structures *******************************/
struct feature;
/** holds feature data relevant to ransac */
struct ransac_data
{
void* orig_feat_data;
int sampled;
};
/******************************* Defs and macros *****************************/
/* RANSAC error tolerance in pixels */
#define RANSAC_ERR_TOL 3
/** pessimistic estimate of fraction of inlers for RANSAC */
#define RANSAC_INLIER_FRAC_EST 0.25
/** estimate of the probability that a correspondence supports a bad model */
#define RANSAC_PROB_BAD_SUPP 0.10
/* extracts a feature's RANSAC data */
#define feat_ransac_data( feat ) ( (struct ransac_data*) (feat)->feature_data )
/**
Prototype for transformation functions passed to ransac_xform(). Functions
of this type should compute a transformation matrix given a set of point
correspondences.
@param pts array of points
@param mpts array of corresponding points; each \a pts[\a i], \a i=0..\a n-1,
corresponds to \a mpts[\a i]
@param n number of points in both \a pts and \a mpts
@return Should return a transformation matrix that transforms each point in
\a pts to the corresponding point in \a mpts or NULL on failure.
*/
typedef CvMat* (*ransac_xform_fn)( CvPoint2D64f* pts, CvPoint2D64f* mpts,
int n );
/**
Prototype for error functions passed to ransac_xform(). For a given
point, its correspondence, and a transform, functions of this type should
compute a measure of error between the correspondence and the point after
the point has been transformed by the transform.
@param pt a point
@param mpt \a pt's correspondence
@param T a transform
@return Should return a measure of error between \a mpt and \a pt after
\a pt has been transformed by the transform \a T.
*/
typedef double (*ransac_err_fn)( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* M );
/***************************** Function Prototypes ***************************/
/**
Calculates a best-fit image transform from image feature correspondences
using RANSAC.
For more information refer to:
Fischler, M. A. and Bolles, R. C. Random sample consensus: a paradigm for
model fitting with applications to image analysis and automated cartography.
<EM>Communications of the ACM, 24</EM>, 6 (1981), pp. 381--395.
@param features an array of features; only features with a non-NULL match
of type \a mtype are used in homography computation
@param n number of features in \a feat
@param mtype determines which of each feature's match fields to use
for transform computation; should be one of FEATURE_FWD_MATCH,
FEATURE_BCK_MATCH, or FEATURE_MDL_MATCH; if this is FEATURE_MDL_MATCH,
correspondences are assumed to be between a feature's img_pt field
and its match's mdl_pt field, otherwise correspondences are assumed to
be between the the feature's img_pt field and its match's img_pt field
@param xform_fn pointer to the function used to compute the desired
transformation from feature correspondences
@param m minimum number of correspondences necessary to instantiate the
transform computed by \a xform_fn
@param p_badxform desired probability that the final transformation
returned by RANSAC is corrupted by outliers (i.e. the probability that
no samples of all inliers were drawn)
@param err_fn pointer to the function used to compute a measure of error
between putative correspondences for a given transform
@param err_tol correspondences within this distance of each other are
considered as inliers for a given transform
@param inliers if not NULL, output as an array of pointers to the final
set of inliers; memory for this array is allocated by this function and
must be freed by the caller using free(*inliers)
@param n_in if not NULL, output as the final number of inliers
@return Returns a transformation matrix computed using RANSAC or NULL
on error or if an acceptable transform could not be computed.
*/
extern CvMat* ransac_xform( struct feature* features, int n, int mtype,
ransac_xform_fn xform_fn, int m,
double p_badxform, ransac_err_fn err_fn,
double err_tol, struct feature*** inliers,
int* n_in );
/**
Calculates a planar homography from point correspondeces using the direct
linear transform. Intended for use as a ransac_xform_fn.
@param pts array of points
@param mpts array of corresponding points; each \a pts[\a i], \a i=0..\a
n-1, corresponds to \a mpts[\a i]
@param n number of points in both \a pts and \a mpts; must be at least 4
@return Returns the \f$3 \times 3\f$ planar homography matrix that
transforms points in \a pts to their corresponding points in \a mpts
or NULL if fewer than 4 correspondences were provided
*/
extern CvMat* dlt_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n );
/**
Calculates a least-squares planar homography from point correspondeces.
Intended for use as a ransac_xform_fn.
@param pts array of points
@param mpts array of corresponding points; each \a pts[\a i], \a i=0..\a n-1,
corresponds to \a mpts[\a i]
@param n number of points in both \a pts and \a mpts; must be at least 4
@return Returns the \f$3 \times 3\f$ least-squares planar homography
matrix that transforms points in \a pts to their corresponding points
in \a mpts or NULL if fewer than 4 correspondences were provided
*/
extern CvMat* lsq_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n );
/**
Calculates the transfer error between a point and its correspondence for
a given homography, i.e. for a point \f$x\f$, it's correspondence \f$x'\f$,
and homography \f$H\f$, computes \f$d(x', Hx)^2\f$. Intended for use as a
ransac_err_fn.
@param pt a point
@param mpt \a pt's correspondence
@param H a homography matrix
@return Returns the transfer error between \a pt and \a mpt given \a H
*/
extern double homog_xfer_err( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* H );
/**
Performs a perspective transformation on a single point. That is, for a
point \f$(x, y)\f$ and a \f$3 \times 3\f$ matrix \f$T\f$ this function
returns the point \f$(u, v)\f$, where<BR>
\f$[x' \ y' \ w']^T = T \times [x \ y \ 1]^T\f$,<BR>
and<BR>
\f$(u, v) = (x'/w', y'/w')\f$.
Note that affine transforms are a subset of perspective transforms.
@param pt a 2D point
@param T a perspective transformation matrix
@return Returns the point \f$(u, v)\f$ as above.
*/
extern CvPoint2D64f persp_xform_pt( CvPoint2D64f pt, CvMat* T );
#endif
xform.c
/*
This file contains definitions for functions to compute transforms from
image feature correspondences
Copyright (C) 2006-2010 Rob Hess <hess@eecs.oregonstate.edu>
@version 1.1.2-20100521
*/
#include "xform.h"
#include "imgfeatures.h"
#include "utils.h"
#include <cxcore.h>
#include <stdlib.h>
#include <time.h>
/************************* Local Function Prototypes *************************/
static __inline struct feature* get_match( struct feature*, int );
static int get_matched_features( struct feature*, int, int, struct feature*** );
static int calc_min_inliers( int, int, double, double );
static __inline double log_factorial( int );
static struct feature** draw_ransac_sample( struct feature**, int, int );
static void extract_corresp_pts( struct feature**, int, int, CvPoint2D64f**,
CvPoint2D64f** );
static int find_consensus( struct feature**, int, int, CvMat*, ransac_err_fn,
double, struct feature*** );
static __inline void release_mem( CvPoint2D64f*, CvPoint2D64f*,
static struct feature** );
/********************** Functions prototyped in xform.h **********************/
/*
Calculates a best-fit image transform from image feature correspondences
using RANSAC.
For more information refer to:
Fischler, M. A. and Bolles, R. C. Random sample consensus: a paradigm for
model fitting with applications to image analysis and automated cartography.
<EM>Communications of the ACM, 24</EM>, 6 (1981), pp. 381--395.
@param features an array of features; only features with a non-NULL match
of type mtype are used in homography computation
@param n number of features in feat
@param mtype determines which of each feature's match fields to use
for model computation; should be one of FEATURE_FWD_MATCH,
FEATURE_BCK_MATCH, or FEATURE_MDL_MATCH; if this is FEATURE_MDL_MATCH,
correspondences are assumed to be between a feature's img_pt field
and its match's mdl_pt field, otherwise correspondences are assumed to
be between the the feature's img_pt field and its match's img_pt field
@param xform_fn pointer to the function used to compute the desired
transformation from feature correspondences
@param m minimum number of correspondences necessary to instantiate the
model computed by xform_fn
@param p_badxform desired probability that the final transformation
returned by RANSAC is corrupted by outliers (i.e. the probability that
no samples of all inliers were drawn)
@param err_fn pointer to the function used to compute a measure of error
between putative correspondences and a computed model
@param err_tol correspondences within this distance of a computed model are
considered as inliers
@param inliers if not NULL, output as an array of pointers to the final
set of inliers
@param n_in if not NULL and \a inliers is not NULL, output as the final
number of inliers
@return Returns a transformation matrix computed using RANSAC or NULL
on error or if an acceptable transform could not be computed.
*/
CvMat* ransac_xform( struct feature* features, int n, int mtype,
ransac_xform_fn xform_fn, int m, double p_badxform,
ransac_err_fn err_fn, double err_tol,
struct feature*** inliers, int* n_in )
{
struct feature** matched, ** sample, ** consensus, ** consensus_max = NULL;
struct ransac_data* rdata;
CvPoint2D64f* pts, * mpts;
CvMat* M = NULL;
double p, in_frac = RANSAC_INLIER_FRAC_EST;
int i, nm, in, in_min, in_max = 0, k = 0;
nm = get_matched_features( features, n, mtype, &matched );
if( nm < m )
{
fprintf( stderr, "Warning: not enough matches to compute xform, %s" \
" line %d\n", __FILE__, __LINE__ );
goto end;
}
/* initialize random number generator */
srand( time(NULL) );
in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform );
p = pow( 1.0 - pow( in_frac, m ), k );
i = 0;
while( p > p_badxform )
{
sample = draw_ransac_sample( matched, nm, m );
extract_corresp_pts( sample, m, mtype, &pts, &mpts );
M = xform_fn( pts, mpts, m );
if( ! M )
goto iteration_end;
in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
if( in > in_max )
{
if( consensus_max )
free( consensus_max );
consensus_max = consensus;
in_max = in;
in_frac = (double)in_max / nm;
}
else
free( consensus );
cvReleaseMat( &M );
iteration_end:
release_mem( pts, mpts, sample );
p = pow( 1.0 - pow( in_frac, m ), ++k );
}
/* calculate final transform based on best consensus set */
if( in_max >= in_min )
{
extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts );
M = xform_fn( pts, mpts, in_max );
in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
cvReleaseMat( &M );
release_mem( pts, mpts, consensus_max );
extract_corresp_pts( consensus, in, mtype, &pts, &mpts );
M = xform_fn( pts, mpts, in );
if( inliers )
{
*inliers = consensus;
consensus = NULL;
}
if( n_in )
*n_in = in;
release_mem( pts, mpts, consensus );
}
else if( consensus_max )
{
if( inliers )
*inliers = NULL;
if( n_in )
*n_in = 0;
free( consensus_max );
}
end:
for( i = 0; i < nm; i++ )
{
rdata = feat_ransac_data( matched[i] );
matched[i]->feature_data = rdata->orig_feat_data;
free( rdata );
}
free( matched );
return M;
}
/*
Calculates a planar homography from point correspondeces using the direct
linear transform. Intended for use as a ransac_xform_fn.
@param pts array of points
@param mpts array of corresponding points; each pts[i], i=0..n-1,
corresponds to mpts[i]
@param n number of points in both pts and mpts; must be at least 4
@return Returns the 3x3 planar homography matrix that transforms points
in pts to their corresponding points in mpts or NULL if fewer than 4
correspondences were provided
*/
CvMat* dlt_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n )
{
CvMat* H, * A, * VT, * D, h, v9;
double _h[9];
int i;
if( n < 4 )
return NULL;
/* set up matrices so we can unstack homography into h; Ah = 0 */
A = cvCreateMat( 2*n, 9, CV_64FC1 );
cvZero( A );
for( i = 0; i < n; i++ )
{
cvmSet( A, 2*i, 3, -pts[i].x );
cvmSet( A, 2*i, 4, -pts[i].y );
cvmSet( A, 2*i, 5, -1.0 );
cvmSet( A, 2*i, 6, mpts[i].y * pts[i].x );
cvmSet( A, 2*i, 7, mpts[i].y * pts[i].y );
cvmSet( A, 2*i, 8, mpts[i].y );
cvmSet( A, 2*i+1, 0, pts[i].x );
cvmSet( A, 2*i+1, 1, pts[i].y );
cvmSet( A, 2*i+1, 2, 1.0 );
cvmSet( A, 2*i+1, 6, -mpts[i].x * pts[i].x );
cvmSet( A, 2*i+1, 7, -mpts[i].x * pts[i].y );
cvmSet( A, 2*i+1, 8, -mpts[i].x );
}
D = cvCreateMat( 9, 9, CV_64FC1 );
VT = cvCreateMat( 9, 9, CV_64FC1 );
cvSVD( A, D, NULL, VT, CV_SVD_MODIFY_A + CV_SVD_V_T );
v9 = cvMat( 1, 9, CV_64FC1, NULL );
cvGetRow( VT, &v9, 8 );
h = cvMat( 1, 9, CV_64FC1, _h );
cvCopy( &v9, &h, NULL );
h = cvMat( 3, 3, CV_64FC1, _h );
H = cvCreateMat( 3, 3, CV_64FC1 );
cvConvert( &h, H );
cvReleaseMat( &A );
cvReleaseMat( &D );
cvReleaseMat( &VT );
return H;
}
/*
Calculates a least-squares planar homography from point correspondeces.
@param pts array of points
@param mpts array of corresponding points; each pts[i], i=0..n-1, corresponds
to mpts[i]
@param n number of points in both pts and mpts; must be at least 4
@return Returns the 3 x 3 least-squares planar homography matrix that
transforms points in pts to their corresponding points in mpts or NULL if
fewer than 4 correspondences were provided
*/
CvMat* lsq_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n )
{
CvMat* H, * A, * B, X;
double x[9];
int i;
if( n < 4 )
{
fprintf( stderr, "Warning: too few points in lsq_homog(), %s line %d\n",
__FILE__, __LINE__ );
return NULL;
}
/* set up matrices so we can unstack homography into X; AX = B */
A = cvCreateMat( 2*n, 8, CV_64FC1 );
B = cvCreateMat( 2*n, 1, CV_64FC1 );
X = cvMat( 8, 1, CV_64FC1, x );
H = cvCreateMat(3, 3, CV_64FC1);
cvZero( A );
for( i = 0; i < n; i++ )
{
cvmSet( A, i, 0, pts[i].x );
cvmSet( A, i+n, 3, pts[i].x );
cvmSet( A, i, 1, pts[i].y );
cvmSet( A, i+n, 4, pts[i].y );
cvmSet( A, i, 2, 1.0 );
cvmSet( A, i+n, 5, 1.0 );
cvmSet( A, i, 6, -pts[i].x * mpts[i].x );
cvmSet( A, i, 7, -pts[i].y * mpts[i].x );
cvmSet( A, i+n, 6, -pts[i].x * mpts[i].y );
cvmSet( A, i+n, 7, -pts[i].y * mpts[i].y );
cvmSet( B, i, 0, mpts[i].x );
cvmSet( B, i+n, 0, mpts[i].y );
}
cvSolve( A, B, &X, CV_SVD );
x[8] = 1.0;
X = cvMat( 3, 3, CV_64FC1, x );
cvConvert( &X, H );
cvReleaseMat( &A );
cvReleaseMat( &B );
return H;
}
/*
Calculates the transfer error between a point and its correspondence for
a given homography, i.e. for a point x, it's correspondence x', and
homography H, computes d(x', Hx)^2.
@param pt a point
@param mpt pt's correspondence
@param H a homography matrix
@return Returns the transfer error between pt and mpt given H
*/
double homog_xfer_err( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* H )
{
CvPoint2D64f xpt = persp_xform_pt( pt, H );
return sqrt( dist_sq_2D( xpt, mpt ) );
}
/*
Performs a perspective transformation on a single point. That is, for a
point (x, y) and a 3 x 3 matrix T this function returns the point
(u, v), where
[x' y' w']^T = T * [x y 1]^T,
and
(u, v) = (x'/w', y'/w').
Note that affine transforms are a subset of perspective transforms.
@param pt a 2D point
@param T a perspective transformation matrix
@return Returns the point (u, v) as above.
*/
CvPoint2D64f persp_xform_pt( CvPoint2D64f pt, CvMat* T )
{
CvMat XY, UV;
double xy[3] = { pt.x, pt.y, 1.0 }, uv[3] = { 0 };
CvPoint2D64f rslt;
cvInitMatHeader( &XY, 3, 1, CV_64FC1, xy, CV_AUTOSTEP );
cvInitMatHeader( &UV, 3, 1, CV_64FC1, uv, CV_AUTOSTEP );
cvMatMul( T, &XY, &UV );
rslt = cvPoint2D64f( uv[0] / uv[2], uv[1] / uv[2] );
return rslt;
}
/************************ Local funciton definitions *************************/
/*
Returns a feature's match according to a specified match type
@param feat feature
@param mtype match type, one of FEATURE_FWD_MATCH, FEATURE_BCK_MATCH, or
FEATURE_MDL_MATCH
@return Returns feat's match corresponding to mtype or NULL for bad mtype
*/
static __inline struct feature* get_match( struct feature* feat, int mtype )
{
if( mtype == FEATURE_MDL_MATCH )
return feat->mdl_match;
if( mtype == FEATURE_BCK_MATCH )
return feat->bck_match;
if( mtype == FEATURE_FWD_MATCH )
return feat->fwd_match;
return NULL;
}
/*
Finds all features with a match of a specified type and stores pointers
to them in an array. Additionally initializes each matched feature's
feature_data field with a ransac_data structure.
@param features array of features
@param n number of features in features
@param mtype match type, one of FEATURE_{FWD,BCK,MDL}_MATCH
@param matched output as an array of pointers to features with a match of
the specified type
@return Returns the number of features output in matched.
*/
static int get_matched_features( struct feature* features, int n, int mtype,
struct feature*** matched )
{
struct feature** _matched;
struct ransac_data* rdata;
int i, m = 0;
_matched = calloc( n, sizeof( struct feature* ) );
for( i = 0; i < n; i++ )
if( get_match( features + i, mtype ) )
{
rdata = malloc( sizeof( struct ransac_data ) );
memset( rdata, 0, sizeof( struct ransac_data ) );
rdata->orig_feat_data = features[i].feature_data;
_matched[m] = features + i;
_matched[m]->feature_data = rdata;
m++;
}
*matched = _matched;
return m;
}
/*
Calculates the minimum number of inliers as a function of the number of
putative correspondences. Based on equation (7) in
Chum, O. and Matas, J. Matching with PROSAC -- Progressive Sample Consensus.
In <EM>Conference on Computer Vision and Pattern Recognition (CVPR)</EM>,
(2005), pp. 220--226.
@param n number of putative correspondences
@param m min number of correspondences to compute the model in question
@param p_badsupp prob. that a bad model is supported by a correspondence
@param p_badxform desired prob. that the final transformation returned is bad
@return Returns the minimum number of inliers required to guarantee, based
on p_badsupp, that the probability that the final transformation returned
by RANSAC is less than p_badxform
*/
static int calc_min_inliers( int n, int m, double p_badsupp, double p_badxform )
{
double pi, sum;
int i, j;
for( j = m+1; j <= n; j++ )
{
sum = 0;
for( i = j; i <= n; i++ )
{
pi = ( i - m ) * log( p_badsupp ) + ( n - i + m ) * log( 1.0 - p_badsupp ) +
log_factorial( n - m ) - log_factorial( i - m ) -
log_factorial( n - i );
/*
* Last three terms above are equivalent to log( n-m choose i-m )
*/
sum += exp( pi );
}
if( sum < p_badxform )
break;
}
return j;
}
/*
Calculates the natural log of the factorial of a number
@param n number
@return Returns log( n! )
*/
static __inline double log_factorial( int n )
{
double f = 0;
int i;
for( i = 1; i <= n; i++ )
f += log( i );
return f;
}
/*
Draws a RANSAC sample from a set of features.
@param features array of pointers to features from which to sample
@param n number of features in features
@param m size of the sample
@return Returns an array of pointers to the sampled features; the sampled
field of each sampled feature's ransac_data is set to 1
*/
static struct feature** draw_ransac_sample( struct feature** features, int n, int m )
{
struct feature** sample, * feat;
struct ransac_data* rdata;
int i, x;
for( i = 0; i < n; i++ )
{
rdata = feat_ransac_data( features[i] );
rdata->sampled = 0;
}
sample = calloc( m, sizeof( struct feature* ) );
for( i = 0; i < m; i++ )
{
do
{
x = rand() % n;
feat = features[x];
rdata = feat_ransac_data( feat );
}
while( rdata->sampled );
sample[i] = feat;
rdata->sampled = 1;
}
return sample;
}
/*
Extrancs raw point correspondence locations from a set of features
@param features array of features from which to extract points and match
points; each of these is assumed to have a match of type mtype
@param n number of features
@param mtype match type; if FEATURE_MDL_MATCH correspondences are assumed
to be between each feature's img_pt field and it's match's mdl_pt field,
otherwise, correspondences are assumed to be between img_pt and img_pt
@param pts output as an array of raw point locations from features
@param mpts output as an array of raw point locations from features' matches
*/
static void extract_corresp_pts( struct feature** features, int n, int mtype,
CvPoint2D64f** pts, CvPoint2D64f** mpts )
{
struct feature* match;
CvPoint2D64f* _pts, * _mpts;
int i;
_pts = calloc( n, sizeof( CvPoint2D64f ) );
_mpts = calloc( n, sizeof( CvPoint2D64f ) );
if( mtype == FEATURE_MDL_MATCH )
for( i = 0; i < n; i++ )
{
match = get_match( features[i], mtype );
if( ! match )
fatal_error( "feature does not have match of type %d, %s line %d",
mtype, __FILE__, __LINE__ );
_pts[i] = features[i]->img_pt;
_mpts[i] = match->mdl_pt;
}
else
for( i = 0; i < n; i++ )
{
match = get_match( features[i], mtype );
if( ! match )
fatal_error( "feature does not have match of type %d, %s line %d",
mtype, __FILE__, __LINE__ );
_pts[i] = features[i]->img_pt;
_mpts[i] = match->img_pt;
}
*pts = _pts;
*mpts = _mpts;
}
/*
For a given model and error function, finds a consensus from a set of
feature correspondences.
@param features set of pointers to features; every feature is assumed to
have a match of type mtype
@param n number of features in features
@param mtype determines the match field of each feature against which to
measure error; if this is FEATURE_MDL_MATCH, correspondences are assumed
to be between the feature's img_pt field and the match's mdl_pt field;
otherwise matches are assumed to be between img_pt and img_pt
@param M model for which a consensus set is being found
@param err_fn error function used to measure distance from M
@param err_tol correspondences within this distance of M are added to the
consensus set
@param consensus output as an array of pointers to features in the
consensus set
@return Returns the number of points in the consensus set
*/
static int find_consensus( struct feature** features, int n, int mtype,
CvMat* M, ransac_err_fn err_fn, double err_tol,
struct feature*** consensus )
{
struct feature** _consensus;
struct feature* match;
CvPoint2D64f pt, mpt;
double err;
int i, in = 0;
_consensus = calloc( n, sizeof( struct feature* ) );
if( mtype == FEATURE_MDL_MATCH )
for( i = 0; i < n; i++ )
{
match = get_match( features[i], mtype );
if( ! match )
fatal_error( "feature does not have match of type %d, %s line %d",
mtype, __FILE__, __LINE__ );
pt = features[i]->img_pt;
mpt = match->mdl_pt;
err = err_fn( pt, mpt, M );
if( err <= err_tol )
_consensus[in++] = features[i];
}
else
for( i = 0; i < n; i++ )
{
match = get_match( features[i], mtype );
if( ! match )
fatal_error( "feature does not have match of type %d, %s line %d",
mtype, __FILE__, __LINE__ );
pt = features[i]->img_pt;
mpt = match->img_pt;
err = err_fn( pt, mpt, M );
if( err <= err_tol )
_consensus[in++] = features[i];
}
*consensus = _consensus;
return in;
}
/*
Releases memory and reduces code size above
@param pts1 an array of points
@param pts2 an array of points
@param features an array of pointers to features; can be NULL
*/
static __inline void release_mem( CvPoint2D64f* pts1, CvPoint2D64f* pts2,
struct feature** features )
{
free( pts1 );
free( pts2 );
if( features )
free( features );
}
4、min_pq.h构建一个队列
struct min_pq
{
struct pq_node* pq_array; /* array containing priority queue */
int nallocd; /* number of elementsallocated */
int n; /**< number ofelements in pq */
};
实质是一个队列，先入先出
基本功能是插入删除，获取队列首位的元素。
队列的数据类型是struct min_pq
min_pq.h
/**@file
Functions and structures for implementing a minimizing priority queue.
Copyright (C) 2006-2010 Rob Hess <hess@eecs.oregonstate.edu>
@version 1.1.2-20100521
*/
#ifndef MINPQ_H
#define MINPQ_H
#include <stdlib.h>
/******************************* Defs and macros *****************************/
/* initial # of priority queue elements for which to allocate space */
#define MINPQ_INIT_NALLOCD 512
/********************************** Structures *******************************/
/** an element in a minimizing priority queue */
struct pq_node
{
void* data;
int key;
};
/** a minimizing priority queue */
struct min_pq
{
struct pq_node* pq_array; /* array containing priority queue */
int nallocd; /* number of elements allocated */
int n; /**< number of elements in pq */
};
/*************************** Function Prototypes *****************************/
/**
Creates a new minimizing priority queue.
*/
extern struct min_pq* minpq_init();
/**
Inserts an element into a minimizing priority queue.
@param min_pq a minimizing priority queue
@param data the data to be inserted
@param key the key to be associated with \a data
@return Returns 0 on success or 1 on failure.
*/
extern int minpq_insert( struct min_pq* min_pq, void* data, int key );
/**
Returns the element of a minimizing priority queue with the smallest key
without removing it from the queue.
@param min_pq a minimizing priority queue
@return Returns the element of \a min_pq with the smallest key or NULL
if \a min_pq is empty
*/
extern void* minpq_get_min( struct min_pq* min_pq );
/**
Removes and returns the element of a minimizing priority queue with the
smallest key.
@param min_pq a minimizing priority queue
@return Returns the element of \a min_pq with the smallest key of NULL
if \a min_pq is empty
*/
extern void* minpq_extract_min( struct min_pq* min_pq );
/**
De-allocates the memory held by a minimizing priorioty queue
@param min_pq pointer to a minimizing priority queue
*/
extern void minpq_release( struct min_pq** min_pq );
#endif
min_pq.c
/*
Functions and structures for implementing a minimizing priority queue.
Copyright (C) 2006-2010 Rob Hess <hess@eecs.oregonstate.edu>
@version 1.1.2-20100521
*/
#include "minpq.h"
#include "utils.h"
#include <limits.h>
/************************* Local Function Prototypes *************************/
static void restore_minpq_order( struct pq_node*, int, int );
static void decrease_pq_node_key( struct pq_node*, int, int );
/************************** Local Inline Functions ***************************/
/* returns the array index of element i's parent */
static __inline int parent( int i )
{
return ( i - 1 ) / 2;
}
/* returns the array index of element i's right child */
static __inline int right( int i )
{
return 2 * i + 2;
}
/* returns the array index of element i's left child */
static __inline int left( int i )
{
return 2 * i + 1;
}
/********************** Functions prototyped in minpq.h **********************/
/*
Creates a new minimizing priority queue.
*/
struct min_pq* minpq_init()
{
struct min_pq* min_pq;
min_pq = malloc( sizeof( struct min_pq ) );
min_pq->pq_array = calloc( MINPQ_INIT_NALLOCD, sizeof( struct pq_node ) );
min_pq->nallocd = MINPQ_INIT_NALLOCD;
min_pq->n = 0;
return min_pq;
}
/**
Inserts an element into a minimizing priority queue.
@param min_pq a minimizing priority queue
@param data the data to be inserted
@param key the key to be associated with \a data
@return Returns 0 on success or 1 on failure.
*/
int minpq_insert( struct min_pq* min_pq, void* data, int key )
{
int n = min_pq->n;
/* double array allocation if necessary */
if( min_pq->nallocd == n )
{
min_pq->nallocd = array_double( &min_pq->pq_array, min_pq->nallocd,
sizeof( struct pq_node ) );
if( ! min_pq->nallocd )
{
fprintf( stderr, "Warning: unable to allocate memory, %s, line %d\n",
__FILE__, __LINE__ );
return 1;
}
}
min_pq->pq_array[n].data = data;
min_pq->pq_array[n].key = INT_MAX;
decrease_pq_node_key( min_pq->pq_array, min_pq->n, key );
min_pq->n++;
return 0;
}
/*
Returns the element of a minimizing priority queue with the smallest key
without removing it from the queue.
@param min_pq a minimizing priority queue
@return Returns the element of \a min_pq with the smallest key or NULL
if \a min_pq is empty
*/
void* minpq_get_min( struct min_pq* min_pq )
{
if( min_pq->n < 1 )
{
fprintf( stderr, "Warning: PQ empty, %s line %d\n", __FILE__, __LINE__ );
return NULL;
}
return min_pq->pq_array[0].data;
}
/*
Removes and returns the element of a minimizing priority queue with the
smallest key.
@param min_pq a minimizing priority queue
@return Returns the element of \a min_pq with the smallest key of NULL
if \a min_pq is empty
*/
void* minpq_extract_min( struct min_pq* min_pq )
{
void* data;
if( min_pq->n < 1 )
{
fprintf( stderr, "Warning: PQ empty, %s line %d\n", __FILE__, __LINE__ );
return NULL;
}
data = min_pq->pq_array[0].data;
min_pq->n--;
min_pq->pq_array[0] = min_pq->pq_array[min_pq->n];
restore_minpq_order( min_pq->pq_array, 0, min_pq->n );
return data;
}
/*
De-allocates the memory held by a minimizing priorioty queue
@param min_pq pointer to a minimizing priority queue
*/
void minpq_release( struct min_pq** min_pq )
{
if( ! min_pq )
{
fprintf( stderr, "Warning: NULL pointer error, %s line %d\n", __FILE__,
__LINE__ );
return;
}
if( *min_pq && (*min_pq)->pq_array )
{
free( (*min_pq)->pq_array );
free( *min_pq );
*min_pq = NULL;
}
}
/************************ Functions prototyped here **************************/
/*
Decrease a minimizing pq element's key, rearranging the pq if necessary
@param pq_array minimizing priority queue array
@param i index of the element whose key is to be decreased
@param key new value of element <EM>i</EM>'s key; if greater than current
key, no action is taken
*/
static void decrease_pq_node_key( struct pq_node* pq_array, int i, int key )
{
struct pq_node tmp;
if( key > pq_array[i].key )
return;
pq_array[i].key = key;
while( i > 0 && pq_array[i].key < pq_array[parent(i)].key )
{
tmp = pq_array[parent(i)];
pq_array[parent(i)] = pq_array[i];
pq_array[i] = tmp;
i = parent(i);
}
}
/*
Recursively restores correct priority queue order to a minimizing pq array
@param pq_array a minimizing priority queue array
@param i index at which to start reordering
@param n number of elements in \a pq_array
*/
static void restore_minpq_order( struct pq_node* pq_array, int i, int n )
{
struct pq_node tmp;
int l, r, min = i;
l = left( i );
r = right( i );
if( l < n )
if( pq_array[l].key < pq_array[i].key )
min = l;
if( r < n )
if( pq_array[r].key < pq_array[min].key )
min = r;
if( min != i )
{
tmp = pq_array[min];
pq_array[min] = pq_array[i];
pq_array[i] = tmp;
restore_minpq_order( pq_array, min, n );
}
}
5、sift.h
int sift_features(IplImage * img,structfeature ** feat)
找到所有的特征点放在feat当中并返回特征点的个数。
特征点存储在strcutfeature * feat中。
draw_features(img,feat,featCount)
绘制出所有的特征点
int _sift_features (IplImage *img, structfeature **feat, int intvls, double sigma, double contr_thr, int curv_thr, intimg_dbl, int descr_width, int descr_hist_bins)
找到所有的特征点放在feat当中并返回特征点的个数。
与intsift_features(IplImage * img,struct feature ** feat)的区别：sift_features实质上只是 调用了一下_sift_features，而在_sift_features用回使用很多用户设置的参数
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 *****************************/
/**
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 );
#endif
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;
}