#include "opencv2/highgui/highgui.hpp"

#include "opencv2/imgproc/imgproc.hpp"

#include <iostream>

#include <stdio.h>

#include <stdlib.h>

using namespace cv;

using namespace std;

/// Global variables

Mat src, src_gray;

int thresh = 200;

int max_thresh = 255;

char* source_window = "Source image";

char* corners_window = "Corners detected";

/// Function header

void cornerHarris_demo( int, void* );

/** @function main */

int main( int argc, char** argv )

{

  /// Load source image and convert it to gray

  src = imread( argv[1], 1 );

  cvtColor( src, src_gray, CV_BGR2GRAY );

  /// Create a window and a trackbar

  namedWindow( source_window, CV_WINDOW_AUTOSIZE );

  createTrackbar( "Threshold: ", source_window, &thresh, max_thresh, cornerHarris_demo );

  imshow( source_window, src );

  cornerHarris_demo( 0, 0 );

  waitKey(0);

  return(0);

}

/** @function cornerHarris_demo */

void cornerHarris_demo( int, void* )

{

  Mat dst, dst_norm, dst_norm_scaled;

  dst = Mat::zeros( src.size(), CV_32FC1 );

  /// Detector parameters

  int blockSize = 2;

  int apertureSize = 3;

  double k = 0.04;

  /// Detecting corners

  cornerHarris( src_gray, dst, blockSize, apertureSize, k, BORDER_DEFAULT );

  /// Normalizing

  normalize( dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1, Mat() );

  convertScaleAbs( dst_norm, dst_norm_scaled );

  /// Drawing a circle around corners

  for( int j = 0; j < dst_norm.rows ; j++ )

     { for( int i = 0; i < dst_norm.cols; i++ )

          {

            if( (int) dst_norm.at<float>(j,i) > thresh )

              {

               circle( dst_norm_scaled, Point( i, j ), 5,  Scalar(0), 2, 8, 0 );

              }

          }

     }

  /// Showing the result

  namedWindow( corners_window, CV_WINDOW_AUTOSIZE );

  imshow( corners_window, dst_norm_scaled );

}

#include "opencv2/highgui/highgui.hpp"

#include "opencv2/imgproc/imgproc.hpp"

#include <cmath>

#include <iostream>

using namespace cv;   

Mat harris(Mat &im, double sigma, int thresh, int radius){

    Mat dx,dy,Ix,Iy,Ix2,Iy2,Ixy,cim;

    Sobel( im, Ix, CV_64F, 1, 0, 3);   //算法第一步。计算水平垂直差分

    Sobel( im, Iy, CV_64F, 0, 1, 3);

    int ksize = max(1, (int)(6*sigma));

    if(ksize % 2 == 0)

    	ksize++;

	GaussianBlur(Ix.mul(Ix), Ix2, Size(ksize, ksize), sigma);   //算法第二步。计算二阶高斯差分矩阵

	GaussianBlur(Iy.mul(Iy), Iy2, Size(ksize, ksize), sigma);

	GaussianBlur(Ix.mul(Iy), Ixy, Size(ksize, ksize), sigma);

	//Harris corner measure

	//cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2);

	cim = (Ix2.mul(Iy2) - Ixy.mul(Ixy)) / (Ix2+Iy2);            //算法第三步。计算响应函数，我使用了第二种

    Mat structedElement(radius, radius, CV_8U, Scalar(1));

    Mat mx,norm_cim;

    normalize( cim, norm_cim, 0, 255, NORM_MINMAX, CV_8U, Mat() );

	dilate(norm_cim, mx, structedElement);

	norm_cim = ( norm_cim == mx) & (norm_cim>thresh);           //算法第4第5步融合，非极大值抑制和阀值检測	

	return norm_cim;

}

int main( int, char** argv )

{

	Mat src,gray;

    src = imread( argv[1] );

    cvtColor( src, gray, CV_RGB2GRAY );

	Mat corners = harris(gray, 1.5, 30, 2);

	for( int j = 0; j < corners.rows ; j++ ) {

		for( int i = 0; i < corners.cols; i++ ) {

            if( corners.at<unsigned char>(j,i) > 0)

			{

				circle( gray, Point( i, j ), 3,  Scalar(0), 2, 8, 0 );

			}

		}

	}

    namedWindow("result", 1);

    imshow("result", gray);

    waitKey();

    return 0;

}