Harris 角点

角点

图像中物体边缘上的角点往往反映了物体局部的空间、颜色、亮度等信息，是一种非常重要的特征。

根据角点其特殊的形态，主要有两种角点检测的方法：基于梯度的方法和基于模板的方法。

其中基于模板的方法主要考虑了邻域内像素点的灰度变化，将变化足够大的点视为角点。

Harris角点

Harris特征检测器是一种经典的基于模板的角点检测方法，通过对目标像素点在各方向上移动模板来观察变化，根据其变化规律将图像分为：平坦区域，

边缘和角点。

平坦区域：在各方向上的变化都不大

边缘：在边缘方向上无变化，其余方向变化大

角点：在各方向上的变化都大

对于图像I(x,y)在点(x,y)处平移(u,v)，得到的模板内平均灰度变化为：

其中w(x,y)为窗口函数，可以为常数，也可以为高斯函数：

通过对E(u,v)进行变换：

将E(u,v)化为以下形式：

，

其中矩阵M

。

至此，E可以被写成一个椭圆函数：

其中

      。

对于一个椭圆函数来说，其两个特征值显示了其变化最快与最慢的两个方向：

其中黄色为变化最快的方向（最大的特征值），红色为变化最慢的方向（最小的特征值）。

特征值与图像中的平坦区域、边缘和角点的关系：

1. 当两个特征值都很小时，为平坦区域；
2. 当一个特征值远大于另一个特征值时，为边缘；
3. 当两个特征值都很大且相近时，为角点。

举个栗子，

为了避免计算特征值，Harris提出了一个角点响应值R：

，

其中Det，Tr分别为矩阵M的行列式和迹：

，

k为Harris系数，取0.04~0.06.

角点响应值与图像中的平坦区域、边缘和角点的关系：

1. R为正且较大时，为角点；
2. R为负且较大时，为边缘；
3. |R|较小时，为平坦区域。

当R>阈值T时，视为角点。在改进算法中，又加入了局部非极大值抑制。

Harris角点检测步骤：

1. 利用水平、垂直差分算子对图像进行滤波，对于每个像素得到Ix、Iy，进而求得M矩阵；
2. 利用高斯平滑对图像滤波；
3. 计算每个像素的角点响应值R；
4. 角点响应值R大于阈值T且满足局部最大的点被视为角点。

C++代码

实现版本

void detectHarrisCorners(const cv::Mat& imgSrc, cv::Mat& imgDst, double k){// Color space transformcv::Mat gimg;if(imgSrc.channels()==3){cvtColor(imgSrc,gimg,CV_BGR2GRAY);}else{gimg=imgSrc.clone();}// Image type transformgimg.convertTo(gimg,CV_64F);// x,y direction difference kernelcv::Mat xKernel(1,3,CV_64F);xKernel.at<double>(0,0)=-1;xKernel.at<double>(0,1)=0;xKernel.at<double>(0,2)=1;cv::Mat yKernel=xKernel.t();// Ix Iycv::Mat Ix,Iy;cv::filter2D(gimg,Ix,CV_64F,xKernel);cv::filter2D(gimg,Iy,CV_64F,yKernel);// Ix^2 Iy^2 Ixycv::Mat Ix2,Iy2,Ixy;Ix2=Ix.mul(Ix);Iy2=Iy.mul(Iy);Ixy=Ix.mul(Iy);// Gaussian weightedcv::Mat gKernel=cv::getGaussianKernel(5,1);cv::filter2D(Ix2,Ix2,CV_64F,gKernel);cv::filter2D(Iy2,Iy2,CV_64F,gKernel);cv::filter2D(Ixy,Ixy,CV_64F,gKernel);// Corner responds// M=[A  C//  C  B]// A=Ix^2// B=Iy^2// C=Ixy// Tr(M)=A+B// Det(M)=AB-C^2// R=Det-k*Tr^2cv::Mat cornerResponds(gimg.size(),gimg.type());for(int i=0;i<gimg.rows;i++){for(int j=0;j<gimg.cols;j++){double trace_M=Ix2.at<double>(i,j)+Iy2.at<double>(i,j);double det_M=Ix2.at<double>(i,j)*Iy2.at<double>(i,j)-cv::pow(Ixy.at<double>(i,j),2.0);cornerResponds.at<double>(i,j)=det_M-k*cv::pow(trace_M,2.0);}}// local maximum corner strengthdouble maxStrength;cv::minMaxLoc(cornerResponds,NULL,&maxStrength);cv::Mat dilated;cv::Mat dKernel=cv::getStructuringElement(cv::MORPH_RECT,cv::Size(5,5));dKernel.setTo(cv::Scalar(1));cv::dilate(cornerResponds,dilated,dKernel);cv::Mat localMax;cv::compare(cornerResponds,dilated,localMax,CV_CMP_EQ);// threshold & non-maximum suppressioncv::Mat cornerThre;double qualityLevel=0.01;double thre=maxStrength*qualityLevel;cv::Mat cornerMap;cornerMap=cornerResponds>thre;//cv::threshold(cornerMap,cornerMap,0,255.0,CV_THRESH_BINARY);cv::bitwise_and(cornerMap,localMax,cornerMap);imgDst=cornerMap.clone();}

使用OpenCV内置Harris角点的改进算法

// opencv refined versionclass HarrisDetection{private:    cv::Mat cornerStrength; // 角点强度图像    cv::Mat cornerThre; // 二值化后的角点强度图    int windowSize; // 方向平均强度变化临域窗口    int sobelSize; // Sobel 内核大小    double k;   // Harris 系数    double maxStrength; // 最大角点强度    cv::Mat localMax; // 局部最大强度图像    double thre; // 阈值public:    HarrisDetection():windowSize(3),sobelSize(3),k(0.05),maxStrength(0.0),thre(0.01){}    // 角点强度＋局部最大值    void detect(cv::Mat& image)    {        cv::cornerHarris(image,cornerStrength,windowSize,sobelSize,k);        double minStrength;        cv::minMaxLoc(cornerStrength,&minStrength,&maxStrength);        cv::Mat dilated;        cv::Mat kernel = cv::getStructuringElement(cv::MORPH_RECT,cv::Size(5,5));        kernel.setTo(cv::Scalar(1));        cv::dilate(cornerStrength,dilated,kernel);        cv::compare(cornerStrength,dilated,localMax,CV_CMP_EQ);    }    // 求角点    void getCorners(std::vector<cv::Point>& points,double level)    {        cv::Mat cornerMap;        cornerMap = getCornerMap(level);        getPoints(points,cornerMap);    }    // 二值化＋非最大值抑制    cv::Mat getCornerMap(double level)    {        cv::Mat cornerMap;        thre=maxStrength*level;        cv::threshold(cornerStrength,cornerThre,thre,255,CV_THRESH_BINARY);        cornerThre.convertTo(cornerMap,CV_8U);        cv::bitwise_and(localMax,cornerMap,cornerMap);        return cornerMap;    }    // 提取交点    void getPoints(std::vector<cv::Point>& points,cv::Mat& corners) {        for (int y = 0; y < corners.rows; y++)        {            const uchar* ptr= corners.ptr<uchar>(y);            for (int x=0;x<corners.cols;x++)            {                if(ptr[x])                {                    points.push_back(cv::Point(x,y));                }            }        }    }    // 绘图    void DrawOnImage(cv::Mat& image,std::vector<cv::Point>& points,cv::Scalar color = cv::Scalar(0,0,255),int radius=2,int thickness=1)    {        std::vector<cv::Point>::const_iterator it = points.begin();        while(it!=points.end())        {            cv::circle(image,*it,radius,color,thickness);            ++it;        }    }};

各版本间的效果比较

#include <opencv2/opencv.hpp>#include "Harris.h"int main(){// our versioncv::Mat imgSrc=cv::imread("building.jpg");cv::Mat imgCopy1=imgSrc.clone();cv::Mat imgCopy2=imgSrc.clone();cv::Mat cornerMap;detectHarrisCorners(imgSrc,cornerMap,0.05);drawHarrisCorners(imgSrc,cornerMap);cv::imshow("Harris Corners",imgSrc);// opencv build-in versioncv::Mat corners;cv::Mat imgG;cv::cvtColor(imgCopy1,imgG,CV_BGR2GRAY);cv::cornerHarris(imgG,corners,3,3,0.05);cv::threshold(corners,corners,0.0001,255,CV_THRESH_BINARY);corners.convertTo(corners,CV_8U);drawHarrisCorners(imgCopy1,corners);cv::imshow("opencv harris",imgCopy1);// opencv refined versionHarrisDetection h;    h.detect(imgG);    std::vector<cv::Point> points;    double level=0.01;    h.getCorners(points,level);    h.DrawOnImage(imgCopy2,points);    cv::imshow("Refined version",imgCopy2);cv::waitKey(0);return 0;}

效果展示：

a. 实现版本

b. OpenCV的Harris版本

c. OpenCV改进版本