QNAP NAS購入-1

自分のPCのハードディスク容量がそもそも256GBしかなく、iPhoneで撮った写真(子供の写真)が増えてきたので、NASを買ってみました。

PCデータとiPhone写真のバックアップが目的。

NAS型番、仕様

購入したのはQNAP社のTS-230。
購入の決め手としては、

• RAID1構成にしたい、HDDの交換もできること
  • HDDが死んでもある程度大丈夫なようにしておきたいので
  • バックアップ用ポータブルHDDが壊れてデータ失ったことがあるので…
• iPhoneからいい感じに写真のアップロードができること
  • IODATA、BUFFALOなどでHDD込みで安いNASがありますが、スマホアプリがいまいちのようなので
  • QNAP、Synologyの2社ではアプリが充実しているよう
• QNAP、Synologyで極力安いものを考えて、以下が候補に。
  • Synology ds220j
  • QNAP TS-230
  • 本当はもう少し高くて高スペックなものがいいけど、経済的な事情から…
• ds220jは色々見ているうちに楽天やamazonで売り切れになった
• au PAY marketで、TS-230とHDDのセット品でお得なものがあったので
  • Western Digital RED 3TB x2とのセットで￥38384
  • TS-230本体は安いところで見ると2万ぐらいなので、HDD1つあたり9000ぐらい。いい感じ。
  • 21/8/2現在もう売り切れになってました。

[パーツセット] QNAP TS-230 + Western Digital WD30EFAX-RT 3TB×2台 セットの通販はau PAY マーケット - パソコン工房｜商品ロットナンバー：508173813

参考にしたサイト色々。

https://botti-bk.com/qnap_qumagie_review/

NASと外付けHDDはどっちが安全？最強の写真バックアップ環境について考えてみた！ - studio9

【2021年2月】買うべきSynologyのNASが分かるフローチャートを作りました！ - 人生総合ライフハック図書館

スペックは、

• CPU: Realtek RTD1296 ARM Cortex クアッドコア 1.4GHz
• RAM 2GB DDR4 (拡張不可)
• ギガビットイーサネット x1
• USB2.0 x1
• USB3.2 Gen1 x2
• 188.6(H) x 90.2(W) x 156.2 (D)

といったところ。

セットアップ

セットアップ手順を書いてみます。

ハードウェアのセットアップ

まず届いた状態。

開封します。
内容はTS-230本体、電源、LANケーブル、HDD固定用ねじ、あとはHDD2つです。

TS-230を開けてみます。
CPUのヒートシンク、SATAコネクタ、ファンなどなど見えます。

HDDを取り付けます。
ねじはなしでもいいみたいですが、一応取り付けておきました。

電源を接続、また、LANケーブルをWi-Fiルータに接続します。
普段PCやスマホからこのルータにつないでいます。

電源ボタンを押して起動。
ビープ音が鳴るらしいですが、ちゃんと聞いてませんでした。
LEDが点灯します。

NAS初期化

• PCにQfinderをインストール
  ユーティリティ | QNAP
  特記することはなし、普通にインストールします。
• Qfinderを起動
  早速見つかりました。

• そのまま置いておくと「スマートインストールガイド」なる画面が出てきました。
  が、「いいえ」をクリックすると、それ以降Qfinderを起動し直しても出なくなりました…

• と思ったら、「リフレッシュ」ボタンを押すとまた出てきました。

• 「スマートインストールガイド」で「はい」を押すと、ブラウザで設定画面が表示されます。
  ファームウェアインストールを開始しておきます。
  HDDも初期化されるそうなので要注意かも。
  本体交換するときはどうするんだろ？

• インストールが進んで、ときどきビープ音が聞こえました。
  その後も何回か鳴っています。

• カウントが終了すると、縦に長い画面が出て、「スマートインストレーションの開始」と出ています。取説に書いてるのと違いますが、スマートインストレーションとやらを開始してみます。

• またファームウェアのインストールとな？
  言われたとおりにやっておきます。

• 何やらぐるぐるしています…
  結構長い。

• 終わらなかったので、ブラウザの更新ボタンを押すと、どうも進めた感じ。
  保証サービスの画面が出ますが、QNAPのサイトに飛ぶだけなので閉じます。

• 今度こそちゃんとスマートインストレーションかぁ？

• ファームウェアバージョンは現状のままで次へ。

• NASの名前と管理者パスワード。ばれないように、忘れないように、慎重に。

• 時刻設定、そのままでいいかな。

一応テストだけしてみる。
ちょっと時間がかかって成功しました。

• IPはDHCPで。

• クロスプラットフォームファイル転送サービス、全OS入れとけばいいか。

• 最後に設定を確認して適用。
  取説にはディスク構成の選択があると書いてあるんだけど…出てこなかった…

• 初期化の警告、了解です。

• 適用してます…
  なかなか99%から進まない…
  ときどきビープ音が聞こえます。

• いけました！おめでとう。

NAS設定

前からの続きですが、ここからはNASの設定をしていく手順になります。

• ログインします。セキュリティ保護されたログインしようとしたらChromeに警告出されたので普通のログインにしておきます。

• ベータプログラムのお誘いがありましたが、いいえで。

• まずストレージ＆スナップショットのようこそ画面が出ました。よくわからんけど設定しとこかな。新規ストレージプールで。

• データ収集契約か… やめとこかな。いいえで。

• ストレージプールの作成ウィザードです。

• SEDセキュアストレージプールというのがある、どうしよう。
  と思いましたが、チェックを入れるとディスクが表示されませんでした。
  今回のHDDはSED機能なしということで。再度チェックを外して、ドライブを選択、RAID1に設定した状態で次に行きます。

• アラートの閾値？そのままにしておきます。

• 設定を確認して作成。 数分かかるそう。

• 作成完了。続いて新規ボリュームの作成。

• シックボリュームにしておきます。

• ボリュームの設定、名前はそのまま、サイズは「最大に設定」を押して最大に設定。
  メッセージが出るがどうしよう？スナップショット機能は使いたいけど…後で設定できるかな？ ちょっと心配なのでボリュームサイズ2TBにしておきます。
  拡張は後でできるでしょう。たぶん。

• 設定を確認して完了。
  初期化してます。

• 完了したかな。

• 左上のバッテンでストレージ＆スナップショットを閉じるとデスクトップ的な画面になります。
  adminというアカウント名は危ないそうなので、別のアカウントを作っておこう。
  「administrators」グループに所属させておきます。ユーザグループ、共有フォルダ権限などは後から設定し直しもできるので焦らず。
  最後に一旦ログアウト、作ったアカウントでログインし直してadminアカウントを無効にしておきます。

（NAS設定編6）ユーザー登録と共有フォルダ設定: 優柔不断流の買物道

• 2段階認証の設定もあるみたいなので、使っておきたい。
  モバイルデバイスが必要になります。私はiPhoneを使っているので、Google Authenticatorアプリをインストールします。

• 次の画面が出るので、Google AuthenticatorでQRコードを読み取って「次へ」を押します。

• Google Authenticatorに表示される6桁の数字を入力。次へ。
  数字の時間制限があるので早めに入力します。
  時間が過ぎたら別の数字になる。

• 代替検証方法の設定を設定。メールがいいかな。
  と設定したら、メール通知サービスの構成が必要とのこと。

• SMTPサービスの設定。自分のGmailからメッセージが送られるようにしました。
  Googleアカウントの承認が必要でした。

• 設定が完了すると、一旦ログアウトされます。再度のログイン時、ユーザ名とパスワードを入力した後に2段階認証の画面が出ます。
  「別の方法で認証」を選択すると、メールでセキュリティコードを送ることができます。メールで届いたセキュリティコードでログインできます。

ひとまず以上

取説の図と違う部分が色々ありました。ファームバージョンによるんだろうか。

なかなか設定項目が多くて進まない…
次回からはNAS運用方法の内容になります。

OpenCVやってみる-14. SIFTによる特徴点検出

今回はSIFTを扱います。

SIFT (Scale-Invariant Feature Transform)の導入 — OpenCV-Python Tutorials 1 documentation

SIFTは、2004年にブリティッシュコロンビア大学のD.Loweが発表したもので、 キーポイントの検出、特徴量の計算を行うアルゴリズムとのことです。

おおざっぱに自分なりにまとめると、

• ステレオ撮影画像や動画での物体追跡では、複数画像間で対応する点を見つける必要があるが、 これに使える画像上の点(特徴点)を見つけ、また、その点の特徴ベクトルを生成するアルゴリズム
• Harrisのコーナー検出でも特徴点検出はできるが、スケール不変でない(拡大・縮小されたり、 撮影時の距離が違う、撮影したカメラの解像度が違う、等に影響される)
• SIFTでは複数の空間スケールで特徴点検出を行うので、スケール不変性を持つ
• 回転不変性も持たせるため、特徴点近傍の角度ごとのヒストグラムを計算して、 回転角を計算する
• 特徴点の特徴量としては、特徴点周囲の4x4ブロック16個のヒストグラムを使う、 各々ビンの数8で、結果1特徴点あたり128要素のベクトルを持つ

以下も参考になりました。

https://qiita.com/icoxfog417/items/adbbf445d357c924b8fc

特許が取られていたようですが、2020年3月に特許権終了になったようです。
OpenCVでフリーに使えます。

SIFTの特許が切れてOpenCV v4.4.0から普通に使えるようになってた話 - Qiita

OpenCVの商用利用 - ari23の研究ノート

実践

やはり前回と同じ画像を使います。

まずは1つめの画像で様子見。

sift = cv2.SIFT_create()
kp1 = sift.detect(img1_gray, None)
img1_kp = cv2.drawKeypoints(img1, kp1, None, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow('KeyPoints1', img1_kp)

こんな形の特徴点が検出されます。いまいち規則性が分かりませんが…
sift.detect()での計算結果を少し確認。

>>> type(kp1)
<class 'list'>
>>> len(kp1)
1334
>>> type(kp1[0])
<class 'cv2.KeyPoint'>
>>> kp1[0].pt
(25.444393157958984, 6.472270965576172)
>>> kp1[0].size
2.4101293087005615
>>> kp1[0].angle
169.22996520996094
>>> kp1[0].response
0.015413731336593628
>>> kp1[0].octave
4588287
>>> kp1[0].class_id
-1

cv2.KeyPoint型のリストで、座標、サイズ(キーポイントの直径)、角度、キーポイントの強さ、 といった情報を含みます。

OpenCV: cv::KeyPoint Class Reference

以下気づき点。

• チュートリアルではcv2.SIFT()を使っていましたが、それでやるとなぜかその後のsift.detect()を実行したときにpythonが落ちてしまった… cv2.SIFT_create()なら大丈夫でした。
• cv2.drawKeypoints()で特徴点検出しますが、第3引数が必要でした。 colorという引数ですが、いまいち何をする引数なのか分からず…

参考

画像の特徴点を抽出する - Qiita

OpenCV: Drawing Function of Keypoints and Matches

colorを変更した結果

img1_kp = cv2.drawKeypoints(img1, kp1, (255,0,0), flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow('KeyPoints1', img1_kp)

描画の色は変わってるけども…

残りの2画像でもやってみます。
同じシール台紙を違う角度から撮っているので、同じような特徴点が出ることを期待。

kp2 = sift.detect(img2_gray, None)
img2_kp = cv2.drawKeypoints(img2, kp2, None, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
kp3 = sift.detect(img3_gray, None)
img3_kp = cv2.drawKeypoints(img3, kp3, None, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
imgs_kp = np.hstack((img1_kp, img2_kp, img3_kp))
cv2.imshow('KeyPoints', imgs_kp)

こっちのほうならまだ規則性が見えます。

• シールとシールの隙間に特徴点があることが多い
  特徴点が検出されていない箇所もあります。
• 右の2つの画像でシール台紙上の同じ点に特徴点が現れていることもある

一応cv2.detect()のflagsを変えてみました。

img1_kp = cv2.drawKeypoints(img1, kp1, None, flags=cv2.DRAW_MATCHES_FLAGS_DEFAULT)
img2_kp = cv2.drawKeypoints(img2, kp2, None, flags=cv2.DRAW_MATCHES_FLAGS_DEFAULT)
img3_kp = cv2.drawKeypoints(img3, kp3, None, flags=cv2.DRAW_MATCHES_FLAGS_DEFAULT)
imgs_kp = np.hstack((img1_kp, img2_kp, img3_kp))
cv2.imshow('KeyPoints', imgs_kp)

DEFAULT設定だと、大きさや方向の表現がなくなるようです。
分かりにくくなったし、あんまり面白くない…

特徴量

sift.compute()で特徴量が計算できます。

>>> kp1, des1 = sift.compute(img1_gray, kp1)
>>> type(des1)
<class 'numpy.ndarray'>
>>> des1.shape
(1334, 128)
>>> des1[0]
array([  0.,   1.,   0.,   4., 108.,  25.,   1.,   0.,  62.,  17.,   0.,
         3., 146.,  11.,   0.,   0., 146.,  49.,   1.,   4.,  18.,   5.,
         1.,   4.,   8.,   5.,   1.,   6.,  36.,  10.,   1.,   2.,   5.,
         6.,   1.,  12., 146.,  12.,   0.,   0.,  71.,  11.,   0.,   9.,
       146.,  12.,   0.,   1., 146.,  21.,   0.,   1.,  31.,  10.,   2.,
        15.,  20.,   1.,   1.,   2.,  63.,  17.,   1.,   7.,  17.,   2.,
         0.,   1., 138.,  42.,   1.,   4.,  58.,   1.,   0.,   0., 146.,
        64.,   2.,   8., 146.,  10.,   0.,   0.,  40.,  15.,   2.,  17.,
        35.,   5.,   1.,   4.,  87.,  19.,   1.,   3.,   7.,   0.,   0.,
         3.,  68.,   9.,   0.,   2.,  27.,   0.,   0.,   3.,  85.,  19.,
         1.,   7., 146.,   1.,   0.,   0.,  21.,  12.,   2.,  26.,  21.,
         1.,   0.,   0.,  60.,  37.,   3.,   4.], dtype=float32)

確かに1点あたり128次元のベクトルとなっています。

以上

まだ有用性は見えませんが、後々の特徴点マッチングで使えるかと思います。
次はSURFでの特徴点検出です。

OpenCVやってみる-13. Harrisコーナー検出

今回はHarrisコーナー検出をやります。
今更ですが、チュートリアルサイトの順番通りでなく、多少飛ばしている部分もあります。

Harrisコーナー検出 — OpenCV-Python Tutorials 1 documentation

Harrisコーナー検出

理論はチュートリアルサイトに書いてある通りです。
今回はcv2.cornerHarris()関数をさくっと試してみます。
画像は前回と同じ。

import cv2
img1 = cv2.imread('harupan_200317_1.jpg')
img2 = cv2.imread('harupan_210402_1.jpg')
img3 = cv2.imread('harupan_210402_2.jpg')
img1 = cv2.resize(img1, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img2 = cv2.resize(img2, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img3 = cv2.resize(img3, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)

import numpy as np
img123 = np.hstack((img1,img2,img3))
cv2.imshow('harupan', img123)

右の2つは同じ台紙を違う角度から撮影しているので、もしかしたらカメラの位置、角度の推定等できるかも。
いや、同じパターンがたくさんあるから難しいのか？
同じシールの検出等できるだろうか。

cv2.cornerHarris()ではいくつかパラメータがありますが、チュートリアルサイトと同じ値でやってみます。

img1_gray = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)
img2_gray = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)
img3_gray = cv2.cvtColor(img3, cv2.COLOR_BGR2GRAY)

img1_harris = cv2.cornerHarris(img1_gray, 2, 3, 0.04)
img2_harris = cv2.cornerHarris(img2_gray, 2, 3, 0.04)
img3_harris = cv2.cornerHarris(img3_gray, 2, 3, 0.04)

結果は、元画像と同じサイズのndaray配列となります。

>>> type(img1_harris)
<class 'numpy.ndarray'>
>>> img1_harris.shape
(403, 302)

特徴点を画像上に描画するとこんな感じ。

img1_cp = img1.copy()
img2_cp = img2.copy()
img3_cp = img3.copy()

img1_cp[img1_harris > 0.01*img1_harris.max()] = [0,255,0]
img2_cp[img2_harris > 0.01*img2_harris.max()] = [0,255,0]
img3_cp[img3_harris > 0.01*img3_harris.max()] = [0,255,0]

img123 = np.hstack((img1_cp, img2_cp, img3_cp));
cv2.imshow("Images", img123)

チュートリアルでは、検出点画像に膨張処理をして見やすくしていましたが、 今回はまずそのまま見てみました。細かい特徴点があるかもしれないので。

シールや台紙の細かい模様が特徴点として抽出されています。
よく見ると、台紙のコーナーも検出されているようです。

膨張処理をして表示してみます。

img1_harris = cv2.dilate(img1_harris, None)
img2_harris = cv2.dilate(img2_harris, None)
img3_harris = cv2.dilate(img3_harris, None)
img1_cp[img1_harris > 0.01*img1_harris.max()] = [0,255,0]
img2_cp[img2_harris > 0.01*img2_harris.max()] = [0,255,0]
img3_cp[img3_harris > 0.01*img3_harris.max()] = [0,255,0]
img123 = np.hstack((img1_cp, img2_cp, img3_cp));
cv2.imshow("Images", img123)

台紙のコーナーがちょっと分かりやすくなりました。

パラメータ変えてみる

シールの細かい文字しか特徴として出て来ませんでしたが、シールの点数のほうを 特徴として抽出できないだろうか。

cv2.cornerHarrisの引数の詳細を見てみましたが、blockSizeは使えなさそう。

特徴検出 — opencv 2.2 documentation

画像フィルタリング — opencv 2.2 documentation

ksizeを7にしてみたらどうだろう？
これはSobelフィルタのカーネルサイズになります。
カーネルサイズ3以外でのカーネルがどうなるのかというと、
以下のサイトにカーネルの取得方法が書かれていました。
cv2.getDerivKernels関数を使えばいいようです。結果も書かれています。

python - What are the kernel coefficients for OpenCV's Sobel filter for sizes larger than 3 x 3? - Stack Overflow

このカーネルの感じからすると、細かい変化はならされるように思います。 カーネルサイズ変更は有効なのでは？

img1_harris = cv2.cornerHarris(img1_gray, 2, 7, 0.04)
img2_harris = cv2.cornerHarris(img2_gray, 2, 7, 0.04)
img3_harris = cv2.cornerHarris(img3_gray, 2, 7, 0.04)
img1_harris = cv2.dilate(img1_harris, None)
img2_harris = cv2.dilate(img2_harris, None)
img3_harris = cv2.dilate(img3_harris, None)
img1_cp = img1.copy()
img2_cp = img2.copy()
img3_cp = img3.copy()
img1_cp[img1_harris > 0.01*img1_harris.max()] = [0,255,0]
img2_cp[img2_harris > 0.01*img2_harris.max()] = [0,255,0]
img3_cp[img3_harris > 0.01*img3_harris.max()] = [0,255,0]
img123 = np.hstack((img1_cp, img2_cp, img3_cp));
cv2.imshow("Images", img123)

細か過ぎるパターンは特徴点としては検出されないようになっています。
一応期待通り。

ここまで

Harrisコーナー検出はこのへんにしておきます。
他にも特徴点抽出の方法が色々紹介されていたので、引き続きやってみようと思います。

OpenCVやってみる-12. ハフ変換(円検出)

どんどん続きを進めます。
次はハフ変換による円検出。

ハフ変換による円検出 — OpenCV-Python Tutorials 1 documentation

ハフ変換による円検出

円上の点は、以下の式を満たします。

直線検出のときと同様、この式のパラメータ()の空間に、 その円上に乗った点の数をマッピングする感じかと。

直線検出と比べてパラメータが1つ多いため、単純に全パラメータ空間を 走査するのは非効率ということで、何がしかの手法を使っているようです。

前回と同じ画像を対象とします。
シールの領域を取得できるといいなと。

import cv2
img1 = cv2.imread('harupan_200317_1.jpg')
img2 = cv2.imread('harupan_210402_1.jpg')
img3 = cv2.imread('harupan_210402_2.jpg')
img1 = cv2.resize(img1, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img2 = cv2.resize(img2, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img3 = cv2.resize(img3, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)

import numpy as np
img123 = np.hstack((img1,img2,img3))
cv2.imshow('harupan', img123)

OpenCVでのハフ変換(円検出)関数

処理の内容はさておき、ハフ変換(円検出)の関数の使い方を確認します。

特徴検出 — opencv 2.2 documentation

【OpenCV】cv2.HoughCircles()の使い方【円を検出する】 | 資格マフィア

ハフ変換(円検出)の中でCannyエッジ検出が行われているようです。
前回は自分でグレースケール化、エッジ検出をやりましたが、今回はグレースケール化するだけになります。

引数が色々ありますが、以下のように設定します。

• dp : 投票空間(パラメータ空間？)の比率の逆数とのこと。参考サイトでは0.8～1.2がよさそうと書いてありますが、1と1.5ぐらいを試してみるかな。
• minDist : 円の中心同士の最小距離。今回は、シールが横に5個並んでいて、台紙が画像(横約300pixel)の横半分ぐらいに写っているものもあるので、シール間は約30pixel、値は20ぐらいにしておくかな。
• param1 : Cannyエッジ検出の閾値の大きいほうということで、前回使った 200 を与えます。
• param2 : 円検出の投票数閾値とのこと。シールが一番小さく写っている画像で直径30pixelぐらいなので、円周100pixelぐらい、閾値は50ぐらいにしておきます。 → やってみると、円が全く検出されなかったので、結局25にしました。
• minRadius : 最小半径、10にします。
• maxRadius : 最大半径、30にします。

img1_gray = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)
img2_gray = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)
img3_gray = cv2.cvtColor(img3, cv2.COLOR_BGR2GRAY)

circles1 = cv2.HoughCircles(img1_gray, cv2.HOUGH_GRADIENT, dp=1, minDist=20, param1=200, param2=25, minRadius=10, maxRadius=30)
circles2 = cv2.HoughCircles(img2_gray, cv2.HOUGH_GRADIENT, dp=1, minDist=20, param1=200, param2=25, minRadius=10, maxRadius=30)
circles3 = cv2.HoughCircles(img3_gray, cv2.HOUGH_GRADIENT, dp=1, minDist=20, param1=200, param2=25, minRadius=10, maxRadius=30)

計算結果を確認すると、

>>> circles1.shape
(1, 40, 3)
>>> circles2.shape
(1, 33, 3)
>>> circles3.shape
(1, 23, 3)
>>> circles1[0,0,:]
array([221.5, 329.5,  27.2], dtype=float32)

まあ妥当か。
結果の次元の先頭に謎の次元が付いています。
描画します。

circles1 = np.uint16(np.around(circles1))
circles2 = np.uint16(np.around(circles2))
circles3 = np.uint16(np.around(circles3))
imgs = [img1.copy(), img2.copy(), img3.copy()]
circles = [circles1, circles2, circles3]

for i in range(3):
    for c in circles[i][0]:
        imgs[i] = cv2.circle(imgs[i], (c[0],c[1]),c[2],(0,255,0),2)

imgs123 = np.hstack((imgs[0], imgs[1], imgs[2]))
cv2.imshow('Images', imgs123)

まあまあといったところでしょうか。

• 一番左の画像では全シールが検出されています。形がきちんとした丸ではないにも関わらず。
  ただし中心、半径は少しずれていたりします。
• 他の画像では、検出できていないシールがあります。
• 全画像で、台紙下側で円を誤検出しています。細かい模様があって、エッジ検出したときにパターンが発生してしまうからか。

パラメータを変えてやってみます。

• dp : 少し増やして1.5
• param2 : パラメータ空間が小さくなるということは、1つのパラメータ組み合わせに対する投票数が増えるのでは、ということで50

circles1 = cv2.HoughCircles(img1_gray, cv2.HOUGH_GRADIENT, dp=1.5, minDist=20, param1=200, param2=50, minRadius=10, maxRadius=30)
circles2 = cv2.HoughCircles(img2_gray, cv2.HOUGH_GRADIENT, dp=1.5, minDist=20, param1=200, param2=50, minRadius=10, maxRadius=30)
circles3 = cv2.HoughCircles(img3_gray, cv2.HOUGH_GRADIENT, dp=1.5, minDist=20, param1=200, param2=50, minRadius=10, maxRadius=30)

circles1 = np.uint16(np.around(circles1))
circles2 = np.uint16(np.around(circles2))
circles3 = np.uint16(np.around(circles3))
imgs = [img1.copy(), img2.copy(), img3.copy()]
circles = [circles1, circles2, circles3]

for i in range(3):
    for c in circles[i][0]:
        imgs[i] = cv2.circle(imgs[i], (c[0],c[1]),c[2],(0,255,0),2)

imgs123 = np.hstack((imgs[0], imgs[1], imgs[2]))
cv2.imshow('Images', imgs123)

だいたい同じような、ただ一番右の画像の結果は悪くなっています。
元の画像シール領域のサイズが小さいので、粗いパラメータ刻みではうまくいかなかった、ということかと。

以上

シール点数集計には使いにくそうなので、このあたりにしておきます。

まだまだチュートリアルは続きます。

OpenCVやってみる-11. ハフ変換

今回はハフ変換をやってみます。

ハフ変換による直線検出 — OpenCV-Python Tutorials 1 documentation

ハフ変換

ハフ変換は2値画像を直線のパラメータ空間に持っていくような変換、 ということでいいのかな。
空間中の直線は、以下の式で表すことができます。

変換後の空間の各点には、元画像の中でその直線に乗った点の数が入ります。
これを使って、閾値以上の値を得たパラメータの直線が検出された直線、とすることができます。

実践

春のパン祭りの台紙外枠を検出してみたいと思います。
3種類用意しています。

import cv2
img1 = cv2.imread('harupan_200317_1.jpg')
img2 = cv2.imread('harupan_210402_1.jpg')
img3 = cv2.imread('harupan_210402_2.jpg')
img1 = cv2.resize(img1, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img2 = cv2.resize(img2, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img3 = cv2.resize(img3, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)

import numpy as np
img123 = np.hstack((img1,img2,img3))
cv2.imshow('harupan', img123)

まずは輪郭検出して2値画像にします。
Canny法を使いますが、閾値は前にやったように100, 200とします。

OpenCVやってみる-8. エッジ検出 - 勉強しないとな～blog

img1_gray = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)
img2_gray = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)
img3_gray = cv2.cvtColor(img3, cv2.COLOR_BGR2GRAY)
img1_canny = cv2.Canny(img1_gray, 100, 200)
img2_canny = cv2.Canny(img2_gray, 100, 200)
img3_canny = cv2.Canny(img3_gray, 100, 200)
img123_canny = np.hstack((img1_canny, img2_canny, img3_canny))
cv2.imshow('Canny', img123_canny)

これらに対してハフ変換を実施します。

直線とみなされるのに必要な最低限の投票数を意味するしきい値

というパラメータがありますが、今回は台紙が画像の縦横の半分ぐらい以上には写っているように見えるので、 横サイズの1/3ぐらい=100にしておこうと思います。

>>> img1.shape
(403, 302, 3)

ハフ変換実施。の刻みは1pixel、の刻みは1°です。

lines1 = cv2.HoughLines(img1_canny, 1, np.pi/180, 100)
lines2 = cv2.HoughLines(img2_canny, 1, np.pi/180, 100)
lines3 = cv2.HoughLines(img3_canny, 1, np.pi/180, 100)

結果は、float32型、3次元のndarrayになっていて、各次元の意味は

     lines[直線インデックス, ?, (0) or (1)]

となっています。 (2つ目の次元の意味は分かりません…)

>>> type(lines1)
<class 'numpy.ndarray'>
>>> lines1.shape
(1155, 1, 2)
>>> lines2.shape
(60, 1, 2)
>>> lines3.shape
(15, 1, 2)
>>> lines1[0]
array([[387.       ,   1.5882496]], dtype=float32)

検出した直線を元の画像に描画します。

imgs = [img1.copy(), img2.copy(), img3.copy()]
lines = [lines1, lines2, lines3]
for i in range(3):
    for line in lines[i]:
        rho = line[0,0]
        theta = line[0,1]
        a = np.cos(theta)
        b = np.sin(theta)
        x0 = a*rho
        y0 = b*rho
        x1 = int(x0 + 1000*(-b))
        y1 = int(y0 + 1000*(a))
        x2 = int(x0 - 1000*(-b))
        y2 = int(y0 - 1000*(a))
        imgs[i] = cv2.line(imgs[i], (x1,y1), (x2,y2), (0,0,255), 2)

imgs123 = np.hstack((imgs[0], imgs[1], imgs[2]))
cv2.imshow('Images', imgs123)

これはあんまりですねー。
1つ目の画像で大量の線が出てしまっていますが、エッジ点が多過ぎるからか？
斜めに見たときに、通過する点の数が閾値(100)を超えてしまったものかと。
3つ目の画像では台紙の直線が検出されていません。ちょっと曲がっているから？

cv2.HoughLines()の3つ目の引数は角度の刻み(rad)ですが、これを粗くしたらどうだろう。
"直線に乗っている"の基準も緩めにならないか？

lines1 = cv2.HoughLines(img1_canny, 1, np.pi/90, 100)
lines2 = cv2.HoughLines(img2_canny, 1, np.pi/90, 100)
lines3 = cv2.HoughLines(img3_canny, 1, np.pi/90, 100)
lines = [lines1, lines2, lines3]
imgs = [img1.copy(), img2.copy(), img3.copy()]
for i in range(3):
    for line in lines[i]:
        rho = line[0,0]
        theta = line[0,1]
        a = np.cos(theta)
        b = np.sin(theta)
        x0 = a*rho
        y0 = b*rho
        x1 = int(x0 + 1000*(-b))
        y1 = int(y0 + 1000*(a))
        x2 = int(x0 - 1000*(-b))
        y2 = int(y0 - 1000*(a))
        imgs[i] = cv2.line(imgs[i], (x1,y1), (x2,y2), (0,0,255), 2)

imgs123 = np.hstack((imgs[0], imgs[1], imgs[2]))
cv2.imshow('Images', imgs123)

そうはなりませんでした。
単純に検出される直線が減っています。

確率的ハフ変換

もう1つのハフ変換の関数cv2.HoughLinesP()(確率的ハフ変換)では 最小の直線長さ、最大の直線ギャップといったパラメータがあるので、うまくやってくれそうな気がします。
1枚目の画像は特にバラバラの点の集まりが直線と認識されてしまっている感じなので、 ギャップの制約で識別できるかも。

最小の直線長さ : 50
最大の直線ギャップ : 5

ぐらいでやってみようと思います。
※以下でcv2.HoughLinesP()関数を使ったとき、minLineLengthとmaxLineGapのパラメータは指定しないと思ったようになりませんでした。実は順番が違うとか？他のパラメータが途中に入っているとか？

linesP1 = cv2.HoughLinesP(img1_canny, 1, np.pi/180, 100, minLineLength=50, maxLineGap=5)
linesP2 = cv2.HoughLinesP(img2_canny, 1, np.pi/180, 100, minLineLength=50, maxLineGap=5)
linesP3 = cv2.HoughLinesP(img3_canny, 1, np.pi/180, 100, minLineLength=50, maxLineGap=5)
lines = [linesP1, linesP2, linesP3]
imgs = [img1.copy(), img2.copy(), img3.copy()]

for i in range(3):
    for line in lines[i]:
        imgs[i] = cv2.line(imgs[i], (line[0,0],line[0,1]), (line[0,2],line[0,3]), (0,0,255), 2)
imgs123 = np.hstack((imgs[0], imgs[1], imgs[2]))
cv2.imshow('Images', imgs123)

うーん、
欲しくない直線は結構消えましたが、欲しい外枠直線があまり検出されていない…

一応、cv2.HoughLinesP()の結果のサイズも確認しておきます。

>>> linesP1.shape
(58, 1, 4)
>>> linesP2.shape
(21, 1, 4)
>>> linesP3.shape
(7, 1, 4)

思うに、"直線に乗っている"の判定がシビアなのでは？
エッジ画像を見るとかなり細い線になっているので、"乗っている"判定になりにくいのでは。

もしかしてモルフォロジー変換が使えないかな？？

モルフォロジー変換併用

エッジ画像を膨張処理で太くしてみます。

kernel = np.ones((3,3), np.uint8)
img1_canny_dlt = cv2.dilate(img1_canny, kernel, iterations=1)
img2_canny_dlt = cv2.dilate(img2_canny, kernel, iterations=1)
img3_canny_dlt = cv2.dilate(img3_canny, kernel, iterations=1)
img123_canny = np.hstack((img1_canny_dlt, img2_canny_dlt, img3_canny_dlt))
cv2.imshow('Images_canny', img123_canny)

さあこれでどうだぁ

linesP1 = cv2.HoughLinesP(img1_canny_dlt, 1, np.pi/180, 100, minLineLength=50, maxLineGap=5)
linesP2 = cv2.HoughLinesP(img2_canny_dlt, 1, np.pi/180, 100, minLineLength=50, maxLineGap=5)
linesP3 = cv2.HoughLinesP(img3_canny_dlt, 1, np.pi/180, 100, minLineLength=50, maxLineGap=5)
lines = [linesP1, linesP2, linesP3]
imgs = [img1.copy(), img2.copy(), img3.copy()]
for i in range(3):
    for line in lines[i]:
        imgs[i] = cv2.line(imgs[i], (line[0,0],line[0,1]), (line[0,2],line[0,3]), (0,0,255), 2)

imgs123 = np.hstack((imgs[0], imgs[1], imgs[2]))
cv2.imshow('Images', imgs123)

ありゃー、台紙内部の点が直線として認識されてしまったー

最大ライン間ギャップを短くしてみます。

linesP1 = cv2.HoughLinesP(img1_canny_dlt, 1, np.pi/180, 100, minLineLength=50, maxLineGap=2)
linesP2 = cv2.HoughLinesP(img2_canny_dlt, 1, np.pi/180, 100, minLineLength=50, maxLineGap=2)
linesP3 = cv2.HoughLinesP(img3_canny_dlt, 1, np.pi/180, 100, minLineLength=50, maxLineGap=2)
...

これもダメ。
内部点はだいぶ減りましたがまだあって、外枠の直線も結局ちゃんと取れていません。

今回の画像は欲しい線が少し曲がっていて、かつ2値化した後に 直線以外のものも残ってしまっているのでいい結果にならなかったのかなと。

これをやるなら前々回の輪郭検出をやったほうがよさそうです。

まとめ

今回はうまくいかなかったパターンですね。
また使える機会があれば使ってみたい。

次回はハフ変換による円検出をやってみます。

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

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

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

テンプレートマッチング

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

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

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

物体検出 — opencv 2.2 documentation

色々な比較方法

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

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

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

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

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

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

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

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

• CV_TM_CCOEFF_NORMED
  
  
  
  CV_TM_CCOEFFの結果を、、それぞれの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)

全輪郭を描画

0番の輪郭

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,:]

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

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

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

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()

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

• 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)

スケーリングでピクセル値と表示色の対応が変わっているので、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)

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

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

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

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)

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

• 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)

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

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)

同じようにできています。

まとめ

テンプレートマッチングでシール位置の検出ができるようになりました。
ついでに2次元配列の極大位置を探す方法も作ることができました。

欲を言えば点数識別までしたいところ。
回転に対するロバスト性も検討しないといけない。
そんなにシールをまっすぐ貼る人ばかりではないと思うので。

まだOpenCVチュートリアルの内容があるので、引き続き進めていきたいと思います。

OpenCVやってみる-9. 輪郭検出

OpenCVチュートリアルの続きです。

OpenCVにおける輪郭(領域) — OpenCV-Python Tutorials 1 documentation

とりあえず輪郭検出、表示

cv2.findContours()関数で、輪郭を検出することができます。
詳細はともかく、やってみます。

2値画像に対しての処理となること、白い物体と黒い背景があり、物体の輪郭を検出するということで、2値化のチュートリアルのところでやった画像を使ってみます。

HSV値からシール領域、台紙領域(および点数数字)を検出したものです。

import cv2
img1 = cv2.imread('harupan_200317_1.jpg')
img1 = cv2.resize(img1, None, fx=0.1, fy=0.1, interpolation=cv2.INTER_AREA)
img1_hsv = cv2.cvtColor(img1, cv2.COLOR_BGR2HSV)
ret, img1_seal = cv2.threshold(img1_hsv[:,:,0], 160, 255, cv2.THRESH_BINARY)
ret, img1_card = cv2.threshold(img1_hsv[:,:,1], 50, 255, cv2.THRESH_BINARY_INV)

from matplotlib import pyplot as plt
plt.subplot(121), plt.imshow(img1_seal, cmap='gray'), plt.title('Seal area'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img1_card, cmap='gray'), plt.title('Card area'), plt.xticks([]), plt.yticks([])
plt.show()

この画像でcv2.findContours()関数を実行、cv2.drawContours()関数で輪郭を表示します。

img1_seal_cp = img1_seal
img1_card_cp = img1_card
img1_seal_contours, img1_seal_hierarchy = cv2.findContours(img1_seal_cp, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
img1_card_contours, img1_card_hierarchy = cv2.findContours(img1_card_cp, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)

img1_seal_ctrs = cv2.drawContours(img1.copy(), img1_seal_contours, -1, (0,255,0), 2)
img1_card_ctrs = cv2.drawContours(img1.copy(), img1_card_contours, -1, (255,0,0), 3)
img1_seal_ctrs = cv2.cvtColor(img1_seal_ctrs, cv2.COLOR_BGR2RGB)
img1_card_ctrs = cv2.cvtColor(img1_card_ctrs, cv2.COLOR_BGR2RGB)
plt.subplot(121), plt.imshow(img1_seal_ctrs), plt.title('Seal contours'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img1_card_ctrs), plt.title('Card contours'), plt.xticks([]), plt.yticks([])
plt.show()

シールや台紙の輪郭はとれていそうですが、細かい余計な輪郭も出ているようです。

いくつか気づき点等。

• 今回の環境ではcv2.findContours()関数の戻り値はcontours, hierarchyの2つだけでした。OpenCVのバージョンによるのか？
  以下のサイトではこの戻り値2つのパターンで記載されていました。
  OpenCV - findContours で画像から輪郭を抽出する方法 - pystyle
• cv2.drawContours()関数では、引数の画像に上書きしてしまうようです。輪郭画像を元の画像に重ねるとき、元の画像をコピーしています。画像のコピーは、単純な代入だと参照になるようだったので、copy()関数を使います。
• カラー画像表示の際、OpenCVではBGRフォーマット、matplotlibではRGBフォーマットで扱われるので、変換しています。

輪郭の特徴

輪郭の面積や周長を出せるようなので、これでほしい輪郭だけフィルタリングすることができるのでは？

まず面積

シール領域画像、台紙領域画像からの輪郭の面積の分布がどうなっているかを調べてみます。

img1_seal_areas = []
for ctr in img1_seal_contours:
    img1_seal_areas += [cv2.contourArea(ctr)]

img1_card_areas = []
for ctr in img1_card_contours:
    img1_card_areas += [cv2.contourArea(ctr)]

以下のように、輪郭面積の範囲を確認しながら、ヒストグラムを見てみました。

>>> max(img1_seal_areas)
1496.0
>>> plt.hist(img1_seal_areas, 150, [0,1500])
(array([819.,  96.,  34.,  30.,  12.,   2.,   5.,   1.,   0.,   0.,   0.,
         0.,   1.,   1.,   4.,   7.,   3.,   0.,   1.,   1.,   0.,   0.,
         2.,   1.,   1.,   0.,   1.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   1.,   0.,   1.,   0.,   0.,   0.,   0.,   0.,   2.,
         0.,   0.,   0.,   0.,   0.,   0.,   1.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   1.,   1.,   0.,   1.,
         0.,   0.,   1.,   0.,   0.,   0.,   0.,   2.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   1.,
         0.,   0.,   0.,   0.,   1.,   0.,   1.,   2.,   0.,   1.,   1.,
         0.,   2.,   0.,   1.,   2.,   1.,   1.,   0.,   0.,   1.,   0.,
         0.,   0.,   1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   1.]), array([   0.,   10.,   20.,   30.,   40.,   50.,   60.,   70.,   80.,
         90.,  100.,  110.,  120.,  130.,  140.,  150.,  160.,  170.,
        180.,  190.,  200.,  210.,  220.,  230.,  240.,  250.,  260.,
        270.,  280.,  290.,  300.,  310.,  320.,  330.,  340.,  350.,
        360.,  370.,  380.,  390.,  400.,  410.,  420.,  430.,  440.,
        450.,  460.,  470.,  480.,  490.,  500.,  510.,  520.,  530.,
        540.,  550.,  560.,  570.,  580.,  590.,  600.,  610.,  620.,
        630.,  640.,  650.,  660.,  670.,  680.,  690.,  700.,  710.,
        720.,  730.,  740.,  750.,  760.,  770.,  780.,  790.,  800.,
        810.,  820.,  830.,  840.,  850.,  860.,  870.,  880.,  890.,
        900.,  910.,  920.,  930.,  940.,  950.,  960.,  970.,  980.,
        990., 1000., 1010., 1020., 1030., 1040., 1050., 1060., 1070.,
       1080., 1090., 1100., 1110., 1120., 1130., 1140., 1150., 1160.,
       1170., 1180., 1190., 1200., 1210., 1220., 1230., 1240., 1250.,
       1260., 1270., 1280., 1290., 1300., 1310., 1320., 1330., 1340.,
       1350., 1360., 1370., 1380., 1390., 1400., 1410., 1420., 1430.,
       1440., 1450., 1460., 1470., 1480., 1490., 1500.]), <BarContainer object of 150 artists>)
>>> 
>>> max(img1_card_areas)
79053.5
>>> plt.hist(img1_card_areas, 80, [0, 80000])
(array([971.,  26.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   1.]), array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,
        8000.,  9000., 10000., 11000., 12000., 13000., 14000., 15000.,
       16000., 17000., 18000., 19000., 20000., 21000., 22000., 23000.,
       24000., 25000., 26000., 27000., 28000., 29000., 30000., 31000.,
       32000., 33000., 34000., 35000., 36000., 37000., 38000., 39000.,
       40000., 41000., 42000., 43000., 44000., 45000., 46000., 47000.,
       48000., 49000., 50000., 51000., 52000., 53000., 54000., 55000.,
       56000., 57000., 58000., 59000., 60000., 61000., 62000., 63000.,
       64000., 65000., 66000., 67000., 68000., 69000., 70000., 71000.,
       72000., 73000., 74000., 75000., 76000., 77000., 78000., 79000.,
       80000.]), <BarContainer object of 80 artists>)

台紙領域画像からの輪郭で、大きな輪郭が1つ見られます。
これが台紙全体の輪郭ではないかと期待されるので、見てみます。

for ctr in img1_card_contours:
    if cv2.contourArea(ctr) > 79000:
        img1_card_contour_largest += [ctr]
img1_card_ctrs = cv2.drawContours(img1.copy(), img1_card_contour_largest, 1, (255,0,0), 3)
cv2.imshow('Card contours', img1_card_ctrs)

予想通りです。
射影変換用の角の4点検出に使えそうな感じがします。

後はシール領域、また、数字の輪郭を検出したい。
ヒストグラムを今度はもう少し細かく見てみます。

>>> plt.hist(img1_seal_areas, 100, [0,100])
(array([408.,  51., 144.,  24.,  69.,  26.,  29.,  21.,  33.,  14.,  18.,
        21.,  11.,  12.,   6.,   7.,   4.,   6.,   7.,   4.,   6.,   5.,
         4.,   5.,   2.,   1.,   2.,   5.,   3.,   1.,   3.,   6.,   5.,
         2.,   1.,   2.,   4.,   1.,   2.,   4.,   1.,   2.,   2.,   2.,
         2.,   1.,   0.,   1.,   0.,   1.,   0.,   0.,   1.,   0.,   0.,
         0.,   0.,   0.,   1.,   0.,   1.,   1.,   0.,   0.,   2.,   1.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.]), array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,
        11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,
        22.,  23.,  24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,
        33.,  34.,  35.,  36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,
        44.,  45.,  46.,  47.,  48.,  49.,  50.,  51.,  52.,  53.,  54.,
        55.,  56.,  57.,  58.,  59.,  60.,  61.,  62.,  63.,  64.,  65.,
        66.,  67.,  68.,  69.,  70.,  71.,  72.,  73.,  74.,  75.,  76.,
        77.,  78.,  79.,  80.,  81.,  82.,  83.,  84.,  85.,  86.,  87.,
        88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,  96.,  97.,  98.,
        99., 100.]), <BarContainer object of 100 artists>)
>>> plt.hist(img1_card_areas, 100, [0,2000])
(array([869.,  49.,  11.,  20.,   6.,   2.,   3.,   1.,   0.,   4.,   1.,
         0.,   0.,   1.,   0.,   0.,   1.,   1.,   0.,   0.,   0.,   0.,
         0.,   1.,   1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   4.,   4.,   5.,   6.,   4.,   1.,   1.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   1.,   0.,
         0.]), array([   0.,   20.,   40.,   60.,   80.,  100.,  120.,  140.,  160.,
        180.,  200.,  220.,  240.,  260.,  280.,  300.,  320.,  340.,
        360.,  380.,  400.,  420.,  440.,  460.,  480.,  500.,  520.,
        540.,  560.,  580.,  600.,  620.,  640.,  660.,  680.,  700.,
        720.,  740.,  760.,  780.,  800.,  820.,  840.,  860.,  880.,
        900.,  920.,  940.,  960.,  980., 1000., 1020., 1040., 1060.,
       1080., 1100., 1120., 1140., 1160., 1180., 1200., 1220., 1240.,
       1260., 1280., 1300., 1320., 1340., 1360., 1380., 1400., 1420.,
       1440., 1460., 1480., 1500., 1520., 1540., 1560., 1580., 1600.,
       1620., 1640., 1660., 1680., 1700., 1720., 1740., 1760., 1780.,
       1800., 1820., 1840., 1860., 1880., 1900., 1920., 1940., 1960.,
       1980., 2000.]), <BarContainer object of 100 artists>)

やっぱり台紙領域画像からの輪郭のほうが分かりやすそうで、面積1000以上に分布している25個の輪郭がシール輪郭に当たるのではと考えられます。(実際にはシールは23個しかないので、2個は違うものです。)

img1_card_contours_seal = []
for ctr in img1_card_contours:
    if 2000 > cv2.contourArea(ctr) > 1000:
        img1_card_contours_seal += [ctr]

img1_card_ctrs = cv2.drawContours(img1.copy(), img1_card_contours_seal, -1, (255,0,0), 3)
cv2.imshow('Card contours', img1_card_ctrs)

これもおよそ予想通り、いい感じです。
シールではない2つの輪郭は、台紙の下のほうの四角い領域でした。

最小外接円で円形領域を探す

シール領域画像の輪郭があまり使えていないので、なんとかならないかと。
シールの輪郭はだいたい円形をしているので、最小外接円との面積の差が小さいもの、で探すことができるのでは？と考えました。

img1_seal_ratios = []
for ctr in img1_seal_contours:
    (x, y), radius = cv2.minEnclosingCircle(ctr)
    img1_seal_ratios += [cv2.contourArea(ctr) / (radius*radius*np.pi)]

>>> max(img1_seal_ratios)
0.8911880248803259
>>> plt.hist(img1_seal_ratios, 100, [0, 1.0])
(array([367.,   0.,   0.,   1.,   1.,   1.,   8.,   6.,   4.,   3.,   6.,
         9.,  12.,   4.,   6.,   6.,   1.,   4.,   5.,   5.,   3.,   3.,
         4.,   4.,   4.,   7.,   6.,   4.,   3.,   3.,   6.,  39.,   3.,
         0.,   4.,   6.,   8.,   7.,  18.,   8.,  13.,   1.,   3.,   0.,
         7.,   6.,   3.,  24.,   3.,   2.,  11.,   2.,   5.,   5.,   6.,
        10.,  74.,   8.,   1.,   3.,   5.,  15.,   3., 160.,   2.,   2.,
         2.,   9.,   1.,   1.,  18.,   3.,   1.,   9.,  12.,   3.,   0.,
         1.,   3.,   5.,   3.,   3.,   4.,   2.,   0.,   1.,   2.,   0.,
         0.,   4.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         0.]), array([0.  , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 ,
       0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21,
       0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3 , 0.31, 0.32,
       0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43,
       0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53, 0.54,
       0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62, 0.63, 0.64, 0.65,
       0.66, 0.67, 0.68, 0.69, 0.7 , 0.71, 0.72, 0.73, 0.74, 0.75, 0.76,
       0.77, 0.78, 0.79, 0.8 , 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87,
       0.88, 0.89, 0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,
       0.99, 1.  ]), <BarContainer object of 100 artists>)

これまた微妙な分布…
とりあえず面積比率の大きい(外接円との面積の差が小さい)ものをいくつかを取ってみたいと思います。

img1_seal_contours_seal = []
for i in range(len(img1_seal_contours)):
    if img1_seal_ratios[i] > 0.77:
        img1_seal_contours_seal += [img1_seal_contours[i]]

img1_seal_ctrs = cv2.drawContours(img1.copy(), img1_seal_contours_seal, -1, (0,255,0), 3)
cv2.imshow("Seal contours", img1_seal_ctrs)

やっぱりいまいちですが、このやり方も使い道はありそうです。
そもそもシール領域画像にもう少し前処理が必要だったかもしれません。

以上

今回はここまでにします。
あっという間に春のパン祭りが終わってしまった…
来年に生かせればいいな。