//#include <iostream>
//#include <opencv.hpp>
//#include <string>
//#include "harris.h"

#ifndef HARRIS_H //先测试 预处理器变量HARRIS_H 是否被宏定义过 
#define HARRIS_H // #define HARRIS_H
                 // 程序段 1 //如果HARRIS_H没有被宏定义过，定义HARRIS_H，并编译程序段 1
                 // #endif
                 // 程序段 2 //如果x已经定义过了则编译程序段2的语句，“忽视”程序段 1。
#include "opencv2/opencv.hpp"  
  
class harris //在C++ 语言中class是定义类的关键字，C++中也可以使用struct定义类。
	//两者区别是，用class定义的类，如果数据成员或成员函数没有说明则默认为private（私有）的，而用struct定义的，默认为public（公共）的。 
{  
private:  
    cv::Mat cornerStrength;  //opencv harris函数检测结果，也就是每个像素的角点响应函数值  
    cv::Mat cornerTh; //cornerStrength阈值化的结果  
    cv::Mat localMax; //局部最大值结果  
    int neighbourhood; //邻域窗口大小  
    int aperture;//sobel边缘检测窗口大小（sobel获取各像素点x，y方向的灰度导数）  
    double k;  
    double maxStrength;//角点响应函数最大值  
    double threshold;//阈值除去响应小的值  
    int nonMaxSize;//这里采用默认的3，就是最大值抑制的邻域窗口大小  
    cv::Mat kernel;//最大值抑制的核，这里也就是膨胀用到的核  
public:  
    harris():neighbourhood(3),aperture(3),k(0.01),maxStrength(0.0),threshold(0.01),nonMaxSize(3){  
  
    };  
  
    void setLocalMaxWindowsize(int nonMaxSize){  
        this->nonMaxSize = nonMaxSize;  
    };  
  
    //计算角点响应函数以及非最大值抑制  
    void detect(const cv::Mat &image){  //const限定变量&image为不可修改的常量，且默认为文件的局部变量
            //opencv自带的角点响应函数计算函数  
            cv::cornerHarris (image,cornerStrength,neighbourhood,aperture,k);  
            double minStrength;  
            //计算最大最小响应值  
            cv::minMaxLoc (cornerStrength,&minStrength,&maxStrength);  
  
            cv::Mat dilated;  
            //默认3*3核膨胀，膨胀之后，除了局部最大值点和原来相同，其它非局部最大值点被  
            //3*3邻域内的最大值点取代  
            cv::dilate (cornerStrength,dilated,cv::Mat());  
            //与原图相比，只剩下和原图值相同的点，这些点都是局部最大值点，保存到localMax  
            cv::compare(cornerStrength,dilated,localMax,cv::CMP_EQ);  
    }  
  
    //获取角点图  
    cv::Mat getCornerMap(double qualityLevel) {  
            cv::Mat cornerMap;  
            // 根据角点响应最大值计算阈值  
            threshold= qualityLevel*maxStrength;  
            cv::threshold(cornerStrength,cornerTh,  
            threshold,255,cv::THRESH_BINARY);  
            // 转为8-bit图  
            cornerTh.convertTo(cornerMap,CV_8U);  
            // 和局部最大值图与，剩下角点局部最大值图，即：完成非最大值抑制  
            cv::bitwise_and(cornerMap,localMax,cornerMap);  
            return cornerMap;  
    }  
  
    void getCorners(std::vector<cv::Point> &points, //表示空类型，它跟int，float是同地位的，一般用在没有返回值的函数中，比如你写void main ()，主函数完了不用写return 语句，但是如果是int main ()或者是main ()，你不写return 语句它就会有warning 
            double qualityLevel) {  
            //获取角点图  
            cv::Mat cornerMap= getCornerMap(qualityLevel);  
            // 获取角点  
            getCorners(points, cornerMap);  
    }  
  
    // 遍历全图，获得角点  
    void getCorners(std::vector<cv::Point> &points,  
    const cv::Mat& cornerMap) {  
  
            for( int y = 0; y < cornerMap.rows; y++ ) {  
                    const uchar* rowPtr = cornerMap.ptr<uchar>(y);  //const定义为常量
                    for( int x = 0; x < cornerMap.cols; x++ ) {  
                    // 非零点就是角点  
                          if (rowPtr[x]) {  
                                points.push_back(cv::Point(x,y));  
                          }  
                     }  
                }  
          }  
  
    //用圈圈标记角点  
    void drawOnImage(cv::Mat &image,  
    const std::vector<cv::Point> &points,  
            cv::Scalar color= cv::Scalar(255,255,255),  
            int radius=3, int thickness=2) {  
                    std::vector<cv::Point>::const_iterator it=points.begin();  
                    while (it!=points.end()) {  
                    // 角点处画圈  
                    cv::circle(image,*it,radius,color,thickness);  
                    ++it;  
            }  
    }  
  
};  
  
#endif // HARRIS_H 
测试文件
int main()
{
   cv::Mat  image, image1 = cv::imread ("D:/13.jpg");  
   //灰度变换  
   cv::cvtColor (image1,image,CV_BGR2GRAY);  
  
  
   // 经典的harris角点方法  
   harris Harris;  
   // 计算角点  
   Harris.detect(image);  
   //获得角点  
   std::vector<cv::Point> pts;  
   Harris.getCorners(pts,0.01);  
   // 标记角点  
   Harris.drawOnImage(image,pts);  
  
   cv::namedWindow ("harris");  
   cv::imshow ("harris",image);  
   cv::waitKey (0);  
   return 0;  

}

http://blog.csdn.net/linnol/article/details/6862005
转载于:https://www.cnblogs.com/CVGNAY/p/3467966.html

from numpy import *
from scipy.ndimage.filters import gaussian_filter
from pylab import *
def get_quo(imarray,sigma=3):#sigma为方差
    imx=zeros(imarray.shape)
    imy=zeros(imarray.shape)
    ##计算在x，y方向上的方向上的导数
    gaussian_filter(imarray,(sigma,sigma),(0,1),imx)
    gaussian_filter(imarray,(sigma,sigma),(1,0),imy)
    #计算MI矩阵，----harris矩阵的分量
    wxx=gaussian_filter(imx*imx,sigma)
    wyy=gaussian_filter(imy*imy,sigma)
    wxy=gaussian_filter(imx*imy,sigma)
    #计算特征值和迹
    det=wxx*wyy-wxy**2#特征值
    trace1=wxx+wyy#迹
    quo=det/trace1#商数
    return quo
def get_harris_point(quo,min_dis,threshold=0.1):
    '''quo为商数，min_dis为两个角点之间最短的距离，threshold为阈值比'''
    conthro=threshold*quo.max()#确定阈值的具体值
    harris=1*(quo>conthro)
    print(harris.nonzero())
    #获取候选坐标
    harris_T=array(harris.nonzero()).T
    # print(harris_T)
    #获取harris的响应值
    harris_value=[quo[c[0],c[1]] for c in harris_T]
    # print(harris_value)
    # 将响应值按从小到大输出其索引值
    index=argsort(harris_value)


    #将可行点的位置保存在数组中
    allow_point=zeros(quo.shape)
    allow_point[min_dis:-min_dis,min_dis:-min_dis]=1
    #按照给出的min_dis-最短距离原则选择最佳的harris点
    filter_point=[]
    for i in index:
        if allow_point[harris_T[i,0],harris_T[i,1]]==1:
            filter_point.append(harris_T[i])
            allow_point[harris_T[i,0]-min_dis:harris_T[i,0]+min_dis,harris_T[i,1]-min_dis:harris_T[i,1]+min_dis]=0
    return filter_point

def plot_harris_point(image,filter_point):
    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filter_point],[p[0] for p in filter_point],'r*')
    axis('off')
    show()



from PIL import Image
from numpy import *
from pylab import *
import harris_fun
im=array(Image.open(r'F:\img_test\t8.jpg').convert('L'))
quo=harris_fun.get_quo(im)
print(quo)
filter_point=harris_fun.get_harris_point(quo,10)
harris_fun.plot_harris_point(im,filter_point)
结果显示：







在做图像匹配时，常需要对两幅图像中的特征点进行匹配。为了保证匹配的准确性，所选择的特征必须有其独特性，角点可以作为一种不错的特征。
那么为什么角点有其独特性呢？角点往往是两条边缘的交点，它是两条边缘方向变换的一种表示，因此其两个方向的梯度变换通常都比较大并且容易检测到。
这里我们理解一下Harris Corner 一种角点检测的算法
角点检测基本原理：
人们通常通过在一个小的窗口区域内观察点的灰度值大小来识别角点，如果往任何方向移动窗口都会引起比较大的灰度变换那么往往这就是我们要找的角点。如下图右
下面我们看一下Harris的数学公式，对于[x,y]平移[u,v]个单位后强度的变换有下式,I(x+u,y+v)是平移后的强度，I(x,y)是原图像像素。对于括号里面的值，如果是强度恒定的区域，那么它就接近于零，反之如果强度变化剧烈那么其值将非常大，所以我们期望E(u,v)很大。
其中w是窗函数，它可以是加权函数，也可以是高斯函数
利用二维泰勒展开式我们有
所以其中一阶可以近似为
于是我们可以给出Harris Corner的如下推导，其中Ix，Iy是x，y方向的Gradient模，乘以位移得到位移后的量
对于小的位移，我们可以用双线性插值方法近似：
其中M为2*2矩阵如下
在本质上我们可以把二次项看成一个椭圆函数，我们对M进行特征值分析有λ1，λ2
根据λ1，λ2的值我们可以把其分为三类：
1.λ1，λ2都很小且近似，E在所以方向接近于常数；
2.λ1>>λ2,或者λ2>>λ1, E将在某一方向上很大；
3.λ1，λ2都很大且近似，E将在所以方向上很大；
如图所示：
最后我们通过计算角点响应值R来判断其属于哪个区间
其中k一般为常数取在0.04-0.06间。
算法步骤：
1.计算图像x,y方向的梯度Ix，Iy
2.计算每个像素点的梯度平方
3.计算梯度在每个像素点的和
4.定义在每个像素点的矩阵H，也就是前面的M
5.计算每个像素的角点响应
6.设置阈值找出可能点并进行非极大值抑制
代码：
close all
clear all
I = imread('empire.jpg');
I = rgb2gray(I);
I = imresize(I,[500,300]);
imshow(I);
sigma = 1;
halfwid = sigma * 3;
[xx, yy] = meshgrid(-halfwid:halfwid, -halfwid:halfwid);
Gxy = exp(-(xx .^ 2 + yy .^ 2) / (2 * sigma ^ 2));
Gx = xx .* exp(-(xx .^ 2 + yy .^ 2) / (2 * sigma ^ 2));
Gy = yy .* exp(-(xx .^ 2 + yy .^ 2) / (2 * sigma ^ 2));
%%apply sobel in herizontal direction and vertical direction compute the
%%gradient
%fx = [-1 0 1;-1 0 1;-1 0 1];
%fy = [1 1 1;0 0 0;-1 -1 -1];
Ix = conv2(I,Gx,'same');
Iy = conv2(I,Gy,'same');
%%compute Ix2, Iy2,Ixy
Ix2 = Ix.*Ix;
Iy2 = Iy.*Iy;
Ixy = Ix.*Iy;
%%apply gaussian filter
h = fspecial('gaussian',[6,6],1);
Ix2 = conv2(Ix2,h,'same');
Iy2 = conv2(Iy2,h,'same');
Ixy = conv2(Ixy,h,'same');
height = size(I,1);
width = size(I,2);
result = zeros(height,width);
R = zeros(height,width);
Rmax = 0;
%% compute M matrix and corner response
for i = 1:height
for j =1:width
M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy(i,j)];
R(i,j) = det(M) - 0.04*(trace(M)^2);
if R(i,j)> Rmax
Rmax = R(i,j);
end
end
end
%% compare whith threshold
count = 0;
for i = 2:height-1
for j = 2:width-1
if R(i,j) > 0.01*Rmax
result(i,j) = 1;
count = count +1;
end
end
end
%non-maxima suppression
result = imdilate(result, [1 1 1; 1 0 1; 1 1 1]);
[posc,posr] = find(result == 1);
imshow(I);
hold on;
plot(posr,posc,'r.');
本文原创，转载请注明出处