Harris Corner Point Detection algorithm detailed

Harris Corner Point algorithm

Feature point detection is widely used in target matching, target tracking, three-dimensional reconstruction and other applications, in the target modeling will be the image of the target feature extraction, commonly used color, corner point, feature points, contour, texture and other characteristics. Now we begin to explain the common feature point detection, in which Harris corner detection is the basis of feature point detection, and the application of the concept of gray difference of neighboring pixel points is put forward to judge whether it is corner, edge and smooth region. Harris Angle Point Detection principle is the use of moving windows in the image to calculate the gray change value, wherein the key process includes converting to grayscale image, calculating the difference image, Gaussian smoothing, calculating local extremum, confirming corner point. first, the basic knowledge

Type of image change:

In the feature point detection, it is often suggested that the scale is constant, the rotation is invariable, and the anti-noise effect is the index to judge whether the characteristic points are stable.

Good corner point: Detection of "real" corner points in the image accurate positioning performance high repeatability detection rate noise robust computational efficiency

Types of corner points:

Based on the gray-level image method, the corner points are detected by calculating the curvature and gradient of the points, which avoids the defects of the first kind of methods, such as Moravec operator, Forstner operator, Harris operator and Susan operator.

Second, Harris Angle point principle 1. Algorithm Idea

The angle point principle derives from the perceptual judgment of the diagonal point of the person, that is, the image has obvious change in each direction. The core of the algorithm is to use the local window to move the image to determine the gray scale changes, so this window is used to calculate the gray change of the image is: [ -1,0,1;-1,0,1;-1,0,1][-1,-1,-1;0,0,0;1,1,1]. The small window that moves this feature in all directions, such as the gray scale of the area within the window in Figure 3, is considered to have encountered a corner point in the window. As shown in Figure 1, the gray level of the image in the window does not change, then there is no corner point in the window, if the window in a certain direction, the image of the gray scale has changed greatly, and in other directions has not changed, then the image inside the window may be a line segment.

2. Mathematical model

According to the algorithm thought, constructs the mathematical model, calculates the moving window the gray scale difference value.

To reduce the computational complexity, the Taylor series is used to simplify the formula:

The W function in the image above represents a window function, and the M matrix is a partial derivative matrix. For the matrix to be able to change the symmetry matrix, assuming that the use of two eigenvalues to replace, its geometric meaning is similar to the expression in the following figure. In the geometry

The properties of a pixel are judged by judging the size of the two eigenvalues in the model.

M is the covariance matrix of gradients, in order to be able to apply better programming in practice, the corner response function r is defined, and the size of R is determined by determining whether the pixel is a corner point.

R depends on the characteristic value of M, for corner points | r| large, flat area | The r| is small and the edges R is negative.

3. Algorithm Flow

1. Using horizontal, vertical difference operators to filter each pixel of the image to obtain Ix,iy, and then obtain the value of four elements in M.

Algorithm Source code:

If array is -1,0,1,-1,0,1,-1,0,1} in the code, the X-direction is solved, and if { -1,-1,-1,0,0,0,1,1,1} is y-direction, the IX and iy solve ends

//c code//The SIZE/2 here is to not count the bounds of the image.

The window function of array is 3*3: Double dx[9]={-1,0,1,-1,0,1,-1,0,1}, or//define the vertical differential operator and ask for 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);                         
        To remove the boundary, start traversing for (i=size1/2;i<ywidth-size1/2;i++) for (j=size2/2;j<xwidth-size2/2;j++) from half the frame body
            {m=0; for (i1=0;i1<size1;i1++) for (j1=0;j1<size2;j1++) {Px=i-size1/2+i
                    1;
                    py=j-size2/2+j1; Cv_mat_elem access matrix elements//each element and form function traversal add M+=cv_mat_elem (*mat,double,px,py) *array[i1*size1            
                +J1];
        }//Assigned value Cv_mat_elem (*mat1,double,i,j) =m;
} return MAT1; }

Solving IX2 is relatively simple and can be multiplied by pixels. The following is a OPENCV-based C + + version, very simple

Build Template
Mat Xkernel = (mat_<double> (1,3) <<-1, 0, 1);
Mat Ykernel = xkernel.t ();
Calculate Ix and Iy
Mat Ix,iy;
can refer to the definition of filter2d
filter2d (Gray, Ix, cv_64f, Xkernel);
Filter2d (Gray, Iy, cv_64f, Ykernel);
Calculates its square
Mat Ix2,iy2,ixy;
Ix2 = Ix.mul (Ix);
Iy2 = Iy.mul (Iy);
Ixy = Ix.mul (Iy);

2. Gaussian smoothing filters the four elements of M to eliminate some unnecessary outliers and bumps to get a new matrix M.

In this example, the use of 5x5 Gaussian template
    //calculation template parameters
    //int gausswidth=5;
    Double sigma=0.8;
    Double *g=new double[gausswidth*gausswidth];
    for (i=0;i<gausswidth;i++)//define template 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));

    Normalization: The sum of the template parameters is 1 (in fact, this step can be omitted)
    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;

    Perform Gaussian smoothing
    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);

Using the OpenCV interface is very simple

Construct the Gaussian function
Mat gausskernel = Getgaussiankernel (7, 1);
Convolution calculation
filter2d (Ix2, Ix2, cv_64f, Gausskernel);
Filter2d (Iy2, Iy2, cv_64f, Gausskernel);
Filter2d (Ixy, Ixy, cv_64f, Gausskernel);

3. Next, use M to calculate the corner response function R for each pixel, i.e.:

You can also use the improved R:
r=[ix^2*iy^2-(Ix*iy) ^2]/(ix^2+iy^2); there is no arbitrarily given parameter k, the value should be more satisfactory than the first one.

Compute CIM: cornerness of image, we call it the ' Corner amount '
    cvmat *mat_cim;
    MAT_CIM=MBCIM (MAT_IX2,MAT_IY2,MAT_IXY,CXDIB,CYDIB);
Used to obtain the response degree
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 + +)
        {
            //Note: to include a small amount in the denominator to prevent a divisor from being zero overflow
            cv_ma T_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 Code

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< ;d ouble> (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, local maximal value suppression, and select its maximum value

--------------------------------------------------------------------------//Fourth step: local non-maximal value suppression//-
    -------------------------------------------------------------------------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;
    Used to derive local maxima 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; }

Applying the Maxminloc function in OpenCV

Threshold
double maxstrength;
Minmaxloc (cornerstrength, NULL, &maxstrength, NULL, NULL);
Mat dilated;
Mat Localmax;
Default 3*3 core expansion, after expansion, in addition to the local maximum point and the original same, the other non-local maximum points are  
 //3*3 adjacent to the largest point in the neighborhood
dilate (cornerstrength, dilated, Mat ());
Compared with the original image, only the same points as the original values, these points are local maximum points, saved to Localmax 
compare (cornerstrength, dilated, Localmax, cmp_eq);

5. In the Matrix R, while satisfying R (I,J) is greater than a certain threshold threshold and R (I,J) is a local maximum value in a domain, it is considered a corner point.

Cv::mat Cornermap;  
  The threshold value is calculated based on the angular response maximum  
  threshold= qualitylevel*maxstrength;  
  Cv::threshold (Cornerstrength,cornerth,  
   threshold,255,cv::thresh_binary);  
  Convert to 8-bit diagram 
 Cornerth.convertto (cornermap,cv_8u);//and local maximum graph with, left corner local maximum graph, that is: complete non-maximum suppression  
  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; + +) 
    {  
               //Not 0 point is corner point  
              if (rowptr[x]) 
    {  
                        points.push_back (CV::P oint (x, y));}}}  
 

third, the algorithm source code