不同旋转的哈里斯角检测器

Question

我正在 python 中使用哈里斯角检测的 opencv 实现。我的问题是关于下面 gif 中显示的行为 - 当图像旋转时，角停止被检测到（在各种旋转中）。完整代码：

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()

基于this。

使用的图片可以在here找到。从 gif 中可能不太清楚（如果您运行代码会更清楚），但是当网格线水平/垂直时会检测到角。

如果增加blocksize 参数，我们可以获得所需的行为（检测所有旋转的角点）。

问题 - 如何解释 gif 中显示的行为？

Answer 1

更新：

使用 goodFeaturesToTrack() 而不是使用手动代码从 Harris 响应图中检索角点。

如果您想使用 OpenCV 获得正确的结果，您必须调整 blockSize 和 qualityLevel 以适应您的用例。

左边是 DIPlib，右边是 OpenCV 结果。

---

DIPlib 代码：

from __future__ import print_function
from timeit import default_timer as timer
import PyDIP as dip
import numpy as np
import imageio
import os

img = imageio.imread('https://i.stack.imgur.com/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 代码：

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)

---

旧：

算法说明

This paper分析Harris角点检测器的实现。

@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},

}

算法如下：

第 1 步未包含在 original paper 中。

---

OpenCV Harris 角点检测实现

• 与some specific conditions，英特尔IPP库用于提取哈里斯角
• ippiHarrisCorner_8u32f_C1R() documentation
• 否则，OpenCV的实现大致是：
  • 使用 Sobel 计算图像梯度
  • 计算梯度的乘积
  • 盒式过滤器
  • 计算角点响应

---

Scikit-image 角点检测实现

• 文档是here
• 代码是here

代码很简单：

    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

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

和_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 角点检测规范

可以找到here。

根据这些规范，供应商可以提供针对其平台优化的自定义实现。例如，ARM compute library。

---

OpenCV 与 Scikit-image 的比较。

• 使用的版本：

OpenCV: 4.2.0-dev
Numpy: 1.18.1
scikit-image: 0.16.2

• 哈里斯参数k=0.04
• OpenCV Harris函数参数修改，实验默认值为：blockSize=3和ksize=1
• 修改的Scikit-image Harris函数参数，实验默认值：sigma=1
• 左：OpenCV 结果；右图：Scikit-image 结果

数独图片

• 默认参数：

• OpenCV：blockSize=7，apertureSize=3； Scikit-image：sigma=5：

可以计算重复率，但不确定我的代码是否完全正确：

距离阈值为 5。

从我的实验来看，Scikit-image 的角点定位精度似乎更高。

Blox 图片

• 默认参数：

• OpenCV：blockSize=3，apertureSize=3； Scikit-image:sigma=3:

OP 数独图片

• 默认参数：

• OpenCV：blockSize=7，apertureSize=3； Scikit-image:sigma=7:

校准器图像

• 默认参数：

---

代码是here。

Answer 2

这个答案只是为了表明哈里斯角检测器的正确实现应该是完全旋转不变的。

如动画图像所示，有一些检测（在背景中和密集簇内）以随机角度出现或消失。这些是由不完美的插值和不可避免的浮点舍入误差引起的。但大多数检测在所有角度都是一致的。 [注意：GIF 太大，无法上传到 SO，我在这里显示帧，请随时运行脚本来查看动画。]

使用的实现是DIPlib 中的实现。它使用高斯梯度来计算导数，并使用高斯低通滤波进行局部平均。 The Gaussian is a perfectly isotropic filter。该实现始终使用单精度浮点运算，以最大限度地减少舍入误差。你可以找到它here。 DIPlib 旨在精确量化，因此试图避免不必要的数值错误。 [披露：我是作者。]

import PyDIP as dip
import numpy as np
import imageio

img = imageio.imread('https://i.stack.imgur.com/hkbFD.jpg')
img = dip.Image(img)
img.SetColorSpace('rgb')
img = dip.ColorSpaceManager.Convert(img,'gray')

animation = []

for angle in range(0,180,4):
   img_rot = dip.Rotation2D(img, angle/180*3.14159, '3-cubic', 'add zeros')
   img_rot = img_rot.Pad(np.maximum(img.Sizes(), img_rot.Sizes()))
   img_rot.Crop(img.Sizes())
   harris = dip.HarrisCornerDetector(img_rot, sigmas=[3.0])
   harris *= dip.Maxima(harris, connectivity=2)
   harris = dip.Dilation(harris > 1, 5)
   harris = dip.Overlay(img_rot, harris)
   harris.Convert('UINT8')
   animation.append(harris)

imageio.mimsave('./so.gif', animation, fps=4) # This file is too large for SO

imageio.imsave('./so_00.png', animation[0]) # Frame at 0 degrees
imageio.imsave('./so_20.png', animation[5]) # Frame at 20 degrees

Answer 3

首先了解哈里斯算法很重要。角点是其局部邻域位于两个主要和不同边缘方向的点。换言之，角点可以理解为两条边的交汇点，其中边是图像亮度的突然变化。

为了找到边缘，Harris 算法获取图像并应用高斯滤波器以消除任何噪声。然后它应用 Sobel 算子来找到灰度图像中每个像素的 x 和 y 梯度值。对于灰度图像中的每个像素p，考虑围绕它的3×3 窗口并计算角强度函数。将此称为 Harris 值。它找到所有超过某个阈值并且是某个窗口内的局部最大值的像素。最后对每个满足5中条件的像素，计算一个特征描述符。

如您所见，Harris 算法高度依赖图像中的X 和Y 梯度。是否可以在 GIF 图片中看到，当您的图片是水平或垂直时，网格表会生成许多高度梯度值。