test_rotation_invariance3.py

#Copyright [2020] [catree]
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.
from __future__ import print_function
from __future__ import division
import cv2 as cv
import numpy as np
import skimage as sk
from skimage import data
from skimage.feature import corner_harris, corner_subpix, corner_peaks
from skimage.transform import warp, AffineTransform
from skimage.draw import ellipse
import matplotlib.pyplot as plt
import os
import argparse

print('OpenCV:', cv.__version__)
print('Numpy:', np.__version__)
print('scikit-image:', sk.__version__)

# https://stackoverflow.com/questions/59473193/harris-corner-detector-at-different-rotations

RADIUS = 5
HARRIS_K = 0.04
OPENCV_HARRIS_BLOCK_SIZE = 3
OPENCV_HARRIS_APERTURE_SIZE = 1
SKIMAGE_HARRIS_SIGMA = 1

def inverse(M):
    D = M[0,0]*M[1,1] - M[0,1]*M[1,0]
    if D != 0:
        D = 1/D
    else:
        D = 0

    Minv = np.ndarray((2,3))

    A11 = M[1,1] * D
    A22 = M[0,0] * D
    Minv[0,0] = A11
    Minv[0,1] = -M[0,1] * D
    Minv[1,0] = -M[1,0] * D
    Minv[1,1] = A22
    b1 = -Minv[0,0] * Minv[0,2] - Minv[0,1] * Minv[1,2]
    b2 = -Minv[1,0] * Minv[0,2] - Minv[1,1] * Minv[1,2]
    Minv[0,2] = b1
    Minv[1,2] = b2

    return Minv

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])), M)

def proj(pt, M):
    x = pt[1]
    y = pt[0]
    x_proj = M[0,0]*x + M[0,1]*y + M[0,2]
    y_proj = M[1,0]*x + M[1,1]*y + M[1,2]
    return np.array([y_proj, x_proj])

def projection(coords, M, width, height):
    coords_proj = np.empty((0,2))
    for c in coords:
        c_proj = proj(c, M)
        if c_proj[1] > RADIUS and c_proj[0] > RADIUS and c_proj[1] < width-RADIUS and c_proj[0] < height-RADIUS:
            coords_proj = np.append(coords_proj, [c_proj], axis=0)

    return coords_proj

def dist(c1, c2):
    return np.linalg.norm(c1 - c2)

def compute_repetability_rate(coords_ref, coords_cur, M, Minv, width, height, tol):
    coords_ref_proj = projection(coords_ref, M, width, height)
    coords_cur_proj = projection(coords_cur, Minv, width, height)

    correct_matches = 0
    min_size = 0;
    if coords_ref_proj.shape[0] < coords_cur_proj.shape[0]:
        min_size = coords_ref_proj.shape[0]
        for c1 in coords_ref_proj:
            for c2 in coords_cur:
                if dist(c1, c2) < tol:
                    correct_matches += 1
                    break
    else:
        min_size = coords_cur_proj.shape[0]
        for c1 in coords_cur_proj:
            for c2 in coords_ref:
                if dist(c1, c2) < tol:
                    correct_matches += 1
                    break

    ratio = 0
    if min_size > 0:
        ratio = correct_matches / min_size
    return ratio

def rotate_detect(image, angle_lim, dist_tol, results_cv, results_sk, debug):
    coords_cv_ref = np.empty((0,0))
    coords_sk_ref = np.empty((0,0))
    width = image.shape[1]
    height = image.shape[0]

    for theta in range(angle_lim[0], angle_lim[1]):
        (img_cv, M) = rotate(image, theta, (height / 2, width / 2))
        Minv = inverse(M)
        img_sk = np.copy(img_cv)

        gray_cv = cv.cvtColor(img_cv, cv.COLOR_BGR2GRAY)
        gray_sk = np.copy(gray_cv)

        # OpenCV Harris
        # response_cv = cv.cornerHarris(gray_cv.astype(np.float32), blockSize=OPENCV_HARRIS_BLOCK_SIZE, ksize=OPENCV_HARRIS_APERTURE_SIZE, k=HARRIS_K)
        response_cv = cv.cornerHarris(gray_cv, blockSize=OPENCV_HARRIS_BLOCK_SIZE, ksize=OPENCV_HARRIS_APERTURE_SIZE, k=HARRIS_K)
        coords_cv = corner_peaks(response_cv, min_distance=5)

        # subpixel detection
        # coords_cv_sub = np.copy(np.expand_dims(coords_cv, axis=1).astype(np.float32))
        coords_cv_sub = np.ndarray((coords_cv.shape[0], 1, 2), dtype=np.float32)
        for i in range(coords_cv.shape[0]):
            coords_cv_sub[i,0,0] = coords_cv[i,1]
            coords_cv_sub[i,0,1] = coords_cv[i,0]

        coords_cv_sub = cv.cornerSubPix(gray_cv, coords_cv_sub, winSize=(5,5), zeroZone=(-1,-1), criteria=(cv.TERM_CRITERIA_EPS + cv.TermCriteria_COUNT, 40, 0.001))
        coords_cv = np.ndarray((coords_cv_sub.shape[0],2), dtype=np.float32)
        for i in range(coords_cv_sub.shape[0]):
            coords_cv[i,0] = coords_cv_sub[i,0,1]
            coords_cv[i,1] = coords_cv_sub[i,0,0]

        if theta == 0:
            coords_cv_ref = np.copy(coords_cv)
        ratio = compute_repetability_rate(coords_cv_ref, coords_cv, M, Minv, width, height, dist_tol)
        results_cv.append(ratio)

        # debug
        image_debug = np.copy(image)
        if debug:
            coords_ref_proj = np.rint(projection(coords_cv_ref, M, width, height)).astype(int)
            for c in coords_ref_proj:
                cv.circle(img_cv, (c[1], c[0]), 8, (0,255,0))
            coords_cur_proj = np.rint(projection(coords_cv, Minv, width, height)).astype(int)
            for c in coords_cv_ref:
                cv.circle(image_debug, (c[1], c[0]), 8, (0,0,255))
            for c in coords_cur_proj:
                cv.circle(image_debug, (c[1], c[0]), 8, (0,255,0))

        for c in coords_cv:
            if not np.isnan(c[1]) and not np.isnan(c[0]):
                cv.circle(img_cv, (int(c[1]), int(c[0])), 8, (0,0,255))

        # skimage Harris
        response_sk = corner_harris(gray_sk, method='k', k=HARRIS_K, sigma=SKIMAGE_HARRIS_SIGMA)
        coords_sk = corner_peaks(response_sk, min_distance=5)

        # subpixel detection
        coords_sk_sub = np.ndarray((coords_sk.shape[0], 1, 2), dtype=np.float32)
        for i in range(coords_sk.shape[0]):
            coords_sk_sub[i,0,0] = coords_sk[i,1]
            coords_sk_sub[i,0,1] = coords_sk[i,0]

        coords_sk_sub = cv.cornerSubPix(gray_sk, coords_sk_sub, winSize=(5,5), zeroZone=(-1,-1), criteria=(cv.TERM_CRITERIA_EPS + cv.TermCriteria_COUNT, 40, 0.001))
        coords_sk = np.ndarray((coords_sk_sub.shape[0],2), dtype=np.float32)
        for i in range(coords_sk_sub.shape[0]):
            coords_sk[i,0] = coords_sk_sub[i,0,1]
            coords_sk[i,1] = coords_sk_sub[i,0,0]

        if theta == 0:
            coords_sk_ref = np.copy(coords_sk)
        ratio = compute_repetability_rate(coords_sk_ref, coords_sk, M, Minv, width, height, dist_tol)
        results_sk.append(ratio)
        for c in coords_sk:
            if not np.isnan(c[1]) and not np.isnan(c[0]):
                cv.circle(img_sk, (int(c[1]), int(c[0])), 8, (0,255,0))

        yield (img_cv, img_sk, image_debug)

def play(video, debug, save=False, wait=10, key='q'):
    idx = 0
    directory = 'video'
    if save and not os.path.exists(directory):
        os.makedirs(directory)

    for img_cv, img_sk, img in video:
        legend_cv = 'CV' + ' block=' + str(OPENCV_HARRIS_BLOCK_SIZE) + ' aper=' + str(OPENCV_HARRIS_APERTURE_SIZE)
        legend_sk = 'SK' + ' sigma=' + str(SKIMAGE_HARRIS_SIGMA)
        img_cv = cv.putText(img_cv, legend_cv, (20,20), cv.FONT_HERSHEY_SIMPLEX, 0.7, (255,0,0))
        img_sk = cv.putText(img_sk, legend_sk, (20,20), cv.FONT_HERSHEY_SIMPLEX, 0.7, (255,0,0))
        cv.imshow('OpenCV', img_cv)
        cv.imshow('skimage', img_sk)
        if debug:
            cv.imshow('OpenCV reference', img)

        if save:
            filename = directory + '/harris_results_%04d.png' % idx
            results = cv.hconcat([img_cv, img_sk])
            cv.imwrite(filename, results)

            idx += 1

        kb = cv.waitKey(wait)
        if kb == ord(key) or kb == 27:
            return

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('--k', default=0.04, type=float, help='Harris k')
parser.add_argument('--block_size', default=3, type=int, help='OpenCV Harris block size')
parser.add_argument('--aperture_size', default=1, type=int, help='OpenCV Harris aperture size')
parser.add_argument('--sigma', default=1.0, type=float, help='scikit-image Harris sigma')
parser.add_argument('--dist_tol', default=5.0, type=float, help='Distance tolerance for repeatability rate')

args = parser.parse_args()

HARRIS_K = args.k
OPENCV_HARRIS_BLOCK_SIZE = args.block_size
OPENCV_HARRIS_APERTURE_SIZE = args.aperture_size
SKIMAGE_HARRIS_SIGMA = args.sigma

print('HARRIS_K', HARRIS_K)
print('OPENCV_HARRIS_BLOCK_SIZE', OPENCV_HARRIS_BLOCK_SIZE)
print('OPENCV_HARRIS_APERTURE_SIZE', OPENCV_HARRIS_APERTURE_SIZE)
print('SKIMAGE_HARRIS_SIGMA', SKIMAGE_HARRIS_SIGMA)

if args.input:
    original_image = cv.imread(args.input)
else:
    # Sheared checkerboard
    tform = AffineTransform(scale=(1.3, 1.1), rotation=1, shear=0.7,
                            translation=(110, 30))
    image = warp(data.checkerboard()[:90, :90], tform.inverse,
                 output_shape=(200, 310))
    # Ellipse
    rr, cc = ellipse(160, 175, 10, 100)
    image[rr, cc] = 1
    # Two squares
    image[30:80, 200:250] = 1
    image[80:130, 250:300] = 1
    image = image*255
    image = cv.cvtColor(image.astype(np.uint8), cv.COLOR_GRAY2BGR)
    original_image = np.copy(image)

results_cv = []
results_sk = []
angle_lim = [0, 180]
dist_tol = args.dist_tol
debug = False
play(rotate_detect(original_image, angle_lim, dist_tol, results_cv, results_sk, debug), debug, save=args.save, wait=10)

plt.plot(results_cv)
plt.plot(results_sk)
plt.xlabel('rotation angle in degrees')
plt.ylabel('repeatability rate')
plt.xlim(0, 180)
plt.legend(['OpenCV', 'scikit-image'], loc='best')
plt.show()