SURF C++代码 详细阅读（二）—— 极值点检测 确定极值点精确位置

阅读（一）进行到了 2.2.2 获取特征点 buildResponseMap();

2.2.3 极值点检测

//! Find the image features and write into vector of features
void FastHessian::getIpoints()
{
    // filter index map
    static const int filter_map [OCTAVES][INTERVALS] = {{0,1,2,3}, {1,3,4,5}, {3,5,6,7}, {5,7,8,9}, {7,9,10,11}};

    // Clear the vector of exisiting ipts
    ipts.clear();

    std::clock_t start;
    start = std::clock();

    // Build the response map
    buildResponseMap();

    std::cout << "buildResponseMap took: " << (std::clock() - start) / (double)(CLOCKS_PER_SEC / 1000) << " ms" << std::endl;

    start = std::clock();

    // Get the response layers 判断是否是极值点
    ResponseLayer *b, *m, *t;

    for (int o = 0; o < octaves; ++o) 
    {
        for (int i = 0; i <= 1; ++i)
        {
            b = responseMap.at(filter_map[o][i]);
            m = responseMap.at(filter_map[o][i+1]);
            t = responseMap.at(filter_map[o][i+2]);

            for (int r = 0; r < t->height; ++r)
            {
                for (int c = 0; c < t->width; ++c)
                {
                    //图像每一点是否极值 参数：行 列 第三张图 第二张图 第一张图
                    if (isExtremum(r, c, t, m, b))
                    {
                        interpolateExtremum(r, c, t, m, b);//, &ipts_tmp[r][c], &changed[r][c]);
                    }
                }
            }
        }
    }

     std::cout << "Ipoint vector took: " << (std::clock() - start) / (double)(CLOCKS_PER_SEC / 1000) << " ms" << std::endl;
}

判断极值点，需要在每一Octave层的四个间隔层中，对于每个像素点，在相邻三个层的领域26个邻居中判断行列式近似值响应是否是极值。
代码中需要判断三项：

• 图像滤波边缘上的点 不是极值点
• 行列式近似值小于阈值的点不是极值点
• 3x3x3 邻域存在更大值的点不是极值点

//! Non Maximal Suppression function
//! 图像每一点是否极值 参数：行 列 第三张图 第二张图 第一张图
int FastHessian::isExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)
{
    // bounds check
    //图像滤波边缘上的点 不是极值点
    int layerBorder = (t->filter + 1) / (2 * t->step);
    if (r <= layerBorder || r >= t->height - layerBorder || c <= layerBorder || c >= t->width - layerBorder)
        return 0;

    // check the candidate point in the middle layer is above thresh
    // 对于中间的图（海森矩阵响应 行列式近似值）小于阈值的点不是极值点
    float candidate = m->getResponse(r, c, t);
    if (candidate < thresh) 
        return 0; 

    int ret_val = 1;//1是极值 0不是极值
    for (int rr = -1; rr <=1; ++rr)
    {
        for (int cc = -1; cc <=1; ++cc)
        {
            // if any response in 3x3x3 is greater candidate not maximum
            if (
                t->getResponse(r+rr, c+cc) >= candidate ||
                ((rr != 0 || cc != 0) && m->getResponse(r+rr, c+cc, t) >= candidate) ||
                b->getResponse(r+rr, c+cc, t) >= candidate
                ) 
                ret_val = 0;
        }
    }

    return ret_val;
}

这里对于邻域极值判断需要还原到第三幅图（大尺度）的比例上，通过重载getResponse函数来实现。

inline float getResponse(unsigned int row, unsigned int column)
    {
        return responses[row * width + column];
    }
inline float getResponse(unsigned int row, unsigned int column, ResponseLayer *src)
    {
        int scale = this->width / src->width;
        return responses[(scale * row) * width + (scale * column)];
    }

2.2.4 确定极值点精确位置

这一步和SIFT基本一致。26个邻域的极值点检测，得到的值是“近似”的，因为最大值/最小值几乎从不精确地位于像素上。它位于像素之间。但我们无法访问像素之间的数据。因此，我们必须从数学上定位子像素位置。

//! Interpolate scale-space extrema to subpixel accuracy to form an image feature.     
void FastHessian::interpolateExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)//, Ipoint* ipts_tmp, int* changed)
{
    // get the step distance between filters
    // check the middle filter is mid way between top and bottom
    int filterStep = (m->filter - b->filter);
    // Get the offsets to the actual location of the extremum
    double xi = 0, xr = 0, xc = 0;
    interpolateStep(r, c, t, m, b, &xi, &xr, &xc );

    // If point is sufficiently close to the actual extremum
	  //如果点距离极值位置足够近（<0.5）获取点的信息
    if( fabs( xi ) < 0.5f    &&    fabs( xr ) < 0.5f    &&    fabs( xc ) < 0.5f )
    {
        Ipoint ipt;
        ipt.x = static_cast<float>((c + xc) * t->step);
        ipt.y = static_cast<float>((r + xr) * t->step);
        ipt.scale = static_cast<float>((0.1333f) * (m->filter + xi * filterStep));
        ipt.laplacian = static_cast<int>(m->getLaplacian(r,c,t));
        ipts.push_back(ipt);
        // *ipts_tmp = ipt;
        // *changed = 1;
    }
}

通过泰勒二阶展开式，近似极值点精确位置，代码在interpolateStep中

对x求导并让方程等于0得到x相对于插值中心偏移量

（其中需要用到二次型的导数，推导过程见https://blog.csdn.net/weixin_43342476/article/details/121562958）：

当它在任一维度上的偏移量小于0.5时（绝对值<0.5）被认为是特征点。

2.3 新建自定义Ipoint类获得特征点

上面的过程中我们得到了特征点的位置、尺度、拉普拉斯（用来加速特征点匹配）数据，
Ipont类中定义了属性和若干方法，包括了主方向，还定义了特征点之间的欧几里得距离计算（64维描述符差的平方和）

class Ipoint {

public:

    //! Destructor
    ~Ipoint() {};

    //! Constructor
    Ipoint() : orientation(0) {};

    //! Gets the distance in descriptor space between Ipoints
	// 两个特征点之间的欧几里得距离
    float operator-(const Ipoint &rhs)
    {
        float sum=0.f;
        for(int i=0; i < 64; ++i)
            sum += (this->descriptor[i] - rhs.descriptor[i])*(this->descriptor[i] - rhs.descriptor[i]);
        return sqrt(sum);
    };

    //! Coordinates of the detected interest point
    float x, y;

    //! Detected scale
    float scale;

    //! Orientation measured anti-clockwise from +ve x-axis
    float orientation;

    //! Sign of laplacian for fast matching purposes
    int laplacian;

    //! Vector of descriptor components
    float descriptor[64];

    //! Placeholds for point motion (can be used for frame to frame motion analysis)
    float dx, dy;

    //! Used to store cluster index
    int clusterIndex;
};

如果用在图像拼接场景，还需在Ipoint类中定义特征点的匹配，扭曲变换、图像混合等函数。

void getMatches(IpVec &ipts1, IpVec &ipts2, IpPairVec &matches);
int translateCorners(IpPairVec &matches, const CvPoint src_corners[4], CvPoint dst_corners[4]);

cv::Mat getCvStitch(IplImage *src, cv::Mat warpped);
cv::Mat getWarpped(IplImage *original, cv::Mat H);
cv::Mat getCvWarpped(IpPairVec &matches, IplImage *original);
cv::Mat getWarppedReMap(IpPairVec &matches, IplImage *original);

cv::Mat findHom(IpPairVec &matches);
cv::Mat getWarppedAcc(IplImage *original, cv::Mat H);
std::vector<cv::Mat> getWarpped_blend(IplImage *original, cv::Mat H);
cv::Mat getHomography(IpPairVec &matches, IplImage *original);
cv::Mat getBlended(IplImage *img1, IplImage *img2, IpPairVec &matches, cv::Mat &warpped, cv::Mat &mask2);

接着来构造SURF描述符，未完待续。