Search code examples
pythonopencvcomputer-visioncorner-detection

Harris corner detector at different rotations


I am using the opencv implementation of Harris Corner detection in python. My question is regarding the behaviour shown in the gif below - as the image is rotated the corners stop being detected (at various rotations). The full code:

import cv2

image_path = 'image1.jpg'

original_image = cv2.imread(image_path)

def play(video, name='video', wait=60, key='q'):
    for f in video:
        cv2.imshow(name, f)
        if cv2.waitKey(wait) == ord(key):
            return

def rotate(image, theta, point=(0,0)):
    M = cv2.getRotationMatrix2D((point[1], point[0]), theta, 1)
    return cv2.warpAffine(image, M, (image.shape[1], image.shape[0]))

def rotate_detect(image):
    for theta in range(0, 360):
        img = rotate(image, theta, (original_image.shape[0] / 2, original_image.shape[1] / 2))

        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

        dst = cv2.cornerHarris(gray,9,13,0.04)

        #result is dilated for marking the corners, not important
        dst = cv2.dilate(dst,None)

        # Threshold for an optimal value, it may vary depending on the image.
        threshold = 0.005
        img[dst>threshold*dst.max()]=[0,0,255]
        yield img

play(rotate_detect(original_image), wait=60)

cv2.destroyAllWindows()

based on this.

enter image description here

The image used can be found here. It is perhaps not so clear from the gif (more clear if you run the code), but the corners are detected when the grid lines are horizontal/vertical.

If the blocksize parameter is increased, we can get the desired behaviour (to detect corners at all rotations).

Question - How can the behaviour shown in the gif be explained?


Solution

  • Updated:

    Use goodFeaturesToTrack() instead of using the manual code to retrieve the corners from the Harris response map.

    You will have to adapt the blockSize and qualityLevel to your use case if you want to have correct results with OpenCV.

    Left is DIPlib, right is OpenCV results.

    enter image description here


    DIPlib code:

    from __future__ import print_function
    from timeit import default_timer as timer
    import diplib as dip
    import numpy as np
    import imageio
    import os
    
    img = imageio.imread('https://i.sstatic.net/hkbFD.jpg')
    img = dip.Image(img)
    img.SetColorSpace('rgb')
    img = dip.ColorSpaceManager.Convert(img,'gray')
    
    animation = []
    times = []
    
    for angle in range(0,360):
        img_rot = dip.Rotation2D(img, np.radians(angle), '3-cubic', 'add zeros')
        img_rot = img_rot.Pad(np.maximum(img.Sizes(), img_rot.Sizes()))
        img_rot.Crop(img.Sizes())
    
        start = timer()
        harris = dip.HarrisCornerDetector(img_rot, sigmas=[3.0])
        end = timer()
        times.append(end - start)
        harris *= dip.Maxima(harris, connectivity=2)
        harris = dip.Dilation(harris > 1, 5)
        harris = dip.Overlay(img_rot, harris)
        harris.Convert('UINT8')
        animation.append(harris)
    
    print('Mean computation time: {:f}s'.format(np.mean(times)))
    print('Median computation time: {:f}s'.format(np.median(times)))
    print('Std computation time: {:f}s'.format(np.std(times)))
    
    save_folder = 'DIPlib'
    if not os.path.exists(save_folder):
        os.mkdir(save_folder)
    for idx, img in enumerate(animation):
        imageio.imsave('{}/Harris_DIPlib_{:03d}.png'.format(save_folder, idx), img)
    

    OpenCV code:

    from __future__ import print_function
    from __future__ import division
    from timeit import default_timer as timer
    import argparse
    import numpy as np
    import cv2 as cv
    import os
    
    parser = argparse.ArgumentParser(description='Test Harris corners rotation invariance.')
    parser.add_argument('--input', default='', type=str, help='Input image path')
    parser.add_argument('--save', default=False, type=bool, help='Save results')
    parser.add_argument('--maxCorners', default=500, type=int, help='Maximum number of corners')
    parser.add_argument('--qualityLevel', default=0.03, type=float, help='Minimal accepted quality of image corners')
    parser.add_argument('--minDistance', default=11, type=int, help='Minimum possible Euclidean distance between the returned corners')
    parser.add_argument('--blockSize', default=11, type=int, help='Size of an average block for computing a derivative covariation matrix over each pixel neighborhood')
    args = parser.parse_args()
    
    harris_params = dict(maxCorners = args.maxCorners,
                         qualityLevel = args.qualityLevel,
                         minDistance = args.minDistance,
                         blockSize = args.blockSize,
                         useHarrisDetector = True)
    print('harris_params:\n', harris_params)
    
    image_path = 'hkbFD.jpg'
    original_image = cv.imread(image_path)
    
    def play(video, name='video', wait=60, key='q'):
        idx = 0
        directory = 'OpenCV'
        if args.save and not os.path.exists(directory):
            os.makedirs(directory)
    
        times = []
        for f, elapsed_time in video:
            times.append(elapsed_time)
            cv.imshow(name, f)
    
            if args.save:
                filename = directory + '/Harris_OpenCV_%03d.png' % idx
                cv.imwrite(filename, f)
                idx += 1
    
            if cv.waitKey(wait) == ord(key):
                return
    
        print('Mean computation time: {:f}s'.format(np.mean(times)))
        print('Median computation time: {:f}s'.format(np.median(times)))
        print('Std computation time: {:f}s'.format(np.std(times)))
    
    def rotate(image, theta, point=(0,0)):
        M = cv.getRotationMatrix2D((point[1], point[0]), theta, 1)
        return cv.warpAffine(image, M, (image.shape[1], image.shape[0]))
    
    def rotate_detect(image):
        for theta in range(0, 360):
            img = rotate(image, -theta, (original_image.shape[0] / 2, original_image.shape[1] / 2))
    
            gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
    
            start = timer()
            harris_cv = cv.goodFeaturesToTrack(gray, **harris_params)
            elapsed_time = timer() - start
    
            for c in harris_cv:
                cv.circle(img, (int(c[0,0]), int(c[0,1])), 8, (0,0,255))
    
            yield (img, elapsed_time)
    
    play(rotate_detect(original_image), wait=60)
    

    Old:

    Description of the algorithm

    This paper analyses the implementation of the Harris corner detector.

    @article{ipol.2018.229,

    title   = {{An Analysis and Implementation of the Harris Corner Detector}},
    author  = {Sánchez, Javier and Monzón, Nelson and Salgado, Agustín},
    journal = {{Image Processing On Line}},
    volume  = {8},
    pages   = {305--328},
    year    = {2018},
    doi     = {10.5201/ipol.2018.229},
    

    }

    The algorithm is the following:

    enter image description here

    Step 1 is not included in the original paper.


    OpenCV Harris corner detection implementation


    Scikit-image corner detection implementation

    The code is straightforward:

        Axx, Axy, Ayy = structure_tensor(image, sigma)
    
        # determinant
        detA = Axx * Ayy - Axy ** 2
        # trace
        traceA = Axx + Ayy
    
        if method == 'k':
            response = detA - k * traceA ** 2
        else:
            response = 2 * detA / (traceA + eps)
    
        return response
    

    with structure_tensor():

        image = _prepare_grayscale_input_2D(image)
    
        imx, imy = _compute_derivatives(image, mode=mode, cval=cval)
    
        # structure tensore
        Axx = ndi.gaussian_filter(imx * imx, sigma, mode=mode, cval=cval)
        Axy = ndi.gaussian_filter(imx * imy, sigma, mode=mode, cval=cval)
        Ayy = ndi.gaussian_filter(imy * imy, sigma, mode=mode, cval=cval)
    
        return Axx, Axy, Ayy
    

    and _compute_derivatives():

        imy = ndi.sobel(image, axis=0, mode=mode, cval=cval)
        imx = ndi.sobel(image, axis=1, mode=mode, cval=cval)
    
        return imx, imy
    

    OpenVX Harris corner detection specifications

    It can be found here.

    From these specifications, vendors can ship custom implementation optimized for their platform. For instance, the ARM compute library.


    Comparison between OpenCV and Scikit-image.

    • Versions used:
    OpenCV: 4.2.0-dev
    Numpy: 1.18.1
    scikit-image: 0.16.2
    
    • Harris parameter k=0.04
    • OpenCV Harris function parameters that are modified, default values for the experiments are: blockSize=3 and ksize=1
    • Scikit-image Harris function parameter that is modified, default value for the experiment: sigma=1
    • Left: OpenCV results; Right: Scikit-image results

    Sudoku image

    • default parameters:

    enter image description here

    • OpenCV: blockSize=7, apertureSize=3 ; Scikit-image: sigma=5:

    enter image description here

    Repeatability rate can be computed, but not sure that my code is totally correct:

    enter image description here

    Distance threshold is 5.

    From my experiments, it looks like Scikit-image yields better corner location accuracy.

    Blox image

    • default parameters:

    enter image description here

    • OpenCV: blockSize=3, apertureSize=3 ; Scikit-image: sigma=3:

    enter image description here

    OP Sudoku image

    • default parameters:

    enter image description here

    • OpenCV: blockSize=7, apertureSize=3 ; Scikit-image: sigma=7:

    enter image description here

    Calibrator image

    • default parameters:

    enter image description here


    Code is here.