概要
位置合わせというといまいちピンとこないですが、画像内のある部分の座標を合わせちゃおう！と言えばイメージがつくかもしれません。
この記事は、その位置合わせをC++で実装してみたよ、というものです。
やはりかなり多くの人がPythonでやっていて、C++ではそのままのコードがあまりなかったので情報共有がてら執筆します。
最終的には、Pythonでも同様の位置合わせを行っていき、速度や精度を観察していきたいと思います。
皆様の位置合わせの引き出しが少しでも増えたら幸いです。
位置合わせって？
こんな大々的に言っているのですが、こーゆージャンルが存在するのかはわかりません。
私も半年くらい前に知りました(笑)
それでは、代表的な位置合わせ手法を１つ紹介しますので雰囲気をつかんでください。
環境
windows10
Python
Python 3.7.4
OpenCV 3.4.2
C++
C++14
OpenCV 4.5.1
テンプレートマッチング
はい、とても有名ですね。
OpenCVで関数化しちゃってます。
なにをしているかを軽く説明すると、テンプレートマッチングは画像を左上から探索していき、輝度情報の差の二乗などを使って類似度を計算しています。
しかし、この方法では物体が回転したり、物体のサイズが違ったりすると全く違う座標として計算されてしまいます。
初耳で「確かに！」とはならないと思うので百聞は一見に如かず。実際に見てみましょう！
実験
まず今回使用する画像です。
適当にお金を並べたものと1円玉をトリミングしたものです。
マークダウンのプレビューだと1円玉単品の画像が大きくなっちゃってますが、実際は同じサイズです。
さっそくこれを1円玉でテンプレートマッチングしてみましょう！
はい。1円玉のところは特に中心が白くなっています。
一番類似度が高かった場所をプログラムで囲ってみます。
左上らしいです。ちなみにトリミングした1円玉は左上のだったのでそれもわかったんですかね(笑)
それでは、実際にこの左上の1円玉を回転させたり大きくさせたりしてみます。
同様にテンプレートマッチングをすると、
[回転]
[拡大]
今回は確実に例が悪かったので回転による影響を受ける様子があまり見られませんでしたが、拡大の影響は一目瞭然です。
最も類似度が高かった座標も他の1円玉の座標になっています。
一応ソースです
# 画像読み込み
ref_img = cv2.imread("./ref_coin.png")
tgt_img = cv2.imread("./tgt_coin.png")
# テンプレートマッチング
match_result = cv2.matchTemplate(ref_img, tgt_img, cv2.TM_CCOEFF_NORMED)
min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(match_result)
# 四角形を描画
h, w = tgt_img.shape[:2]
top_left = max_loc
bottom_right = (top_left[0] + h, top_left[1] + w)
cv2.rectangle(result,top_left, bottom_right, (0,0,255), 3)
// 画像読み込み
Mat ref_img = imread("./ref_coin.png");
Mat tgt_img = imread("./tgt_coin.png");
// 出力先変数定義
Mat match_result;
Point pos;
// テンプレートマッチング
matchTemplate(ref_img, tgt_img, match_result, CV_TM_CCOEFF_NORMED);
minMaxLoc(match_result, NULL, NULL, NULL, &pos);
// 四角形を描画
rectangle(ref_img, Rect(pos.x, pos.y, tgt_img.cols, tgt_img.rows), Scalar(0, 0, 255), 3);
これでテンプレートマッチングの元画像への依存度の高さがわかったのではないのでしょうか。
長ったらしい前置きはここまでとして本題に入っていきます！
アプローチ１
Phase Only Correlation (POC)、日本語で位相限定相関法というものがあります。
どんな手法かというと、それぞれの画像をフーリエ変換して、パワースペクトル画像を生成します。それらの相関をとった新たなパワースペクトル画像を生成し、最後に逆フーリエ変換することで相関の強かった座標のみがハイライトされたような画像が生成されます。
そのピーク値の座標を求めることで一番相関の強い座標を得ることができるという手法です。
このようにフーリエ変換することで、現画像は振幅(輝度)と位相(像の輪郭)に分解されます。こうすることで輝度情報のみを扱っているテンプレートマッチングよりも判断材料が１つ増え、外乱の影響を受けにくくなったことがわかります。
このサイトで可視化もしながらとても丁寧に説明してくれているので、ぜひ参考にしてください。
そして！長々と原理を説明したから１から実装するのかと思いきや、OpenCVで関数化されています。(https://docs.opencv.org/master/d7/df3/group__imgproc__motion.html)
ですので、それを使って楽に実装していきましょう！
実装
今回使用する画像です。
ランドルト環を位置合わせしていきましょう！
２枚目のように右に100ピクセル、下に100ピクセル平行移動させた画像を用意しました。
また、平行移動させる前に左に90度回転させた画像も用意しました (以降、平行移動画像、回転画像と呼ぶ)。
仕様上、OpenCVのPOCはグレースケール画像を渡さないと動いてくれません。IMREAD_GRAYSCALEをそっと添えてあげましょう。
また、画像の要素も小数である必要があります。そっとどうにかしてあげましょう。
それではコードを書いていきます！
# グレースケールとして画像読み込み
ref_img = cv2.imread("./Landolt.png", cv2.IMREAD_GRAYSCALE)
tgt_img = cv2.imread("./ShiftLandolt.png", cv2.IMREAD_GRAYSCALE)
# 位相限定相関法
(x, y), response = cv2.phaseCorrelate(ref_img .astype(np.float32), tgt_img .astype(np.float32))
# affine変換による平行移動
M = np.float32([[1, 0, -x], [0, 1, -y]])
aligned_img= cv2.warpAffine(tgt_img, M, (w, h))
print(f'Shift X: {x}\nShift Y: {y}')
>>Shift X: 99.99693842207353
Shift Y: 99.99310716347935
// グレースケールとして画像読み込み
Mat ref_img = imread("./Landolt.png", IMREAD_GRAYSCALE);
Mat tgt_img = imread("./ShiftLandolt.png", IMREAD_GRAYSCALE);
// 出力先変数定義
Mat aligned_img;
double response;
// Float型に変換したため、matの範囲は[0:1]となり、255で割る必要がある。
ref_img.convertTo(ref_img, CV_64FC1, 1.0 / 255.0);
tgt_img.convertTo(tgt_img, CV_64FC1, 1.0 / 255.0);
// 位相限定相関法
Point2d shift = phaseCorrelate(ref_img, tgt_img);
// affine変換による平行移動
Mat M = (Mat_<double>(2, 3) << 1.0, 0.0, -shift.x, 0.0, 1.0, -shift.y);
warpAffine(tgt_img, aligned_img, M, aligned_img.size());
printf("Shift X: %g\nShift Y: %g\n", shift.x, shift.y);
>>Shift X: 99.997
Shift Y: 99.9932
aligned_imgを出力したらこんな感じ
回転画像に対してPOCを適用したらこんな感じ
意外と位置合いましたね(笑)
アプローチ２
Rotation Invariant Phase Only Correlation (RIPOC)、日本語で回転不変位相限定相関法というものがあります。
字面からわかるように先ほどのPOCを改良した手法です。
これもかなりかいつまんで説明すると、フーリエ変換した後に、その画像を対数極座標に変換します。（対数極座標について）
こうすることで２次元座標(x, y)が対数極座標(ρ, θ) (ρ：ある点と原点との距離の対数、θ：x軸とのなす角)に変換されます。ですので、この極座標系では角度を2つの変数で表現できるようになります。
そこに先ほどのPOCを使って画像の移動量を取得することで、角度を補正するという手法です。目からうろこですね🥺
角度を補正できたら、後はまた平行移動の位置合わせを行って無事完了です。
それでは確実に簡単であろうPythonから実装していきます(笑)
実装（Python）
対象画像を回転画像として実装していきます。
それでは段階に分けて見ていきましょう！
まずはフーリエ変換まで
# 画像読み込み
ref_img = cv2.imread("./Landolt.png", cv2.IMREAD_GRAYSCALE)
tgt_img = cv2.imread("./RotateLandolt.png", cv2.IMREAD_GRAYSCALE)
# Numpy配列に変換
ref_img_gray = np.array(ref_img,dtype=np.float64)
tgt_img_gray = np.array(tgt_img,dtype=np.float64)
# 画像サイズ定義
h, w = ref_img_gray.shape
center = (w/2, h/2)
# 窓関数定義
hunning_y = np.hanning(h)
hunning_x = np.hanning(h)
hunning_w = hunning_y.reshape(h, 1)*hunning_x
# フーリエ変換
ref_img_fft = np.fft.fftshift(np.log(np.abs(np.fft.fft2(ref_img_gray*hunning_w))))
tgt_img_fft = np.fft.fftshift(np.log(np.abs(np.fft.fft2(tgt_img_gray*hunning_w))))
tgt_img_fftを出力してみると、
真っ白です。原因は、フーリエ変換後の画像はMat形式のFloat64なんて無視して普通に１以上の値を返してきたからです。
正規化して確認すると、
美しい画像が出力されました！
この画像の中心を削ると低周波数成分がカットされるハイパスフィルタの挙動になったりします。
興味がある方はこのサイトで詳細に説明されてたので、見てみてください。
では、対数極座標変換して、回転とスケールを導出してみます！
# 対数極座標変換 (lanczos法補間)
l = np.sqrt(w*w + h*h)
m = l/np.log(l)
flags = cv2.INTER_LANCZOS4 + cv2.WARP_POLAR_LOG
ref_img_pol = cv2.warpPolar(ref_img_fft, (w, h), center, m, flags)
tgt_img_pol = cv2.warpPolar(tgt_img_fft, (w, h), center, m, flags)
# 位相限定相関法
(x, y), response = cv2.phaseCorrelate(ref_img_pol, tgt_img_pol, hunning_w)
# 角度と大きさを導出
angle = y*360/h
scale = (np.e)**(x/m)
print(f'angle: {angle}\nscale:{scale}')
>>angle: 0.017308490013098775
scale:0.9986876396874065
tgt_img_polを出力してみると、
人間にはもうなにがなんだかわからないですが、printの出力結果から、角度と大きさが取れているのは確認できます。
求められた角度と大きさからaffine変換行列を求めて、変換させると、
後はまたPOCとかして位置合わせをすれば合いそう！
それではこのコードを平行移動画像と回転画像に対して実行していきます！
一応コード全文
# 画像読み込み
ref_img = cv2.imread("./Landolt.png", cv2.IMREAD_GRAYSCALE)
tgt_img = cv2.imread("./RotateLandolt.png", cv2.IMREAD_GRAYSCALE)
# Numpy配列に変換
ref_img_gray = np.array(ref_img,dtype=np.float64)
tgt_img_gray = np.array(tgt_img,dtype=np.float64)
# 画像サイズ定義
h, w = ref_img_gray.shape
center = (w/2, h/2)
# 窓関数定義
hunning_y = np.hanning(h)
hunning_x = np.hanning(h)
hunning_w = hunning_y.reshape(h, 1)*hunning_x
# フーリエ変換
ref_img_fft = np.fft.fftshift(np.log(np.abs(np.fft.fft2(ref_img_gray*hunning_w))))
tgt_img_fft = np.fft.fftshift(np.log(np.abs(np.fft.fft2(tgt_img_gray*hunning_w))))
# 対数極座標変換 (lanczos法補間)
l = np.sqrt(w*w + h*h)
m = l/np.log(l)
flags = cv2.INTER_LANCZOS4 + cv2.WARP_POLAR_LOG
ref_img_pol = cv2.warpPolar(ref_img_fft, (w, h), center, m, flags)
tgt_img_pol = cv2.warpPolar(tgt_img_fft, (w, h), center, m, flags)
# 位相限定相関法
(x, y), response = cv2.phaseCorrelate(ref_img_pol, tgt_img_pol, hunning_w)
# affine変換による平行移動
angle = y*360/h
scale = (np.e)**(x/m)
M = cv2.getRotationMatrix2D(center, angle, scale)
tgt_img_affine = cv2.warpAffine((tgt_img_gray), M, (w, h))
# 位相限定相関法
(x, y), response = cv2.phaseCorrelate(ref_img_gray, tgt_img_affine)
#位相限定相関法の返り値から画像を生成
M[0][2] -= x
M[1][2] -= y
dst = cv2.warpAffine(tgt_img, M, (w, h))
平行移動画像の出力結果
print(f'Shift X: {x}\nShift Y: {y}\nangle: {angle}\nscale:{scale}')
>>Shift X: 99.93655082583143
Shift Y: 99.88039360791896
angle: 0.017308490013098775
scale:0.9986876396874065
回転画像の出力結果
print(f'Shift X: {x}\nShift Y: {y}\nangle: {angle}\nscale:{scale}')
>>Shift X: -98.93914068317139
Shift Y: 100.10505647091378
angle: -89.94978000131137
scale:1.0002553379376475
両画像ともにかなりの精度で合っています（主観）。回転量もいい具合に補正できていますね。
これでみんなの視力が良くなるはずです。
また、この主観をカバーするために、ある座標の輝度値の差分をとったりして精度を評価するといった手段など、様々なものがあります。
実装（C++）
上記のPythonコードを参考に、C++で実装していきます。
まず、準備としてフーリエ変換をする関数を用意します。
フーリエ変換は、このサイトを参考に実装しました。
Mat fft(Mat src) {
// 変数定義
Mat padded;
int m = getOptimalDFTSize(src.rows);
int n = getOptimalDFTSize(src.cols);
// 入力画像を中央に配置。周囲は0でパディングする。
copyMakeBorder(src, padded, 0, m - src.rows, 0, n - src.cols, BORDER_CONSTANT, Scalar::all(0));
// フーリエ変換画像の入れ子定義
Mat planes[] = { Mat_<double>(padded), Mat::zeros(padded.size(), CV_64F) };
Mat complexI;
merge(planes, 2, complexI);
// フーリエ変換
dft(complexI, complexI);
// 対数に変換
split(complexI, planes);
magnitude(planes[0], planes[1], planes[0]);
Mat magI = planes[0];
log(magI, magI);
// 行・列が奇数の時に、偶数に変換するためトリミングする
magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
// 以下で象限の入れ替え
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy)); // 左上（第二象限）
Mat q1(magI, Rect(cx, 0, cx, cy)); // 右上（第一象限）
Mat q2(magI, Rect(0, cy, cx, cy)); // 左下（第三象限）
Mat q3(magI, Rect(cx, cy, cx, cy)); // 右下（第四象限）
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
// 正規化
normalize(magI, magI, 0, 1, NORM_MINMAX, -1);
return magI;
}
大体のやっていることはコメントに書いておきました。
さて、フーリエ変換された画像を返す関数が作成できたのでRIPOCを実装していきます！
// 変数定義
const int height = ref_img.rows;
const int width = ref_img.cols;
Point center = Point(width / 2, height / 2);
const Size window_size(height, height);
Mat aligned_img = Mat::zeros(tgt_img.size(), CV_64FC1);
// グレースケールに変換＋Folat型に変換
ref_img.convertTo(ref_img, CV_64FC1, 1.0 / 255.0);
tgt_img.convertTo(tgt_img, CV_64FC1, 1.0 / 255.0);
// フーリエ変換
Mat ref_img_fft = fft(ref_img);
Mat tgt_img_fft = fft(tgt_img);
// 対数極座標変換された画像の入れ子
Mat bg_ref = Mat::zeros(ref_img.size(), CV_64FC1);
Mat bg_tgt = Mat::zeros(tgt_img.size(), CV_64FC1);
// 対数極座標変換 (lanczos法補間)
double l = sqrt(height * height + width * width);
double m = l / log(l);
int flags = INTER_LANCZOS4 + WARP_POLAR_LOG;
warpPolar(ref_img_fft, bg_ref, window_size, center, m, flags);
warpPolar(tgt_img_fft, bg_tgt, window_size, center, m, flags);
// 位相限定相関法
Point2d pt = phaseCorrelate(bg_ref, bg_tgt);
// 角度と大きさを導出
float rotate = pt.y * 360 / width;
float scale = exp(pt.x / m);
// アフィン変換係数行列取得
Mat M = getRotationMatrix2D(center, rotate, scale);
Mat rotated_img = Mat::zeros(tgt_img.size(), CV_64FC1);
// 角度と大きさを補正
warpAffine(tgt_img, rotated_img, M, tgt_img.size());
rotated_img.convertTo(rotated_img, CV_64FC1, 1.0 / 255.0);
// 位相限定相関法
pt = phaseCorrelate(ref_img, rotated_img);
// affine変換行列を調整してaffine変換
M.at<double>(0, 2) -= pt.x;
M.at<double>(1, 2) -= pt.y;
warpAffine(tgt_img, aligned_img, M, tgt_img.size());
これで完成です！aligned_imgの出力を見てみましょう！
平行移動画像の出力結果
// ripoc.cpp内
printf("Rotation: %g\nScale: %g\n", rotate, scale);
// poc.cpp内
printf("Shift X: %g\nShift Y: %g\n", shift.x, shift.y);
>>Rotation: 0.011399
Scale: 1.00014
Shift X: 98.1662
Shift Y: 100.914
回転画像の出力結果
// ripoc.cpp内
printf("Rotation: %g\nScale: %g\n", rotate, scale);
// poc.cpp内
printf("Shift X: %g\nShift Y: %g\n", shift.x, shift.y);
>>Rotation: -89.9924
Scale: 1.00039
Shift X: -99.0014
Shift Y: 100.003
C++でも両画像ともにかなりの精度で合っています。
それでは、どうせなので速度比較を行っていきたいと思います！！
C++とPythonの速度比較
今回使用した画像サイズ：(800[px], 800[px], 3[ch])
PCスペック：メモリ：16[GB]
CPU：i5-9300H
ベースクロック：2.4[GHz]
計算時間はfor文で100回実行してその平均としました。
| C++計算時間[ms] | Python計算時間[ms] | |
|---|---|---|
| POC（位相限定相関法） | 34.036662 | 18.470275 |
| RIPOC（回転不変位相限定相関法） | 135.769934 | 133.652456 |
これは意外です。
差が顕著だったPOCをステップ実行して計算時間を調査したところ、C++ではcv::phaseCorrelate関数の計算に30[ms]前後かかっていました。この時点でPythonを下回っています。
私の実装ではRIPOCで2回cv::phaseCorrelate関数を使っているので、それ以外の計算はC++の方が早いと単純には考えられます。
まぁ最大のボトルネックは私の書き方にあると思いますが←
さいごに
テンプレートマッチング以外にもこんなマッチングの手法があるよという紹介でした
他にもcv::findtransformECCという関数とかもありますので、気になる方はそちらもチェック！
ちなみに、ランドルト環の画像で上記関数を実行したところ、無相関過ぎてダメというエラーメッセージが表示され、他の画像ではうまくいったので、本当に選んだ画像が良くなかったです💦
今回は同じサイズのランドルト環に対して実行しましたが、cv::findtransformECCは回転以外にもスケール変化に対応しているので、精度が気になる方は是非実験してみてください
Pythonの方がOpenCV関数の実行速度が速いという意外な結果になりました
これからは脳死でC++の方が早いよというのはやめようと思います
それでは、楽しい位相限定相関法・回転不変位相限定相関法ライフをenjoyしてください！