图像处理-sift算法

不念不忘少年蓝@ 2022-06-18 11:10 386阅读 0赞

sift是图像匹配的非常经典的算法，但是很复杂，要想自己拿C或C++实现很麻烦，如果只是使用的话，有国外某高人维护的sift库，前期只要自己能够调用即可，关键是要熟悉大致的流程，对sift库有个了解，具体的工作只要调用其中的函数即可。匹配效果：

![0_1320230911PPmR.gif][]

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 FEATURE_LOWE */
 enum feature_type
 {
 	FEATURE_OXFD,
 	FEATURE_LOWE,
 };
 
 /** FEATURE_FWD_MATCH FEATURE_BCK_MATCH 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
 	 
 	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 );
 }

![0_1320231793fV4x.gif][]

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.
 Communications of the ACM, 24, 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 
 
 \f$[x' \ y' \ w']^T = T \times [x \ y \ 1]^T\f$, 
 
 and 
 
 \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.
 Communications of the ACM, 24, 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 Conference on Computer Vision and Pattern Recognition (CVPR),
 (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 i'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.
 International Journal of Computer Vision, 60, 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.
 International Journal of Computer Vision, 60, 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 \ 
 	| Ixy Iyy Iys | 
 	\ 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;
 }

[0_1320230911PPmR.gif]: /images/20220618/717fbdf38e5048109abff38df71a50f8.png
[0_1320231793fV4x.gif]: /images/20220618/1cd373b6af7e44ef9e614fd3293a8569.png