目的

SIFT (Scale Invariant Feature Transform) を Python で numpy を用いて実装することを目標に、画像処理の基礎である、画像中の物体のエッジとコーナーを検出する代表的な Moravec と Harris のアルゴリズムを習得する。本記事は [1] と [2] を主に参考として書かれている。

使用する環境

ソフトウェア名	バージョン
Python	3.4 or 3.5
Pillow	3.1.0
Numpy	1.10
Scipy	0.16.1
Matplotlib	1.5.0

注意

本記事は、Morning Project Samurai 2016/9/17 ミーティングのために書かれたものである。
本記事は今後も加筆修正されていく。
訂正したほうが良い事項などありましたら、ご指摘いただけると幸いです。

Moravec's corner detector (MCD)

MCD では、まず元画像を右、右下、下、左下にずらし、元画像に重ねることを考える。次に、元画像中に小さなウィンドウを考え、そのウィンドウ内において、元の画像とずらした後の画像のピクセルごとの輝度差の2乗を合計する。上記4パターンの動かし方の中で得た4つの合計値の最小値がある基準以上であれば、そのウィンドウ内に物体の角がある、と判定する。

ここまでの動きを図で描くと、下図のようになる。下図では、左上から右に向かって、元画像、右、右下にずらされた画像、左下から右に向かって、下、左下にずらされた画像となっている。各方向へのズレは1画素である。赤枠の格子と青枠の格子はそれぞれウィンドウを表している。それぞれのセルの大きさは1画素である。青枠の格子は元画像上で赤枠の格子に対応する位置にあるものとする。差分は、赤枠の格子中で黒 (または白) だった部分が青枠の格子中で白 (または黒) に変化した画素の個数である。閾値を3と設定しておけば、下図の状況において赤く塗りつぶされた部分はコーナーとして認識される。

式で表せば次式のようになる。

E_{x, y} = \sum_{u, v} w_{u, v} |I_{x+u, y+v} - I_{u, v}|^{2}

ここで、 $I$ は元画像、$u, v$ は元画像のピクセル、$x, y$ は元画像をずらす方向を表し $(1, 0), (1, 1), (0, 1), (-1, 1)$ のいずれかである。$w$ はウィンドウを表し、その内部で $1$、外部で $0$ をとる元画像と同サイズの長方形の binary フィルタである。上図では、赤枠の格子で表される部分が 1 となる部分である。

今、我々は元画像 $I$ の全ての画素 $(s, t)$ に対して上式を適用したい。つまり、ある画素 $(s, t)$ 周りで $1$ となるような binary フィルタ $w^{s, t}$ を用意し、$(s, t)$ がコーナーであるか否かを判定するという作業を、すべての $(s, t) \in I$ に対して適用する。

E_{x, y}^{s, t} = \sum_{u, v} w_{u, v}^{s, t} |I_{x+u, y+v} - I_{u, v}|^{2}

上式は、$u, v$ に関する convolution を用いて下記のように簡潔に書くことができる。

E_{x, y}^{s, t} = w_{u, v}^{s, t} * |I_{x+u, y+v} - I_{u, v}|^{2}

下記は MCD の Python コードである。

from copy import deepcopy
import numpy as np
from scipy.misc import imread, imsave
from scipy.signal import convolve2d


def moravec(orig_img, window_shape, threshold):
    def extend_img(_orig_img):
        _extended_img = np.zeros(shape=(_orig_img.shape[0] + 2, _orig_img.shape[1] + 2))
        _extended_img[1:-1, 1:-1] = _orig_img
        return _extended_img

    def ssd(_orig_img, _extended_img, _moving_direction, window):
        row, col = _moving_direction
        _diff = deepcopy(_extended_img)
        _diff[1 + row:_orig_img.shape[0] + 1 + row, 1 + col:_orig_img.shape[1] + 1 + col] -= _orig_img
        return convolve2d(np.power(_diff[1:-1, 1:-1], 2), window, boundary='symm', mode='same')

    def corner_map(_ssds, _threshold, window_shape):
        _map = np.zeros(shape=_ssds[0].shape)
        for i in range(window_shape[0], _ssds[0].shape[0] - window_shape[0]):
            for j in range(window_shape[1], _ssds[0].shape[1] - window_shape[1]):
                _map[i, j] = 1 if np.min(_ssds[:, i, j]) > _threshold else 0
        return _map

    moving_directions = [(1, 0), (1, 1), (0, 1), (-1, 1)]
    extended_img = extend_img(orig_img)
    ssds = np.array([ssd(orig_img, extended_img, md, np.ones(shape=window_shape)) for md in moving_directions])
    return corner_map(ssds, threshold, window_shape)


if __name__ == '__main__':
    img = imread('img/kizuna.png', flatten=True)/255
    corners = moravec(img, (5, 5), 0.8)
    imsave('img/kizuna_moravec.png', 1 - (img + corners) / np.max(img + corners))

下は上記オペレータを画像に適用した結果である。黒い点の部分がコーナーとして認識された場所である。

(出典: http://getnews.jp/img/archives/imp/and_157244.jpg)

MDC には、(1) パラメータの値によって精度が変わる、(2) 画像のずれが $\pi/4$ 刻であるので、$\pi n/4$ ($n$ は自然数) 以外の傾きを持つ直線上に謝ってコーナーを認識してしまう、(3) ノイズに弱い、という問題点が存在する。 (2) の例として、$\pi /4$ の傾きを持つ直線が描かれた画像と、$\pi /4$ よりも小さな傾きを持つ直線が描かれた画像にたいして MDC を適用した結果を、それぞれ順に示す。今度は、白く描かれた部分がコーナーとして認識された部分であることに注意する。

($\pi /4$ の傾きを持つ直線に対する結果)

($\pi /4$ よりも小さな傾きを持つ直線に対する結果)

下記は上記結果を再現するためのコードである。

from copy import deepcopy
import numpy as np
from scipy.misc import imread, imsave
from scipy.signal import convolve2d


def moravec(orig_img, window_shape, threshold):
    def extend_img(_orig_img):
        _extended_img = np.zeros(shape=(_orig_img.shape[0] + 2, _orig_img.shape[1] + 2))
        _extended_img[1:-1, 1:-1] = _orig_img
        return _extended_img

    def ssd(_orig_img, _extended_img, _moving_direction, window):
        row, col = _moving_direction
        _diff = deepcopy(_extended_img)
        _diff[1 + row:_orig_img.shape[0] + 1 + row, 1 + col:_orig_img.shape[1] + 1 + col] -= _orig_img
        return convolve2d(np.power(_diff[1:-1, 1:-1], 2), window, boundary='symm', mode='same')

    def corner_map(_ssds, _threshold, window_shape):
        _map = np.zeros(shape=_ssds[0].shape)
        for i in range(window_shape[0], _ssds[0].shape[0] - window_shape[0]):
            for j in range(window_shape[1], _ssds[0].shape[1] - window_shape[1]):
                _map[i, j] = 1 if np.min(_ssds[:, i, j]) > _threshold else 0
        return _map

    moving_directions = [(1, 0), (1, 1), (0, 1), (-1, 1)]
    extended_img = extend_img(orig_img)
    ssds = np.array([ssd(orig_img, extended_img, md, np.ones(shape=window_shape)) for md in moving_directions])
    return corner_map(ssds, threshold, window_shape)


if __name__ == '__main__':
    img = np.zeros(shape=(256, 256))
    for i in range(256):
        img[i, i] = 1.0

    corners = moravec(img, (3, 3), 0.8)
    imsave('img/pi_4.png', corners)

    img = np.zeros(shape=(256, 256))
    row = 0
    for j in range(256):
        if j % 3 == 0:
            row += 1
        img[row, j] = 1.0

    corners = moravec(img, (5, 5), 0.8)
    imsave('img/not_pi_4.png', corners)

Harris' Corner Detector (HCD)

HCD は解析的な手法を用いて MCD を改良したものである。まず、画像を連続信号と捉え、1次の Teylor 展開を用いて、$E^{s, t}$ を次式で近似する。ここで $(s, t)$ は注目画素である。

\begin{eqnarray}
E^{s, t}(x, y) &=& \sum_{u, v} w^{s, t}_{u, v} |I(x+u, y+v) - I(u, v)|^{2}\\
&\approx& \sum_{u, v} w^{s, t}_{u, v} |\frac{\partial{I}}{\partial{x}}x + \frac{\partial{I}}{\partial{y}}y|^{2}\\
&=& \sum_{u, v} w^{s, t}_{u, v} \{(\frac{\partial{I}}{\partial{x}}x)^{2} + (\frac{\partial{I}}{\partial{y}}y)^{2} + (\frac{\partial{I}}{\partial{x}}x)(\frac{\partial{I}}{\partial{y}}y)\}
\end{eqnarray}

画像の列方向の微分は画像と $(-1, 0, 1)$ というフィルタの convolution で表され、列方向の微分は $(-1, 0, 1)^{T}$ というフィルタで表されるので、上式は次式のように書き直すことができる。ここで $*$ は convolution を表す。

\begin{eqnarray}
X &=& I * (-1, 0, 1)\\
Y &=& I * (-1, 0, 1)^{T}
\end{eqnarray}

\begin{eqnarray}
E^{s, t}(x, y) &=& \sum_{u, v} w^{s, t}_{u, v} \{(X(u, v)x)^{2} + (Y(u, v)y)^{2} + (X(u, v)x)(Y(u, v)y)\}\\
&=& \sum_{u, v} w^{s, t}_{u, v} \{(X(u, v)x)^{2} + (Y(u, v)y)^{2} + X(u, v)Y(u, v)xy\}\\
\end{eqnarray}

さらに、上式は $w$ と $(X(u,v)x)^{2}+(Y(u,v)y)^{2}+(X(u,v)x)(Y(u,v)y)$ の convolution であるので、次式のように書ける。

\begin{eqnarray}
A^{s, t} &=& X^{2} * w^{s, t}\\
B^{s, t} &=& Y^{2} * w^{s, t}\\
C^{s, t} &=& (XY) * w^{s, t}
\end{eqnarray}

\begin{eqnarray}
E^{s, t}(x, y) &=& A^{s, t}x^{2} + B^{s, t}y^{2} + C^{s, t}xy \\
&=& (x, y)M^{s, t}(x, y)^{T}
\end{eqnarray}

ここで、

m = 
\begin{matrix}
A^{s, t} & C^{s, t} \\
C^{s, t} & B^{s, t}
\end{matrix}

である。

さらに、ノイズの影響を受けにくくするため、$w^{s, t}$ を長方形の binary フィルタから gaussian フィルタに変更する。 gaussian フィルタは、次のように表される。

w^{s, t}_{u, v} = exp(-\frac{(u - s)^{2} + (v - t)^{2}}{2 \sigma^{2}})

下記は、空白画像とコーナーを持つ画像の二つの画像に対して $M$ を計算し、それによって表現される2次曲面をグラフ化するコードである。

from copy import deepcopy
import numpy as np
from scipy.ndimage import convolve, gaussian_filter
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt


class Harris:
    def __init__(self, orig_img, sigma):
        diff_filter = np.array([[-1.0, 0, 1.0],])
        dIdx = convolve(orig_img, diff_filter)
        dIdy = convolve(orig_img, diff_filter.T)
        self._A = gaussian_filter(np.power(dIdx, 2), sigma=sigma)
        self._B = gaussian_filter(np.power(dIdy, 2), sigma=sigma)
        self._C = gaussian_filter(dIdx * dIdy, sigma=sigma)

    def M(self, row, col):
        return np.array([[self._A[row, col], self._C[row, col]], [self._C[row, col], self._B[row, col]]])


def calc_quadratic_form(M, start, stop, num):
    rows = np.linspace(start[0], stop[0], num[0])
    cols = np.linspace(start[1], stop[1], num[1])
    values = []
    for row in rows:
        for col in cols:
            values.append((np.array([[row, col]]) @ M @ np.array([[row, col]]).T)[0, 0])
    return rows, cols, np.array(values)


def plot(rows, cols, values):
    fig = plt.figure()
    ax = Axes3D(fig)
    X, Y = np.meshgrid(cols, rows)
    ax.plot_wireframe(X, Y, values.reshape((len(rows), len(cols))))


if __name__ == '__main__':

    # Create an empty image
    img = np.zeros(shape=(256, 256))
    harris = Harris(img, 1.0)

    # Draw a graph of the quadratic form of the empty image
    rows, cols, values = calc_quadratic_form(harris.M(10, 10), (-1, -1), (1, 1), (10, 10))
    plot(rows, cols, values)

    # Create an image having a corner
    img = np.zeros(shape=(256, 256))
    col = 0
    for i in range(254):
        img[i, col] = img[i + 1, col] = img[i + 2, col] = 1.0
        if i < 128:
            col += 1
        else:
            col -= 1
    harris = Harris(img, 1.0)

    img *= 0.5

    # Draw a graph of the quadratic form on a line
    rows, cols, values = calc_quadratic_form(harris.M(120, 120), (-1, -1), (1, 1), (10, 10))
    plot(rows, cols, values)

    plt.figure()
    line_img = deepcopy(img)
    line_img[115:125, 115:125] = 1.0
    plt.imshow(line_img, cmap='Greys_r')

    # Draw a graph of the quadratic form on the corner
    rows, cols, values = calc_quadratic_form(harris.M(129, 129), (-1, -1), (1, 1), (10, 10))
    plot(rows, cols, values)

    plt.figure()
    corner_img = deepcopy(img)
    corner_img[125:135, 125:135] = 1.0
    plt.imshow(corner_img, cmap='Greys_r')

    plt.show()

まず、面 (輝度が一定) 上の点における2次形式のグラフを示す。

画像のライン上にある点 (下図1枚目の白い点の部分) における2次形式のグラフ (下図2枚目) を示す。

画像のコーナー上にある点(下図1枚目の白い点の部分) における2次形式のグラフ (下図2枚目) を示す。

上記3つの2次形式のグラフを見れば、面上の点について $M$ で表現される2次形式は平坦であり、ライン上の点についてはライン方向には平坦でラインから外れる方向については傾斜がついている。コーナー上の点についてはどの方向を見ても傾斜がついていることがわかる。これらより、$M$ で表される2次形式の傾斜の性質を数値化することができれば、それを用いてコーナー判定ができるのではないかと考えられる。この性質を数値化したものを曲率と呼び、その中でも特に主曲率と呼ばれるものは $M$の固有値として求められ、コーナー判定にはこの主曲率を間接的に用いる。

ここからは、主曲率を数学的に定義し、$M^{s, t}$ との関係について少し詳しく見てゆく。

我々が考えている曲面は、$z = E^{s, t}(x, y) = (x, y)M^{s, t}(x, y)^T$ である。簡単のため、しばらく $(s, t)$ は固定し明記しないこととする。

$x = x(p, q), y = y(p, q), z = z(p, q)$ とパラメータ表示することを考える。そうすると例えば、次式のように表される。

\begin{array}
x &=& p\\
y &=& q\\
z &=& Ap^{2} + 2Cpq + Bq^{2}
\end{array}

上記をまとめて $P = P(p, q) = (x(p, q), y(p, q), z(p, q))$ と書く。

$p$ と $q$ に対して共通の $r$ というパラメータを導入し $p(r)$、$q(r)$ とかくことにより、曲面上の曲線を表現することができる。

曲面を表現する方法として、第I基本形式と第II基本形式と呼ばれるものがある。第I基本形式は、$p$ 方向、 $q$ 方向にそれぞれ微小変化したときの曲面上での接平面に沿った変化の大きさを表現する。第II基本形式は、$p$ 方向、 $q$ 方向にそれぞれ微小変化したときの法線方向における曲面の変化 (凹凸度合い) を表す。

我々は、まず $z = E(x, y)$ の凹凸度合いを知りたいので、第II基本形式から見てゆく。単位法線ベクトル $e(p, q)$ を下記のように定義する。

e(p, q) = \frac{P_{p} \times P_{q}}{|P_{p} \times P_{q}|}

ある点 $(p_{0}, q_{0})$ における法線ベクトルを $a = e(p_{0}, q_{0})$ を固定する。$P$ の $a$ 方向に関する (符号付き) 長さは $f(p, q) = a P(p, q)^{T}$ で表される。これを微分すれば、次式が得られる。

df(p, q) = a P_{p}^{T}(p, q) dp + a P_{q}^{T}(p, q) dq

$a P_{p}^{T}(p_{0}, q_{0}) = 0$ かつ $a P_{q}^{T}(p_{0}, q_{0}) = 0$ であるので、$(p, q) = (p_{0}, q_{0})$ において、 $df(p_{0}, q_{0}) = 0$ となる。これより、$f$ は $(p_{0}, q_{0})$ において critical value をとることがわかる。

我々の扱う曲面は今、2次のテーラー展開で十分に近似可能なものであるので、$(p_{0}, q_{0})$ における Hesse 行列 $H$ を用いてその近傍の $f$ の動きを検証することができる。これを詳しく見て行く。$H$ は次のように定義される。

H = 
\begin{matrix}
f_{pp} & f_{pq} \\
f_{pq} & f_{qq}
\end{matrix}

$f(p, q)$ を $(p_{0}, q_{0})$ 周りで $H$ を用いて2次近似したものを $(p_{0}, q_{0})$ が原点に来るように並行移動すると下記のよう書ける。

\tilde{f}(p, q) = (p, q)H(p, q)^{T}

$H$ は実対称行列であるので必ず実数の固有値が存在する。それらを $\lambda_{0}$、$\lambda_{1}$ とおき、それらに対応する単位固有ベクトルを $v_{0}$、$v_{1}$ とする。固有値　$\lambda_{0}$、$\lambda_{1}$ を対角成分に順に並べた行列を $A$、固有ベクトル $v_{0}$、$v_{1}$ を左から順に並べた行列を $V$ と置けば、次式が得られる。

\begin{array}
\tilde{f}(p, q) &=& (p, q)H(p, q)^{T}\\
&=& (p, q)V^{-1} A V(p, q)^{T}
\end{array}

$H$ が対称行列であるので $V$ は直行行列である。これを用いると次式が得られる。

\tilde{f}(p, q) = (p, q)V^{T} A V(p, q)^{T}

$(\bar{p}, \bar{q}) = (p, q)V^{T}$ と置くと、次式が得られる。

\begin{array}
\tilde{f}(\bar{p}, \bar{q}) &=& (\bar{p}, \bar{q}) A (\bar{p}, \bar{q})^{T}\\
&=& \lambda_{0} (\bar{p})^{2} + \lambda_{1} (\bar{q})^{2}
\end{array}

以上より、$\lambda_{0}$ と $\lambda_{1}$ が $0$ より大きい時、$f$ は原点から離れるに従って大きくなるので、 $P$ は $a$とは逆方向に凸である。$\lambda_{0}$ と $\lambda_{1}$ が $0$ より小さい時、$P$ は $a$ 方向に凸である。$\lambda_{0}$ と $\lambda_{1}$ の符号が一致していなければ鞍点となる。

このように $f$ を近似する Hesse 行列 $H$ を用いれば、曲面 $P$ の局所的な形についての情報が得られる。Hesse 行列は、$P$ と $a$ を用いれば、次のように表される。

H = \begin{matrix}
aP_{pp}^{T} & aP_{pq}^{T}\\
aP_{qp}^{T} & aP_{qq}^{T}
\end{matrix}

$a = e(p_{0}, q_{0})$ であったので、より一般的にするために、次式のように書き直す。

H(p_{0}, q_{0}) = \begin{matrix}
e(p_{0}, q_{0})P_{pp}^{T} & e(p_{0}, q_{0})P_{pq}^{T}\\
e(p_{0}, q_{0})P_{qp}^{T} & e(p_{0}, q_{0})P_{qq}^{T}
\end{matrix}

簡単のために、次式のように書く。

H(p_{0}, q_{0}) = \begin{matrix}
L(p_{0}, q_{0}) & M(p_{0}, q_{0})\\
M(p_{0}, q_{0}) & N(p_{0}, q_{0})
\end{matrix}

第II基本形式は、次のように定義される。

\begin{array}
II(p_{0}, q_{0}) &=& (dp, dq)H(p_{0}, q_{0})(dp, dq)^{T}\\
&=& L(p_{0}, q_{0})dpdp + 2M(p_{0}, q_{0})dpdq + N(p_{0}, q_{0})dqdq
\end{array}

上述までの議論を第II基本形式を用いて言い表すと、次のようになる。第II基本形式が正値2次形式の場合、$P$ は $a$とは逆方向に凸である。第II基本形式が負値2次形式の場合、$P$ は $a$とは逆方向に凸である。

第II基本形式は、ある点 $p_{0}, q_{0}$ まわりの法線方向の2次までの情報となっている。

ここまでで第II基本形式を用いて、ある点 $(p_{0}, q_{0})$ において曲面が凸であるか否かの判定ができるようになったが、我々は、ある点 $(p_{0}, q_{0})$ において曲面がどの程度曲がっているのかについて定量的に知りたい。この定量的な指標である主曲率を、第II基本形式と第I基本形式 (後に定義する) を用いて導いてゆく。

我々は、これまで $p$、$q$ という二つのパラメータを使って曲面を表現してきた。$p$ と $q$ に共通のパラメータ $r$ を導入すれば $(p(r), q(r))$ は $p, q$ 平面上に曲線を描き、$P(r) = P(p(r), q(r))$ は曲面上に描かれる一本の曲線となる。ここで $|dP/dr| = 1$ とする。曲面には無数にこの曲線を描くことができ、その描き方によって同じ点を通る場合でも、通過する時の曲面の曲り度合いが変化する。そこで、我々は曲面上のある点を通る任意の曲面上の曲線に対してその点における曲り度合いを表す曲率を定義し、その最大と最小のものを主曲率として定義する。

曲線に対する曲率は $P''(r)$ として定義される。

$P''(r)$ を曲面に対する接ベクトル方向 $k_{n}(r)$ とその法線方向ベクトル $k_{g}(r)$ へ分解する。法線方向ベクトルは、接ベクトル方向からどの程度の速さで離脱していくか (どの程度曲線が曲がっているか) を表す。

法線方向ベクトル $k_{n}(r)$ は法線方向の単位ベクトル $e(r)$ を用いれば、スカラー量 $\kappa_{n}$ を用いて $\kappa_{n}(r)e(r)$ として表すことができる。スカラー量 $k_{n}$ を法曲率と呼ぶ。この法曲率の最大値、最小値を主曲率と呼び、それぞれ値が得られた方向を主方向と呼ぶ。我々の目的はこの主曲率を求めることである。

$\kappa_{n}(r)$ に単位法ベクトル同士の内積 $e(r)e(r)^{T}$ を掛けてもその値は変わらないので、次式を得る。

\begin{array}
\kappa_{n}(r) &=& \kappa_{n}(r) e(s)e(r)^{T}\\
&=& k_{n}(r)e(r)^{T} \\
&=& (P''(r) - k_{g}(r))e(r)^{T}\\
&=& P''(r)e(r)^{T}
\end{array}

ここで、

(P'(r)e(r)^{T})' = P''(r)e(r)^{T} + P'(r)e'(r)^{T}

であるので、

P''(r)e(r)^{T} = (P'(r)e(r)^{T})' - P'(r)e'(r)^{T}

となり、$P'(r)e(r)^{T} = 0$ であるので、結局次式を得る。

P''(r)e(r)^{T} = - P'(r)e'(r)^{T}

これより、$\kappa_{n}(r)$ は次のようになる。

\kappa_{n}(r) = - P'(r)e'(r)^{T}

上式に chain rule を適用すれば、次式を得る。

\begin{array}
\kappa_{n}(r) &=& - P'(r)e'(r)^{T}\\
&=& -(P_{p}(p(r), q(r))p'(r) + P_{q}(p(r), q(r))q'(r))\\
&& (e_{p}(p(r), q(r))p'(r) + e{q}(p(r), q(r))q'(r))^{T}
\end{array}

さらに、上式を展開し、$L = e(p, q)P_{pp}^{T}(p, q) = - e_{p}(p, q)P_{p}^{T}(p, q)$ といった関係を用いれば (TODO: もう少し丁寧に書く) 、最終的に次式を得る。

\begin{array}
\kappa_{n}(r) &=& L(p(r), q(r))p'(r)^{2} + 2M(p(r), q(r))p'(r)q'(r) + N(p(r), q(r))q'(r)^{2}\\
&=& (p'(r), q'(r)) \begin{matrix}
L(p(r), q(r)) & M(p(r), q(r)) \\
M(p(r), q(r)) & N(p(r), q(r))
\end{matrix}
(p'(r), q'(r))^{T}
\end{array}

今 $r_{0}$ で $(p_{0}, q_{0})$ を通る任意の曲線を考え、$\kappa_{n}(r_{0})$ を考える。そうすれば上式の右辺で変化するのは $p'(r_{0})$ と $q'(r_{0})$ のみになり、これらは任意の実数をとり得る。これに基づき、$\xi = p'(r_{0})$、$\zeta = q'(r_{0})$ とおき上式を次のように書く。

\kappa_{n}(\xi, \zeta) = (\xi, \zeta) \begin{matrix}
L & M \\
M & N
\end{matrix}
(\xi, \zeta)^{T}

我々は、上式を最大化または最小化したいのだが、その前に、1つ考えておかなければならない条件がある。それは、少し前に設定した制約 $|dP/dr| = 1$ である。

$|dP/dr| = 1$ を計算すると

\begin{array}
1 &=&  \sqrt{(P_{p}p' + P_{q}q')(P_{p}p' + P_{q}q')^{T}} \\
&=& \sqrt{P_{p}P_{p}^{T}p'^{2} + 2P_{p}P_{q}^{T}p'q' + P_{q}P_{q}^{T}q'^{2}}\\
&=& \sqrt{Ep'^{2} + 2Fp'q' + Gq'^{2}}
\end{array}

である。ここで、$E=P_{p}P_{p}^{T}$、$F=P_{p}P_{q}^{T}$、$G=P_{q}P_{q}^{T}$ と置いた。さらに計算すれば、下記のように表される。

Ep'^{2} + 2Fp'q' + Gq'^{2} = 1

この左辺を第I基本形式と呼ぶ。今、$\xi = p'(r_{0})$、$\zeta = q'(r_{0})$ であるので、さらに上式を書き換えて、次式を得る。

E\xi^{2} + 2F\xi\zeta + G\zeta^{2} = (\xi, \zeta)J(\xi, \zeta)^{T} = 1

ここで、J は下記のような行列である。

J = \begin{matrix}
E & F\\
F & G
\end{matrix}

まとめると、我々は、上記制約 $E\xi^{2} + 2F\xi\zeta + G\zeta^{2} = 1$ のもとで、$\kappa_{n}(\xi, \zeta)$ の最大値、最小値を求めることとなる。

この最大値、最小値は、次式の $\lambda$ を最大値、最小値を求めることと等しい。

\lambda = \frac{Lp'^{2} + 2Mp'q' + Nq'^{2}}{Ep'^{2} + 2Fp'q' + Gq'^{2}}

次のように式を変形して、

L\xi^{2} + 2M\xi\zeta + N\zeta^{2} - \lambda (E\xi^{2} + 2F\xi\zeta + G\zeta^{2}) = 0

$\xi$ と $\zeta$ について偏微分したものを 0 としておき、それを満たす $(\xi, \zeta) \neq (0, 0)$ について $\lambda$ を計算すれば良い。$\xi$ と $\zeta$ について偏微分したものを 0 としておいたものを次に示す。

  \left\{
    \begin{array}{l}
      (L - \lambda E)\xi + (M - \lambda F)\zeta = 0 \\
      (M - \lambda F)\xi + (N - \lambda G)\zeta = 0
    \end{array}
  \right.

これは得られた2つの偏微分に関する式を連立方程式として見た時に、$(\xi, \theta) \neq (0, 0)$ であるような解が存在するような $\lambda$ を見つけることに等しい。これは、ベクトル $(L - \lambda E, M - \lambda F)$ と $(M - \lambda F, N - \lambda G)$ が1次従属であることを意味し、それは下記行列式が $0$ となることに等しい。

\left|
\begin{array}{cc}
L - \lambda E & M - \lambda F\\
M - \lambda F & N - \lambda G\\
\end{array}
\right| = 0

上式を $H$ と $J$ を用いて書き直せば、下記のようになる。

|H - \lambda J| = |(H - \lambda J)J^{-1}J| 
= |(H - \lambda J)J^{-1}J|
= |HJ^{1-} - \lambda I||J| = 0

両辺を $|J|$ で割れば次式を得る。

|HJ^{-1} - \lambda I| = 0

よって、$\lambda$ は $HJ^{-1}$ の固有値となる。上式に基づいて求めた曲率の最大値と最小値をそれぞれ $k_{0}$、$k_{1}$ とする。

さて、我々の考えている $P$ の具体的な形は次式で与えられた。

\begin{array}
x &=& p\\
y &=& q\\
z &=& Ap^{2} + 2Cpq + Bq^{2}
\end{array}

これをこれまでの議論に当てはめて計算すると、$HJ^{-1}$ の固有値を求めることと、$M$ の固有値を求めることが同じであることに気づく (TODO: 計算過程を詳しく書く)。

下記は、各画素における $M$ の最大固有値と最小固有値をプロットするコードである。

from copy import deepcopy
import numpy as np
from scipy.ndimage import convolve, gaussian_filter
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt


class Harris:
    def __init__(self, orig_img, sigma):
        diff_filter = np.array([[-1.0, 0, 1.0],])
        dIdx = convolve(orig_img, diff_filter)
        dIdy = convolve(orig_img, diff_filter.T)
        self._A = gaussian_filter(np.power(dIdx, 2), sigma=sigma)
        self._B = gaussian_filter(np.power(dIdy, 2), sigma=sigma)
        self._C = gaussian_filter(dIdx * dIdy, sigma=sigma)

    def M(self, row, col):
        return np.array([[self._A[row, col], self._C[row, col]], [self._C[row, col], self._B[row, col]]])


def calc_quadratic_form(M, start, stop, num):
    rows = np.linspace(start[0], stop[0], num[0])
    cols = np.linspace(start[1], stop[1], num[1])
    values = []
    for row in rows:
        for col in cols:
            values.append((np.array([[row, col]]) @ M @ np.array([[row, col]]).T)[0, 0])
    return rows, cols, np.array(values)


def plot(rows, cols, values, fname):
    fig = plt.figure()
    ax = Axes3D(fig)
    X, Y = np.meshgrid(cols, rows)
    ax.plot_wireframe(X, Y, values.reshape((len(rows), len(cols))))
    plt.savefig(fname)


if __name__ == '__main__':

    # Create an empty image
    img = np.zeros(shape=(256, 256))
    harris = Harris(img, 1.0)

    # Draw a graph of the quadratic form of the empty image
    rows, cols, values = calc_quadratic_form(harris.M(10, 10), (-1, -1), (1, 1), (10, 10))
    plot(rows, cols, values, 'img/harris_quadratic_empty.png')

    # Create an image having a corner
    img = np.zeros(shape=(256, 256))
    col = 0
    for i in range(254):
        img[i, col] = img[i + 1, col] = img[i + 2, col] = 1.0
        if i < 128:
            col += 1
        else:
            col -= 1
    harris = Harris(img, 1.0)

    # Draw eigen values
    k0_list = []
    k1_list = []
    plot_range = np.arange(100, 140)
    for i in plot_range:
        for j in plot_range:
            ks, vs = np.linalg.eig(harris.M(i, j))
            k0_list.append(ks[0])
            k1_list.append(ks[1])
    plot(plot_range, plot_range, np.array(k0_list), 'img/k0_list.png')
    plot(plot_range, plot_range, np.array(k1_list), 'img/k1_list.png')

    plt.show()

最大固有値のグラフを次に示す。

最小固有値のグラフを次に示す。

上のグラフを見ると、エッジ付近では最大固有値は大きく最小固有値は0に近いことがわかる。またコーナー付近では、最大固有値と最小固有値がともに、そこそこの大きさを持っていることがわかる。面 (輝度値が一様な部分) では、最大固有値、最小固有値ともに $0$ であることがわかる。

このように最大固有値と最小固有値の大きさとその関係を見ることによって、コーナー検出ができる可能性があることがわかった。ただし、固有値を求めるためには行列の固有値を求めるという計算量の多い操作が必要となるという問題がある。

これを回避するために、ガウスの曲率と平均曲率という概念を導入する。

(TODO: 続きを書く)
(TODO: Qiita の数式エディタの関係で数式が正しく表示できていない部分があるので修正する)

画像中のエッジとコーナーの検出 (Morning Project Samurai Image Processing Series)

目的

関連して得られる知識

使用する環境

注意

Moravec's corner detector (MCD)

Harris' Corner Detector (HCD)

参考文献