SIFT算法步骤二(在DOG空间中查找极值特征点)

void __scaleSpaceExtrema(const vector<Mat> &dog_pyramid, vector<Feature> &feats, int octaves, int intervals,
    double contrast_thres, int curvature_thres)
{
    double prelim_constrast_thres = 0.5*contrast_thres / intervals;

    int layer_per_octave_dog = intervals + 2;
    for (int oct = 0; oct < octaves; oct++)
    {
        Size s = dog_pyramid[oct*layer_per_octave_dog].size();//获取每一层图片大小
        for (int lay = 1; lay <= intervals; lay++)            //DOG金字塔每组有intervals+2层，首尾两张没法求极值，所以每组求极值的图片有intervals张
        {
            int idx = oct*layer_per_octave_dog + lay;
            const Mat& dog = dog_pyramid[idx];
            
            for (int r = SIFT_IMG_BORDER; r < s.height - SIFT_IMG_BORDER; r++)
            {
                for (int c = SIFT_IMG_BORDER; c < s.width - SIFT_IMG_BORDER; c++)
                {
                    float pixel_val = dog.at<float>(r, c);
                    if (abs(pixel_val) <= prelim_constrast_thres)  //去除低对比度图片
                        continue;
                    if (!__isExtremum(dog_pyramid, idx, r, c))        //去除非极值点
                        continue;

                    Feature feat;
                    //检查初始极值点修正后的点是否仍为合法极值点
                    if (!__interpExtremum(dog_pyramid, feat, idx, r, c, intervals, contrast_thres))
                        continue;
                    //检查该点是否在边缘上，去除边缘点
                    if (__isTooEdgeLike(dog_pyramid[feat.__idx], feat.__r, feat.__c, curvature_thres))
                        continue;
                    feats.push_back(feat);
                }
            }
        }
    }
}

1.去除非极值点

bool __isExtremum(const vector<Mat>& dog_pyramid,int idx,int r,int c)
{
    float pixel_val = dog_pyramid[idx].at<float>(r, c);
    int flag = pixel_val>0 ? 1 : -1;
    //如果pixel_val大于0，则该点必须比周围点都大;否则就必须比周围点都小
    //这样才满足是一个极值点
    //（但为什么不能是pixel_val>0而且pixel_val比周围点都小呢？？）
    for (int i = -1; i <= 1; i++)
    {
        for (int j = -1; j <= 1; j++)
        {
            for (int k = -1; k <= 1; k++)
            {
                if (pixel_val * flag < dog_pyramid[idx + i].at<float>(r + j, c + k) * flag)
                {
                    return false;
                }
            }
        }
    }
}

2.通过插值更新极值点

bool __interpExtremum(const vector<Mat> &dog_pyramid, Feature &feat, int idx, int r, int c, int intervals, double contrast_thres)
{
    int layer_per_octave_dog = intervals + 2;
    Size s = dog_pyramid[idx].size();

    int i = 0;
    double xi, xr, xc;
    while (i < SIFT_MAX_INTERP_STEPS) //如果超过了最大迭代步数还不收敛，则说明该点不是极值点
    {
        __interpStep(dog_pyramid, idx, r, c, xi, xr, xc);    //根据导数=0所求出的方向进行一次迭代
        if (abs(xi) < 0.5&&abs(xr) < 0.5&&abs(xc) < 0.5)    //如果这几个偏移都小于0.5，则认为收敛到了极值点附近
        {
            break;
        }

        int lay = idx%layer_per_octave_dog + round(xi);        //lay用来检测是否超出了某一组范围

        //迭代一次后更新系数
        c += round(xc);
        r += round(xr);
        idx += round(xi);        //更新索引其实是在更新尺度

        //lay用来检测是否超出了某一组范围
        //r和c分别不能超过[SIFT_IMG_BORDER,s.height - SIFT_IMG_BORDER]和[SIFT_IMG_BORDER,s.width - SIFT_IMG_BORDER]
        if (lay < 1 || lay > intervals || c < SIFT_IMG_BORDER || r < SIFT_IMG_BORDER ||
            c >= s.width - SIFT_IMG_BORDER || r >= s.height - SIFT_IMG_BORDER) {
            return false;
        }
        i++;
    }

    if (i >= SIFT_MAX_INTERP_STEPS)
        return false;

    //求插值后所得到点的对比度
    double contrast = __interpContrast(dog_pyramid, idx, r, c, xi, xr, xc);

    if (abs(contrast) < contrast_thres / intervals)
        return false;
    
    int octave = idx / layer_per_octave_dog;

    //The final offset is added to the location
    //of its sample point to get the interpolated estimate for the location of the extremum.
    //但为啥要乘以pow(2.0,octave)呢？因为该图片已经缩小了这个倍数，而__x和__y表示在原图片尺寸中，该特征点的坐标

    feat.__x = (c + xc) * pow(2.0, octave);
    feat.__y = (r + xr) * pow(2.0, octave);
    feat.__contrast = contrast;

    feat.__r = r;
    feat.__c = c;
    feat.__octave = octave; //组数
    feat.__interval = idx % layer_per_octave_dog; //组内层数
    feat.__idx = idx; //全局索引号
    feat.__sub_interval = xi; //final offset中xi的值

    return true;
}

 插值公式是将尺度空间函数在点(r,c,idx)处泰勒展开至二次项，然后对其求导，另导数等于0, 求出更新一次后的δx

void __interpStep(const vector<Mat> &dog_pyramid, int idx, int r, int c, double &xi, double &xr, double &xc) 
{
    //在点(idx,r,c)处展开
    Mat dD = __derivative(dog_pyramid, idx, r, c);
    Mat H = __hessian(dog_pyramid, idx, r, c);
    Mat H_inv = H.inv(DECOMP_SVD);
    Mat X = -H_inv*dD;
    xi = X.at<double>(2, 0);
    xr = X.at<double>(1, 0);
    xc = X.at<double>(0, 0);
}

接下来是把所求的这个δx代入

得到

求出更新后新的点的值(也就是对比度)

double __interpContrast(const vector<Mat> &dog_pyramid, int idx, int r, int c, double xi, double xr, double xc)
{
    Mat dD = __derivative(dog_pyramid, idx, r, c);
    Mat X(3, 1, CV_64FC1);
    X.at<double>(2, 0) = xi;
    X.at<double>(1, 0) = xr;
    X.at<double>(0, 0) = xc;

    Mat t = dD.t()*X;
    return (dog_pyramid[idx].at<float>(r, c) + 0.5 * t.at<double>(0, 0));
}

若所求值小于阈值，则去除该点。

 

3.判断该点是否位于边缘部分

bool __isTooEdgeLike(const Mat &dog, int r, int c, int curvature_thres)
{
    //这个其实就是用Harris角点相关知识判断是否在边上，论文中curvature_thres=10
    double d = dog.at<float>(r, c);
    double dxx = dog.at<float>(r, c + 1) + dog.at<float>(r, c - 1) - 2 * d;
    double dyy = dog.at<float>(r + 1, c) + dog.at<float>(r - 1, c) - 2 * d;
    double dxy = (dog.at<float>(r + 1, c + 1) - dog.at<float>(r + 1, c - 1) -
        dog.at<float>(r - 1, c + 1) + dog.at<float>(r - 1, c - 1)) / 4.0;
    double tr = dxx + dyy;
    double det = dxx*dyy - dxy*dxy;

    if (det <= 0)
        return true;
    if (tr*tr / det < (curvature_thres + 1.0)*(curvature_thres + 1.0) / curvature_thres)
        return false;
    return true;
}

 

至此，求出极值特征点步骤结束

 

 =========================================================================

之后就是对所求极值特征点进行一些后处理

void __calcFeatureScales(vector<Feature> &feats,double sigma,int intervals)
{
    for (auto it = feats.begin(); it != feats.end(); it++)
    {
        float interval = (*it).__interval + (*it).__sub_interval;
        //全局而言，这个特征点对应的sigma
        (*it).__scl = sigma*pow(2.0, (*it).__octave + interval / intervals);
        //组内而言，这个特征点对应的sigma
        (*it).__scl_octave = sigma*pow(2.0, interval / intervals);
    }
}

 

如果曾经将原始图片扩大一倍，那么相应的x，y，sigma都要除以2

void __adjustForImgDbl(vector<Feature> &feats) {
    for (auto it = feats.begin(); it != feats.end(); it++) {
        (*it).__x /= 2.0;
        (*it).__y /= 2.0;
        (*it).__scl /= 2.0;
    }
}