sift算法特征描述子构建代码实现---尺度金字塔原理及代码

0.引言

sift针对局部特征进行特征提取，在尺度空间寻找极值点，提取位置，尺度，旋转不变量，生成特征描述子。

总共分四个步骤：

• 尺度金字塔生成
• 关键点/极值点提取
• 生成梯度直方图
• 特征描述子构建

step1

---

1.1 尺度空间理论

目的：模拟图像数据的多尺度特征
基本思想：引入一个被视为尺度的参数sigma, 通过连续变化尺度参数获得不同尺度下的视觉处理信息, 然后综合这些信息以深入地挖掘图像的本质特征。

1.2 构建尺度空间步骤

• 生成-1层灰度层
• guass金字塔
• 差分高斯金字塔

1.2.1生成-1层灰度层

在最开始建立高斯金字塔的时候，要预先模糊输入的图像来作为第0个组的第0层图像，但是，这相当于丢弃了最高的空域的信息。

因此，通常在生成金字塔前，先将图像的尺度扩大一倍来生成第-1组。以便保留原始图像的信息(特别是图像的高频信息，比如边缘，角点)，增加特征点的数量。

在 Lowe的论文中 ，将第0层的初始尺度 σ定为1.6（最模糊），图片的初始尺度定为0.5（最清晰）。如果输入图像的尺寸用双线性插值扩大一倍，那么相当于σ=1.0(放大一倍变得模糊)。
第0层σ为1.6，第-1层sqrt(1.6^2-(2*0.5)^2).
（虽然论文中尺度空间是kσ倍数，但是认为相机拍摄也对图片进行了模糊且服从平方关系，所以此处及后面金字塔代码实现中的σ都按照开方计算
即：第一层本身σ为1，经过-1层高斯模糊得到1.6，则-1层的σ为sqrt(1.6^2-(2*0.5)^2)）

void CreateInitSmoothGray(const Mat& src, Mat &dst, double sigma)
{
    Mat gray;
    Mat up;
    ConvertToGray(src, gray);               //转为灰度
    resize(gray, up, Size(), 2, 2);         //插值放大
                                            //高斯金字塔的-1组的sigma的计算公式  
    double sigma_init = sqrt(sigma*sigma - (INIT_SIGMA * 2)*(INIT_SIGMA * 2));
    GaussianBlur(up, dst, Size(0, 0), sigma_init, sigma_init);
}

1.2.2 高斯金字塔

得到-1层灰度图像，以其为最初图通过下采样生成不同组（octave）的图像，此时构成一个简单的金字塔。但是尺度不具备连续性。因此还需要经过不同σ对某组图像进行高斯模糊，生成一组内多层（layor）尺度图像，构成高斯金字塔。

高斯卷积核是实现尺度变换的唯一线性核，对一副某组（octave）某层(layor)二维图像进行高斯模糊，其尺度空间定义为：

G(x,y,σ) 是尺度可变高斯函数。

首先，需要根据图片尺寸确定组数和层数：

//计算高斯金字塔的组数 
    int octaves = log((double)min(init_gray.rows, init_gray.cols)) / log(2.0) - 2;

每组层数一般设为3~5层。

#define INTERVALS 3

注：为了保持尺度的连续性，往往每组内多增加3层；
下一组的第一张图为上一组的倒数第三张图。
取倒数第三层和多加3层尺度的目的实际上就是为了在生成差分金字塔最后，各层之间的尺度都是一个固定的倍数：2^(1/s)。否则，经过作差后，尺度是不连续的。
其次，计算各层尺度

s为总共的第s层，S为不加3前的层数。
可以看出：
组内相邻层的尺度相差2^(1/s)
组间同一层对应位置尺度相差2.

void GaussianPyramid(const Mat& src, vector<Mat>&guass_pyr, int octaves, int intervals, double sigma)
{
    double *sigmas = new double[intervals + 3];
    //高斯金字塔需要 intervals+3层图像得到nOctaveLayers+2层DoG金字塔  
    guass_pyr.resize(octaves*(intervals + 3));

    //  \sigma_{total}^2 = \sigma_{i}^2 - \sigma_{i-1}^2、  
    // 计算对图像做不同尺度高斯模糊的尺度因子  
    double k = pow(2, 1.0 / intervals);
    sigmas[0] = sigma;
    for (int i = 1; i < intervals + 3; i++)
    {
        double sig_prev = pow(k, i - 1)*sigma;
        double sig_total = sig_prev*k;
        sigmas[i] = sqrt(sig_total*sig_total - sig_prev*sig_prev);//不是成
        printf("sigma[%d]=%lf\n", i, sigmas[i]);
    }

    /**/
    for (int o = 0; o < octaves; o++)
    {
        //为保持尺度连续，每组多3层
        for (int i = 0; i < intervals + 3; i++)
        {
            // dst为第o组（Octave）金字塔  
            Mat &dst = guass_pyr[o*(intervals + 3) + i];
            // 第0组第0层为原始图像  
            if (o == 0 && i == 0)
            {
                dst = src;
            }
            // 每一组第0副图像时上一组倒数第三幅图像隔点采样得到  
            else if (i == 0)
            {
                const Mat base = guass_pyr[(o - 1)*(intervals + 3) + intervals];
                resize(base, dst, Size(base.rows / 2, base.cols / 2), 0, 0);
            }
            else
            {
                // 每一组第i副图像是由第i-1副图像进行sig[i]的高斯模糊得到  
                // 也就是本组图像在sig[i]的尺度空间下的图像 
                const Mat base = guass_pyr[o*(intervals + 3) + i - 1];
                GaussianBlur(base, dst, Size(), sigmas[i], sigmas[i]);//sigmas没有组的差别？？不是隔组乘2吗
            }
        }

    }
    delete[] sigmas;
}

1.2.3 DOG金字塔构建

使用高斯差分金字塔近似LoG算子，在尺度空间检测稳定的关键点。

差分金字塔顾名思义，就是对上述高斯金字塔每相邻两层作差。所以组内层数为INTERNALS+2.比高斯少一层。
为什么使用差分呢？目的是为了找到更好地表征图像特征的极值点。
这里有两个前提：
1 而Lindberg早在1994年发现高斯差分函数（Difference of Gaussian,简称DOG算子）与尺度归一化的拉普拉斯函数非常相似，因此，利用图像间的差分代替微分。
2 2002年，Mikolajczyk在详细的实验比较中发现尺度归一化的高斯拉普拉斯函数的极大值和极小值同其他的特征提取函数如：梯度，Hessian或者Harris角点特征比较，能够产生稳定的图像特征。

void DogPyramid(const vector<Mat>&guass_pyr, vector<Mat>& dog_pyr, int intervals)
{
    int octaves = (int)guass_pyr.size() / (intervals + 3);
    dog_pyr.resize(octaves*(intervals + 2));
    for (int o = 0; o < octaves; o++)
    {
        for (int i = 0; i < intervals + 2; i++)
        {
            // 第o组第i副图像为高斯金字塔中第o组第i+1和i组图像相减得到 
            const Mat& src1 = guass_pyr[o*(intervals + 3) + i];
            const Mat& src2 = guass_pyr[o*(intervals + 3) + i + 1];
            Mat& dst = dog_pyr[o*(intervals + 2) + i];
            subtract(src1, src2, dst, noArray(), CV_64F);
        }
    }
}