end0tknr's kipple - 新web写経開発

http://d.hatena.ne.jp/end0tknr/ から移転しました

正規化相互相関を用いた Harrisコーナー検出結果の対応点探索 (python)

前回と以前のエントリの続き。というか、これまで理解不足で、さかのぼっていた感じ

Harrisコーナー検出(再び)と、pythonによる実装 - end0tknr's kipple - 新web写経開発

自己相関と相互相関の違いや、正規化相互相関 - end0tknr's kipple - 新web写経開発

#!/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

im1 = numpy.array(Image.open('sf_view1.jpg').convert('L'))
im2 = numpy.array(Image.open('sf_view2.jpg').convert('L'))
wid = 5

def main():
    # 画像1のHariisコーナー検出
    harrisim = compute_harris_response(im1,5)
    filtered_coords1 = get_harris_points(harrisim,wid+1)
    d1 = get_descriptors(im1,filtered_coords1,wid)
    # 画像2のHariisコーナー検出
    harrisim = compute_harris_response(im2,5)
    filtered_coords2 = get_harris_points(harrisim,wid+1)
    d2 = get_descriptors(im2,filtered_coords2,wid)

    # 画像1,2のコーナーの対応点探索
    matches = match_twosided(d1,d2)

    plb.figure()
    plb.gray()
    plot_matches(im1,im2,filtered_coords1,filtered_coords2,matches[:100])
    plb.savefig( '2_1_1.png' )

# グレースケール画像の画素に対応する 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

# 各点について、点の周辺で幅 2*wid+1 の近傍ピクセル値を返す。
#(点の最小距離 min_distance > wid を仮定する)
def get_descriptors(image,filtered_coords,wid=5):
    desc = []
    for coords in filtered_coords:
        patch = image[coords[0]-wid:coords[0]+wid+1,
                      coords[1]-wid:coords[1]+wid+1].flatten()
        desc.append(patch)

    return desc


# match()の双方向で一致を調べる
def match_twosided(desc1,desc2,threshold=0.5):

    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)

    ndx_12 = numpy.where(matches_12 >= 0)[0]

    # 非対称の場合を除去する
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1

    return matches_12

# 【正規化相互相関】を用い、第1の画像の各コーナー点記述子について
#  第2の画像の対応点を選択
def match(desc1,desc2,threshold=0.5):
    n = len(desc1[0])

    # 対応点ごとの距離
    d = -numpy.ones((len(desc1),len(desc2)))
    for i in range(len(desc1)):
        for j in range(len(desc2)):
            d1 = (desc1[i] - numpy.mean(desc1[i])) / numpy.std(desc1[i])
            d2 = (desc2[j] - numpy.mean(desc2[j])) / numpy.std(desc2[j])
            ncc_value = sum(d1 * d2) / (n-1)
            if ncc_value > threshold:
                d[i,j] = ncc_value

    ndx = numpy.argsort(-d)
    matchscores = ndx[:,0]
    return matchscores


# 対応点を線で結んで画像を表示する
# 入力: im1,im2(配列形式の画像)、locs1,locs2(特徴点座標)
# machescores(match()の出力)、
# show_below(対応の下に画像を表示するならTrue)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):

    im3 = appendimages(im1,im2)
    if show_below:
        im3 = numpy.vstack((im3,im3))

    plb.imshow(im3)

    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0: plb.plot([locs1[i][1],locs2[m][1]+cols1],
                         [locs1[i][0],locs2[m][0]],'c')
    plb.axis('off')

# 2つの画像を左右に並べる
def appendimages(im1,im2):

    # 行の少ない方を選び、空行を0で埋める
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]

    if rows1 < rows2:
        im1 = numpy.concatenate((im1,numpy.zeros((rows2-rows1,im1.shape[1]))),
                                axis=0)
    elif rows1 > rows2:
        im2 = numpy.concatenate((im2,numpy.zeros((rows1-rows2,im2.shape[1]))),
                                axis=0)
    # 行が同じなら、0で埋める必要なし

    return numpy.concatenate((im1,im2), axis=1)

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます。数分待ちますが... f:id:end0tknr:20171003214833p:plain