Harris-Laplace与Harris-Affine原理

yunleaf

已于 2023-06-20 07:54:27 修改

阅读量230

点赞数

分类专栏：摄影测量与遥感文章标签： c++ 计算机视觉 opencv

于 2023-06-20 07:43:11 首次发布

本文链接：https://blog.csdn.net/yunleaf/article/details/131298984

版权

摄影测量与遥感专栏收录该内容

3 篇文章 0 订阅

订阅专栏

文章详细阐述了Harris-Laplace检测器的原理，它通过引入多尺度空间增强了Harris角点检测的尺度不变性。同时，介绍了自适应尺度选择方法，以降低误检率和计算复杂度。接着，概述了Harris-Affine方法，该方法通过迭代优化达到仿射不变性，包括特征点检测、归一化参考帧生成和迭代更新过程。这两种技术旨在提高图像匹配的稳健性和准确性。

摘要由CSDN通过智能技术生成

Harris-Laplace

Harris-Affine

Harris-Laplace

Harris-Laplace [1]检测原理及实现

Harris 角点检测虽然对于光照强度、旋转角度改变具有较好的检测不变性，但是却不具有尺度不变性及仿射不变性，然后在现实生活中，两张图片中目标物体发生尺度变化，或由视点变化而引起仿射变化是非常常见的。为了获得尺度不变性，比较直观的方法就是建立多尺度空间，对于每个特征位置都有在不同尺度下的表示，那么在匹配时只要找到对应尺度空间下的特征点就可以了。所以我们只需要在经典的方法里引入多尺度空间，在原特征点空间里增加了多个其他尺度空间的特征点，这些增加的特征点对应于不同的尺度空间的图像，增加了目标尺度变化的鲁棒性，使其具有了一定程度的尺度不变性。由此的新的二阶矩 $M\left( x,y \right)$

可以表示为：

$M\left( x,y \right) =\sigma _{D}^{2}g\left( \sigma _I \right) *\left[ \begin{matrix} L_{x}^{2}\left( x,\sigma _D \right)& L_xL_y\left( x,\sigma _D \right)\\ L_xL_y\left( x,\sigma _D \right)& L_{y}^{2}\left( x,\sigma _D \right)\\\end{matrix} \right]$

新的 $M\left( x,y \right)$ 同原来相比增加了两个高斯尺度参数（ $\sigma _I$ 称为积分尺度，它是决定Harris 角点当前尺度的变量，其越大说明对于尺度越大，有放大效果， $\sigma _D$ 为微分尺度或局部尺度，它是决定角点附近微分值变化的变量，实际上可以认为是一个高斯平滑参数，滤除 Laplace 变换后形成的细小噪声点）， $\sigma _{D}^{2}$ 为归一化系数。

然而应用多尺度空间后，会增加大量的特征点，然而特征点的增加会提高误检测的概率及计算复杂度，另外尺度空间的变化并不是连续的，尺度空间变化间隔越少，那么最后匹配的尺度参数与实际尺度的误差越少。基于这个原因，论文作者提出了自适应尺度选择的 Harris 角点检测方法。

自适应尺度选择的 Harris 角点检测方法主要是找到局部的结构最显著的尺度，只需要确定一个显著尺度的特征点，而不需要要原来多个尺度来表示，那么在匹配时只需要将两副图片里目标点分别变换到其显著尺度下，进行匹配就可以了。这里的最显著尺度的衡量是根据所选择的点对应不同尺度 LoG 响应程度，选择出响应程度最大的尺度（这里的尺度实际上指的是LoG的模板半径，由 $\sigma$ 来确定，窗口半径为 $3\sigma$ ，即窗口参数）。衡量公式如下：

$LoG\left( x,\sigma _n \right) =\sigma _{n}^{2}\left| L_{xx}\left( x,\sigma _n \right) +L_{yy}\left( x,\sigma _n \right) \right|$

我们需要找到具有最显著尺度的点和其显著尺度，原文给我们介绍了两种方法，一种是精确的迭代方法，另一种权衡了精确度及计算效率的简化方法：

精确的迭代方法：

第一步：积分尺度 $\sigma _I$ 从 $\sigma _n=k^n\sigma _0\left( k=1.4,n=1,2,3,...,N \right)$ 中遍历选取，并且令 $\sigma _D=s\sigma _n\left( s=0.7 \right)$ ，对于每一个积分尺度和微分尺度，进行harris角点检测，获得初始的点集 $X$ 。

第二步：针对每个初始点 $x\in X$ ，在尺度空间 $\sigma _{n}^{\left( k+1 \right)}=t\sigma _{n}^{\left( k \right)},t\in \left[ 0.7,1.4 \right]$ 寻找的 $LoG\left( x,\sigma _n \right)$ 极值点所对应的尺度 $\sigma _{n}^{\left( k+1 \right)}$ 。如果不存在极值点，则舍弃该点，退出循环，继续评估点集X 中的下一个点。若存在LOG极值点，则在此时对应的尺度为 $\sigma _{n}^{\left( k+1 \right)}$ 的尺度层按8领域寻找最大的Harris响应极值点。

第三步：当该点的尺度和空间位置都不再变化，即 $\sigma _{n}^{\left( k+1 \right)}=\sigma _{n}^{\left( k \right)},x^{\left( k+1 \right)}=x^{\left( k \right)}$ 时就可以停止搜索，否则转到第二步继续迭代。

简化方法：首先同精确迭代算法第一步一样，建立 Harris 函数的多尺度空间，然后计算每个尺度里的局部 cornerness 角程度最大值，作为初始点；然后验证每一个初始点在LoG在尺度空间是否是极值点，不是极值点则删除。

简化方法的相关opencv代码如下

void HarrisLaplaceFeatureDetector_Impl::detect(InputArray img, std::vector<KeyPoint>& keypoints, InputArray msk )

{

    Mat image = img.getMat();

    if( image.empty() )

    {

        keypoints.clear();

        return;

    }

    Mat mask = msk.getMat();

    if( !mask.empty() )

    {

        CV_Assert(mask.type() == CV_8UC1);

        CV_Assert(mask.size == image.size);

    }

    Mat_<float> dx2, dy2, dxy;

    Mat Lx, Ly;

    float si, sd;

    int gsize;

    Mat fimage;

    image.convertTo(fimage, CV_32F, 1.f/255);

    /*Build gaussian pyramid*/高斯尺度金字塔

    Pyramid pyr(fimage, numOctaves, num_layers, 1, -1, true);

    keypoints = std::vector<KeyPoint> (0);

//opencv中的算法和论文中的描述有点区别，opencv是在尺度空间中提取Harris点，原文和wiki中，是直接构造不同的积分尺度进行检索。

    /*Find Harris corners on each layer*/在尺度空间进行harris角点检测

    for (int octave = 0; octave <= numOctaves; octave++)

    {

        for (int layer = 1; layer <= num_layers; layer++)

        {

            if (octave == 0)

                layer = num_layers;



            Mat Lxm2smooth, Lxmysmooth, Lym2smooth;



            si = powf(2.f, layer / (float) num_layers);//积分尺度，组内的layer尺度

            sd = si * 0.7f;//微分尺度



            Mat curr_layer;

            if (num_layers == 4)

            {

                if (layer == 1)

                {

                    Mat tmp = pyr.getLayer(octave - 1, num_layers - 1);

                    resize(tmp, curr_layer, Size(0, 0), 0.5, 0.5, INTER_AREA);



                } else

                    curr_layer = pyr.getLayer(octave, layer - 2);

            } else /*if num_layer==2*/

            {



                curr_layer = pyr.getLayer(octave, layer - 1);

            }



            /*Calculates second moment matrix二阶矩*/



            /*Derivatives*/

            Sobel(curr_layer, Lx, CV_32F, 1, 0, 1);

            Sobel(curr_layer, Ly, CV_32F, 0, 1, 1);



            /*Normalization*/微分尺度越大，导致图像越来越模糊，因此梯度越来越小，这里可以消除这一部分影响，在DOG和LOG中，也有类似操作

            Lx = Lx * sd;

            Ly = Ly * sd;



            Mat Lxm2 = Lx.mul(Lx);

            Mat Lym2 = Ly.mul(Ly);

            Mat Lxmy = Lx.mul(Ly);



            gsize = int(ceil(si * 3)) * 2 + 1;



            /*Convolution*/积分尺度

            GaussianBlur(Lxm2, Lxm2smooth, Size(gsize, gsize), si, si, BORDER_REPLICATE);

            GaussianBlur(Lym2, Lym2smooth, Size(gsize, gsize), si, si, BORDER_REPLICATE);

            GaussianBlur(Lxmy, Lxmysmooth, Size(gsize, gsize), si, si, BORDER_REPLICATE);



            Mat cornern_mat(curr_layer.size(), CV_32F);



            /*Calculates cornerness in each pixel of the image*/计算角点响应值

            for (int row = 0; row < curr_layer.rows; row++)

            {

                for (int col = 0; col < curr_layer.cols; col++)

                {

                    float dx2f = Lxm2smooth.at<float> (row, col);

                    float dy2f = Lym2smooth.at<float> (row, col);

                    float dxyf = Lxmysmooth.at<float> (row, col);

                    float det = dx2f * dy2f - dxyf * dxyf;

                    float tr = dx2f + dy2f;

                    float cornerness = det - (0.04f * tr * tr);

                    cornern_mat.at<float> (row, col) = cornerness;

                }

            }



            double maxVal = 0;

            Mat corn_dilate;



            /*Find max cornerness value and rejects all corners that are lower than a threshold*/

            minMaxLoc(cornern_mat, 0, &maxVal, 0, 0);

            threshold(cornern_mat, cornern_mat, maxVal * corn_thresh, 0, THRESH_TOZERO);

            dilate(cornern_mat, corn_dilate, Mat());//非极大值抑制，只有cornern_mat, corn_dilate相等且不为0的位置，才是检测到的Harris角点。



            Size imgsize = curr_layer.size();

//验证harris点对应的DOG响应是否为尺度极值

            /*Verify for each of the initial points whether the DoG attains a maximum at the scale of the point*/

            Mat prevDOG, curDOG, succDOG;

            prevDOG = pyr.getDOGLayer(octave, layer - 1);

            curDOG = pyr.getDOGLayer(octave, layer);

            succDOG = pyr.getDOGLayer(octave, layer + 1);



            for (int y = 1; y < imgsize.height - 1; y++)

            {

                for (int x = 1; x < imgsize.width - 1; x++)

                {

                    float val = cornern_mat.at<float> (y, x);

                    if (val != 0 && val == corn_dilate.at<float> (y, x))

                    {

//DOG确定尺度，值比较了相邻尺度，并没有比较8领域

                        float curVal = curDOG.at<float> (y, x);

                        float prevVal =  prevDOG.at<float> (y, x);

                        float succVal = succDOG.at<float> (y, x);

//将在尺度空间检测到的点换算到原始分辨率s

                        KeyPoint kp(

                                Point2f(x * powf(2.0f, (float) octave - 1) + powf(2.0f, (float) octave - 1) / 2,

                                        y * powf(2.0f, (float) octave - 1) + powf(2.0f, (float) octave - 1) / 2),

                                3 * powf(2.0f, (float) octave - 1) * si * 2, 0, val, octave);



                        if(!mask.empty() && mask.at<unsigned char>(int(kp.pt.y), int(kp.pt.x)) == 0)

                        {

                            // ignore keypoints where mask is zero

                            continue;

                        }



                        /*Check whether keypoint size is inside the image*/

                        float start_kp_x = kp.pt.x - kp.size / 2;

                        float start_kp_y = kp.pt.y - kp.size / 2;

                        float end_kp_x = start_kp_x + kp.size;

                        float end_kp_y = start_kp_y + kp.size;

//判断当前harris对应的DOG，是不是局部尺度极值

                        if (curVal > prevVal && curVal > succVal && curVal >= DOG_thresh

                                && start_kp_x > 0 && start_kp_y > 0 && end_kp_x < image.cols

                                && end_kp_y < image.rows)

                            keypoints.push_back(kp);



                    }

                }

            }



        }



    }



    /*Sort keypoints in decreasing cornerness order*/

    sort(keypoints.begin(), keypoints.end(), sort_func);

    for (size_t i = 1; i < keypoints.size(); i++)

    {

        float max_diff = powf(2, keypoints[i].octave + 1.f / 2);

//去除在不同尺度上检测到的重复点

        if (keypoints[i].response == keypoints[i - 1].response && norm(

                keypoints[i].pt - keypoints[i - 1].pt) <= max_diff)

        {



            float x = (keypoints[i].pt.x + keypoints[i - 1].pt.x) / 2;

            float y = (keypoints[i].pt.y + keypoints[i - 1].pt.y) / 2;



            keypoints[i].pt = Point2f(x, y);

            --i;

            keypoints.erase(keypoints.begin() + i);



        }

    }



    /*Select strongest keypoints*/只保留检测到的角点响应前N的harris点

    if (maxCorners > 0 && maxCorners < (int) keypoints.size())

        keypoints.resize(maxCorners);

}

Harris-Affine

Harris-Affine[2–5]的主要计算步骤如下所示：

Harris-laplace进行特征点检测，并且 $U^{\left( 0 \right)}=identity$ ， $x^{\left( 0 \right)},\sigma _{D}^{\left( 0 \right)},\sigma _{I}^{\left( 0 \right)}$ 由Harris-laplace确定。
应用前面迭代的形状自适应矩阵 $U^{\left( k-1 \right)}$ 生成归一化参考帧， $U^{\left( k-1 \right)}\mathbf{x}_{w}^{\left( k-1 \right)}=\mathbf{x}^{\left( k-1 \right)}$ 。第一次迭代中，使用 $U^{\left( 0 \right)}=identity$ 。//这里存在一个假设，即归一化后二阶矩为单位矩阵，但实际过程中归一化后，邻域范围发生了变化，因此新邻域的二阶矩并不是单位矩阵。这也是为何后面的迭代过程需要不断地求解微分和积分尺度。
应用一个类似于harris-laplace检测，选择积分尺度 $\sigma _{I}^{\left( k \right)}$ 。即选择LOG极值对应的尺度。

$\sigma _{I}^{\left( k \right)}=arg\max \sigma _{I}^{2}\det \left( L_{xx}\left( x,\sigma _I \right) +L_{yy}\left( x,\sigma _I \right) \right) \,\, \sigma _I=t\sigma _{I}^{\left( k-1 \right)}\,\, t\in \left[ 0.7,....,1.4 \right]$

在U归一化空间的积分尺度和非归一化空间时非常不一样的，因此相对于非归一化空间的积分尺度，搜索归一化空间的积分尺度是非常有必要的。

4.选择微分尺度，假设微分尺度 $\sigma _{D}^{k}=s\sigma _{I}^{k}$ ，显然，常数s应该小于1，假设 $s\in \left[ 0.5,0.75 \right]$ ，选择微分尺度，能够最大化

$Q=\frac{\lambda _{\min}\left( u \right)}{\lambda _{\max}\left( u \right)}$

即最大化

$\sigma _{D}^{\left( k \right)}=arg\max \frac{\lambda _{\min}\left( \mu \left( x_{w}^{\left( k-1 \right)},\sigma _{I}^{k},\sigma _D \right) \right)}{\lambda _{\max}\left( \mu \left( x_{w}^{\left( k-1 \right)},\sigma _{I}^{k},\sigma _D \right) \right)}\,\, \sigma _D=s\sigma _{I}^{\left( k \right)}\,\, s\in \left[ 0.5,...0.75 \right]$ $\mu \left( x_{w}^{\left( k \right)},\sigma _{I}^{k},\sigma _D \right)$ 是归一化参考帧的二阶矩。//积分和微分尺度确定后，由于在第2步中假设当前邻域的二阶矩是单位矩阵，因此在下一步求解二阶矩时的卷积核心是各向均匀的高斯核心。

5.空间定位，选择在 $x_{w}^{\left( k-1 \right)}$ 的8邻域内，harris角点响应最大的点 $x_{w}^{\left( k \right)}$

$x_{w}^{\left( k \right)}=\mathop {arg\max} \limits_{x_w\in W\left( x_{w}^{\left( k-1 \right)} \right)}\,\,\det \left( \mu \left( x_w,\sigma _{I}^{k},\sigma _{D}^{\left( k \right)} \right) -\alpha trace^2\left( \mu \left( x_w,\sigma _{I}^{k},\sigma _{D}^{\left( k \right)} \right) \right) \right)$

$\mu \left( x_w,\sigma _{I}^{k},\sigma _{D}^{\left( k \right)} \right)$ 是二阶矩， $W\left( x_{w}^{\left( k-1 \right)} \right)$ 是上一次在归一化参考帧的迭代点的8邻域。因为新的空间定位实在归一化参考帧上进行的，因此新的定位点必须变换回原始的参考帧上面。计算方法如下：

$\mathbf{x}^{\left( k \right)}=\mathbf{x}^{\left( k-1 \right)}+U^{\left( k-1 \right)}\cdot \left( \mathbf{x}_{w}^{\left( k \right)}-\mathbf{x}_{w}^{\left( k-1 \right)} \right)$

6.更新U矩阵，正如上面提到的，二阶矩的平方根定义了生成归一化参考帧的变换矩阵，因此我们需要保存矩阵， $\mu _{i}^{\left( k \right)}=\mu ^{-\frac{1}{2}}\left( \mathbf{x}_{w}^{\left( k \right)},\sigma _{I}^{\left( k \right)},\sigma _{D}^{\left( k \right)} \right)$ ，变换矩阵U被更新为 $U^{\left( k \right)}=\mu _{i}^{\left( k \right)}\cdot U^{\left( k-1 \right)}$ 。为了确保图像采样正确，我们在最小变化方向扩展图像，即最小特征值方向，我们固定最大特征值为 $\lambda _{\max}\left( U^{\left( k \right)} \right) =1$ .使用这种更新方法，能够非常简单的得到最终的U矩阵，以如下形式

$U=\prod_k{\mu _{i}^{\left( k \right)}}\cdot U^{\left( 0 \right)}=\prod_k{\left( \mu ^{-\frac{1}{2}} \right) ^{\left( k \right)}}\cdot U^{\left( 0 \right)}$

7.如果中止准则没有满足，则继续从步骤2开始下一次迭代。因为算法迭代求解U归一化矩阵，转换一个各向异性的区域到各向同性的区域，非常显然的，当 $Q=\frac{\lambda _{\min}\left( u \right)}{\lambda _{\max}\left( u \right)}$ 足够接近1，则迭代终止.中止条件如下所示

$1-\frac{\lambda _{\min}\left( \mu _{i}^{\left( k \right)} \right)}{\lambda _{\max}\left( \mu _{i}^{\left( k \right)} \right)}<\varepsilon _C$

在论文 [3]中， $\varepsilon _C=0.05$

代码参考opencv的affine2D类或者vl-feat相关代码。

References

1. Mikolajczyk, K.; Schmid, C. Indexing based on scale invariant interest points. In Eighth IEEE International Conference on Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on Computer Vision, Vancouver, BC, Canada, 7-14 July 2001; IEEE / Institute of Electrical and Electronics Engineers Incorporated, 2001; pp 525–531, ISBN 0-7695-1143-0.

2. Baumberg, A. Reliable feature matching across widely separated views. In IEEE conference on computer vision and pattern recognition. IEEE Conference on Computer Vision and Pattern Recognition. CVPR 2000, Hilton Head Island, SC, USA, 13-15 June 2000; IEEE Comput. Soc, 2000; pp 774–781, ISBN 0-7695-0662-3.

3. Mikolajczyk, K. Scale & Affine Invariant Interest Point Detectors. International Journal of Computer Vision 2004, 60, 63–86, doi:10.1023/B:VISI.0000027790.02288.f2.

4. Mikolajczyk, K.; Schmid, C. An Affine Invariant Interest Point Detector. In Computer Vision — ECCV 2002, Berlin, Heidelberg, 2002; Heyden, A., Sparr, G., Nielsen, M., Johansen, P., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2002; pp 128–142, ISBN 978-3-540-47969-7.

5. PEREZ-DANIEL, K.; ESCAMILLA, E.; NAKANO, M.; PEREZ-MEANA, H. An Improved Harris-Affine Invariant Interest Point Detector.