SIFT算法步骤

void extractSiftFeatures(const Mat &img, vector<KeyPoint> &keypoints, Mat &descriptor, int intervals, double sigma,
    double contrast_thres, int curvature_thres, bool img_dbl, int descr_width, int descr_hist_bins)
{
    Mat init_img = __createInitImg(img, img_dbl, sigma); //根据img_dbl决定是否需要将原来的图片扩大作为第0层

    Size s = init_img.size();
    int octaves = log(min(s.width, s.height)) / log(2) - 2;   //确定金字塔的组数

    vector<Mat> gaussian_pyramid;
    __buildGaussPyramid(init_img, gaussian_pyramid, octaves, intervals, sigma); //建立高斯金字塔

    vector<Mat> dog_pyramid;
    __buildDogPyramid(gaussian_pyramid, dog_pyramid, octaves, intervals); //建立DOG金字塔

    vector<Feature> feats;
    __scaleSpaceExtrema(dog_pyramid, feats, octaves, intervals, contrast_thres, curvature_thres); //在DOG空间中查找极值，并记录为特征

    __calcFeatureScales(feats, sigma, intervals); //计算特征点的尺度

    if (img_dbl)
        __adjustForImgDbl(feats);

    __calcFeatureOris(feats, gaussian_pyramid, intervals + 3); //计算特征点的方向

    __computeDescriptors(feats, gaussian_pyramid, intervals + 3, descr_width, descr_hist_bins);//计算描述子

    __feats2KeyPoints(feats, keypoints);

    __featsVec2Mat(feats, descriptor);
}

1.根据img_dbl确定是否需要扩大原来的图像

Mat __createInitImg(const Mat &img, bool img_dbl, double sigma) {
    if (img_dbl) {
        Mat init_img(img.size() * 2, CV_32FC1);
        resize(img, init_img, init_img.size(), 0, 0, INTER_CUBIC);

        //这一句的理解是，图片初始尺度已经是SIFT_INIT_SIGMA(0.5)了，因为第0层尺度是1.6，所以还需要再把
        //图片进行一次高斯模糊，使得尺度变成1.6，那么这次高斯模糊所需要的尺度就是sqrt(1.6^2-(0.5*2)^2)
        //这里0.5*2是因为图片大小扩大了一倍，相当于现在图片的尺度是1.0了
        double sig_diff = sqrt(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4);
        GaussianBlur(init_img, init_img, Size(), sig_diff, sig_diff);
        return init_img;
    }
    else {
        Mat init_img = Mat(img.size(), CV_32FC1);

        double sig_diff = sqrt(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA);
        GaussianBlur(img, init_img, Size(), sig_diff, sig_diff);
        return init_img;
    }
}

 2.有了初始图像，就可以建立高斯金字塔了

void __buildGaussPyramid(const Mat& base, vector<Mat>& gaussian_pyramid, int octaves, int intervals, double sigma)
{
    //octaves是组数，intervals是组内层数，sigma默认1.6
    int layer_per_octave = intervals + 3;    //每一组中层数加3
    vector<double> sigmas(layer_per_octave);
    double k = pow(2.0f, 1.0f / intervals);

    //计算sigma序列，注意因为第0层图片已经经过高斯模糊了，之后的模糊是在第0层的基础上做的
    //那第1层尺度要变成k*sigma，需要在已经模糊化的第0层的基础上再模糊多少，才能让第1层尺度变成k*sigma呢？
    //根据高斯模糊性质，sigmas[1]=sqrt((k*sigma)^2-sigma^2)，这其实和__createInitImg函数中sig_diff道理是一样的
    //之后的sigmas[2]=sqrt((k*k*sigma)^2-sigma^2)=k*sigmas[1]，依次类推
    sigmas[0] = sigma;
    sigmas[1] = sigma*sqrt(k*k - 1); //sigma默认值是1.6，不过这里sigmas[1]的计算方式有点疑问
    for (int i = 2; i < layer_per_octave; i++)
    {
        sigmas[i] = sigmas[i - 1] * k;
    }
    int layers = octaves*layer_per_octave;
    gaussian_pyramid.reserve(layers);

    for (int oct = 0; oct < octaves; oct++)
    {
        for (int lay = 0; lay < layer_per_octave; lay++)
        {
            if (oct == 0 && lay == 0)
            {
                gaussian_pyramid.push_back(base);
                continue;
            }

            if (lay == 0)
            {
                //因为这一层尺度大小就是2*sigma，所以采样后不需要再高斯模糊。这里也可以看出取倒数第三张可以保证尺度连续
                const Mat &last = gaussian_pyramid[gaussian_pyramid.size() - 3]; 
                Size s = last.size();
                Mat down_img(s.height / 2, s.width / 2, CV_32FC1);
                resize(last, down_img, down_img.size(), 0, 0, CV_INTER_NN);
                gaussian_pyramid.push_back(down_img);
                continue;
            }

            const Mat &last = gaussian_pyramid.back();
            Mat smooth_img(last.size(), CV_32FC1);
            GaussianBlur(last, smooth_img, Size(), sigmas[lay], sigmas[lay]);
            gaussian_pyramid.push_back(smooth_img);
        }
    }
}

3.通过高斯金字塔，计算DOG金字塔，就是同一组中高斯金字塔上一层减去下一层

void __buildDogPyramid(const vector<Mat>& gaussian_pyramid, vector<Mat>& dog_pyramid, int octaves, int intervals)
{
    int layer_per_octave_dog = intervals + 2;
    int layer_per_octave_gaussian = intervals + 3;
    dog_pyramid.reserve(octaves*layer_per_octave_dog);
    for (int oct = 0; oct < octaves; oct++)
    {
        for (int lay = 0; lay < layer_per_octave_dog; lay++)
        {
            int idx = oct*layer_per_octave_gaussian + lay;
            dog_pyramid.push_back(gaussian_pyramid[idx + 1] - gaussian_pyramid[idx]);
        }
    }
}

 

这样，SIFT算法初始化工作算是做好了。