前回と以前のエントリの続き。というか、これまで理解不足で、さかのぼっていた感じ
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()
↑こう書くと↓こう表示されます。数分待ちますが...
