先日までのエントリで「Harrisコーナー検出は理解できた」と思い実践コンピュータビジョンのサンプルプログラムを写経しようとしたら、理解度が不十分だった為、やり直し。

Harrisコーナー検出の基礎

まず、ある画像上の点xについて、行列 Miを定義します。

$\Large M_I = \nabla I \ \nabla I^{T} = \left[ \begin{array}{c} I_x \\ I_y \end{array} \right] \left[ \begin{array}{cc} I_x & I_y \end{array} \right] = \left[ \begin{array}{cc} I_x^2 & I_x I_y \\\ I_x I_y & I_y^2 \end{array} \right]$

ここで、∇については以前のエントリに記載した通りです。

ベクトル偏微分演算子 : ∇ (ナブラ) - end0tknr's kipple - 新web写経開発

また、∇ I は微分係数Ix, Iyから構成され、次のように求められます。

$\Large I_x = I \ast G_{σx} \ 、\ I_y = I \ast G_{σy}$

( 上記において、Ｇσx, Ｇσy, ＊はそれぞれ、標準偏差σのガウス関数と、畳込み記号 )

また、先程のMiにおいては、階数(rank)=1 で、固有値は $\lambda_1 = | \nabla I |^2 と \lambda_2 = 0$ 。

次にガウス分布Ｇσ 等の重み行列を用いて $\Large \overline{M_{I}}$ を定義します ( Ｈａｒｒｉｓ行列 ) 。

$\Large \overline{M_{I}} = W \ast M_I$

$\overline{M_{I}}$ は、成分ごとの畳み込みで周辺ピクセルの局所平均となります。

( ガウス関数や畳み込みについては、以前のエントリ参照

画像処理 - エッジ抽出フィルタ(sobel, prewitt)と、平滑化フィルタ(gaussian) - end0tknr's kipple - 新web写経開発 )

$\nabla I$ の値により $\overline{M_{I}}$ の固有値は変化し、次のケースに分けられます。

固有値	判定
λ1, λ2 が大きな正値	コーナーあり
λ1 が大きく, λ2 ≒ 0	エッジあり
λ1 ≒0, λ2 ≒ 0	平面

ただし、コーナー判定において固有値を算出すると、計算コスト:大の為、次のような判定式を使用します。

$\Large \frac { det( \overline{M_{I}} ) }{ trace( \overline{M_{I}} ) }$

( det():行列式、trace():転置 )

pythonによるHarrisコーナー検出

実践コンピュータビジョンサンプルプログラム

相変わらずの上記の写経ですが...

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
from scipy.ndimage import filters

def main():
    im = numpy.array(Image.open('empire.jpg').convert('L'))
    # Harrisコーナー検出器の応答関数
    harrisim = compute_harris_response(im)
    #
    filtered_coords = get_harris_points(harrisim,6)
    # 検出されたコーナーの表示
    plot_harris_points(im, filtered_coords)

# グレースケール画像の画素に対応する Harrisコーナー検出器の応答関数
def compute_harris_response(im,sigma=3):
    # 微分係数
    imx = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    # Harris行列の成分
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    Wdet = Wxx*Wyy - Wxy**2 # 判別式   det()
    Wtr = Wxx + Wyy         # 対角成分 trace()

    return Wdet / Wtr

# Harris応答画像からコーナーを返す。
# min_distはコーナーや画像境界から分離する最小ピクセル数
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    # 閾値(threshold)を超えるコーナー候補を探す
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # 候補の座標を取得. 「T」は転置method
    coords = plb.array(harrisim_t.nonzero()).T
    # 候補の値を得る
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # 候補をsort
    index = plb.argsort(candidate_values)
    # 許容する点の座標を配列に格納
    allowed_locations = plb.zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # 最小距離を考慮しながら、最良の点を得る
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
                              (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0

    return filtered_coords

def plot_harris_points(image,filtered_coords):
    plb.figure()
    plb.gray()
    plb.imshow(image)
    plb.plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    plb.axis('off')
    plb.savefig( '2_1.png' )

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます f:id:end0tknr:20171003212441p:plain

end0tknr's kipple - web写経開発

太宰府天満宮の狛犬って、妙にカワイイ

Harrisコーナー検出(再び)と、pythonによる実装

Harrisコーナー検出の基礎

pythonによるHarrisコーナー検出