• 热门专题

图像处理-sift算法

作者:renshengrumenglibing的专栏  发布日期:2011-11-02 19:16:22
Tag标签:图像处理  
  • sift是图像匹配的非常经典的算法,但是很复杂,要想自己拿C或C++实现很麻烦,如果只是使用的话,有国外某高人维护的sift库,前期只要自己能够调用即可,关键是要熟悉大致的流程,对sift库有个了解,具体的工作只要调用其中的函数即可。匹配效果:

     

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

    一、sift简介

    1、sift算法应用典型场合:

      物体识别

      机器人定位与导航

      图像拼接

      三维建模

      手势识别

      视频跟踪

      笔记鉴定

      指纹与人脸识别

      犯罪现场特征提取

    2、sift算法解决的问题

      目标的旋转、缩放、平移(RST)

      图像仿射/投影变换(视点viewpoint)

      光照影响(illumination)

      目标遮挡(occlusion)

      杂物场景(clutter)

      噪声

    二、sift基本概念

    1、SIFT匹配的定义

      SIFT匹配(Scale-invariant feature transform,尺度不变特征转换)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征,它在空间尺度中寻找极值点,并提取出其位置、尺度、旋转不变量,此算法由 David Lowe 在1999年所发表,2004年完善总结。其应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对。

      局部影像特征的描述与侦测可以帮助辨识物体,SIFT 特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、些微视角改变的容忍度也相当高。基于这些特性,它们是高度显著而且相对容易撷取,在母数庞大的特征数据库中,很容易辨识物体而且鲜有误认。使用 SIFT特征描述对于部分物体遮蔽的侦测率也相当高,甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下,辨识速度可接近即时运算。SIFT特征的信息量大,适合在海量数据库中快速准确匹配。

    2、SIFT特征的主要特点

      从理论上说,SIFT是一种相似不变量,即对图像尺度变化和旋转是不变量。然而,由于构造SIFT特征时,在很多细节上进行了特殊处理,使得SIFT对图像的复杂变形和光照变化具有了较强的适应性,同时运算速度比较快,定位精度比较高。如:

      在多尺度空间采用DOG算子检测关键点,相比传统的基于LOG算子的检测方法,运算速度大大加快;

      关键点的精确定位不仅提高了精度,而且大大提高了关键点的稳定性;

      在构造描述子时,以子区域的统计特性,而不是以单个像素作为研究对象,提高了对图像局部变形的适应能力;

      对于16*16的关键点邻域和4*4的子区域,在处理梯度幅度时都进行了类似于高斯函数的加权处理,强化了中心区域,淡化了边缘区域的影响,从而提高了算法对几何变形的适应性;

      该方法不仅对通用的线性光照模型具有不变性,而且对复杂的光照变化亦具有一定的适应性。关于这部分内容的细节,可参看文献“Distinctive Image Features from Scale-Invariant Keypoints”

      SIFT算法的特点:

             a.SIFT特征是图像的局部特征,其对旋转、尺度缩放、亮度变化保持不变性,对视角变化、仿射变换、噪声也保持一定程     度的稳定性;

             b.独特性(Distinctiveness)好,信息量丰富,适用于在海量特征数据库中进行快速、准确的匹配;

             c.多量性,即使少数的几    个物体也可以产生大量的SIFT特征向量;

             d.高速性,经优化的SIFT匹配算法甚至可以达到实时的要求;

             e.可扩展性,可以很方便的与其他形式   的特征向量进行联合。

    三、sift算法的应用

    1、sift的main函数

     

    #include "sift.h"
    #include "imgfeatures.h"
    #include "kdtree.h"
    #include "utils.h"
    #include "xform.h"
    
    #include <cv.h>
    #include <cxcore.h>
    #include <highgui.h>
    
    #include <stdio.h>
    
    
    /* the maximum number of keypoint NN candidates to check during BBF search */
    #define KDTREE_BBF_MAX_NN_CHKS 200
    
    /* threshold on squared ratio of distances between NN and 2nd NN */
    #define NN_SQ_DIST_RATIO_THR 0.49
    
    /******************************** Globals ************************************/
    
    char img1_file[] = "beaver.png";
    char img2_file[] = "beaver_xform.png";
    
    /********************************** Main *************************************/
    
    
    int main( int argc, char** argv )
    {
    	IplImage* img1, * img2, * stacked;
    	struct feature* feat1, * feat2, * feat;
    	struct feature** nbrs;
    	struct kd_node* kd_root;
    	CvPoint pt1, pt2;
    	double d0, d1;
    	int n1, n2, k, i, m = 0;
    
    	img1 = cvLoadImage( img1_file, 1 );
    	if( ! img1 )
    		fatal_error( "unable to load image from %s", img1_file );
    	img2 = cvLoadImage( img2_file, 1 );
    	if( ! img2 )
    		fatal_error( "unable to load image from %s", img2_file );
    	stacked = stack_imgs( img1, img2 );
    
    	fprintf( stderr, "Finding features in %s...\n", img1_file );
    	n1 = sift_features( img1, &feat1 );
    	fprintf( stderr, "Finding features in %s...\n", img2_file );
    	n2 = sift_features( img2, &feat2 );
    
    
    	kd_root = kdtree_build( feat2, n2 );
    	for( i = 0; i < n1; i++ )
    	{
    		feat = feat1 + i;
    		k = kdtree_bbf_knn( kd_root, feat, 2, &nbrs, KDTREE_BBF_MAX_NN_CHKS );
    		if( k == 2 )
    		{
    			d0 = descr_dist_sq( feat, nbrs[0] );
    			d1 = descr_dist_sq( feat, nbrs[1] );
    			if( d0 < d1 * NN_SQ_DIST_RATIO_THR )
    			{
    				pt1 = cvPoint( cvRound( feat->x ), cvRound( feat->y ) );
    				pt2 = cvPoint( cvRound( nbrs[0]->x ), cvRound( nbrs[0]->y ) );
    				pt2.y += img1->height;
    				cvLine( stacked, pt1, pt2, CV_RGB(255,0,255), 1, 8, 0 );
    				m++;
    				feat1[i].fwd_match = nbrs[0];
    			}
    		}
    		free( nbrs );
    	}
    
    	fprintf( stderr, "Found %d total matches\n", m );
    	cvNamedWindow( "Matches", 1 );
    	cvShowImage( "Matches", stacked );
    	cvSaveImage("1.bmp",stacked,0);
    	cvWaitKey( 0 );
    
    	cvReleaseImage( &stacked );
    	cvReleaseImage( &img1 );
    	cvReleaseImage( &img2 );
    	kdtree_release( kd_root );
    	free( feat1 );
    	free( feat2 );
    	return 0;
    }

    Sift(Scale InvariantFeature Transform)是一个很好的图像匹配算法,同时能处理亮度、平移、旋转、尺度的变化,利用特征点来提取特征描述符,最后在特征描述符之间寻找匹配。

    sift库分析:

    总共包含6 个文件再加上一个main.c文件。

    untils.c

    minpq.c

    imgfeature.c

    sift.c

    xform.c

    kdtree.c

    以及其对应的.h文件

    最核心的文件是sift.c 和 kdtree.c

    1 sift.c 寻找图片中的特征点

    主要步骤:

    a、构建尺度空间,检测极值点,获得尺度不变性;

    b、特征点过滤并进行精确定位,剔除不稳定的特征点;

    c、在特征点处提取特征描述符,为特征点分配方向值;

    d、生成特征描述子,利用特征描述符寻找匹配点;

    以特征点为中心取16*16的邻域作为采样窗口,

    将采样点与特征点的相对方向通过高斯加权后归入包含8个bin的方向直方图,

    最后获得4*4*8的128维特征描述子。

    2 kdtree.c 对两幅图像进行匹配

        主要步骤:

        当两幅图像的Sift特征向量生成以后,下一步就可以采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。

    取图1的某个关键点,通过遍历找到图像2中的距离最近的两个关键点。在这两个关键点中,如果次近距离除以最近距离小于某个阙值,则判定为一对匹配点。

    3 imgfeature.c在以上图片中画上标记。连接对应的匹配的点。

     

    四、具体文件分析

    1、utils.h

             utils是sift算法的一个比较基础的头文件,主要是对图像的基本操纵:

             获取某位置的像素点

             设置某位置的像素点(8位,32位和64位),

             计算两点之间的距离的平方

             在图片某一点画一个“X”

             将两张图片合成为一个,高是二者之和,宽是二者的较大者。img1 在左上角,img2在右下角。

    函数的解释如下:

    static int pixval8 (IplImage *img, int r,int c)

             从灰度图像中获取某点像素

    static void setpix8 (IplImage *img, int r,int c, uchar val)

             设置灰度图像中某点的像素

     

    static float pixval32f (IplImage *img, intr, int c)  

    static void setpix32f (IplImage *img, intr, int c, float val)     

    static double pixval64f (IplImage *img, intr, int c)

    static void setpix64f (IplImage *img, intr, int c, double val)

             在32位和64位图像中获取和设置某点的像素值

     

    void fatal_error (char *format,...)

             错误处理,用VS2010+opencv2.3报错,直接把内部实现注释掉即可,无影响。

     

    char * replace_extension (const char *file,const char *extn)

            获取一个文件全名,将文件名和扩展名连接到一起如 traffic + jpg => traffic.jpg

    char * prepend_path (const char *path,const char *file)

            获取一个文件的全路径,将路径名添加到文件名之前 C:\\ + traffic.jpg => c:\\traffic.jpg

    char * basename (const char *pathname)

            文件名中去掉路径c:\\traffic.jpg => traffic.jpg

    void  progress (int done)

            显示程序进展,用|/-\表示

    int   array_double(void **array, int n, int size)

            数组长度加倍

    double     dist_sq_2D(CvPoint2D64f p1, CvPoint2D64f p2)

            计算两点的对角线距离

    void          draw_x(IplImage *img, CvPoint pt, int r, int w, CvScalar color)

            在点pt画个叉,本质就是在那个点画四条线

    IplImage * stack_imgs (IplImage *img1,IplImage *img2)

            两张图片生成一张图片,高是二者之和,宽是二者的较大者       

    int   win_closed(char *name)

            查看某个窗口是否已经关闭

    utils.h

     

    /**@file
    Miscellaneous utility functions.
    
    Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>
    
    @version 1.1.2-20100521
    */
    
    #ifndef UTILS_H
    #define UTILS_H
    
    #include "cxcore.h"
    
    #include <stdio.h>
    
    /* absolute value */
    #ifndef ABS
    #define ABS(x) ( ( (x) < 0 )? -(x) : (x) )
    #endif
    
    
    /***************************** Inline Functions ******************************/
    
    
    /**
    A function to get a pixel value from an 8-bit unsigned image.
    
    @param img an image
    @param r row
    @param c column
    @return Returns the value of the pixel at (\a r, \a c) in \a img
    */
    static __inline int pixval8( IplImage* img, int r, int c )
    {
    	return (int)( ( (uchar*)(img->imageData + img->widthStep*r) )[c] );
    }
    
    
    /**
    A function to set a pixel value in an 8-bit unsigned image.
    
    @param img an image
    @param r row
    @param c column
    @param val pixel value
    */
    static __inline void setpix8( IplImage* img, int r, int c, uchar val)
    {
    	( (uchar*)(img->imageData + img->widthStep*r) )[c] = val;
    }
    
    
    /**
    A function to get a pixel value from a 32-bit floating-point image.
    
    @param img an image
    @param r row
    @param c column
    @return Returns the value of the pixel at (\a r, \a c) in \a img
    */
    static __inline float pixval32f( IplImage* img, int r, int c )
    {
    	return ( (float*)(img->imageData + img->widthStep*r) )[c];
    }
    
    
    /**
    A function to set a pixel value in a 32-bit floating-point image.
    
    @param img an image
    @param r row
    @param c column
    @param val pixel value
    */
    static __inline void setpix32f( IplImage* img, int r, int c, float val )
    {
    	( (float*)(img->imageData + img->widthStep*r) )[c] = val;
    }
    
    
    /**
    A function to get a pixel value from a 64-bit floating-point image.
    
    @param img an image
    @param r row
    @param c column
    @return Returns the value of the pixel at (\a r, \a c) in \a img
    */
    static __inline double pixval64f( IplImage* img, int r, int c )
    {
    	return (double)( ( (double*)(img->imageData + img->widthStep*r) )[c] );
    }
    
    
    /**
    A function to set a pixel value in a 64-bit floating-point image.
    
    @param img an image
    @param r row
    @param c column
    @param val pixel value
    */
    static __inline void setpix64f( IplImage* img, int r, int c, double val )
    {
    	( (double*)(img->imageData + img->widthStep*r) )[c] = val;
    }
    
    
    /**************************** Function Prototypes ****************************/
    
    
    /**
    Prints an error message and aborts the program.  The error message is
    of the form "Error: ...", where the ... is specified by the \a format
    argument
    
    @param format an error message format string (as with \c printf(3)).
    */
    extern void fatal_error( char* format, ... );
    
    
    /**
    Replaces a file's extension, which is assumed to be everything after the
    last dot ('.') character.
    
    @param file the name of a file
    
    @param extn a new extension for \a file; should not include a dot (i.e.
    	\c "jpg", not \c ".jpg") unless the new file extension should contain
    	two dots.
    
    @return Returns a new string formed as described above.  If \a file does
    	not have an extension, this function simply adds one.
    */
    extern char* replace_extension( const char* file, const char* extn );
    
    
    /**
    A function that removes the path from a filename.  Similar to the Unix
    basename command.
    
    @param pathname a (full) path name
    
    @return Returns the basename of \a pathname.
    */
    extern char* basename( const char* pathname );
    
    
    /**
    Displays progress in the console with a spinning pinwheel.  Every time this
    function is called, the state of the pinwheel is incremented.  The pinwheel
    has four states that loop indefinitely: '|', '/', '-', '\'.
    
    @param done if 0, this function simply increments the state of the pinwheel;
    	otherwise it prints "done"
    */
    extern void progress( int done );
    
    
    /**
    Erases a specified number of characters from a stream.
    
    @param stream the stream from which to erase characters
    @param n the number of characters to erase
    */
    extern void erase_from_stream( FILE* stream, int n );
    
    
    /**
    Doubles the size of an array with error checking
    
    @param array pointer to an array whose size is to be doubled
    @param n number of elements allocated for \a array
    @param size size in bytes of elements in \a array
    
    @return Returns the new number of elements allocated for \a array.  If no
    	memory is available, returns 0 and frees array.
    */
    extern int array_double( void** array, int n, int size );
    
    
    /**
    Calculates the squared distance between two points.
    
    @param p1 a point
    @param p2 another point
    */
    extern double dist_sq_2D( CvPoint2D64f p1, CvPoint2D64f p2 );
    
    
    /**
    Draws an x on an image.
    
    @param img an image
    @param pt the center point of the x
    @param r the x's radius
    @param w the x's line weight
    @param color the color of the x
    */
    extern void draw_x( IplImage* img, CvPoint pt, int r, int w, CvScalar color );
    
    
    /**
    Combines two images by scacking one on top of the other
    
    @param img1 top image
    @param img2 bottom image
    
    @return Returns the image resulting from stacking \a img1 on top if \a img2
    */
    extern IplImage* stack_imgs( IplImage* img1, IplImage* img2 );
    
    
    /**
    Allows user to view an array of images as a video.  Keyboard controls
    are as follows:
    
    <ul>
    <li>Space - start and pause playback</li>
    <li>Page Up - skip forward 10 frames</li>
    <li>Page Down - jump back 10 frames</li>
    <li>Right Arrow - skip forward 1 frame</li>
    <li>Left Arrow - jump back 1 frame</li>
    <li>Backspace - jump back to beginning</li>
    <li>Esc - exit playback</li>
    <li>Closing the window also exits playback</li>
    </ul>
    
    @param imgs an array of images
    @param n number of images in \a imgs
    @param win_name name of window in which images are displayed
    */
    extern void vid_view( IplImage** imgs, int n, char* win_name );
    
    
    /**
    Checks if a HighGUI window is still open or not
    
    @param name the name of the window we're checking
    
    @return Returns 1 if the window named \a name has been closed or 0 otherwise
    */
    extern int win_closed( char* name );
    
    #endif

     ###utils.c

     

    #include "utils.h"
    
    #include <cv.h>
    #include <cxcore.h>
    #include <highgui.h>
    
    #include <errno.h>
    #include <string.h>
    #include <stdlib.h>
    
    
    
    void fatal_error(char* format, ...)
    {
    	
    }
    
    
    
    
    char* replace_extension( const char* file, const char* extn )
    {
    	char* new_file, * lastdot;
    
    	new_file = calloc( strlen( file ) + strlen( extn ) + 2,  sizeof( char ) );
    	strcpy( new_file, file );
    	lastdot = strrchr( new_file, '.' );
    	if( lastdot )
    		*(lastdot + 1) = '\0';
    	else
    		strcat( new_file, "." );
    	strcat( new_file, extn );
    
    	return new_file;
    }
    
    
    
    
    char* basename( const char* pathname )
    {
    	char* base, * last_slash;
    
    	last_slash = strrchr( pathname, '/' );
    	if( ! last_slash )
    	{
    		base = calloc( strlen( pathname ) + 1, sizeof( char ) );
    		strcpy( base, pathname );
    	}
    	else
    	{
    		base = calloc( strlen( last_slash++ ), sizeof( char ) );
    		strcpy( base, last_slash );
    	}
    
    	return base;
    }
    
    
    
    
    void progress( int done )
    {
    	char state[4] = { '|', '/', '-', '\\' };
    	static int cur = -1;
    
    	if( cur == -1 )
    		fprintf( stderr, "  " );
    
    	if( done )
    	{
    		fprintf( stderr, "\b\bdone\n");
    		cur = -1;
    	}
    	else
    	{
    		cur = ( cur + 1 ) % 4;
    		fprintf( stdout, "\b\b%c ", state[cur] );
    		fflush(stderr);
    	}
    }
    
    
    
    
    void erase_from_stream( FILE* stream, int n )
    {
    	int j;
    	for( j = 0; j < n; j++ )
    		fprintf( stream, "\b" );
    	for( j = 0; j < n; j++ )
    		fprintf( stream, " " );
    	for( j = 0; j < n; j++ )
    		fprintf( stream, "\b" );
    }
    
    
    
    
    int array_double( void** array, int n, int size )
    {
    	void* tmp;
    
    	tmp = realloc( *array, 2 * n * size );
    	if( ! tmp )
    	{
    		fprintf( stderr, "Warning: unable to allocate memory in array_double(),"
    				" %s line %d\n", __FILE__, __LINE__ );
    		if( *array )
    			free( *array );
    		*array = NULL;
    		return 0;
    	}
    	*array = tmp;
    	return n*2;
    }
    
    
    
    
    double dist_sq_2D( CvPoint2D64f p1, CvPoint2D64f p2 )
    {
    	double x_diff = p1.x - p2.x;
    	double y_diff = p1.y - p2.y;
    
    	return x_diff * x_diff + y_diff * y_diff;
    }
    
    
    
    
    void draw_x( IplImage* img, CvPoint pt, int r, int w, CvScalar color )
    {
    	cvLine( img, pt, cvPoint( pt.x + r, pt.y + r), color, w, 8, 0 );
    	cvLine( img, pt, cvPoint( pt.x - r, pt.y + r), color, w, 8, 0 );
    	cvLine( img, pt, cvPoint( pt.x + r, pt.y - r), color, w, 8, 0 );
    	cvLine( img, pt, cvPoint( pt.x - r, pt.y - r), color, w, 8, 0 );
    }
    
    
    
    
    extern IplImage* stack_imgs( IplImage* img1, IplImage* img2 )
    {
    	IplImage* stacked = cvCreateImage( cvSize( MAX(img1->width, img2->width),
    										img1->height + img2->height ),
    										IPL_DEPTH_8U, 3 );
    
    	cvZero( stacked );
    	cvSetImageROI( stacked, cvRect( 0, 0, img1->width, img1->height ) );
    	cvAdd( img1, stacked, stacked, NULL );
    	cvSetImageROI( stacked, cvRect(0, img1->height, img2->width, img2->height) );
    	cvAdd( img2, stacked, stacked, NULL );
    	cvResetImageROI( stacked );
    
    	return stacked;
    }
    
    
    
    
    void vid_view( IplImage** imgs, int n, char* win_name )
    {
    	int k, i = 0, playing = 0;
    
    	cvNamedWindow( win_name, 1 );
    	cvShowImage( win_name, imgs[i] );
    	while( ! win_closed( win_name ) )
    	{
    		/* if already playing, advance frame and check for pause */
    		if( playing )
    		{
    			i = MIN( i + 1, n - 1 );
    			cvNamedWindow( win_name, 1 );
    			cvShowImage( win_name, imgs[i] );
    			k = cvWaitKey( 33 );
    			if( k == ' '  ||  i == n - 1 )
    				playing = 0;
    		}
    
    		else
    
    		{
    			k = cvWaitKey( 0 );
    			switch( k )
    			{
    				/* space */
    			case ' ':
    				playing = 1;
    				break;
    
    				/* esc */
    			case 27:
    			case 1048603:
    				cvDestroyWindow( win_name );
    				break;
    
    				/* backspace */
    			case '\b':
    				i = 0;
    				cvNamedWindow( win_name, 1 );
    				cvShowImage( win_name, imgs[i] );
    				break;
    
    				/* left arrow */
    			case 65288:
    			case 1113937:
    				i = MAX( i - 1, 0 );
    				cvNamedWindow( win_name, 1 );
    				cvShowImage( win_name, imgs[i] );
    				break;
    
    				/* right arrow */
    			case 65363:
    			case 1113939:
    				i = MIN( i + 1, n - 1 );
    				cvNamedWindow( win_name, 1 );
    				cvShowImage( win_name, imgs[i] );
    				break;
    
    				/* page up */
    			case 65365:
    			case 1113941:
    				i = MAX( i - 10, 0 );
    				cvNamedWindow( win_name, 1 );
    				cvShowImage( win_name, imgs[i] );
    				break;
    
    				/* page down */
    			case 65366:
    			case 1113942:
    				i = MIN( i + 10, n - 1 );
    				cvNamedWindow( win_name, 1 );
    				cvShowImage( win_name, imgs[i] );
    				break;
    			}
    		}
    	}
    }
    
    
    
    /*
    Checks if a HighGUI window is still open or not
    
    @param name the name of the window we're checking
    
    @return Returns 1 if the window named \a name has been closed or 0 otherwise
    */
    int win_closed( char* win_name )
    {
    	if( ! cvGetWindowHandle(win_name) )
    		return 1;
    	return 0;
    }

    2、imfeature.h

    int import_features (char *filename, inttype, struct feature **feat)

             从一个文件中(txt)读入几个数据,作为特征的初始化参数

             type指的是特征的类型FEATURE_OXFD或者FEATURE_LOWE

    int export_features (char *filename, structfeature *feat, int n)

             将几个参数存到一个文件中去,参数指的是feature 这个结构,feature**指的是这个结构的数组的指针      ,传递指针可以避免数据的复制

    void draw_features (IplImage *img, structfeature *feat, int n)

            将特征绘制到一个图片上

    double     descr_dist_sq(struct feature *f1, struct feature *f2)

            计算两个特征描绘子之间的欧式距离

     

    feature数据结构分析:

    struct feature

    {

            

    double x;                      /**< x coord */

    double y;                      /**< y coord */

    double a;                      /**< Oxford-typeaffine region parameter */

    double b;                      /**< Oxford-typeaffine region parameter */

    double c;                      /**< Oxford-type affineregion parameter */

    double scl;                    /**< scale of aLowe-style feature */

    double ori;                    /**< orientation of aLowe-style feature */

    int d;                         /**< descriptorlength */

    double descr[FEATURE_MAX_D];   /**< descriptor */

    int type;                      /**< feature type,OXFD or LOWE */

    int category;                  /**< all-purpose featurecategory */

    struct feature* fwd_match;     /**< matching feature from forwardimage */

    struct feature* bck_match;     /**< matching feature from backmwardimage */

    struct feature* mdl_match;     /**< matching feature from model */

    CvPoint2D64f img_pt;           /**< location in image */

    CvPoint2D64f mdl_pt;           /**< location in model */

    void* feature_data;            /**< user-definable data */

    };

    imfeature.h

     

    /**@file
    Functions and structures for dealing with image features
    
    Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>
    
    @version 1.1.2-20100521
    */
    
    #ifndef IMGFEATURES_H
    #define IMGFEATURES_H
    
    #include "cxcore.h"
    
    /** FEATURE_OXFD <BR> FEATURE_LOWE */
    enum feature_type
    {
    	FEATURE_OXFD,
    	FEATURE_LOWE,
    };
    
    /** FEATURE_FWD_MATCH <BR> FEATURE_BCK_MATCH <BR> FEATURE_MDL_MATCH */
    enum feature_match_type
    {
    	FEATURE_FWD_MATCH,
    	FEATURE_BCK_MATCH,
    	FEATURE_MDL_MATCH,
    };
    
    
    /* colors in which to display different feature types */
    #define FEATURE_OXFD_COLOR CV_RGB(255,255,0)
    #define FEATURE_LOWE_COLOR CV_RGB(255,0,255)
    
    /** max feature descriptor length */
    #define FEATURE_MAX_D 128
    
    
    /**
    Structure to represent an affine invariant image feature.  The fields
    x, y, a, b, c represent the affine region around the feature:
    
    a(x-u)(x-u) + 2b(x-u)(y-v) + c(y-v)(y-v) = 1
    */
    struct feature
    {
    	double x;                      /**< x coord */
    	double y;                      /**< y coord */
    	double a;                      /**< Oxford-type affine region parameter */
    	double b;                      /**< Oxford-type affine region parameter */
    	double c;                      /**< Oxford-type affine region parameter */
    	double scl;                    /**< scale of a Lowe-style feature */
    	double ori;                    /**< orientation of a Lowe-style feature */
    	int d;                         /**< descriptor length */
    	double descr[FEATURE_MAX_D];   /**< descriptor */
    	int type;                      /**< feature type, OXFD or LOWE */
    	int category;                  /**< all-purpose feature category */
    	struct feature* fwd_match;     /**< matching feature from forward image */
    	struct feature* bck_match;     /**< matching feature from backmward image */
    	struct feature* mdl_match;     /**< matching feature from model */
    	CvPoint2D64f img_pt;           /**< location in image */
    	CvPoint2D64f mdl_pt;           /**< location in model */
    	void* feature_data;            /**< user-definable data */
    };
    
    
    /**
    Reads image features from file.  The file should be formatted as from
    the code provided by the Visual Geometry Group at Oxford or from the
    code provided by David Lowe.
    
    
    @param filename location of a file containing image features
    @param type determines how features are input.  If \a type is FEATURE_OXFD,
    	the input file is treated as if it is from the code provided by the VGG
    	at Oxford: http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html
    	<BR><BR>
    	If \a type is FEATURE_LOWE, the input file is treated as if it is from
    	David Lowe's SIFT code: http://www.cs.ubc.ca/~lowe/keypoints  
    @param feat pointer to an array in which to store imported features; memory for
    	this array is allocated by this function and must be freed by the caller using
    	free(*feat)
    
    @return Returns the number of features imported from filename or -1 on error
    */
    extern int import_features( char* filename, int type, struct feature** feat );
    
    
    /**
    Exports a feature set to a file formatted depending on the type of
    features, as specified in the feature struct's type field.
    
    @param filename name of file to which to export features
    @param feat feature array
    @param n number of features 
    
    @return Returns 0 on success or 1 on error
    */
    extern int export_features( char* filename, struct feature* feat, int n );
    
    
    /**
    Displays a set of features on an image
    
    @param img image on which to display features
    @param feat array of Oxford-type features
    @param n number of features
    */
    extern void draw_features( IplImage* img, struct feature* feat, int n );
    
    
    /**
    Calculates the squared Euclidian distance between two feature descriptors.
    
    @param f1 first feature
    @param f2 second feature
    
    @return Returns the squared Euclidian distance between the descriptors of
    \a f1 and \a f2.
    */
    extern double descr_dist_sq( struct feature* f1, struct feature* f2 );
    
    
    #endif

    imfeature.c

    /*
    Functions and structures for dealing with image features
    
    Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>
    
    @version 1.1.2-20100521
    */
    
    #include "utils.h"
    #include "imgfeatures.h"
    
    #include <cxcore.h>
    
    #include <math.h>
    
    static int import_oxfd_features( char*, struct feature** );
    static int export_oxfd_features( char*, struct feature*, int );
    static void draw_oxfd_features( IplImage*, struct feature*, int );
    static void draw_oxfd_feature( IplImage*, struct feature*, CvScalar );
    
    static int import_lowe_features( char*, struct feature** );
    static int export_lowe_features( char*, struct feature*, int );
    static void draw_lowe_features( IplImage*, struct feature*, int );
    static void draw_lowe_feature( IplImage*, struct feature*, CvScalar );
    
    
    /*
    Reads image features from file.  The file should be formatted as from
    the code provided by the Visual Geometry Group at Oxford:
    
    
    @param filename location of a file containing image features
    @param type determines how features are input.  If \a type is FEATURE_OXFD,
    	the input file is treated as if it is from the code provided by the VGG
    	at Oxford:
    
    	http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html
    
    	If \a type is FEATURE_LOWE, the input file is treated as if it is from
    	David Lowe's SIFT code:
    
    	http://www.cs.ubc.ca/~lowe/keypoints  
    @param features pointer to an array in which to store features
    
    @return Returns the number of features imported from filename or -1 on error
    */
    int import_features( char* filename, int type, struct feature** feat )
    {
    	int n;
    
    	switch( type )
    	{
    	case FEATURE_OXFD:
    		n = import_oxfd_features( filename, feat );
    		break;
    	case FEATURE_LOWE:
    		n = import_lowe_features( filename, feat );
    		break;
    	default:
    		fprintf( stderr, "Warning: import_features(): unrecognized feature" \
    				"type, %s, line %d\n", __FILE__, __LINE__ );
    		return -1;
    	}
    
    	if( n == -1 )
    		fprintf( stderr, "Warning: unable to import features from %s,"	\
    			" %s, line %d\n", filename, __FILE__, __LINE__ );
    	return n;
    }
    
    
    
    /*
    Exports a feature set to a file formatted depending on the type of
    features, as specified in the feature struct's type field.
    
    @param filename name of file to which to export features
    @param feat feature array
    @param n number of features 
    
    @return Returns 0 on success or 1 on error
    */
    int export_features( char* filename, struct feature* feat, int n )
    {
    	int r, type;
    
    	if( n <= 0  ||  ! feat )
    	{
    		fprintf( stderr, "Warning: no features to export, %s line %d\n",
    				__FILE__, __LINE__ );
    		return 1;
    	}
    	type = feat[0].type;
    	switch( type )
    	{
    	case FEATURE_OXFD:
    		r = export_oxfd_features( filename, feat, n );
    		break;
    	case FEATURE_LOWE:
    		r = export_lowe_features( filename, feat, n );
    		break;
    	default:
    		fprintf( stderr, "Warning: export_features(): unrecognized feature" \
    				"type, %s, line %d\n", __FILE__, __LINE__ );
    		return -1;
    	}
    
    	if( r )
    		fprintf( stderr, "Warning: unable to export features to %s,"	\
    				" %s, line %d\n", filename, __FILE__, __LINE__ );
    	return r;
    }
    
    
    /*
    Draws a set of features on an image
    
    @param img image on which to draw features
    @param feat array of Oxford-type features
    @param n number of features
    */
    void draw_features( IplImage* img, struct feature* feat, int n )
    {
    	int type;
    
    	if( n <= 0  ||  ! feat )
    	{
    		fprintf( stderr, "Warning: no features to draw, %s line %d\n",
    				__FILE__, __LINE__ );
    		return;
    	}
    	type = feat[0].type;
    	switch( type )
    	{
    	case FEATURE_OXFD:
    		draw_oxfd_features( img, feat, n );
    		break;
    	case FEATURE_LOWE:
    		draw_lowe_features( img, feat, n );
    		break;
    	default:
    		fprintf( stderr, "Warning: draw_features(): unrecognized feature" \
    			" type, %s, line %d\n", __FILE__, __LINE__ );
    		break;
    	}	
    }
    
    
    
    /*
    Calculates the squared Euclidian distance between two feature descriptors.
    
    @param f1 first feature	
    @param f2 second feature
    
    @return Returns the squared Euclidian distance between the descriptors of
    f1 and f2.
    */
    double descr_dist_sq( struct feature* f1, struct feature* f2 )
    {
    	double diff, dsq = 0;
    	double* descr1, * descr2;
    	int i, d;
    
    	d = f1->d;
    	if( f2->d != d )
    		return DBL_MAX;
    	descr1 = f1->descr;
    	descr2 = f2->descr;
    
    	for( i = 0; i < d; i++ )
    	{
    		diff = descr1[i] - descr2[i];
    		dsq += diff*diff;
    	}
    	return dsq;
    }
    
    
    
    /***************************** Local Functions *******************************/
    
    
    /*
    Reads image features from file.  The file should be formatted as from
    the code provided by the Visual Geometry Group at Oxford:
    
    http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html
    
    @param filename location of a file containing image features
    @param features pointer to an array in which to store features
    
    @return Returns the number of features imported from filename or -1 on error
    */
    static int import_oxfd_features( char* filename, struct feature** features )
    {
    	struct feature* f;
    	int i, j, n, d;
    	double x, y, a, b, c, dv;
    	FILE* file;
    
    	if( ! features )
    		fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );
    
    	if( ! ( file = fopen( filename, "r" ) ) )
    	{
    		fprintf( stderr, "Warning: error opening %s, %s, line %d\n",
    				filename, __FILE__, __LINE__ );
    		return -1;
    	}
    
    	/* read dimension and number of features */
    	if( fscanf( file, " %d %d ", &d, &n ) != 2 )
    	{
    		fprintf( stderr, "Warning: file read error, %s, line %d\n",
    				__FILE__, __LINE__ );
    		return -1;
    	}
    	if( d > FEATURE_MAX_D )
    	{
    		fprintf( stderr, "Warning: descriptor too long, %s, line %d\n",
    				__FILE__, __LINE__ );
    		return -1;
    	}
    
    
    	f = calloc( n, sizeof(struct feature) );
    	for( i = 0; i < n; i++ )
    	{
    		/* read affine region parameters */
    		if( fscanf( file, " %lf %lf %lf %lf %lf ", &x, &y, &a, &b, &c ) != 5 )
    		{
    			fprintf( stderr, "Warning: error reading feature #%d, %s, line %d\n",
    					i+1, __FILE__, __LINE__ );
    			free( f );
    			return -1;
    		}
    		f[i].img_pt.x = f[i].x = x;
    		f[i].img_pt.y = f[i].y = y;
    		f[i].a = a;
    		f[i].b = b;
    		f[i].c = c;
    		f[i].d = d;
    		f[i].type = FEATURE_OXFD;
    
    		/* read descriptor */
    		for( j = 0; j < d; j++ )
    		{
    			if( ! fscanf( file, " %lf ", &dv ) )
    			{
    				fprintf( stderr, "Warning: error reading feature descriptor" \
    						" #%d, %s, line %d\n", i+1, __FILE__, __LINE__ );
    				free( f );
    				return -1;
    			}
    			f[i].descr[j] = dv;
    		}
    
    		f[i].scl = f[i].ori = 0;
    		f[i].category = 0;
    		f[i].fwd_match = f[i].bck_match = f[i].mdl_match = NULL;
    		f[i].mdl_pt.x = f[i].mdl_pt.y = -1;
    		f[i].feature_data = NULL;
    	}
    
    	if( fclose(file) )
    	{
    		fprintf( stderr, "Warning: file close error, %s, line %d\n",
    				__FILE__, __LINE__ );
    		free( f );
    		return -1;
    	}
    
    	*features = f;
    	return n;
    }
    
    
    
    
    /*
    Exports a feature set to a file formatted as one from the code provided
    by the Visual Geometry Group at Oxford:
    
    http://www.robots.ox.ac.uk:5000/~vgg/research/affine/index.html
    
    @param filename name of file to which to export features
    @param feat feature array
    @param n number of features
    
    @return Returns 0 on success or 1 on error
    */
    static int export_oxfd_features( char* filename, struct feature* feat, int n )
    {
    	FILE* file;
    	int i, j, d;
    
    	if( n <= 0 )
    	{
    		fprintf( stderr, "Warning: feature count %d, %s, line %s\n",
    				n, __FILE__, __LINE__ );
    		return 1;
    	}
    	if( ! ( file = fopen( filename, "w" ) ) )
    	{
    		fprintf( stderr, "Warning: error opening %s, %s, line %d\n",
    				filename, __FILE__, __LINE__ );
    		return 1;
    	}
    
    	d = feat[0].d;
    	fprintf( file, "%d\n%d\n", d, n );
    	for( i = 0; i < n; i++ )
    	{
    		fprintf( file, "%f %f %f %f %f", feat[i].x, feat[i].y, feat[i].a,
    				feat[i].b, feat[i].c );
    		for( j = 0; j < d; j++ )
    			fprintf( file, " %f", feat[i].descr[j] );
    		fprintf( file, "\n" );
    	}
    
    	if( fclose(file) )
    	{
    		fprintf( stderr, "Warning: file close error, %s, line %d\n",
    				__FILE__, __LINE__ );
    		return 1;
    	}
    
    	return 0;
    }
    
    
    
    /*
    Draws Oxford-type affine features
    
    @param img image on which to draw features
    @param feat array of Oxford-type features
    @param n number of features
    */
    static void draw_oxfd_features( IplImage* img, struct feature* feat, int n )
    {
    	CvScalar color = CV_RGB( 255, 255, 255 );
    	int i;
    
    	if( img-> nChannels > 1 )
    		color = FEATURE_OXFD_COLOR;
    	for( i = 0; i < n; i++ )
    		draw_oxfd_feature( img, feat + i, color );
    }
    
    
    
    /*
    Draws a single Oxford-type feature
    
    @param img image on which to draw
    @param feat feature to be drawn
    @param color color in which to draw
    */
    static void draw_oxfd_feature( IplImage* img, struct feature* feat, CvScalar color )
    {
    	double m[4] = { feat->a, feat->b, feat->b, feat->c };
    	double v[4] = { 0 };
    	double e[2] = { 0 };
    	CvMat M, V, E;
    	double alpha, l1, l2;
    
    	/* compute axes and orientation of ellipse surrounding affine region */
    	cvInitMatHeader( &M, 2, 2, CV_64FC1, m, CV_AUTOSTEP );
    	cvInitMatHeader( &V, 2, 2, CV_64FC1, v, CV_AUTOSTEP );
    	cvInitMatHeader( &E, 2, 1, CV_64FC1, e, CV_AUTOSTEP );
    	cvEigenVV( &M, &V, &E, DBL_EPSILON, 0, 0 );
    	l1 = 1 / sqrt( e[1] );
    	l2 = 1 / sqrt( e[0] );
    	alpha = -atan2( v[1], v[0] );
    	alpha *= 180 / CV_PI;
    
    	cvEllipse( img, cvPoint( feat->x, feat->y ), cvSize( l2, l1 ), alpha,
    				0, 360, CV_RGB(0,0,0), 3, 8, 0 );
    	cvEllipse( img, cvPoint( feat->x, feat->y ), cvSize( l2, l1 ), alpha,
    				0, 360, color, 1, 8, 0 );
    	cvLine( img, cvPoint( feat->x+2, feat->y ), cvPoint( feat->x-2, feat->y ),
    			color, 1, 8, 0 );
    	cvLine( img, cvPoint( feat->x, feat->y+2 ), cvPoint( feat->x, feat->y-2 ),
    			color, 1, 8, 0 );
    }
    
    
    
    /*
    Reads image features from file.  The file should be formatted as from
    the code provided by David Lowe:
    
    http://www.cs.ubc.ca/~lowe/keypoints/
    
    @param filename location of a file containing image features
    @param features pointer to an array in which to store features
    
    @return Returns the number of features imported from filename or -1 on error
    */
    static int import_lowe_features( char* filename, struct feature** features )
    {
    	struct feature* f;
    	int i, j, n, d;
    	double x, y, s, o, dv;
    	FILE* file;
    
    	if( ! features )
    		fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );
    
    	if( ! ( file = fopen( filename, "r" ) ) )
    	{
    		fprintf( stderr, "Warning: error opening %s, %s, line %d\n",
    			filename, __FILE__, __LINE__ );
    		return -1;
    	}
    
    	/* read number of features and dimension */
    	if( fscanf( file, " %d %d ", &n, &d ) != 2 )
    	{
    		fprintf( stderr, "Warning: file read error, %s, line %d\n",
    				__FILE__, __LINE__ );
    		return -1;
    	}
    	if( d > FEATURE_MAX_D )
    	{
    		fprintf( stderr, "Warning: descriptor too long, %s, line %d\n",
    				__FILE__, __LINE__ );
    		return -1;
    	}
    
    	f = calloc( n, sizeof(struct feature) );
    	for( i = 0; i < n; i++ )
    	{
    		/* read affine region parameters */
    		if( fscanf( file, " %lf %lf %lf %lf ", &y, &x, &s, &o ) != 4 )
    		{
    			fprintf( stderr, "Warning: error reading feature #%d, %s, line %d\n",
    					i+1, __FILE__, __LINE__ );
    			free( f );
    			return -1;
    		}
    		f[i].img_pt.x = f[i].x = x;
    		f[i].img_pt.y = f[i].y = y;
    		f[i].scl = s;
    		f[i].ori = o;
    		f[i].d = d;
    		f[i].type = FEATURE_LOWE;
    
    		/* read descriptor */
    		for( j = 0; j < d; j++ )
    		{
    			if( ! fscanf( file, " %lf ", &dv ) )
    			{
    				fprintf( stderr, "Warning: error reading feature descriptor" \
    						" #%d, %s, line %d\n", i+1, __FILE__, __LINE__ );
    				free( f );
    				return -1;
    			}
    			f[i].descr[j] = dv;
    		}
    
    		f[i].a = f[i].b = f[i].c = 0;
    		f[i].category = 0;
    		f[i].fwd_match = f[i].bck_match = f[i].mdl_match = NULL;
    		f[i].mdl_pt.x = f[i].mdl_pt.y = -1;
    	}
    
    	if( fclose(file) )
    	{
    		fprintf( stderr, "Warning: file close error, %s, line %d\n",
    				__FILE__, __LINE__ );
    		free( f );
    		return -1;
    	}
    
    	*features = f;
    	return n;
    }
    
    
    
    /*
    Exports a feature set to a file formatted as one from the code provided
    by David Lowe:
    
    http://www.cs.ubc.ca/~lowe/keypoints/
    
    @param filename name of file to which to export features
    @param feat feature array
    @param n number of features
    
    @return Returns 0 on success or 1 on error
    */
    static int export_lowe_features( char* filename, struct feature* feat, int n )
    {
    	FILE* file;
    	int i, j, d;
    
    	if( n <= 0 )
    	{
    		fprintf( stderr, "Warning: feature count %d, %s, line %s\n",
    				n, __FILE__, __LINE__ );
    		return 1;
    	}
    	if( ! ( file = fopen( filename, "w" ) ) )
    	{
    		fprintf( stderr, "Warning: error opening %s, %s, line %d\n",
    				filename, __FILE__, __LINE__ );
    		return 1;
    	}
    
    	d = feat[0].d;
    	fprintf( file, "%d %d\n", n, d );
    	for( i = 0; i < n; i++ )
    	{
    		fprintf( file, "%f %f %f %f", feat[i].y, feat[i].x,
    				feat[i].scl, feat[i].ori );
    		for( j = 0; j < d; j++ )
    		{
    			/* write 20 descriptor values per line */
    			if( j % 20 == 0 )
    				fprintf( file, "\n" );
    			fprintf( file, " %d", (int)(feat[i].descr[j]) );
    		}
    		fprintf( file, "\n" );
    	}
    
    	if( fclose(file) )
    	{
    		fprintf( stderr, "Warning: file close error, %s, line %d\n",
    				__FILE__, __LINE__ );
    		return 1;
    	}
    
    	return 0;
    }
    
    
    /*
    Draws Lowe-type features
    
    @param img image on which to draw features
    @param feat array of Oxford-type features
    @param n number of features
    */
    static void draw_lowe_features( IplImage* img, struct feature* feat, int n )
    {
    	CvScalar color = CV_RGB( 255, 255, 255 );
    	int i;
    
    	if( img-> nChannels > 1 )
    		color = FEATURE_LOWE_COLOR;
    	for( i = 0; i < n; i++ )
    		draw_lowe_feature( img, feat + i, color );
    }
    
    
    
    /*
    Draws a single Lowe-type feature
    
    @param img image on which to draw
    @param feat feature to be drawn
    @param color color in which to draw
    */
    static void draw_lowe_feature( IplImage* img, struct feature* feat, CvScalar color )
    {
    	int len, hlen, blen, start_x, start_y, end_x, end_y, h1_x, h1_y, h2_x, h2_y;
    	double scl, ori;
    	double scale = 5.0;
    	double hscale = 0.75;
    	CvPoint start, end, h1, h2;
    
    	/* compute points for an arrow scaled and rotated by feat's scl and ori */
    	start_x = cvRound( feat->x );
    	start_y = cvRound( feat->y );
    	scl = feat->scl;
    	ori = feat->ori;
    	len = cvRound( scl * scale );
    	hlen = cvRound( scl * hscale );
    	blen = len - hlen;
    	end_x = cvRound( len *  cos( ori ) ) + start_x;
    	end_y = cvRound( len * -sin( ori ) ) + start_y;
    	h1_x = cvRound( blen *  cos( ori + CV_PI / 18.0 ) ) + start_x;
    	h1_y = cvRound( blen * -sin( ori + CV_PI / 18.0 ) ) + start_y;
    	h2_x = cvRound( blen *  cos( ori - CV_PI / 18.0 ) ) + start_x;
    	h2_y = cvRound( blen * -sin( ori - CV_PI / 18.0 ) ) + start_y;
    	start = cvPoint( start_x, start_y );
    	end = cvPoint( end_x, end_y );
    	h1 = cvPoint( h1_x, h1_y );
    	h2 = cvPoint( h2_x, h2_y );
    
    	cvLine( img, start, end, color, 1, 8, 0 );
    	cvLine( img, end, h1, color, 1, 8, 0 );
    	cvLine( img, end, h2, color, 1, 8, 0 );
    }

    3、xforms.h

    CvMat * ransac_xform (struct feature*features, int n, int mtype, ransac_xform_fn xform_fn, int m, doublep_badxform, ransac_err_fn err_fn, double err_tol, struct feature ***inliers,int *n_in)

            运用RANSAC方法进行特征匹配,计算出匹配最佳的图像转化矩阵

    CvMat * dlt_homog (CvPoint2D64f *pts,CvPoint2D64f *mpts, int n)

             运用线性变换,进行点匹配计算平面单应性

    CvMat * lsq_homog (CvPoint2D64f *pts,CvPoint2D64f *mpts, int n)

             进行点匹配,计算对角线最短的平面单应性

    double     homog_xfer_err(CvPoint2D64f pt, CvPoint2D64f mpt, CvMat *H)

             对于给定的单应性,计算某点及其匹配的转换误差

    CvPoint2D64f persp_xform_pt (CvPoint2D64fpt, CvMat *T)

            透视变换

    xforms.h

     

    /**@file
    Functions for computing transforms from image feature correspondences
    
    Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>
    @version 1.1.2-20100521
    */
    
    #ifndef XFORM_H
    #define XFORM_H
    
    #include <cxcore.h>
    
    
    /********************************** Structures *******************************/
    
    struct feature;
    
    /** holds feature data relevant to ransac */
    struct ransac_data
    {
    	void* orig_feat_data;
    	int sampled;
    };
    
    /******************************* Defs and macros *****************************/
    
    /* RANSAC error tolerance in pixels */
    #define RANSAC_ERR_TOL 3
    
    /** pessimistic estimate of fraction of inlers for RANSAC */
    #define RANSAC_INLIER_FRAC_EST 0.25
    
    /** estimate of the probability that a correspondence supports a bad model */
    #define RANSAC_PROB_BAD_SUPP 0.10
    
    /* extracts a feature's RANSAC data */
    #define feat_ransac_data( feat ) ( (struct ransac_data*) (feat)->feature_data )
    
    
    /**
    Prototype for transformation functions passed to ransac_xform().  Functions
    of this type should compute a transformation matrix given a set of point
    correspondences.
    
    @param pts array of points
    @param mpts array of corresponding points; each \a pts[\a i], \a i=0..\a n-1,
    	corresponds to \a mpts[\a i]
    @param n number of points in both \a pts and \a mpts
    
    @return Should return a transformation matrix that transforms each point in
    	\a pts to the corresponding point in \a mpts or NULL on failure.
    */
    typedef CvMat* (*ransac_xform_fn)( CvPoint2D64f* pts, CvPoint2D64f* mpts,
    								  int n );
    
    
    /**
    Prototype for error functions passed to ransac_xform().  For a given
    point, its correspondence, and a transform, functions of this type should
    compute a measure of error between the correspondence and the point after
    the point has been transformed by the transform.
    
    @param pt a point
    @param mpt \a pt's correspondence
    @param T a transform
    
    @return Should return a measure of error between \a mpt and \a pt after
    	\a pt has been transformed by the transform \a T.
    */
    typedef double (*ransac_err_fn)( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* M );
    
    
    /***************************** Function Prototypes ***************************/
    
    
    /**
    Calculates a best-fit image transform from image feature correspondences
    using RANSAC.
    
    For more information refer to:
    
    Fischler, M. A. and Bolles, R. C.  Random sample consensus: a paradigm for
    model fitting with applications to image analysis and automated cartography.
    <EM>Communications of the ACM, 24</EM>, 6 (1981), pp. 381--395.
    
    @param features an array of features; only features with a non-NULL match
    	of type \a mtype are used in homography computation
    @param n number of features in \a feat
    @param mtype determines which of each feature's match fields to use
    	for transform computation; should be one of FEATURE_FWD_MATCH,
    	FEATURE_BCK_MATCH, or FEATURE_MDL_MATCH; if this is FEATURE_MDL_MATCH,
    	correspondences are assumed to be between a feature's img_pt field
    	and its match's mdl_pt field, otherwise correspondences are assumed to
    	be between the the feature's img_pt field and its match's img_pt field
    @param xform_fn pointer to the function used to compute the desired
    	transformation from feature correspondences
    @param m minimum number of correspondences necessary to instantiate the
    	transform computed by \a xform_fn
    @param p_badxform desired probability that the final transformation
    	returned by RANSAC is corrupted by outliers (i.e. the probability that
    	no samples of all inliers were drawn)
    @param err_fn pointer to the function used to compute a measure of error
    	between putative correspondences for a given transform
    @param err_tol correspondences within this distance of each other are
    	considered as inliers for a given transform
    @param inliers if not NULL, output as an array of pointers to the final
    	set of inliers; memory for this array is allocated by this function and
    	must be freed by the caller using free(*inliers)
    @param n_in if not NULL, output as the final number of inliers
    
    @return Returns a transformation matrix computed using RANSAC or NULL
    	on error or if an acceptable transform could not be computed.
    */
    extern CvMat* ransac_xform( struct feature* features, int n, int mtype,
    						   ransac_xform_fn xform_fn, int m,
    						   double p_badxform, ransac_err_fn err_fn,
    						   double err_tol, struct feature*** inliers,
    						   int* n_in );
    
    
    /**
    Calculates a planar homography from point correspondeces using the direct
    linear transform.  Intended for use as a ransac_xform_fn.
    
    @param pts array of points
    @param mpts array of corresponding points; each \a pts[\a i], \a i=0..\a
    	n-1, corresponds to \a mpts[\a i]
    @param n number of points in both \a pts and \a mpts; must be at least 4
    
    @return Returns the \f$3 \times 3\f$ planar homography matrix that
    	transforms points in \a pts to their corresponding points in \a mpts
    	or NULL if fewer than 4 correspondences were provided
    */
    extern CvMat* dlt_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n );
    
    
    /**
    Calculates a least-squares planar homography from point correspondeces.
    Intended for use as a ransac_xform_fn.
    
    @param pts array of points
    @param mpts array of corresponding points; each \a pts[\a i], \a i=0..\a n-1,
    	corresponds to \a mpts[\a i]
    @param n number of points in both \a pts and \a mpts; must be at least 4
    
    @return Returns the \f$3 \times 3\f$ least-squares planar homography
    	matrix that transforms points in \a pts to their corresponding points
    	in \a mpts or NULL if fewer than 4 correspondences were provided
    */
    extern CvMat* lsq_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n );
    
    
    /**
    Calculates the transfer error between a point and its correspondence for
    a given homography, i.e. for a point \f$x\f$, it's correspondence \f$x'\f$,
    and homography \f$H\f$, computes \f$d(x', Hx)^2\f$.  Intended for use as a
    ransac_err_fn.
    
    @param pt a point
    @param mpt \a pt's correspondence
    @param H a homography matrix
    
    @return Returns the transfer error between \a pt and \a mpt given \a H
    */
    extern double homog_xfer_err( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* H );
    
    
    /**
    Performs a perspective transformation on a single point.  That is, for a
    point \f$(x, y)\f$ and a \f$3 \times 3\f$ matrix \f$T\f$ this function
    returns the point \f$(u, v)\f$, where<BR>
    
    \f$[x' \ y' \ w']^T = T \times [x \ y \ 1]^T\f$,<BR>
    
    and<BR>
    
    \f$(u, v) = (x'/w', y'/w')\f$.
    
    Note that affine transforms are a subset of perspective transforms.
    
    @param pt a 2D point
    @param T a perspective transformation matrix
    
    @return Returns the point \f$(u, v)\f$ as above.
    */
    extern CvPoint2D64f persp_xform_pt( CvPoint2D64f pt, CvMat* T );
    
    
    #endif

    xform.c

    /*
    This file contains definitions for functions to compute transforms from
    image feature correspondences
    
    Copyright (C) 2006-2010  Rob Hess <hess@eecs.oregonstate.edu>
    
    @version 1.1.2-20100521
    */
    
    #include "xform.h"
    #include "imgfeatures.h"
    #include "utils.h"
    
    #include <cxcore.h>
    
    #include <stdlib.h>
    #include <time.h>
    
    /************************* Local Function Prototypes *************************/
    
    static __inline struct feature* get_match( struct feature*, int );
    static int get_matched_features( struct feature*, int, int, struct feature*** );
    static int calc_min_inliers( int, int, double, double );
    static __inline double log_factorial( int );
    static struct feature** draw_ransac_sample( struct feature**, int, int );
    static void extract_corresp_pts( struct feature**, int, int, CvPoint2D64f**,
    								 CvPoint2D64f** );
    static int find_consensus( struct feature**, int, int, CvMat*, ransac_err_fn,
    						   double, struct feature*** );
    static __inline void release_mem( CvPoint2D64f*, CvPoint2D64f*,
    static struct feature** );
    
    /********************** Functions prototyped in xform.h **********************/
    
    
    /*
    Calculates a best-fit image transform from image feature correspondences
    using RANSAC.
    
    For more information refer to:
    
    Fischler, M. A. and Bolles, R. C.  Random sample consensus: a paradigm for
    model fitting with applications to image analysis and automated cartography.
    <EM>Communications of the ACM, 24</EM>, 6 (1981), pp. 381--395.
    
    @param features an array of features; only features with a non-NULL match
    	of type mtype are used in homography computation
    @param n number of features in feat
    @param mtype determines which of each feature's match fields to use
    	for model computation; should be one of FEATURE_FWD_MATCH,
    	FEATURE_BCK_MATCH, or FEATURE_MDL_MATCH; if this is FEATURE_MDL_MATCH,
    	correspondences are assumed to be between a feature's img_pt field
    	and its match's mdl_pt field, otherwise correspondences are assumed to
    	be between the the feature's img_pt field and its match's img_pt field
    @param xform_fn pointer to the function used to compute the desired
    	transformation from feature correspondences
    @param m minimum number of correspondences necessary to instantiate the
    	model computed by xform_fn
    @param p_badxform desired probability that the final transformation
    	returned by RANSAC is corrupted by outliers (i.e. the probability that
    	no samples of all inliers were drawn)
    @param err_fn pointer to the function used to compute a measure of error
    	between putative correspondences and a computed model
    @param err_tol correspondences within this distance of a computed model are
    	considered as inliers
    @param inliers if not NULL, output as an array of pointers to the final
    	set of inliers
    @param n_in if not NULL and \a inliers is not NULL, output as the final
    	number of inliers
    
    @return Returns a transformation matrix computed using RANSAC or NULL
    	on error or if an acceptable transform could not be computed.
    */
    CvMat* ransac_xform( struct feature* features, int n, int mtype,
    					ransac_xform_fn xform_fn, int m, double p_badxform,
    					ransac_err_fn err_fn, double err_tol,
    struct feature*** inliers, int* n_in )
    {
    	struct feature** matched, ** sample, ** consensus, ** consensus_max = NULL;
    	struct ransac_data* rdata;
    	CvPoint2D64f* pts, * mpts;
    	CvMat* M = NULL;
    	double p, in_frac = RANSAC_INLIER_FRAC_EST;
    	int i, nm, in, in_min, in_max = 0, k = 0;
    
    	nm = get_matched_features( features, n, mtype, &matched );
    	if( nm < m )
    	{
    		fprintf( stderr, "Warning: not enough matches to compute xform, %s" \
    			" line %d\n", __FILE__, __LINE__ );
    		goto end;
    	}
    
    	/* initialize random number generator */
    	srand( time(NULL) );
    
    	in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform );
    	p = pow( 1.0 - pow( in_frac, m ), k );
    	i = 0;
    	while( p > p_badxform )
    	{
    		sample = draw_ransac_sample( matched, nm, m );
    		extract_corresp_pts( sample, m, mtype, &pts, &mpts );
    		M = xform_fn( pts, mpts, m );
    		if( ! M )
    			goto iteration_end;
    		in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
    		if( in > in_max )
    		{
    			if( consensus_max )
    				free( consensus_max );
    			consensus_max = consensus;
    			in_max = in;
    			in_frac = (double)in_max / nm;
    		}
    		else
    			free( consensus );
    		cvReleaseMat( &M );
    
    iteration_end:
    		release_mem( pts, mpts, sample );
    		p = pow( 1.0 - pow( in_frac, m ), ++k );
    	}
    
    	/* calculate final transform based on best consensus set */
    	if( in_max >= in_min )
    	{
    		extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts );
    		M = xform_fn( pts, mpts, in_max );
    		in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
    		cvReleaseMat( &M );
    		release_mem( pts, mpts, consensus_max );
    		extract_corresp_pts( consensus, in, mtype, &pts, &mpts );
    		M = xform_fn( pts, mpts, in );
    		if( inliers )
    		{
    			*inliers = consensus;
    			consensus = NULL;
    		}
    		if( n_in )
    			*n_in = in;
    		release_mem( pts, mpts, consensus );
    	}
    	else if( consensus_max )
    	{
    		if( inliers )
    			*inliers = NULL;
    		if( n_in )
    			*n_in = 0;
    		free( consensus_max );
    	}
    
    end:
    	for( i = 0; i < nm; i++ )
    	{
    		rdata = feat_ransac_data( matched[i] );
    		matched[i]->feature_data = rdata->orig_feat_data;
    		free( rdata );
    	}
    	free( matched );
    	return M;
    }
    
    
    
    /*
    Calculates a planar homography from point correspondeces using the direct
    linear transform.  Intended for use as a ransac_xform_fn.
      
    @param pts array of points
    @param mpts array of corresponding points; each pts[i], i=0..n-1,
      corresponds to mpts[i]
    @param n number of points in both pts and mpts; must be at least 4
      
    @return Returns the 3x3 planar homography matrix that transforms points
      in pts to their corresponding points in mpts or NULL if fewer than 4
      correspondences were provided
    */
    CvMat* dlt_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n )
    {
    	CvMat* H, * A, * VT, * D, h, v9;
    	double _h[9];
    	int i;
    
    	if( n < 4 )
    	    return NULL;
    
    	/* set up matrices so we can unstack homography into h; Ah = 0 */
    	A = cvCreateMat( 2*n, 9, CV_64FC1 );
    	cvZero( A );
    	for( i = 0; i < n; i++ )
    	{
    		cvmSet( A, 2*i, 3, -pts[i].x );
    		cvmSet( A, 2*i, 4, -pts[i].y );
    		cvmSet( A, 2*i, 5, -1.0  );
    		cvmSet( A, 2*i, 6, mpts[i].y * pts[i].x );
    		cvmSet( A, 2*i, 7, mpts[i].y * pts[i].y );
    		cvmSet( A, 2*i, 8, mpts[i].y );
    		cvmSet( A, 2*i+1, 0, pts[i].x );
    		cvmSet( A, 2*i+1, 1, pts[i].y );
    		cvmSet( A, 2*i+1, 2, 1.0  );
    		cvmSet( A, 2*i+1, 6, -mpts[i].x * pts[i].x );
    		cvmSet( A, 2*i+1, 7, -mpts[i].x * pts[i].y );
    		cvmSet( A, 2*i+1, 8, -mpts[i].x );
    	}
    	D = cvCreateMat( 9, 9, CV_64FC1 );
    	VT = cvCreateMat( 9, 9, CV_64FC1 );
    	cvSVD( A, D, NULL, VT, CV_SVD_MODIFY_A + CV_SVD_V_T );
    	v9 = cvMat( 1, 9, CV_64FC1, NULL );
    	cvGetRow( VT, &v9, 8 );
    	h = cvMat( 1, 9, CV_64FC1, _h );
    	cvCopy( &v9, &h, NULL );
    	h = cvMat( 3, 3, CV_64FC1, _h );
    	H = cvCreateMat( 3, 3, CV_64FC1 );
    	cvConvert( &h, H );
    
    	cvReleaseMat( &A );
    	cvReleaseMat( &D );
    	cvReleaseMat( &VT );
    	return H;
    }
    
    
    
    /*
    Calculates a least-squares planar homography from point correspondeces.
    
    @param pts array of points
    @param mpts array of corresponding points; each pts[i], i=0..n-1, corresponds
    	to mpts[i]
    @param n number of points in both pts and mpts; must be at least 4
    
    @return Returns the 3 x 3 least-squares planar homography matrix that
    	transforms points in pts to their corresponding points in mpts or NULL if
    	fewer than 4 correspondences were provided
    */
    CvMat* lsq_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n )
    {
    	CvMat* H, * A, * B, X;
    	double x[9];
    	int i;
    
    	if( n < 4 )
    	{
    		fprintf( stderr, "Warning: too few points in lsq_homog(), %s line %d\n",
    			__FILE__, __LINE__ );
    		return NULL;
    	}
    
    	/* set up matrices so we can unstack homography into X; AX = B */
    	A = cvCreateMat( 2*n, 8, CV_64FC1 );
    	B = cvCreateMat( 2*n, 1, CV_64FC1 );
    	X = cvMat( 8, 1, CV_64FC1, x );
    	H = cvCreateMat(3, 3, CV_64FC1);
    	cvZero( A );
    	for( i = 0; i < n; i++ )
    	{
    		cvmSet( A, i, 0, pts[i].x );
    		cvmSet( A, i+n, 3, pts[i].x );
    		cvmSet( A, i, 1, pts[i].y );
    		cvmSet( A, i+n, 4, pts[i].y );
    		cvmSet( A, i, 2, 1.0 );
    		cvmSet( A, i+n, 5, 1.0 );
    		cvmSet( A, i, 6, -pts[i].x * mpts[i].x );
    		cvmSet( A, i, 7, -pts[i].y * mpts[i].x );
    		cvmSet( A, i+n, 6, -pts[i].x * mpts[i].y );
    		cvmSet( A, i+n, 7, -pts[i].y * mpts[i].y );
    		cvmSet( B, i, 0, mpts[i].x );
    		cvmSet( B, i+n, 0, mpts[i].y );
    	}
    	cvSolve( A, B, &X, CV_SVD );
    	x[8] = 1.0;
    	X = cvMat( 3, 3, CV_64FC1, x );
    	cvConvert( &X, H );
    
    	cvReleaseMat( &A );
    	cvReleaseMat( &B );
    	return H;
    }
    
    
    
    /*
    Calculates the transfer error between a point and its correspondence for
    a given homography, i.e. for a point x, it's correspondence x', and
    homography H, computes d(x', Hx)^2.
    
    @param pt a point
    @param mpt pt's correspondence
    @param H a homography matrix
    
    @return Returns the transfer error between pt and mpt given H
    */
    double homog_xfer_err( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* H )
    {
    	CvPoint2D64f xpt = persp_xform_pt( pt, H );
    
    	return sqrt( dist_sq_2D( xpt, mpt ) );
    }
    
    
    
    /*
    Performs a perspective transformation on a single point.  That is, for a
    point (x, y) and a 3 x 3 matrix T this function returns the point
    (u, v), where
    
    [x' y' w']^T = T * [x y 1]^T,
    
    and
    
    (u, v) = (x'/w', y'/w').
    
    Note that affine transforms are a subset of perspective transforms.
    
    @param pt a 2D point
    @param T a perspective transformation matrix
    
    @return Returns the point (u, v) as above.
    */
    CvPoint2D64f persp_xform_pt( CvPoint2D64f pt, CvMat* T )
    {
    	CvMat XY, UV;
    	double xy[3] = { pt.x, pt.y, 1.0 }, uv[3] = { 0 };
    	CvPoint2D64f rslt;
    
    	cvInitMatHeader( &XY, 3, 1, CV_64FC1, xy, CV_AUTOSTEP );
    	cvInitMatHeader( &UV, 3, 1, CV_64FC1, uv, CV_AUTOSTEP );
    	cvMatMul( T, &XY, &UV );
    	rslt = cvPoint2D64f( uv[0] / uv[2], uv[1] / uv[2] );
    
    	return rslt;
    }
    
    
    /************************ Local funciton definitions *************************/
    
    /*
    Returns a feature's match according to a specified match type
    
    @param feat feature
    @param mtype match type, one of FEATURE_FWD_MATCH, FEATURE_BCK_MATCH, or
    FEATURE_MDL_MATCH
    
    @return Returns feat's match corresponding to mtype or NULL for bad mtype
    */
    static __inline struct feature* get_match( struct feature* feat, int mtype )
    {
    	if( mtype == FEATURE_MDL_MATCH )
    		return feat->mdl_match;
    	if( mtype == FEATURE_BCK_MATCH )
    		return feat->bck_match;
    	if( mtype == FEATURE_FWD_MATCH )
    		return feat->fwd_match;
    	return NULL;
    }
    
    
    
    /*
    Finds all features with a match of a specified type and stores pointers
    to them in an array.  Additionally initializes each matched feature's
    feature_data f				

延伸阅读:

About IT165 - 广告服务 - 隐私声明 - 版权申明 - 免责条款 - 网站地图 - 网友投稿 - 联系方式
本站内容来自于互联网,仅供用于网络技术学习,学习中请遵循相关法律法规