harris角点 亚像素 matlab,亚像素级的多尺度Harris角点检测算法的制作方法

本发明属于一种图像检测方法，具体为一种针对传统Harris角点的亚像素级别的检测方法。

技术背景

在机器视觉研究中，角点检测是一个重要环节。在摄像机标定，立体匹配，三维重建等计算机视觉处理任务中起重要作用。角点作为图像的一个重要的直观局部特征，能够有效保证图像的重要信息，同时减少了信息的数据量，运算量。角点检测效果直接影响后期图像处理工作的进行，所以好的检测算法很重要。

目前角点检测算法大致分为两类：一类是通过以链码的形式提取图像的边缘，从而进行角点的寻找，另一类基于图像的灰度信息的进行角点检测。第一类算法需要对图像的边缘进行操作，对边缘线完整性的要求较高，且计算量大，目前算法有K-R(Kitchen-Rosenfeld)算法；第二类算法是直接针对灰度图像进行操作，计算图像灰度分布的曲率，以最大曲率的点作为角点，但是算法精度较差，容易受噪声的影响，此类算法较多，主要有Harris，Moravec，Susan算子，其中Harris算法是最为广泛运用的算法。

一般来说，图像中的角点的尺度信息都是未知且具有尺度上的偏差，然而在检测角点时，利用较小的尺度寻找角点时定位精度精确，但会检测到较多的伪角点。利用大尺度检测检点时，可以有效的检测真实的角点，但是角点的定位不准确。结合大尺度下角点检测的去伪和小尺度的定位精确性，可以实现角点检测的精度和准确性。

Harris算法通过微分运算和自相关矩阵来检测角点。记像素点(x，y)的灰度值为f(x，y)，像素点(x，y)移动(u，v)后灰度强度变化可表示为：

矩阵M的计算公式为：

由于矩阵M的特征值与自相关函数的极值曲率成比例，当矩阵M的两个特征值都很大时，该点的自相关函数在两个正交方向上的极值曲率都比较大，即认为该点为角点。为了避免计算M矩阵的特征值，于是Harris角点检测器被定义为：

R＝det(M)-k×trace(M)2 (3)

k值的选择一般为0.04。检测时，针对像素点为中心的窗口区域进行计算，要求中心点像素的R值为该窗口的极大值，并且需求该点的R值大于设定的阈值，即可认定该中心点为角点。

传统Harris算法虽然在角点检测效果良好，但是存在缺点：(1)角点检测依赖高斯平滑，平滑窗口过大容易造成角点的丢失，窗口过小会提取出伪角点；(2)在单一尺度角点检测时，过分依赖阈值T的设定，当角点探测器R大于T时，然而阈值T是人为设定，T的选值过大会造成角点的丢失，T的选值过小则会导致提取伪角点，只有经过多次试验尝试才能获得比较理想的角点；(3)提取角点为像素级，精度不准确。

技术实现要素：

为了解决传统Harris的问题，本发明结合Harris角点检测算子和多尺度焦点检测方法。使Harris焦点检测算法具有尺度上的变化，两次非极大值抑制，使本发明无需考虑阈值，并且加以亚像素级别的焦点检测，使角点检测的更加精确有效。

本发明提取亚像素级角点的步骤如下：

步骤一、对图片进行预筛选，具体过程如下：

步骤1-1、选择3*3区域为备选角点筛选模块，计算中心点像素灰度值与周边一点的像素灰度值之差的绝对值Δt；

步骤1-2、选择相似度阈值t，t的选值为像素灰度最大值的10％～15％；

步骤1-3、比较Δt值和t值，如果Δt值小于t值，则判定中心点与周边点相似；

步骤1-4、中心点与周边店的相似度m(中心点与相似的周边点的个数)，m为2～6时，认为中心点为备选角点；

步骤二：选择尺度空间核进行尺度变换，具体细分为：

步骤2-1、选择高斯核为尺度变换核，本发明选择尺度变换核模型为

步骤2-2、将Harris算子结合尺度空间，通过式得到Harris尺度表达式；

步骤2-3、计算带尺度变化的自相关矩阵M，计算得到M矩阵为：

其中，Lu(x，σD)和Lv(x，σD)分别表示L(x，σ)在x和y方向上的导数，σ1＝σn是计算特征点时所选择的尺度参数，σD＝sσn是微分尺度。

步骤三：选择合适的角点响应函数，本发明选择R＝det(M)/(trace(M)+ε)为角点响应函数，避免了传统Harris角点响应函数中人为选择k值带着来的误差；

步骤四：进行二次非极大值抑制寻找角点，本发明针对步骤一得到的备选角点进行非极大值抑制，选择一个10*10模板区域，计算出模板区域内的极大值，该极大值点为模板区域的中心点，而此极大值就是候补备选角点，再在此基础上进行第二次非极大值抑制，经过两次非极大值抑制，可以得到较大局部响应极大值点，并且本文认为该极大值点为候选角点。同时避免了角点响应函数计算时阈值的选择。

步骤五：对角点进行Laplace检测，确定角点；

针对在步骤四位置空间上的候选角点，对其进行尺度空间上的搜索。本发明采用Laplace算子验证角点是否在尺度维上获得Laplace算子的局部极大值。如果是，则确定为特征点，并将获得极大值的点所在的尺度作为特征尺度。通过对不同的尺度空间的Laplace响应值进行比较，获得特征点，计算过程如下：

F(x，σn)＞F(x，σm)，m∈{n-1，n+1} (5)

其中：

步骤六：进行高斯曲面拟合，得到角点的亚像素级坐标。

步骤6-1、选择高斯曲面函数，由于高斯函数具有可分离性，所以高斯曲面函数为：

步骤6-2、由于高斯曲面函数是非线性的，直接对上式进行线性拟合很困难，所以本发明针对高斯曲面函数进行双边取对数得到：

通过双边取对数，函数具有线性，所以上式同样可以用函数形式表示：

ln z＝a+bx+cy+d(x2+y2)

步骤6-3、其中x，y为像素点的坐标值，根据最小二乘法，只要求：

Q＝min∑(a+bx+cy+d(x2+y2)-ln z)2

步骤6-4、由于极值条件，可知分别对其进行求偏导得到：

步骤6-5、利用上式可以表示为AX＝B，其中X＝(a，B，c，d)T，线性求得：X＝A-1B计算得到a，b，c，d，最终得到u0，u1，(u0，u1)即为极值点的位置，也是角点亚像素级所在的位置。

综上所述本发明的算法具有很好的实用性，稳定性，精确性。实现Harris角点提取的多尺度亚像素级的交点提取。

1、通过在Harris角点提取算法的基础上添加尺度变化的理论，解决了Harris角点提取过程中单尺度的缺点。

2、采取两次非极大值抑制算法，有效地解决了Harris角点提取过程中人为设定阈值这一环节。

3、使用高斯曲面拟合技术，实现了角点的亚像素提取，提高了角点提取的精确性。

附图说明

图1是本发明的流程图。

图2是多尺度Harris角点的提取流程图。

图3是算法实验结果图。

具体实施方法

下面是结合附图对本发明进一步说明。

本发明提取亚像素级角点的步骤如下：

步骤一、对图片进行预筛选，具体过程如下：

步骤1-1、选择3*3区域为备选角点筛选模块，计算中心点像素灰度值与周边一点的像素灰度值之差的绝对值Δt；

步骤1-2、选择相似度阈值t，t的选值为像素灰度最大值的10％～15％；

步骤1-3、比较Δt值和t值，如果Δt值小于t值，则判定中心点与周边点相似；

步骤1-4、中心点与周边店的相似度m(中心点与相似的周边点的个数)，m为2～6时，认为中心点为备选角点；

步骤二：选择尺度空间核进行尺度变换，具体细分为：

步骤2-1、选择高斯核为尺度变换核，本发明选择尺度变换核模型为

步骤2-2、将Harris算子结合尺度空间，通过式得到Harris尺度表达式；

步骤2-3、计算带尺度变化的自相关矩阵M，计算得到M矩阵为：

其中，Lu(x，σD)＝I(x)*Gu(x，σD)

Lv(x，σD)＝I(x)*Gv(x，σD)

Lu(x，σD)和Lv(x，σD)分别表示L(x，σ)在u和V方向上的导数。

σ1＝σn是计算特征点时所选择的尺度参数，σD＝sσn是微分尺度。

步骤三：选择合适的角点响应函数，本发明选择R＝det(M)/(trace(M)+ε)为角点响应函数，避免了传统Harris角点响应函数中人为选择k值带着来的误差。

其中，det(M)＝AB-C2；trace(M)＝A+B；

步骤四：进行二次非极大值抑制寻找角点，避免了角点响应函数计算时阈值的选择。

步骤4-1、选择非极大值抑制的模板大小，本发明选择大小为10*10的模板；

步骤4-2、针对步骤一得到的备选角点进行第一次非极大值抑制，计算出模板区域内的极大值，该极大值点为模板区域的中心点，而此极大值就是候补备选角点；

步骤4-3、进行第二次非极大值抑制，得到较大局部响应极大值点，对备选角点的筛选，得到候选角点。

步骤五：对候选角点进行Laplace检测，确定角点；

针对在步骤四位置空间上的候选角点，对其进行尺度空间上的搜索。本发明采用Laplace算子验证角点是否在尺度维上获得Laplace算子的局部极大值。如果是，则确定为特征点，并将获得极大值的点所在的尺度作为特征尺度。通过对不同的尺度空间的Laplace响应值进行比较，获得特征点，计算过程如下：

F(x，σn)＞F(x，σm)，m∈{n-1，n+1} (7)

其中：

步骤六：进行高斯曲面拟合，得到角点的亚像素级坐标。

步骤6-1、选择高斯曲面函数，由于高斯函数具有可分离性，所以高斯曲面函数为：

步骤6-2、由于高斯曲面函数是非线性的，直接对上式进行线性拟合很困难，所以本发明针对高斯曲面函数进行双边取对数得到：

通过双边取对数，函数具有线性，所以上式同样可以用函数形式表示：

ln z＝a+bx+cy+d(x2+y2)

步骤6-3、其中x，y为像素点的坐标值，根据最小二乘法，只要求：

Q＝min∑(a+bx+cy+d(x2+y2)-ln z)2

步骤6-4、由于极值条件，可知分别对其进行求偏导得到：

步骤6-5、利用上式可以表示为AX＝B，其中X＝(a，b，c，d)T，线性求解：X＝A-1B。计算得到a，b，c，d，最终得到u0，u1，(u0，u1)即为极值点的位置，也是角点亚像素级所在的位置。