opencv 算子 角点算法 Harris

参考https://blog.csdn.net/poem_qianmo/article/details/29356187

        https://www.cnblogs.com/polly333/p/5416172.html

 

算法流程：

1. 利用水平，竖直差分算子对图像的每个像素进行滤波以求得Ix,Iy,进而求得M中的四个元素的值。

算法源码：

代码中如果array为-1,0,1,-1,0,1,-1,0,1}则是求解X方向的，如果为{-1,-1,-1,0,0,0,1,1,1}为Y方向的，则Ix和Iy求解结束

 

//C代码
//这里的size/2是为了不把图像边界算进去。
//array也就是3*3的窗口函数为：double dx[9]={-1,0,1,-1,0,1,-1,0,1};或者 //定义垂直方向差分算子并求Iy
    double dy[9]={-1,-1,-1,0,0,0,1,1,1};

CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *array,int size1,int size2)//size 
{
    int i,j;
    int i1,j1;
    int px,py;
    int m;
    CvMat *mat1;
    mat1=cvCloneMat(mat);
//为了去除边界，从框体一半开始遍历
    for(i=size1/2;i<ywidth-size1/2;i++)
        for(j=size2/2;j<xwidth-size2/2;j++)                         
        {
            m=0;
            for(i1=0;i1<size1;i1++)
                for(j1=0;j1<size2;j1++)
                {
                    px=i-size1/2+i1;
                    py=j-size2/2+j1;
                    //CV_MAT_ELEM访问矩阵元素
                  //每个元素和窗体函数遍历相加
                    m+=CV_MAT_ELEM(*mat,double,px,py)*array[i1*size1+j1];            
                }
               //赋值
                CV_MAT_ELEM(*mat1,double,i,j)=m;
        }
        return mat1;
}

求解IX2相对比较简单，像素相乘即可。下面为基于opencv的C++版本，很是简单

//构建模板
Mat xKernel = (Mat_<double>(1,3) << -1, 0, 1);
Mat yKernel = xKernel.t();
//计算IX和IY
Mat Ix,Iy;
//可参考filter2d的定义
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);
//计算其平方
Mat Ix2,Iy2,Ixy;
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);

2.对M的四个元素进行高斯平滑滤波，为的是消除一些不必要的孤立点和凸起，得到新的矩阵M。

//本例中使用5×5的高斯模板
    //计算模板参数
    //int gausswidth=5;
    //double sigma=0.8;
    double *g=new double[gausswidth*gausswidth];
    for(i=0;i<gausswidth;i++)//定义模板
        for(j=0;j<gausswidth;j++)
            g[i*gausswidth+j]=exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma));

    //归一化：使模板参数之和为1（其实此步可以省略）
    double total=0;
    for(i=0;i<gausswidth*gausswidth;i++)
        total+=g[i];
    for(i=0;i<gausswidth;i++)
        for(j=0;j<gausswidth;j++)
            g[i*gausswidth+j]/=total;

    //进行高斯平滑
    mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);
    mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);
    mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);

利用opencv接口则是很简单

//构建高斯函数
Mat gaussKernel = getGaussianKernel(7, 1);
//卷积计算
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);

3.接下来利用M计算对应每个像素的角点响应函数R，即：

也可以使用改进的R：
R=[Ix^2*Iy^2-(Ix*Iy)^2]/(Ix^2+Iy^2);里面没有随意给定的参数k，取值应当比第一个令人满意。

//计算cim：即cornerness of image，我们把它称做‘角点量’
    CvMat *mat_cim;
    mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB);
//用来求得响应度
CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,int ywidth)
{
    int i,j;
    CvMat *mat;
    mat=cvCloneMat(mat1);
    for(i = 0; i <ywidth; i++)
    {
        for(j = 0; j < xwidth; j++)
        {
            //注意：要在分母中加入一个极小量以防止除数为零溢出
            CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)-
                CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/
                (CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001);
        }
    }
    return mat;
}

//opencv代码
Mat cornerStrength(gray.size(), gray.type());
    for (int i = 0; i < gray.rows; i++)
    {
        for (int j = 0; j < gray.cols; j++)
        {
            double det_m = Ix2.at<double>(i,j) * Iy2.at<double>(i,j) - Ixy.at<double>(i,j) * Ixy.at<double>(i,j);
            double trace_m = Ix2.at<double>(i,j) + Iy2.at<double>(i,j);
            cornerStrength.at<double>(i,j) = det_m - alpha * trace_m *trace_m;
        }
    }

4、局部极大值抑制，同时选取其极大值

//--------------------------------------------------------------------------
    //                 第四步：进行局部非极大值抑制
    //--------------------------------------------------------------------------
    CvMat *mat_locmax;
    //const int size=7;
    mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);
//     cout<<CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)<<endl;
//用来求得局部极大值
CvMat *mblocmax(CvMat *mat1,int xwidth,int ywidth,int size)
{
    int i,j;
    double max=-1000;
    int i1,j1;
    int px,py;
    CvMat *mat;
    mat=cvCloneMat(mat1);
    for(i=size/2;i<ywidth-size/2;i++)
        for(j=size/2;j<xwidth-size/2;j++)
        {
            max=-10000;
            for(i1=0;i1<size;i1++)
                for(j1=0;j1<size;j1++)
                {
                    px=i-size/2+i1;
                    py=j-size/2+j1;
                    if(CV_MAT_ELEM(*mat1,double,px,py)>max)
                        max=CV_MAT_ELEM(*mat1,double,px,py);

                }
                if(max>0)
                    CV_MAT_ELEM(*mat,double,i,j)=max;
                else
                    CV_MAT_ELEM(*mat,double,i,j)=0;
        }
        return mat;
}

在opencv中应用maxminloc函数

// threshold
double maxStrength;
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat dilated;
Mat localMax;
//默认3*3核膨胀，膨胀之后，除了局部最大值点和原来相同，其它非局部最大值点被  
 //3*3邻域内的最大值点取代
dilate(cornerStrength, dilated, Mat());
//与原图相比，只剩下和原图值相同的点，这些点都是局部最大值点，保存到localMax 
compare(cornerStrength, dilated, localMax, CMP_EQ);

5.在矩阵R中，同时满足R(i,j)大于一定阈值threshold和R(i,j)是某领域内的局部极大值，则被认为是角点。

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);
imgDst = cornerMap.clone();

for( int y = 0; y < cornerMap.rows; y++ ) 
{  
        const uchar* rowPtr = cornerMap.ptr<uchar>(y);  
        for( int x = 0; x < cornerMap.cols; x++ ) 
        {  
               // 非零点就是角点  
            if (rowPtr[x]) 
            {  
                        points.push_back(cv::Point(x,y));  
            }  
        }  
 }

 

基于opencv源码

1 #include "imgFeat.h"
 2 
 3 void feat::detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha)
 4 {
 5     Mat gray;
 6     if (imgSrc.channels() == 3)
 7     {
 8         cvtColor(imgSrc, gray, CV_BGR2GRAY);
 9     }
10     else
11     {
12         gray = imgSrc.clone();
13     }
14     gray.convertTo(gray, CV_64F);
15 
16     Mat xKernel = (Mat_<double>(1,3) << -1, 0, 1);
17     Mat yKernel = xKernel.t();
18 
19     Mat Ix,Iy;
20     filter2D(gray, Ix, CV_64F, xKernel);
21     filter2D(gray, Iy, CV_64F, yKernel);
22 
23     Mat Ix2,Iy2,Ixy;
24     Ix2 = Ix.mul(Ix);
25     Iy2 = Iy.mul(Iy);
26     Ixy = Ix.mul(Iy);
27 
28     Mat gaussKernel = getGaussianKernel(7, 1);
29     filter2D(Ix2, Ix2, CV_64F, gaussKernel);
30     filter2D(Iy2, Iy2, CV_64F, gaussKernel);
31     filter2D(Ixy, Ixy, CV_64F, gaussKernel);
32     
33 
34     Mat cornerStrength(gray.size(), gray.type());
35     for (int i = 0; i < gray.rows; i++)
36     {
37         for (int j = 0; j < gray.cols; j++)
38         {
39             double det_m = Ix2.at<double>(i,j) * Iy2.at<double>(i,j) - Ixy.at<double>(i,j) * Ixy.at<double>(i,j);
40             double trace_m = Ix2.at<double>(i,j) + Iy2.at<double>(i,j);
41             cornerStrength.at<double>(i,j) = det_m - alpha * trace_m *trace_m;
42         }
43     }
44     // threshold
45     double maxStrength;
46     minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
47     Mat dilated;
48     Mat localMax;
49     dilate(cornerStrength, dilated, Mat());
50     compare(cornerStrength, dilated, localMax, CMP_EQ);
51     
52 
53     Mat cornerMap;
54     double qualityLevel = 0.01;
55     double thresh = qualityLevel * maxStrength;
56     cornerMap = cornerStrength > thresh;
57     bitwise_and(cornerMap, localMax, cornerMap);
58     
59     imgDst = cornerMap.clone();
60     
61 }
62 
63 void feat::drawCornerOnImage(Mat& image, const Mat&binary)
64 {
65     Mat_<uchar>::const_iterator it = binary.begin<uchar>();
66     Mat_<uchar>::const_iterator itd = binary.end<uchar>();
67     for (int i = 0; it != itd; it++, i++)
68     {
69         if (*it)
70             circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
71     }
72 }

opencv+C++推荐

opencv中Harris接口

OpenCV的Hairrs角点检测的函数为cornerHairrs()，但是它的输出是一幅浮点值图像，浮点值越高，表明越可能是特征角点，我们需要对图像进行阈值化。

C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);

• src – 输入的单通道8-bit或浮点图像。
• dst – 存储着Harris角点响应的图像矩阵，大小与输入图像大小相同，是一个浮点型矩阵。
• blockSize – 邻域大小。
• apertureSize – 扩展的微分算子大。
• k – 响应公式中的，参数$\alpha$。
• boderType – 边界处理的类型。

 

int main()
{
  Mat image = imread("../buliding.png");
  Mat gray;
  cvtColor(image, gray, CV_BGR2GRAY);

  Mat cornerStrength;
  cornerHarris(gray, cornerStrength, 3, 3, 0.01);
  threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY);
  return 0;
}

从上面上间一幅图像我们可以看到，有很多角点都是粘连在一起的，我们下面通过加入非极大值抑制来进一步去除一些粘在一起的角点。

非极大值抑制原理是，在一个窗口内，如果有多个角点则用值最大的那个角点，其他的角点都删除，窗口大小这里我们用3*3，程序中通过图像的膨胀运算来达到检测极大值的目的，因为默认参数的膨胀运算就是用窗口内的最大值替代当前的灰度值。

int main()
{
  Mat image = imread("buliding.png");
  Mat gray;
  cvtColor(image, gray, CV_BGR2GRAY);

  Mat cornerStrength;
  cornerHarris(gray, cornerStrength, 3, 3, 0.01);

  double maxStrength;
  double minStrength;
  // 找到图像中的最大、最小值
  minMaxLoc(cornerStrength, &minStrength, &maxStrength);

  Mat dilated;
  Mat locaMax;
  // 膨胀图像，最找出图像中全部的局部最大值点
  dilate(cornerStrength, dilated, Mat());
  // compare是一个逻辑比较函数，返回两幅图像中对应点相同的二值图像
  compare(cornerStrength, dilated, locaMax, CMP_EQ);

  Mat cornerMap;
  double qualityLevel = 0.01;
  double th = qualityLevel*maxStrength; // 阈值计算
  threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY);
  cornerMap.convertTo(cornerMap, CV_8U);
  // 逐点的位运算黑色图片显示的结果
  bitwise_and(cornerMap, locaMax, cornerMap);

  drawCornerOnImage(image, cornerMap);
  namedWindow("result");
  imshow("result", image);
  waitKey();

  return 0;
}
void drawCornerOnImage(Mat& image, const Mat&binary)
{
  Mat_<uchar>::const_iterator it = binary.begin<uchar>();
  Mat_<uchar>::const_iterator itd = binary.end<uchar>();
  for (int i = 0; it != itd; it++, i++)
  {
    if (*it)
      circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
  }
}

函数备注

1、compare

功能：两个数组之间或者一个数组和一个常数之间的比较

结构：

void compare(InputArray src1, InputArray src2, OutputArray dst, int cmpop)

src1 ：第一个数组或者标量，如果是数组，必须是单通道数组。
src2 ：第二个数组或者标量，如果是数组，必须是单通道数组。
dst ：输出数组，和输入数组有同样的size和type=CV_8UC1
cmpop ：

标志指明了元素之间的对比关系
CMP_EQ src1 相等 src2.

CMP_GT src1 大于 src2.
CMP_GE src1 大于或等于 src2.
CMP_LT src1 小于 src2.
CMP_LE src1 小于或等于 src2.
CMP_NE src1 不等于 src2.

如果对比结果为true，那么输出数组对应元素的值为255，否则为0

//获取角点图
    Mat getCornerMap(double qualityLevel) {
      Mat cornerMap;
      // 根据角点响应最大值计算阈值
      thresholdvalue= qualityLevel*maxStrength;
      threshold(cornerStrength,cornerTh,
      thresholdvalue,255,cv::THRESH_BINARY);
      // 转为8-bit图
      cornerTh.convertTo(cornerMap,CV_8U);
      // 和局部最大值图与，剩下角点局部最大值图，即：完成非最大值抑制
      bitwise_and(cornerMap,localMax,cornerMap);
      return cornerMap;
    }

2、bitwise_and（位运算函数）

功能：计算两个数组或数组和常量之间与的关系

结构：

void bitwise_and(InputArray src1, InputArray src2, OutputArray dst, InputArray mask=noArray())

src1 ：第一个输入的数组或常量
src2 ：第二个输入的数组或常量
dst ：输出数组，和输入数组有同样的size和type
mask ：可选择的操作掩码，为8位单通道数组，指定了输出数组哪些元素可以被改变，哪些不可以

操作过程为：

如果为多通道数组，每个通道单独处理

 

 

 

3、filter2D

OpenCV中的内部函数filter2D可以直接用来做图像卷积操作

void filter2D(InputArray src, OutputArray dst, int ddepth, InputArray kernel, Point anchor=Point(-1,-1), double delta=0, int borderType=BORDER_DEFAULT )