【Pythonでコーナー検出】ヘシアンコーナー検出器の理論と実装

前回はコーナー検出の概要の勉強と雑実装を行いました。 ayatalog.hatenablog.jp 今回はヘシアンコーナー検出器の理論と実装について勉強したいと思います。
↓目次

1. ヘシアンコーナー検出器の概要
2. ヘシアン行列式の演算部の実装
3. 極大値の計算
4. ソース全体
まとめ

1. ヘシアンコーナー検出器の概要

ヘシアンコーナー検出器では、画像をx,y平面 + 輝度の3次元空間としてとらえたとき、コーナーは曲率が大きくなる点と考えます。
試しに画像をx,y平面 + 輝度の3次元プロットすると以下のようになります。
3次元グラフの書き方は次のサイトを参考にしました。 Python 3：3次元グラフの書き方（matplotlib, pyplot, mplot3d, MPL） - Qiita

f:id:Ayata:20210323175200p:plain:w300 — 元画像

f:id:Ayata:20210323175303p:plain:w300 — 3次元プロット

0か255の2値しかないので極端ですが、確かにコーナー周辺では曲率が大きく変化しているように見えます。
曲率の大きさはガウス曲率の大きさで判定します。
ガウス曲率は以下の式で表されます。
f:id:Ayata:20210323143552p:plain:w250

式中のHはヘッセ行列と呼ばれ、以下の式になります。
f:id:Ayata:20210323143601p:plain:w300

また、det(H)はヘシアンの行列式と呼ばれ、次式で定義されます。
f:id:Ayata:20210323143613p:plain:w250

ヘシアンコーナー検出器では、ガウス曲率の分子であるヘシアンの行列式det(H)が極大となる点をコーナーとして検出します。
ガウス曲率の求め方は以下サイトが参考になります。
画像処理 (5) コーナー検出

2. ヘシアン行列式の演算部の実装

まずヘシアンの行列式の演算部を作っていきます。
離散系の二次微分はエッジ処理のところで勉強したので、これをベースに実装していきます。
ayatalog.hatenablog.jp Ixx, Iyy, Ixyはそれぞれ以下の式で表されます。
f:id:Ayata:20210323155239p:plain:w550 Ixyはこのままフィルタを作ると偏るため、Δxが-1と1の場合の平均を取って以下の式としました。
f:id:Ayata:20210323164548p:plain:w550 これを元にフィルタを作ると以下のようになります。
f:id:Ayata:20210323164828p:plain:w550 このフィルタを用いて微分を演算し、そのあとにヘシアン行列式detHを求める流れで実装していきます。
フィルタおよびヘシアン行列式の演算部は以下のようになりました。

#フィルタ
#Ixx
fil_xx = np.array([[0, 0, 0],
                  [1, -2, 1],
                  [0, 0, 0]])

#Iyy
fil_yy = np.array([[0, 1, 0],
                  [0, -2, 0],
                  [0, 1, 0]])

#Ixy
fil_xy = np.array([[1, -1, 0],
                  [-1, 2, -1],
                  [0, -1, 1]])

#空間フィルタリングの演算部
N = 1
Ixx = np.zeros((height - N, width - N))
Iyy = np.zeros((height - N, width - N))
Ixy = np.zeros((height - N, width - N))

#フィルタの演算
for x in range(N, width - N, 1):
    for y in range(N, height - N, 1):
        xx = 0
        yy = 0
        xy = 0
        for i in range(-N, N+1, 1):
            for j in range(-N, N+1, 1):
                xx = xx + img_gray[y + j, x + i] * fil_xx[j + N, i + N]
                yy = yy + img_gray[y + j, x + i] * fil_yy[j + N, i + N]
                xy = xy + img_gray[y + j, x + i] * fil_xy[j + N, i + N]
        Ixx[y - 1, x - 1] = xx
        Iyy[y - 1, x - 1] = yy
        Ixy[y - 1, x - 1] = xy / 2

#ヘシアンの演算
detH = np.zeros((height - N, width - N))
for x in range(width - N):
    for y in range(height - N):
        detH[x, y] = abs(Ixx[y, x] * Iyy[y, x] - Ixy[x, y] * Ixy[x, y])

3. 極大値の計算

ヘシアン行列式の演算部ができたため、続いて極大値の計算部を作っていきます。
自身(x,y)と周辺{(x-1,y), (x+1,y), (x,y-1), (x,y+1), (x-1, y-1), (x+1, y+1), (x+1, y-1), (x-1, y+1)}を比較し、すべてに対して大きければ極大値と判定するものとしました。
すごくごり押しの判定ですが一応できました。もう少し頭のいい方法を知りたいです。
コーナーとして判定された位置を白で置き換えています。

#極大値を求める
img_corner = np.zeros((height - N-1, width - N-1))
for x in range(1, width - N -1):
    for y in range(1, height - N-1):
        if detH[x, y] > detH[x-1,y] and detH[x, y] > detH[x+1,y]:
            if detH[x, y] > detH[x,y-1] and detH[x, y] > detH[x,y+1]:
                if detH[x, y] > detH[x-1, y - 1] and detH[x, y] > detH[x+1, y + 1]:
                    if detH[x, y] > detH[x+1, y - 1] and detH[x, y] > detH[x-1, y + 1]:
                        img_corner[x-1,y-1] = 255

コーナーが求められ、出力は以下のようになりました。
f:id:Ayata:20210323175348p:plain:w350

4. ソース全体

ソースコード全体を載せておきます。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import cv2 #OpenCVのインポート
import numpy as np

#画像パス
FILE_PATH = 'C:/.../aaa.png'

#画像の読み込み
img = cv2.imread(FILE_PATH)

#画像サイズの取得
height, width, color = img.shape
print("width = " + str(height))
print("height = " + str(height))
print("color channel = " + str(color))

#グレースケールへの変換
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

#2次元メッシュの作成
x = np.arange(0, width, 1)
y = np.arange(0, height, 1)
X,Y = np.meshgrid(x,y)
Z = img_gray[X, Y,]

#3次元プロット
fig = plt.figure()
ax = Axes3D(fig)

ax.plot_wireframe(X,Y,Z)
plt.show()

#フィルタ
#Ixx
fil_xx = np.array([[0, 0, 0],
                  [1, -2, 1],
                  [0, 0, 0]])

#Iyy
fil_yy = np.array([[0, 1, 0],
                  [0, -2, 0],
                  [0, 1, 0]])

#Ixy
fil_xy = np.array([[1, -1, 0],
                  [-1, 2, -1],
                  [0, -1, 1]])

#空間フィルタリングの演算部
N = 1
Ixx = np.zeros((height - N, width - N))
Iyy = np.zeros((height - N, width - N))
Ixy = np.zeros((height - N, width - N))

#フィルタの演算
for x in range(N, width - N, 1):
    for y in range(N, height - N, 1):
        xx = 0
        yy = 0
        xy = 0
        for i in range(-N, N+1, 1):
            for j in range(-N, N+1, 1):
                xx = xx + img_gray[y + j, x + i] * fil_xx[j + N, i + N]
                yy = yy + img_gray[y + j, x + i] * fil_yy[j + N, i + N]
                xy = xy + img_gray[y + j, x + i] * fil_xy[j + N, i + N]
        Ixx[y - 1, x - 1] = xx
        Iyy[y - 1, x - 1] = yy
        Ixy[y - 1, x - 1] = xy / 2

#ヘシアンの演算
detH = np.zeros((height - N, width - N))
for x in range(width - N):
    for y in range(height - N):
        detH[x, y] = abs(Ixx[y, x] * Iyy[y, x] - Ixy[x, y] * Ixy[x, y])

#極大値を求める
img_corner = np.zeros((height - N-1, width - N-1))
for x in range(1, width - N -1):
    for y in range(1, height - N-1):
        if detH[x, y] > detH[x-1,y] and detH[x, y] > detH[x+1,y]:
            if detH[x, y] > detH[x,y-1] and detH[x, y] > detH[x,y+1]:
                if detH[x, y] > detH[x-1, y - 1] and detH[x, y] > detH[x+1, y + 1]:
                    if detH[x, y] > detH[x+1, y - 1] and detH[x, y] > detH[x-1, y + 1]:
                        print(x-1, y-1)
                        img_corner[x-1,y-1] = 255

#画像のプロット
cv2.namedWindow("img", cv2.WINDOW_AUTOSIZE)
cv2.namedWindow("img1", cv2.WINDOW_AUTOSIZE)
cv2.namedWindow("img2", cv2.WINDOW_AUTOSIZE)
cv2.imshow("img",img)
cv2.imshow("img1",detH)
cv2.imshow("img2",img_corner)

cv2.waitKey(0)
cv2.destroyAllWindows()

#出力の保存
cv2.imwrite('img.png',img)
cv2.imwrite('detH.png',detH)
cv2.imwrite('img_corner_hessian.png',img_corner)