This is a tutorial on how to create a Harris detector in Python!

We start by importing OpenCV and numpy.
        
import cv2 as cv
import numpy as np
        
    
Next, we create the function we will call in order to get our covariance matrix. The covariance matrix is the same size as the original image and contains values that will later be used to compute determinants and traces in windows around pixels, in order to get a score on whether the current pixel has "features" (corners).

The function takes an image and the k value used to scale the square of the trace, as seen further below. Two variables are defined: window size (size of the window around the pixel that is used to check for features), sobel aperture size (the size of the Sobel filter kernel). These are hardcoded in this case.
        
def corner_harris(img, k):
    window_size = 3
    sobel_aperture_size = 3
    rows, cols = img.shape

    scale = np.float32((1 << (sobel_aperture_size - 1)) * window_size)
    scale *= 255.0
    scale = 1.0 / scale

    img_gradient_x = cv.Sobel(img, cv.CV_32F, 1, 0, ksize=sobel_aperture_size, scale=scale)
    img_gradient_y = cv.Sobel(img, cv.CV_32F, 0, 1, ksize=sobel_aperture_size, scale=scale)

    covariance_matrix = np.zeros((rows, cols, 3), np.float32)

    for x in range(rows):
        for y in range(cols):
            dx = float(img_gradient_x[x, y])
            dy = float(img_gradient_y[x, y])
            
            covariance_matrix[x, y, 0] = dx*dx
            covariance_matrix[x, y, 1] = dx*dy
            covariance_matrix[x, y, 2] = dy*dy

    covariance_matrix = cv.boxFilter(covariance_matrix, cv.CV_32F, (window_size, window_size))

    return calc_harris(covariance_matrix, k)
        
    
The scale is used in the Sobel filters in order to get nice values in the img_gradients. The Sobel filter is used for edge detection. Gradients of the image is taken in x and y directions in two separate matrices (img_gradient_x and img_gradient_y).

A covariance matrix is made. It is also called structure tensor.

A box-filter is applied to the covariance matrix. If this step is not done, the determinant operation done with the values in the covariance matrix will always be 0! This adds weights to all values in the covariance matrix. This is done with a convolution kernel. In this case the kernel will have the same size as the window we detect features (corners) with.

In the return statement we call another function that returns a finished matrix with the scores of whether or not the current pixel (and window) has a feature (corner).

This is that function:
        
def calc_harris(covariance_matrix, k):
    rows, cols, channels = covariance_matrix.shape

    out_matrix = np.zeros((rows, cols), np.float32)

    for x in range(rows):
        for y in range(cols):
            a = covariance_matrix[x, y, 0]
            b = covariance_matrix[x, y, 1]
            c = covariance_matrix[x, y, 2]

            h_det = a*c - b*b
            h_trace = a + c
            r = h_det - k * h_trace * h_trace
            
            out_matrix[x, y] = r

    return out_matrix
        
    
Then we call the function and mark the places where features are detected.
        
img_raw = cv.imread('pic2_small.jpg')
img_raw = cv.cvtColor(img_raw, cv.COLOR_BGR2GRAY)

rows, cols = img_raw.shape

cv.imshow('Original', img_raw)

img_r = corner_harris(img_raw, 0.06)
img_r = img_r * 81 # No idea why...

max_r = np.amax(img_r)
min_r = np.amin(img_r)

img_markers = np.array(img_raw)
for i in range(rows):
    for j in range(cols):
        if img_r[i, j] > 0.0005:
            cv.circle(img_markers, center=(j, i), radius=9, color=(255, 0, 0), thickness=1)

cv.imshow('Detected features', img_markers)
        
    
This is the result on an image.