OpenCVやってみる-10. テンプレートマッチング

かなり時間が空いてしまいました…
引き続きOpenCVのチュートリアルをやっていきたいと思います。
今回はテンプレートマッチングをやろうかと。
春のパン祭りシール点数集計的には結構使えそうな内容です。

テンプレートマッチング — OpenCV-Python Tutorials 1 documentation

テンプレートマッチング

cv2.matchTemplate()関数を使います。
処理内容は、

この関数はテンプレート画像を入力画像全体にスライド(2D convolutionと同様に)させ，テンプレート画像と画像の注目領域とを比較します．

とのことです。
詳細は以下に記載があります。

物体検出 — opencv 2.2 documentation

色々な比較方法

比較方法がいくつかあるので、ざっと確認してみます。

このあたりのサイトも参考に。

opencv matchTemplateの種類まとめ | 株式会社Ｚｉｎ

CV_TM_SQDIFF

$R(x,y)=\displaystyle \sum_{x',y'}(T(x',y')-I(x+x',y+y'))^{2}$

テンプレート画像と入力画像のピクセル値差分の2乗和のようです。直観的ですね。
値が小さいほど類似度が高いということになります。
CV_TM_SQDIFF_NORMED

$\displaystyle R(x,y)=\frac{\sum_{x',y'}(T(x',y')-I(x+x',y+y'))^{2}}{\sqrt{\sum_{x',y'}T(x',y')^2 \cdot \sum_{x',y'}I(x+x',y+y')^2}}$

CV_TM_SQDIFFでの結果を、テンプレート画像の2乗和、入力画像(のテンプレート画像サイズ分)の2乗和の積の平方根で割って(正規化して)います。
テンプレート画像のサイズ、テンプレート画像と入力画像の平均輝度に影響されにくい値という感じでしょうか。
入力画像の輝度が画像内で大きく違う場合はこれを使ったほうがいいか？
CV_TM_CCORR

$R(x,y)=\displaystyle \sum_{x',y'}(T(x',y') \cdot I(x+x',y+y'))$

テンプレート画像と入力画像の内積に当たるような計算。
輝度の高いエリアの位置がそろっているほど大きな値が出ます。
CV_TM_CCORR_NORMED

$\displaystyle R(x,y)=\frac{\sum_{x',y'}(T(x',y') \cdot I(x+x',y+y'))}{\sqrt{\sum_{x',y'}T(x',y')^2 \cdot \sum_{x',y'}I(x+x',y+y')^2}}$

CV_TM_CCORRの結果を、CVM_TM_SQDIFF_NORMEDと同様のもので割っています。
これも画像サイズ、平均輝度に影響されにくくするためのものと思われます。
CV_TM_CCOEFF

$R(x,y)=\displaystyle \sum_{x',y'}(T'(x',y') \cdot I'(x+x',y+y'))$
$T'(x',y')=T(x',y')-\frac{1}{w\cdot h} \sum_{x'',y''}T(x'',y'')$
$I'(x+x',y+y')=I(x+x',y+y')-\frac{1}{w\cdot h} \sum_{x'',y''}I(x+x'',y+y'')$

これは統計でいう共分散に当たるようです。
テンプレート画像、入力画像それぞれの平均を基準とした値を使用しています。
正負の値が出そうですが、画像の各点で類似度が高ければ同じ符号の値となるので、結局正の大きい値になります。
CV_TM_CCOEFF_NORMED

$\displaystyle R(x,y)=\frac{\sum_{x',y'}(T'(x',y') \cdot I'(x+x',y+y'))}{\sqrt{\sum_{x',y'}T'(x',y')^2 \cdot \sum_{x',y'}I'(x+x',y+y')^2}}$

CV_TM_CCOEFFの結果を、 $T'(x',y')$ 、 $I'(x+x',y+y')$ それぞれの2乗和の積の平方根で割っていますが、これは統計でいう相関係数に近いものになっています。
ただし、微妙に違いがあります。
相関係数の意味と求め方 - 公式と計算例

実際の画像で確認

今まで通り春のパン祭りシール台紙画像で試してみます。
1つのシール画像をテンプレートとしたいですが、前回の輪郭検出を利用してみます。
cv2.boundingRect()関数で外接矩形を得て、これをテンプレート領域とします。

領域(輪郭)の特徴 — OpenCV-Python Tutorials 1 documentation

まずはシールの輪郭を探します。
前回の記事で、台紙領域画像から輪郭検出をするといい具合にシール輪郭が取れていたので、それでやってみます。前回調べた閾値を使用します。

import cv2
img = cv2.imread('harupan_200317_1.jpg')
img = cv2.resize(img, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img_hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
ret, img_card = cv2.threshold(img_hsv[:,:,1], 50, 255, cv2.THRESH_BINARY_INV)

img_card_contours, img_card_hierarchy = cv2.findContours(img_card.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
img_card_contours_seal = []
for ctr in img_card_contours:
    if 2000 > cv2.contourArea(ctr) > 1000:
        img_card_contours_seal += [ctr]

これで検出した輪郭は以下の通りでしたが、シール輪郭でないものも含まれていました。
輪郭を1個ずつ順に見ていったところ、3つ目の輪郭がシール輪郭となっていました。

img_card_ctr = cv2.drawContours(img.copy(), img_card_contours_seal, 0, (255,0,0), 3)
cv2.imshow("img_card_ctr", img_card_ctr)
cv2.waitKey(0)

...

img_card_ctr = cv2.drawContours(img.copy(), img_card_contours_seal, 3, (255,0,0), 3)
cv2.imshow("img_card_ctr", img_card_ctr)
cv2.waitKey(0)

f:id:nokixa:20210507005223p:plain — 全輪郭を描画

f:id:nokixa:20210625005049p:plain — 0番の輪郭

f:id:nokixa:20210625004431p:plain — 3番の輪郭

外接矩形の位置、サイズを取得して、テンプレート画像を用意します。

>>> x,y,w,h = cv2.boundingRect(img_card_contours_seal[3])
>>> x,y,w,h
(96, 236, 41, 41)
>>> template = img[y:y+h,x:x+w,:]

テンプレート画像はこんな感じ。

f:id:nokixa:20210625012632p:plain

テンプレートマッチング適用

各比較方法でテンプレートマッチングしてみます。

res_sqdiff = cv2.matchTemplate(img.copy(), template, cv2.TM_SQDIFF)
res_sqdiff_normed = cv2.matchTemplate(img.copy(), template, cv2.TM_SQDIFF_NORMED)
res_ccorr = cv2.matchTemplate(img.copy(), template, cv2.TM_CCORR)
res_ccorr_normed = cv2.matchTemplate(img.copy(), template, cv2.TM_CCORR_NORMED)
res_ccoeff = cv2.matchTemplate(img.copy(), template, cv2.TM_CCOEFF)
res_ccoeff_normed = cv2.matchTemplate(img.copy(), template, cv2.TM_CCOEFF_NORMED)

from matplotlib import pyplot as plt
plt.subplot(231),plt.imshow(res_sqdiff,cmap='gray'),plt.title('SQDIFF'),plt.xticks([]),plt.yticks([])
plt.subplot(232),plt.imshow(res_sqdiff_normed,cmap='gray'),plt.title('SQDIFF_NORMED'),plt.xticks([]),plt.yticks([])
plt.subplot(233),plt.imshow(res_ccorr,cmap='gray'),plt.title('CCORR'),plt.xticks([]),plt.yticks([])
plt.subplot(234),plt.imshow(res_ccorr_normed,cmap='gray'),plt.title('CCORR_NORMED'),plt.xticks([]),plt.yticks([])
plt.subplot(235),plt.imshow(res_ccoeff,cmap='gray'),plt.title('CCOEFF'),plt.xticks([]),plt.yticks([])
plt.subplot(236),plt.imshow(res_ccoeff_normed,cmap='gray'),plt.title('CCOEFF_NORMED'),plt.xticks([]),plt.yticks([])
plt.show()

f:id:nokixa:20210625014037p:plain

結果をざっと見てみると、

SQDIFF、SQDIFF_NORMEDは値の小さいところ(黒いところ)、その他は値の大きいところが一致度が高いところになる
シールの中心と思われる点で一致度が高くなっている
シールとシールの間の隙間も一致度が高くなってしまっている
CCORRではうまく検出できていない、正規化の必要性は高いか

という感じです。

値の範囲の確認。

>>> res_sqdiff.max()
81624490.0
>>> res_sqdiff.min()
26.0
>>> res_sqdiff_normed.max()
1.0
>>> res_sqdiff_normed.min()
1.9568046e-07
>>> res_ccorr.max()
151010820.0
>>> res_ccorr.min()
30270670.0
>>> res_ccorr_normed.max()
0.9999999
>>> res_ccorr_normed.min()
0.7027236
>>> res_ccoeff.max()
16118729.0
>>> res_ccoeff.min()
-6135498.5
>>> res_ccoeff_normed.max()
0.99999917
>>> res_ccoeff_normed.min()
-0.38832977

正規化すると分かりやすい。

テンプレートマッチング結果データのサイズの確認。

>>> img.shape
(403, 302, 3)
>>> res_sqdiff.shape
(363, 262)

テンプレート画像のサイズは(w,h)=(41,41)でしたが、テンプレートマッチング結果は (403-41+1, 302-41+1)と、元画像サイズ-テンプレートサイズ+1となっています。

cv2.TM_CCOEFF、cv2.TM_CCOEFF_NORMEDで結構はっきりしているように見えるので、これでシール位置を出してみたいと思います。

極大位置算出

シール位置を点で算出したいですが、テンプレートマッチング結果は一致度合いの分布になっているので、このままでは使いづらいです。
2次元分布の極大位置を検出したい。

SciPyのsignal.argrelmax()関数で1次元方向の極大位置を出すことができます。
2次元での極大というのは2軸方向いずれでも極大になっている点、というのを探せばいいかと。

SciPy で離散データのピークを検出 | org-技術

scipy.signal.argrelmax — SciPy v1.7.0 Manual

signal.argrelmax()関数のorder引数で、前後の何個分のデータと比較するか指定できます。
まず垂直方向の極大位置を出してみます。

>>> from scipy import signal
>>> ccoeff_normed_peak_v = signal.argrelmax(res_ccoeff_normed, axis=0)
>>> ccoeff_normed_peak_v
(array([  1,   1,   1, ..., 361, 361, 361], dtype=int64), array([ 6,  8,  9, ...,  2, 26, 27], dtype=int64))
>>> ccoeff_normed_peak_v[0].shape
(10692,)
>>> ccoeff_normed_peak_v = signal.argrelmax(res_ccoeff_normed, axis=0, order=2)
>>> ccoeff_normed_peak_v[0].shape
(7825,)
>>> ccoeff_normed_peak_v = signal.argrelmax(res_ccoeff_normed, axis=0, order=5)
>>> ccoeff_normed_peak_v[0].shape
(3802,)
>>> ccoeff_normed_peak_v = signal.argrelmax(res_ccoeff_normed, axis=0, order=10)
>>> ccoeff_normed_peak_v[0].shape
(2652,)
>>> ccoeff_normed_peak_v = signal.argrelmax(res_ccoeff_normed, axis=0, order=20)
>>> ccoeff_normed_peak_v[0].shape

order引数の値を大きくするとノイズによるピーク値が除去されていきます。
今回はテンプレート画像のサイズが41x41なので、20の値でいいかと。
水平方向も計算します。

ccoeff_normed_peak_h = signal.argrelmax(res_ccoeff_normed, axis=1, order=20)

これらの極大位置をテンプレートマッチング結果画像上に表示してみます。

その前に画像フォーマットを整えておきます。
CCOEFF_NORMEDで計算したテンプレートマッチング結果はfloat32型で-1.0～1.0の範囲になるので、これをint8型の0～255にスケーリング、型変換します。
また、後で極大位置に色を付けて表示したいので、cv2.merge()関数を使ってカラー画像にしておきます。

>>> res_ccoeff_normed.dtype
dtype('float32')
>>> tm_img_gray = (res_ccoeff_normed + 1.0) * 128
>>> tm_img_gray.max()
255.9999
>>> tm_img_gray.min()
78.29379
>>> tm_img_gray = tm_img_gray.astype('uint8')
>>> tm_img_bgr = cv2.merge((tm_img_gray, tm_img_gray, tm_img_gray))
>>> tm_img_bgr.shape
(363, 262, 3)
>>> cv2.imshow('ccoeff_normed', tm_img_bgr)
>>> cv2.waitKey(0)

f:id:nokixa:20210626223412p:plain

スケーリングでピクセル値と表示色の対応が変わっているので、pyplotを使った時と比べて見た目も少し変わっています。

垂直方向、水平方向それぞれの極大位置座標のリストを用意します。

pts_v = []
for i in range(ccoeff_normed_peak_v[0].shape[0]):
    pts_v += [[ccoeff_normed_peak_v[0][i], ccoeff_normed_peak_v[1][i]]]

pts_h = []
for i in range(ccoeff_normed_peak_h[0].shape[0]):
    pts_h += [[ccoeff_normed_peak_h[0][i], ccoeff_normed_peak_h[1][i]]]

※signal.argrelmax()関数で得られる結果は2つのndarrayのタプルになりますが、 1つ目が水平座標、2つ目が垂直座標になるようです。
なので、上記のやり方だと画像上に描画等する際に逆にして座標指定しないといけません。
ちょっと失敗。

これらの点をテンプレートマッチング結果画像に描画すると、

for pt in pts_v:
    tm_img_bgr = cv2.circle(tm_img_bgr, (pt[1],pt[0]), 2, (255,0,0), -1)

for pt in pts_h:
    tm_img_bgr = cv2.circle(tm_img_bgr, (pt[1],pt[0]), 2, (0,255,0), -1)

cv2.imshow('ccoeff_normed', tm_img_bgr)
cv2.waitKey(0)

f:id:nokixa:20210626224839p:plain

やはり垂直、水平方向にいずれも極大を取る点を選べばよさそうです。
期待していない点もありそうですが、ピクセル値で選別すればよさそう。

選別閾値調査のためヒストグラムを表示。

plt.hist(tm_img_gray.ravel(), 256, [0,256]);plt.show()

f:id:nokixa:20210626233230p:plain:w300

200ぐらいでいいかな。
そして2次元極大点を探します。
閾値も適用。

pts_b = []
for pt in pts_v:
    if pt in pts_h:
        pts_b += [pt]

pts_b_th = []
for pt in pts_b:
    if tm_img_gray[pt[0],pt[1]] > 200:
        pts_b_th += [pt]

tm_img_bgr = cv2.merge((tm_img_gray, tm_img_gray, tm_img_gray))
for pt in pts_b_th:
    tm_img_bgr = cv2.circle(tm_img_bgr, (pt[1], pt[0]), 2, (0,0,255), -1)

cv2.imshow('ccoeff_normed', tm_img_bgr)
cv2.waitKey(0)

f:id:nokixa:20210627001411p:plain

思った以上にうまくいっています！

1点、2点のシール位置が取れている
0.5点のシール位置は除去されている
シール間の極大も除去されている

元画像にも表示してみます。

img_cp = img.copy()
for pt in pts_b_th:
    img_cp = cv2.rectangle(img_cp, (pt[1], pt[0]), (pt[1]+w, pt[0]+h), (0,0,255), 2)

cv2.imshow('img', img_cp)
cv2.waitKey(0)

f:id:nokixa:20210627003205p:plain

かなりいい結果と言えるんじゃないでしょうか。

2次元極大位置検出を関数化

他でも使えそうなので、関数化してみました。

def relmax_2d(m, order):
    from scipy import signal
    
    relmax = signal.argrelmax(m, axis=0, order=order)
    pts0 = []
    for i in range(relmax[0].shape[0]):
        pts0 += [[relmax[1][i], relmax[0][i]]]
    
    relmax = signal.argrelmax(m, axis=1, order=order)
    pts1 = []
    for i in range(relmax[0].shape[0]):
        pts1 += [[relmax[1][i], relmax[0][i]]]
    
    pts = []
    for pt in pts0:
        if pt in pts1:
            pts += [pt]
    
    return pts

引数 m : 入力の2次元配列(ndarray)
引数 order : 内部でsignal.argrelmax()関数に渡すorder引数の値
返り値 : 2次元極大位置座標のリスト、座標軸の順は上でやっていたのと逆にしました。

これを使って同じことをやってみると、

tm_pts = relmax_2d(res_ccoeff_normed, 20)
img_cp = img.copy()
for pt in tm_pts:
    if tm_img_gray[pt[1], pt[0]] > 200:
        img_cp = cv2.rectangle(img_cp, (pt[0], pt[1]), (pt[0]+w, pt[1]+h), (255,0,255), 2)
cv2.imshow('with function', img_cp)
cv2.waitKey(0)