RobHess的SIFT代码解析之kd树

刺骨的言语ヽ痛彻心扉 2021-11-11 03:36 225阅读 0赞

**平台：win10 x64 +VS 2015专业版 +opencv-2.4.11 + gtk\_-bundle\_2.24.10\_win32**

主要参考：1.代码：RobHess的SIFT源码：SIFT+KD树+BBF算法+RANSAC算法

2.书：王永明 王贵锦 《图像局部不变性特征与描述》

[**RobHess的SIFT源码中的几个文件说明？**][RobHess_SIFT]

**RobHess****的****SIFT****源码分析：**

(1) minpq.h和minpq.c文件  
这两个文件中实现了最小优先级队列(Minimizing Priority Queue)，也就是小顶堆，在k-d树的建立和搜索过程中要用到。  
  
(2) kdtree.h和kdtree.c文件  
这两个文件中实现了k-d树的建立以及用BBF(Best Bin First)算法搜索匹配点的函数。  
如果你需要对两个图片中的特征点进行匹配，就要用到这两个文件。

**RobHess的SIFT源码中使用的结构体及存储说明？**

（1）minpq.h——结构体min\_pq和pq\_node  
　　  minpq.c——无

struct pq\_node //最小化优先级队列中的元素

void\* data;

int key;

\};

struct min\_pq //最小化优先级队列

struct pq\_node\* pq\_array;    /\*数组包含优先级队列 \*/

int nallocd;                 /\* 分配的元素数 \*/

int n;                       /\*pq中的元素数量\*/

\};

（2）kdtree.h——结构体kd\_node  
　　  kdtree.c文件——结构体bbf\_data

struct kd\_node //`K-D树中的结点结构`

int ki;                      /\*分割位置(枢轴)的维数索引(哪一维是分割位置)，取值为1-128\*/

double kv;                   /\*枢轴的值(所有特征向量在枢轴索引维数上的分量的中值)\*/

int leaf;                    /\*是否叶子结点的标志 1 if node is a leaf, 0 otherwise \*/

struct feature\* features;    /\*此结点对应的特征点集合(数组)\*/

int n;                       /\*特征点的个数\*/

struct kd\_node\* kd\_left;     /\*左子树\*/

struct kd\_node\* kd\_right;    /\*右子树\*/

\};

struct bbf\_data

double d;

void\* old\_data;

\};

**请在看这个之前先看完以下内容前四个：**

SIFT四步骤和特征匹配及筛选：

[步骤一：建立尺度空间，即建立高斯差分(DoG)金字塔dog\_pyr][RobHess_SIFT]

[步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器][Link 1]

[步骤三：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向][Link 2]

[步骤四：计算特征描述子][Link 3]

[SIFT后特征匹配：KD树+BBF算法][SIFT_KD_BBF]

[SIFT后特征匹配后错误点筛选：RANSAC算法][SIFT_RANSAC]

**[SIFT后特征匹配：KD树+BBF算法][SIFT_KD_BBF]**

**问题及解答：**

（1）问题描述：SIFT检测的特征点后，为什么用KD树算法进行特征匹配？

答： **找到问题的本质：**

之前网上有blog内曾经介绍过SIFT特征匹配算法，特征点匹配和数据库查、图像检索本质上是同一个问题，都可以归结为一个通过距离函数在高维矢量之间进行相似性检索的问题，如何快速而准确地找到查询点的近邻，不少人提出了很多高维空间索引结构和近似查询的算法。  
一般说来，索引结构中相似性查询有两种基本的方式：  
一种是范围查询，范围查询时给定查询点和查询距离阈值，从数据集中查找所有与查询点距离小于阈值的数据  
另一种是K近邻查询，就是给定查询点及正整数K，从数据集中找到距离查询点最近的K个数据，当K=1时，它就是最近邻查询。  
同样，针对特征点匹配也有两种方法：  
最容易的办法就是线性扫描，也就是我们常说的穷举搜索，依次计算样本集E中每个样本到输入实例点的距离，然后抽取出计算出来的最小距离的点即为最近邻点。此种办法简单直白，但当样本集或训练集很大时，它的缺点就立马暴露出来了，举个例子，在物体识别的问题中，可能有数千个甚至数万个SIFT特征点，而去一一计算这成千上万的特征点与输入实例点的距离，明显是不足取的。  
另外一种，就是构建数据索引，因为实际数据一般都会呈现簇状的聚类形态，因此我们想到建立数据索引，然后再进行快速匹配。索引树是一种树结构索引方法，其基本思想是对搜索空间进行层次划分。根据划分的空间是否有混叠可以分为Clipping和Overlapping两种。前者划分空间没有重叠，其代表就是k-d树；后者划分空间相互有交叠，其代表为R树。

**什么是k-d树？**

首先必须搞清楚的是，k-d树是一种空间划分树，说白了，就是把整个空间划分为特定的几个部分，然后在特定空间的部分内进行相关搜索操作。想像一个三维(多维有点为难你的想象力了)空间，kd树按照一定的划分规则把这个三维空间划分了多个空间，如下图所示：

![1471528-20190605173635742-531646544.png][]

（2）问题描述：RobHess的源码如何实现RANSAC，大概思路是这样的？

答：（2.1）代码及说明：——match.c的main函数里

//特征匹配

fprintf( stderr, "Building kd tree...\\n" ); //建立kd树

kd\_root = kdtree\_build( feat2, n2 ); //根据图2的特征点集feat2建立k-d树，返回k-d树根给kd\_root

//遍历特征点集feat1，针对feat1中每个特征点feat，选取符合距离比值条件的匹配点，放到feat的fwd\_match域中

for( i = 0; i < n1; i++ )    //逐点匹配

feat = feat1 + i; //第i个特征点的指针

//在kd\_root中搜索目标点feat的2个最近邻点，存放在nbrs中，返回实际找到的近邻点个数

k = kdtree\_bbf\_knn( kd\_root, feat, 2, &nbrs, KDTREE\_BBF\_MAX\_NN\_CHKS );          //找2个最近点

if( k == 2 ) //只有进行2次以上匹配过程,才算是正常匹配过程

d0 = descr\_dist\_sq( feat, nbrs\[0\] ); //feat与最近邻点的距离的平方

d1 = descr\_dist\_sq( feat, nbrs\[1\] ); //feat与次近邻点的距离的平方

//若d0和d1的比值小于阈值NN\_SQ\_DIST\_RATIO\_THR，则接受此匹配，否则剔除

if( d0 < d1 \* NN\_SQ\_DIST\_RATIO\_THR )      //最近点与次最近点距离之比要小才当做正确匹配，然后画一条线

//pt1,pt2为连线的两个端点，将目标点feat和最近邻点作为匹配点对

pt1 = cvPoint( cvRound( feat->x ), cvRound( feat->y ) );

pt2 = cvPoint( cvRound( nbrs\[0\]->x ), cvRound( nbrs\[0\]->y ) );

pt2.y += img1->height; //由于两幅图是上下排列的，pt2的纵坐标加上图1的高度，作为连线的终点

cvLine( stacked, pt1, pt2, CV\_RGB(255,0,255), 1, 8, 0 );  //画出连线

m++; //统计匹配点对的个数

feat1\[i\].fwd\_match = nbrs\[0\]; //使点feat的fwd\_match域指向其对应的匹配点

free( nbrs ); //释放近邻数组

fprintf( stderr, "Found %d total matches\\n", m ); //总共找到多少组匹配

display\_big\_img( stacked, "Matches" );

cvWaitKey( 0 );

（2.2）kdtree\_build代码及说明：

答：

/\*根据给定的特征点集合建立k-d树  
参数：  
features：特征点数组，注意：此函数将会改变features数组中元素的排列顺序  
n：特征点个数  
返回值：建立好的k-d树的树根指针  
\*/  
struct kd\_node\* kdtree\_build( struct feature\* features, int n )  
\{   
    struct kd\_node\* kd\_root;  
   
    //输入参数检查  
    if( ! features  ||  n <= 0 )  
    \{   
        fprintf( stderr, "Warning: kdtree\_build(): no features, %s, line %d\\n",  
                \_\_FILE\_\_, \_\_LINE\_\_ );  
        return NULL;  
    \}  
   
    //调用函数，用给定的特征点集初始化k-d树节点，返回值作为树根  
    kd\_root = kd\_node\_init( features, n );  
    //调用函数，扩展根节点kd\_root及其左右孩子  
    expand\_kd\_node\_subtree( kd\_root );  
   
    return kd\_root;  
\}

（2.2.1）kd\_node\_init代码及说明：

答：

/\*用给定的特征点集初始化k-d树节点  
参数：  
features：特征点集  
n：特征点个数  
返回值：k-d树节点指针  
\*/  
static struct kd\_node\* kd\_node\_init( struct feature\* features, int n )  
\{   
    struct kd\_node\* kd\_node;  
   
    kd\_node = malloc( sizeof( struct kd\_node ) );//分配内存  
    memset( kd\_node, 0, sizeof( struct kd\_node ) ); //用0填充  
    kd\_node->ki = -1;//枢轴索引  
    kd\_node->features = features;//节点对应的特征点集  
    kd\_node->n = n;//特征点的个数  
   
    return kd\_node;  
\}

问题：memset函数说明？

答：void \*memset(void \*s, int ch, size\_t n);  
函数解释：将s中当前位置后面的n个字节 （typedef unsigned int size\_t ）用 ch 替换并返回 s 。  
memset：作用是在一段内存块中填充某个给定的值，它是对较大的结构体或数组进行清零操作的一种最快方法  。  
memset()函数原型是extern void \*memset(void \*buffer, int c, int count) buffer：为指针或是数组,c：是赋给buffer的值,count：是buffer的长度.

参看：百度百科（memset）：[https://baike.baidu.com/item/memset/4747579?fr=aladdin][https_baike.baidu.com_item_memset_4747579_fr_aladdin]

（2.2.2）expand\_kd\_node\_subtree代码及说明：

答：

/\*

递归的扩展指定的k-d树节点及其左右孩子

kd\_node为k-d树中未展开节点

\*/  
static void expand\_kd\_node\_subtree( struct kd\_node\* kd\_node )  
\{   
    //基本情况：叶子节点  
    /\* base case: leaf node \*/  
    if( kd\_node->n == 1  ||  kd\_node->n == 0 )  
    \{   
        kd\_node->leaf = 1;//叶节点标志位设为1  
        return;  
    \}  
   
    //调用函数，确定节点的枢轴索引和枢轴值  
    assign\_part\_key( kd\_node ); //ki：方差最大的那个维度；kv：各特征点在那个维度上的中间值  
    //在指定k-d树节点上划分特征点集(即根据指定节点的ki和kv值来划分特征点集)  
    partition\_features( kd\_node );//特征点ki位置左树比右树模值小,kv作为分界模值

//继续扩展左右孩子  
    if( kd\_node->kd\_left )  
        expand\_kd\_node\_subtree( kd\_node->kd\_left ); //递归，展开左子树  
    if( kd\_node->kd\_right )  
        expand\_kd\_node\_subtree( kd\_node->kd\_right ); //递归，展开右子树  
\}

（2.2.2.1）assign\_part\_key代码及说明：

答：

/\*确定输入节点的枢轴索引和枢轴值  
参数：kd\_node：输入的k-d树节点  
函数执行完后将给kd\_node的ki和kv成员复制  
\*/  
static void assign\_part\_key( struct kd\_node\* kd\_node )  
\{   
    struct feature\* features;  
    //枢轴的值kv，均值mean，方差var，方差最大值var\_max  
    double kv, x, mean, var, var\_max = 0;  
    double\* tmp;  
    int d, n, i, j, ki = 0; //枢轴索引ki  
   
    features = kd\_node->features;  
    n = kd\_node->n;//结点个数  
    d = features\[0\].d;//特征向量的维数  
   
    //枢轴的索引值就是方差最大的那一维的维数,即n个128维的特征向量中，若第k维的方差最大，则k就是枢轴(分割位置)  
    /\* partition key index is that along which descriptors have most variance \*/  
    for( j = 0; j < d; j++ )  
    \{   
        mean = var = 0;  
        //求第j维的均值  
        for( i = 0; i < n; i++ )  
            mean += features\[i\].descr\[j\];  
        mean /= n;  
        //求第j维的方差  
        for( i = 0; i < n; i++ )  
        \{   
            x = features\[i\].descr\[j\] - mean;  
            var += x \* x;  
        \}  
        var /= n;  
        //找到最大方差的维数  
        if( var > var\_max ) //书P154  
        \{   
            ki = j;//最大方差的维数就是枢轴  
            var\_max = var;  
        \}  
    \}  
   
    //枢轴的值就是所有特征向量的ki维的中值(按ki维排序后中间的那个值)  
    tmp = calloc( n, sizeof( double ) ); //分配空间  
    for( i = 0; i < n; i++ )  
        tmp\[i\] = features\[i\].descr\[ki\]; //n个特征点在ki维上的描述子，赋给tmp数组  
    //调用函数，找tmp数组的中值  
    kv = median\_select( tmp, n ); //kv：各特征点在那个维度上的中间值  
    free( tmp );  
   
    kd\_node->ki = ki;//枢轴的维数索引  
    kd\_node->kv = kv;//枢轴的值  
\}

（2.2.2.1.1）median\_select代码及说明：

答：

/\*找到输入数组的中值  
参数：  
array：输入数组，元素顺序将会被改动  
n：元素个数  
返回值：中值  
\*/  
static double median\_select( double\* array, int n )  
\{   
    //调用函数，找array数组中的第(n-1)/2小的数，即中值  
    return rank\_select( array, n, (n - 1) / 2 );  
\}

（2.2.2.1.1.1）rank\_select代码及说明：

答：

/\*找到输入数组中第r小的数  
参数：  
array：输入数组，元素顺序将会被改动  
n：元素个数  
r：找第r小元素，要选择的元素从0开始的等级  
返回值：第r小的元素值  
\*/ //该函数重新排列数组中的元素  
static double rank\_select( double\* array, int n, int r )  
\{   
    double\* tmp, med;  
    int gr\_5, gr\_tot, rem\_elts, i, j;  
   
    /\* base case \*/  
    if( n == 1 ) //1个特征点，返回array\[0\]  
        return array\[0\];  
   
    //将数组分成5个一组，共gr\_tot组  
    /\* divide array into groups of 5 and sort them \*/  
    gr\_5 = n / 5; //组的个数-1，n/5向下取整  
    gr\_tot = cvCeil( n / 5.0 ); //组的个数，n/5向上取整（cvCeil向上取整）  
    rem\_elts = n % 5;//最后一组中的元素个数  
    tmp = array;  
    //对每组进行插入排序  
    for( i = 0; i < gr\_5; i++ )  
    \{   
        insertion\_sort( tmp, 5 );  
        tmp += 5; //数组下标加5  
    \}  
    //最后一组，不足5个  
    insertion\_sort( tmp, rem\_elts );  
   
    //以递归的方式查找5组中位数的中位数，找中值的中值  
    /\* recursively find the median of the medians of the groups of 5 \*/  
    tmp = calloc( gr\_tot, sizeof( double ) );  
    //将每个5元组中的中值(即下标为2,2+5,...的元素)复制到temp数组  
    for( i = 0, j = 2; i < gr\_5; i++, j += 5 )  
        tmp\[i\] = array\[j\];  
    //最后一组的中值  
    if( rem\_elts )  
        tmp\[i++\] = array\[n - 1 - rem\_elts/2\];  
    //找temp中的中值med，即中值的中值  
    med = rank\_select( tmp, i, ( i - 1 ) / 2 );  
    free( tmp );  
   
    //利用中值的中值划分数组，看划分结果是否是第r小的数，若不是则递归调用rank\_select重新选择  
    /\* partition around median of medians and recursively select if necessary \*/  
    j = partition\_array( array, n, med );//划分数组，返回med在新数组中的索引  
    if( r == j )//结果是第r小的数  
        return med;  
    else if( r < j )//第r小的数在前半部分  
        return rank\_select( array, j, r );  
    else//第r小的数在后半部分  
    \{   
        array += j+1;  
        return rank\_select( array, ( n - j - 1 ), ( r - j - 1 ) );  
    \}  
\}

（2.2.2.1.1.1.1）insertion\_sort代码及说明：

答：

/\*用插入法对输入数组进行升序排序  
参数：  
array：输入数组  
n：元素个数  
\*/  
static void insertion\_sort( double\* array, int n )  
\{   
    double k;  
    int i, j;  
   
    for( i = 1; i < n; i++ )  
    \{   
        k = array\[i\];  
        j = i-1;  
        while( j >= 0  &&  array\[j\] > k )  
        \{   
            array\[j+1\] = array\[j\];  
            j -= 1;  
        \}  
        array\[j+1\] = k;  
    \}  
\}

（2.2.2.1.1.1.2）partition\_array代码及说明：

答：

/\*根据给定的枢轴值分割数组，使数组前部分小于pivot，后部分大于pivot  
参数：  
array：输入数组  
n：数组的元素个数  
pivot：枢轴值  
返回值：分割后枢轴的下标  
\*/  
static int partition\_array( double\* array, int n, double pivot )  
\{   
    double tmp;  
    int p, i, j;  
   
    i = -1;  
    for( j = 0; j < n; j++ )  
        if( array\[j\] <= pivot )  
        \{   
            tmp = array\[++i\];  
            array\[i\] = array\[j\];  
            array\[j\] = tmp;  
            if( array\[i\] == pivot )  
                p = i;//p保存枢轴的下标  
        \}  
    //将枢轴和最后一个小于枢轴的数对换  
    array\[p\] = array\[i\];  
    array\[i\] = pivot;  
   
    return i; //在分区后返回数据透视的索引  
\}

（2.2.2.2）partition\_features代码及说明：

答：

/\*在指定的k-d树节点上划分特征点集  
使得特征点集的前半部分是第ki维小于枢轴的点，后半部分是第ki维大于枢轴的点  
\*/  
static void partition\_features( struct kd\_node\* kd\_node )  
\{   
    struct feature\* features, tmp;  
    double kv;  
    int n, ki, p, i, j = -1;  
   
    features = kd\_node->features;//特征点数组  
    n = kd\_node->n;//特征点个数  
    //printf("%d\\n",n);  
    ki = kd\_node->ki;//枢轴的维数索引(哪一维是枢轴)ki：方差最大的那个维度  
    kv = kd\_node->kv;//枢轴的值 kv：各特征点在那个维度上的中间值  
    for( i = 0; i < n; i++ )  
    \{   
        //若第i个特征点的特征向量的第ki维的值小于kv  
        if( features\[i\].descr\[ki\] <= kv ) //左子树特征点的数据  
        \{   
            tmp = features\[++j\];  
            features\[j\] = features\[i\];  
            features\[i\] = tmp;  
            if( features\[j\].descr\[ki\] == kv )  
                p = j;//p保存枢轴所在的位置  
        \}  
    \}  
    //将枢轴features\[p\]和最后一个小于枢轴的点features\[j\]对换  
    tmp = features\[p\];  
    features\[p\] = features\[j\];  
    features\[j\] = tmp;  
    //此后，枢轴的位置下标为j  
   
    //若所有特征点落在同一侧，则此节点成为叶节点  
    /\* if all records fall on same side of partition, make node a leaf \*/  
    if( j == n - 1 )  
    \{   
        kd\_node->leaf = 1;  
        return;  
    \}  
   
    //初始化左孩子的根节点，左孩子共j+1个特征点  
    kd\_node->kd\_left = kd\_node\_init( features, j + 1 );  
    //初始化右孩子的根节点，右孩子共n-j-1个特征点  
    kd\_node->kd\_right = kd\_node\_init( features + ( j + 1 ), ( n - j - 1 ) );  
\}

（2.3）kdtree\_bbf\_knn代码及说明：

答：

/\*用BBF算法在k-d树中查找指定特征点的k个最近邻特征点  
参数：  
kd\_root：图像特征的k-d树的树根  
feat：目标特征点  
k：近邻个数  
nbrs：k个近邻特征点的指针数组，按到目标特征点的距离升序排列  
      此数组的内存将在本函数中被分配，使用完后必须在调用出释放：free(\*nbrs)  
max\_nn\_chks：搜索的最大次数，超过此值不再搜索  
返回值：存储在nbrs中的近邻个数，返回-1表示失败  
\*/  
int kdtree\_bbf\_knn( struct kd\_node\* kd\_root, struct feature\* feat, int k,  
                    struct feature\*\*\* nbrs, int max\_nn\_chks )  
\{   
    struct kd\_node\* expl; //expl是当前搜索节点  
    struct min\_pq\* min\_pq; //优先级队列  
    struct feature\* tree\_feat, \*\* \_nbrs; //tree\_feat是单个SIFT特征，\_nbrs中存放着查找出来的近邻特征节点  
    struct bbf\_data\* bbf\_data; //bbf\_data是一个用来存放临时特征数据和特征间距离的缓存结构  
    int i, t = 0, n = 0; //t是搜索的最大次数，n是当前最近邻数组中的元素个数  
   
    //输入参数检查  
    if( ! nbrs  ||  ! feat  ||  ! kd\_root )  
    \{   
        fprintf( stderr, "Warning: NULL pointer error, %s, line %d\\n", \_\_FILE\_\_, \_\_LINE\_\_ );  
        return -1;  
    \}  
   
    \_nbrs = calloc( k, sizeof( struct feature\* ) ); //给查找结果分配相应大小的内存  
    min\_pq = minpq\_init(); //min\_pq队列初始化，分配默认大小的空间  
    minpq\_insert( min\_pq, kd\_root, 0 ); //将根节点先插入到min\_pq优先级队列中  
   
    //min\_pq队列没有回溯完且未达到搜索最大次数  
    while( min\_pq->n > 0  &&  t < max\_nn\_chks )  
    \{   
        //从min\_pq中提取(并移除)优先级最高的节点，赋值给当前节点expl  
        expl = (struct kd\_node\*)minpq\_extract\_min( min\_pq );  
        if( ! expl )  
        \{   //出错处理  
            fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\\n",\_\_FILE\_\_, \_\_LINE\_\_ );  
            goto fail;  
        \}  
        //从当前搜索节点expl一直搜索到叶子节点，搜索过程中将未搜索的节点根据优先级放入队列，返回值为叶子节点  
        expl = explore\_to\_leaf( expl, feat, min\_pq );  
        if( ! expl )  
        \{   //出错处理  
            fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\\n",\_\_FILE\_\_, \_\_LINE\_\_ );  
            goto fail;  
        \}  
   
        //比较查找最近邻  
        for( i = 0; i < expl->n; i++ )  
        \{   
            tree\_feat = &expl->features\[i\];//第i个特征点的指针  
            bbf\_data = malloc( sizeof( struct bbf\_data ) );//新建bbf结构  
            if( ! bbf\_data )  
            \{   //出错处理  
                fprintf( stderr, "Warning: unable to allocate memory," " %s line %d\\n", \_\_FILE\_\_, \_\_LINE\_\_ );  
                goto fail;  
            \}  
            bbf\_data->old\_data = tree\_feat->feature\_data;//保存第i个特征点的feature\_data域以前的值  
            bbf\_data->d = descr\_dist\_sq(feat, tree\_feat);//当前搜索点和目标点之间的欧氏距离  
            tree\_feat->feature\_data = bbf\_data;//将bbf结构赋给此特征点的feature\_data域  
            //判断并插入符合条件的特征点到最近邻数组\_nbrs中,插入成功返回1  
            //当最近邻数组中元素个数已达到k时，继续插入元素个数不会增加，但会更新元素的值  
            n += insert\_into\_nbr\_array( tree\_feat, \_nbrs, n, k );  
        \}  
        t++;//搜索次数  
    \}  
   
    minpq\_release( &min\_pq );//释放优先队列  
   
    //对于最近邻数组中的特征点，恢复其feature\_data域的值  
    for( i = 0; i < n; i++ )  
    \{   
        bbf\_data = \_nbrs\[i\]->feature\_data;  
        \_nbrs\[i\]->feature\_data = bbf\_data->old\_data;//将之前的数据赋值给feature\_data域  
        free( bbf\_data );  
    \}  
    \*nbrs = \_nbrs;  
    return n;  
   
    //失败处理  
fail:  
    minpq\_release( &min\_pq );  
    //对于最近邻数组中的特征点，恢复其feature\_data域的值  
    for( i = 0; i < n; i++ )  
    \{   
        bbf\_data = \_nbrs\[i\]->feature\_data;  
        \_nbrs\[i\]->feature\_data = bbf\_data->old\_data;  
        free( bbf\_data );  
    \}  
    free( \_nbrs );  
    \*nbrs = NULL;  
    return -1;  
\}

（2.3.1）minpq\_init代码及说明：

答：

/\*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; //分配的元素数，512  
  min\_pq->n = 0; //pq中元素数量  
  
  return min\_pq;  
\}

问题：calloc函数说明？calloc与malloc有什么区别？

答：函数原型：void \*calloc(size\_t n, size\_t size)；  
功 能： 在内存的动态存储区中分配n个长度为size的连续空间，函数返回一个指向分配起始地址的指针；如果分配不成功，返回NULL。

malloc/calloc/realloc/alloca内存分配函数的区别？

在说明它们具体含义之前，先简单从字面上加以认识，前3个函数有个共同的特点，就是都带有字符”alloc”，就是”allocate”，”分配”的意思，也就是给对象分配足够的内存，” calloc()”是”分配内存给多个对象”，” malloc()”是”分配内存给一个对象”，”realloc()”是”重新分配内存”之意。”free()”就比较简单了，”释放”的意思，就是把之前所分配的内存空间给释放出来。

void \*calloc(size\_t nobj, size\_t size);  
  
分配足够的内存给nobj个大小为size的对象组成的数组, 并返回指向所分配区域的第一个字节的指针;  
若内存不够,则返回NULL. 该空间的初始化大小为0字节.  
char \*p = (char \*) calloc(100,sizeof(char));  
  
void \*malloc(size\_t size);  
分配足够的内存给大小为size的对象, 并返回指向所分配区域的第一个字节的指针;  
若内存不够,则返回NULL. 不对分配的空间进行初始化.  
char \*p = (char \*)malloc(sizeof(char));  
void \*realloc(void \*p, size\_t size);  
将p所指向的对象的大小改为size个字节.  
如果新分配的内存比原内存大, 那么原内存的内容保持不变, 增加的空间不进行初始化.  
如果新分配的内存比原内存小, 那么新内存保持原内存的内容, 增加的空间不进行初始化.  
返回指向新分配空间的指针; 若内存不够,则返回NULL, 原p指向的内存区不变.  
char \*p = (char \*)malloc(sizeof(char));  
p= (char \*)realloc(p, 256);  
void free(void \*p);  
释放p所指向的内存空间; 当p为NULL时, 不起作用.  
p必先调用calloc, malloc或realloc.  
值得注意的有以下5点：  
  
　　1）通过malloc函数得到的堆内存必须使用memset函数来初始化  
malloc函数分配得到的内存空间是未初始化的。因此，一般在使用该内存空间时，要调用另一个函数memset来将其初始化为全0，memset函数的声明如下：void \* memset (void \* p,int c,int n) ;  
该函数可以将指定的内存空间按字节单位置为指定的字符c，其中，p为要清零的内存空间的首地址，c为要设定的值，n为被操作的内存空间的字节长度。如果要用memset清0，变量c实参要为0。  
  
malloc函数和memset函数的操作语句一般如下：  
int \* p=NULL;  
p=(int\*)malloc(sizeof(int));  
if(p==NULL)  
printf(“Can’t get memory!\\n”);  
memset(p,0,siezeof(int));  
　　2）使用malloc函数分配的堆空间在程序结束之前必须释放  
从堆上获得的内存空间在程序结束以后，系统不会将其自动释放，需要程序员来自己管理。一个程序结束时，必须保证所有从堆上获得的内存空间已被安全释放，否则，会导致内存泄露。  
我们可以使用free()函数来释放内存空间，但是，free函数只是释放指针指向的内容，而该指针仍然指向原来指向的地方，此时，指针为野指针，如果此时操作该指针会导致不可预期的错误。安全做法是：在使用free函数释放指针指向的空间之后，将指针的值置为NULL。  
  
　　3）calloc函数的分配的内存也需要自行释放  
calloc函数的功能与malloc函数的功能相似，都是从堆分配内存,它与malloc函数的一个显著不同时是，calloc函数得到的内存空间是经过初始化的，其内容全为0。calloc函数适合为数组申请空间，可以将size设置为数组元素的空间长度，将n设置为数组的容量。  
　　4）如果要使用realloc函数分配的内存，必须使用memset函数对其内存初始化  
realloc函数的功能比malloc函数和calloc函数的功能更为丰富，可以实现内存分配和内存释放的功能。realloc 可以对给定的指针所指的空间进行扩大或者缩小，无论是扩张或是缩小，原有内存的中内容将保持不变。当然，对于缩小，则被缩小的那一部分的内容会丢失。realloc 并不保证调整后的内存空间和原来的内存空间保持同一内存地址。相反，realloc 返回的指针很可能指向一个新的地址。  
所以，在代码中，我们必须将realloc返回的值，重新赋值给 p :  
p = (int \*) realloc(p, sizeof(int) \*15);  
甚至，你可以传一个空指针（0）给 realloc ，则此时realloc 作用完全相当于malloc。  
int\* p = (int \*)realloc (0,sizeof(int) \* 10); //分配一个全新的内存空间，  
这一行，作用完全等同于：  
int\* p = (int \*)malloc(sizeof(int) \* 10);  
　　5）关于alloca()函数  
还有一个函数也值得一提，这就是alloca()。其调用序列与malloc相同，但是它是在当前函数的栈帧上分配存储空间，而不是在堆中。其优点是：当 函数返回时，自动释放它所使用的栈帧，所以不必再为释放空间而费心。其缺点是：某些系统在函数已被调用后不能增加栈帧长度，于是也就不能支持alloca 函数。尽管如此，很多软件包还是使用alloca函数，也有很多系统支持它。  
  
总结：应用时候需要记得，只有calloc可以指定个数和大小，而且能够对分配内存进行初始化，其余函数均不会对内存进行初始化工作，需要自行调用memset（）函数.

参看：1）百度百科（calloc）——[https://baike.baidu.com/item/calloc/10931444?fr=aladdin][https_baike.baidu.com_item_calloc_10931444_fr_aladdin]

2）malloc/calloc/realloc/alloca内存分配函数——[https://www.cnblogs.com/3me-linux/p/3962152.html][https_www.cnblogs.com_3me-linux_p_3962152.html]

（2.3.2）minpq\_insert代码及说明：

答：

//将元素插入到最小化优先级队列，成功时返回0或失败返回1

//min\_pq最小化优先级队列,  data要插入的元素,  key与数据关联的键

int minpq\_insert( struct min\_pq\* min\_pq, void\* data, int key )  
\{   
  int n = min\_pq->n;  
  //必要时进行双数组分配  
  
  if( min\_pq->nallocd == n ) //越界检查  
    \{   
      min\_pq->nallocd = array\_double( (void\*\*)&min\_pq->pq\_array,  
                      min\_pq->nallocd,  
                      sizeof( struct pq\_node ) ); //array\_double在utils.c中，函数功能：数组大小加倍，错误检查  
      if( ! min\_pq->nallocd )  
    \{   
      fprintf( stderr, "Warning: unable to allocate memory, %s, line %d\\n",  
           \_\_FILE\_\_, \_\_LINE\_\_ );  
      return 1; //失败返回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->pq\_array为初始根节点, min\_pq->n初始为0, key初始为0  
  min\_pq->n++;  
  
  return 0; //成功返回0  
\}

（2.3.2.1）decrease\_pq\_node\_key代码及说明：

答：

/\*减少最小化pq元素的键，必要时重新排列pq

pq\_array最小化优先级队列数组

i要减少其键的元素的索引

key键的新值; 如果大于当前关键，不采取任何行动\*/  
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; //初始pq\_array\[0\].key=0  
  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);  
    \}  
\}

（2.3.2.1.1）parent代码及说明：

答：

/\*返回元素i的父级的数组索引\*/  
static inline int parent( int i )  
\{   
  return ( i - 1 ) / 2;  
\}

（2.3.3）minpq\_extract\_min代码及说明：

答：

//删除并返回具有最小键的最小化优先级队列的元素，返回最小键或如果为空，返回NULL

//min\_pq最小化优先级队列

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;  
\}

（2.3.3.1）restore\_minpq\_order代码及说明：

答：

/\*递归地将正确的优先级队列顺序恢复为最小化pq数组

pq\_array最小化优先级队列数组  
i索引开始重新排序

n 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 ); //递归  
    \}  
\}

（2.3.3.1）left和right代码及说明：

答：

/\*返回元素i右子的数组索引\*/  
static inline int right( int i )  
\{  
  return 2 \* i + 2;  
\}  
//返回元素i的左子节点的数组索引  
static inline int left( int i )  
\{  
  return 2 \* i + 1;  
\}

（2.3.4）explore\_to\_leaf代码及说明：

答：

/\*从给定结点搜索k-d树直到叶节点，搜索过程中将未搜索的节点根据优先级放入队列  
优先级队列和搜索路径是同时生成的，这也是BBF算法的精髓所在：在二叉搜索的时  
候将搜索路径另一侧的分支加入到优先级队列中，供回溯时查找。而优先级队列的排  
序就是根据目标特征与分割超平面的距离ABS(kv - feat->descr\[ki\])  
参数：  
kd\_node：要搜索的子树的树根  
feat：目标特征点  
min\_pq：优先级队列  
返回值：叶子节点的指针  
\*/  
static struct kd\_node\* explore\_to\_leaf( struct kd\_node\* kd\_node, struct feature\* feat,  
                                        struct min\_pq\* min\_pq )  
\{   
    //unexpl中存放着优先级队列的候选节点(还未搜索的节点)，expl为当前搜索节点  
    struct kd\_node\* unexpl, \* expl = kd\_node;  
    double kv;//分割枢轴的值  
    int ki;//分割枢轴的维数索引  
   
    //一直搜索到叶子节点，搜索过程中将未搜索的节点根据优先级放入队列  
    while( expl  &&  ! expl->leaf )  
    \{   
        ki = expl->ki; //ki：方差最大的那个维度  
        kv = expl->kv; //kv：各特征点在那个维度上的中间值  
   
        //枢轴的维数索引大于特征点的维数，出错  
        if( ki >= feat->d )  
        \{   
            fprintf( stderr, "Warning: comparing imcompatible descriptors, %s" " line %d\\n", \_\_FILE\_\_, \_\_LINE\_\_ );  
            return NULL;  
        \}  
        //目标特征点ki维的数据小于等于kv，进入左子树搜索  
        if( feat->descr\[ki\] <= kv )  
        \{   
            unexpl = expl->kd\_right;//候选搜索节点是expl的右子树的根  
            expl = expl->kd\_left;//当前搜索节点是expl的左子树的根  
        \}  
        else//目标特征点ki维的数据大于kv，进入右子树搜索  
        \{   
            unexpl = expl->kd\_left;//候选搜索节点是expl的左子树  
            expl = expl->kd\_right;//当前搜索节点是expl的右子树  
        \}  
   
        //将候选节点unexpl根据目标特征点ki维与其父节点的距离插入到优先队列中，距离越小，优先级越大  
        if( minpq\_insert( min\_pq, unexpl, ABS( kv - feat->descr\[ki\] ) ) )  
        \{   
            fprintf( stderr, "Warning: unable to insert into PQ, %s, line %d\\n",\_\_FILE\_\_, \_\_LINE\_\_ );  
            return NULL;  
        \}  
    \}  
   
    return expl;//返回叶子节点的指针  
\}

（2.3.5）descr\_dist\_sq代码及说明：

答：

//计算两个特征描述子间的欧式距离的平方，返回值为欧式距离的平方

//f1为第一个特征点，f2为第二个特征点

double descr\_dist\_sq( struct feature\* f1, struct feature\* f2 )  
\{   
  double diff, dsq = 0;  
  double\* descr1, \* descr2;  
  int i, d;  
  
  d = f1->d; //f1的特征描述子的长度，128  
  if( f2->d != d ) //若f1和f2的特征描述子长度不同，返回  
    return DBL\_MAX; //DBL\_MAX其他文档未找到定义  
  descr1 = f1->descr; //f1的特征描述子，一个double数组  
  descr2 = f2->descr; //f2的特征描述子，一个double数组  
　//计算欧式距离的平方，即对应元素的差的平方和  
  for( i = 0; i < d; i++ )  
    \{   
      diff = descr1\[i\] - descr2\[i\];  
      dsq += diff\*diff;  
    \}  
  return dsq; //返回欧式距离的平方  
\}

（2.3.6）insert\_into\_nbr\_array代码及说明：

答：

/\*插入一个特征点到最近邻数组，使数组中的点按到目标点的距离升序排列  
参数：  
feat：要插入的特征点，其feature\_data域应是指向bbf\_data结构的指针，其中的d值时feat和目标点的距离的平方  
nbrs：最近邻数组  
n：已在最近邻数组中的元素个数  
k：最近邻数组元素个数的最大值  
返回值：若feat成功插入，返回1，否则返回0  
\*/  
static int insert\_into\_nbr\_array( struct feature\* feat, struct feature\*\* nbrs, int n, int k )  
\{   
    struct bbf\_data\* fdata, \* ndata;//fdata是要插入的点的bbf结构，ndata是最近邻数组中的点的bbf结构  
    double dn, df; //dn是最近邻数组中特征点的bbf结构中的距离值，df是要插入的特征点的bbf结构中的距离值  
    int i, ret = 0;  
   
    //原最近邻数组为空  
    if( n == 0 )  
    \{   
        nbrs\[0\] = feat;  
        return 1;//插入成功，返回1  
    \}  
   
    /\* check at end of array \*/  
    fdata = (struct bbf\_data\*)feat->feature\_data;//要插入点的bbf结构  
    df = fdata->d;//要插入的特征点的bbf结构中的距离值  
    ndata = (struct bbf\_data\*)nbrs\[n-1\]->feature\_data;//最近邻数组中的点的bbf结构  
    dn = ndata->d;//最近邻数组中最后一个特征点的bbf结构中的距离值  
   
    //df>=dn，说明要插入在最近邻数组的末尾  
    if( df >= dn )  
    \{   
        //最近邻数组中元素个数已达到最大值，无法插入  
        if( n == k )  
        \{   
            feat->feature\_data = fdata->old\_data; //原来的值还原回去

//在SIFT极值点检测中，是detection\_data结构的指针  
                 //在k-d树搜索中，是bbf\_data结构的指针，搜索结束，还原回detection\_data  
                 //在RANSAC算法中，是ransac\_data结构的指针\*/  
            free( fdata );  
            return 0;//插入失败，返回0  
        \}  
        nbrs\[n\] = feat;//插入到末尾  
        return 1;//插入成功，返回1  
    \}  
   
    //运行到此处说明插入位置不在数组末尾  
    /\* find the right place in the array \*/  
    if( n < k )//最近邻数组中元素个数小于最大值，可插入  
    \{   
        nbrs\[n\] = nbrs\[n-1\];//原数组最后一个点后移  
        ret = 1;//插入结果设为1  
    \}  
    else//最近邻数组中元素个数大于或等于最大值，无法插入，插入结果ret还是0  
    \{//其实也不是无法插入，而是最近邻数组中元素个数不变，但值会更新  
        nbrs\[n-1\]->feature\_data = ndata->old\_data; //原来的值还原回去  
        free( ndata );  
    \}  
    i = n-2;  
    //在最近邻数组中查找要插入的位置  
    while( i >= 0 )  
    \{   
        ndata = (struct bbf\_data\*)nbrs\[i\]->feature\_data;  
        dn = ndata->d;  
        if( dn <= df )//找到插入点  
            break;  
        nbrs\[i+1\] = nbrs\[i\];//一次后移  
        i--;  
    \}  
    i++;  
    nbrs\[i\] = feat;//插入  
   
    return ret;//返回结果  
\}

（2.3.7）minpq\_release代码及说明：

答：

/\* 取消分配由最小化的先前队列保持的内存

min\_pq指向最小化优先级队列的指针\*/  
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;  
    \}  
\}

转载于:https://www.cnblogs.com/Alliswell-WP/p/SIFT\_code6.html

[RobHess_SIFT]: https://www.cnblogs.com/Alliswell-WP/p/SIFT_code1.html
[Link 1]: https://www.cnblogs.com/Alliswell-WP/p/SIFT_code2.html
[Link 2]: https://www.cnblogs.com/Alliswell-WP/p/SIFT_code3.html
[Link 3]: https://www.cnblogs.com/Alliswell-WP/p/SIFT_code4.html
[SIFT_KD_BBF]: https://www.cnblogs.com/Alliswell-WP/p/SIFT_code6.html
[SIFT_RANSAC]: https://www.cnblogs.com/Alliswell-WP/p/SIFT_code5.html
[1471528-20190605173635742-531646544.png]: /images/20211109/0ea958050ad6452b91a0a418537db9b4.png
[https_baike.baidu.com_item_memset_4747579_fr_aladdin]: https://baike.baidu.com/item/memset/4747579?fr=aladdin
[https_baike.baidu.com_item_calloc_10931444_fr_aladdin]: https://baike.baidu.com/item/calloc/10931444?fr=aladdin
[https_www.cnblogs.com_3me-linux_p_3962152.html]: https://www.cnblogs.com/3me-linux/p/3962152.html