長きにわたって OpenCV の Stitch モジュールについてコードリーディングを進めています。

こちらの記事の続きです。
https://qiita.com/itoshogo3/items/0b0869ae3a119964f036

全体の概要はこちらです。
https://qiita.com/itoshogo3/items/7a3279668b24008a3761

概要

今回は各画像で検出された特徴点を見比べていきます。もし、ある画像に写った特徴点が別の画像にも写り込んでいたら、それらの画像は隣り合っているとみなすことが出来ます。逆に特徴点を見比べても一致するものがなければ、それらの画像は離れて撮影されたものと考えることが出来ます。こんな感じで重なりを持つ画像のペアを探していきます。

特徴点のマッチング

下図のように異なる画像から抽出した特徴点は、抽出される点の数も違いますし、全ての点が一対一に対応するわけではありません。この中から同じ場所と思われる特徴点のペアを抜き出していくのが、特徴点のマッチング処理になります。

前回、各特徴点につき特徴量と呼ばれる多次元のベクトル値を計算しました。この時の次元数というのが特徴量記述のアルゴリズムによって異なるのですが、だいたい 32 次元くらいは必要で、ものによっては 128 次元とかもっと多次元なベクトルになります。この特徴量ベクトル同士の距離を計算して、近ければ近いほど同じ点であるとみなせるわけです。

２つの画像の特徴点同士を総当たりで比較すれば、確実に最も似ていると思われる点のペア（特徴量ベクトルの距離が近いペア）を見つけることができそうです。この総当たりで比較する方法はブルートフォースと呼ばれていますが、比較するものが多次元のベクトルということもあって、明らかに計算量が多くて時間もかかりそうです。そこで stitch モジュールでは全ての点で確実とまではいかなくても、まあまあな精度でマッチングできれば良いという方針で効率的にマッチングを行っています。

このようなたくさんのベクトルデータの中から、与えられたベクトルに近いと思われる値を探す問題は、近似最近傍探索問題として知られており、これまた奥が深い分野となっております。ここで話しきれるほど簡単な内容ではないので詳細な説明は割愛しますが、ざっくり言うとデータ構造を工夫して捜索範囲を小さくしているという感じです。高速ですが最近傍にたどり着かないことがあります。ちなみに OpenCV には FLANN という近似最近傍探索のモジュールが用意されており、stitch モジュールもこの手法で実装されております。下図はマッチング処理の結果です。だいたい正しくマッチングできているのですが、間違ったペアもあることに注目してください。

間違ったマッチング

特徴量を比較する方法は、特徴点の周りの局所的な情報を見比べただけなので、もし画像の中にたまたま似たような部分があると間違ってマッチングしてしまいます。これは困ったことです。下図は間違ってマッチングしたペアの例です。

ただし、間違ったペアというのは他の正しくマッチングしたペアと比べて、下図のように明らかに変な移動の仕方をするのです。この特性を利用すると間違ってマッチングしたペアを取り除くことが出来ます。ちなみに正しくマッチングしたペアのことをインライア（ inlier ）、間違ってマッチングしたペアをアウトライア ( outlier ) と呼びます。

もともと入力される画像は、その場を移動せずに向きだけを変えて撮影したものである必要がありました。もしそのような制約のもとで撮影された画像であれば、全ての特徴点はおおよそ下図のように規則的な動きをするはずです。

画像 $i$ の $k$ 番目のマッチングペアの特徴点の２次元座標を $u_i^k=( x_i^k, y_i^k ) $ とすると、上図の関係は下記のような数式で表せます。

s\begin{pmatrix}x_i^k\\y_i^k\\1\end{pmatrix}=H_{ij}\begin{pmatrix}x_j^k\\y_j^k\\1\end{pmatrix}

$H$ はホモグラフィー行列と呼ばれる $3 \times 3$ の行列で、ある画像の点を別の画像の点へと移す行列となっています。画像 $j$ 上の特徴点に $H_{ij}$ を作用させると、画像 $i$ 上のマッチングした特徴点の座標が得られるのです。正しくマッチングした特徴点ペアは、みなこちらの式の条件を満たします。ちなみにここでさりげなくパラメータ $s$ を導入しました。こちらは任意のスケール係数となります。ホモグラフィー行列は定数倍の不定性が残るのですが、このことを明示的に示すために、$s$ は慣習的に $H$ からくくり出されます。

ちなみにホモグラフィー行列といえば、基本的には同一平面上の点の画像間の射影変換を表すものです。今回の状況では、球面上の点を無限遠の位置にある平面上の点とみなすことで、画像間の点の移動をホモグラフィー行列で表すことが出来るのです。このことは OpenCV のチュートリアルでも少し触れられています。

ホモグラフィー行列のコンセプト
https://docs.opencv.org/master/d9/dab/tutorial_homography.html

ホモグラフィー行列はマッチングした特徴点ペアが４組以上あれば、それらの座標から次に示す手順で計算することが出来ます。

ホモグラフィー行列の計算

まずは $H_{ij}$ を下記のように要素で表します。

H_{ij}=\begin{pmatrix}\,
h_{11} & h_{12} & h_{13} \\
h_{21} & h_{22} & h_{23} \\
h_{31} & h_{32} & 1
\end{pmatrix}

すると、先ほどのホモグラフィーの条件式は下記のように書き直せます。

s\begin{pmatrix}x_i^k\\y_i^k\\1\end{pmatrix}=\begin{pmatrix}\,
h_{11} & h_{12} & h_{13} \\
h_{21} & h_{22} & h_{23} \\
h_{31} & h_{32} & 1
\end{pmatrix}
\begin{pmatrix}x_j^k\\y_j^k\\1\end{pmatrix}

いったん３つの式にして、

sx_i^k=h_{11}x_j^k+h_{12}y_j^k+h_{13}\\
sy_i^k=h_{21}x_j^k+h_{22}y_j^k+h_{23}\\
s=h_{31}x_j^k+h_{32}y_j^k+1\\

ごにょごにょして $s$ を消します。

h_{11}x_j^k+h_{12}y_j^k+h_{13}-h_{31}x_j^kx_i^k-h_{32}y_j^kx_i^k=x_i^k\\
h_{21}x_j^k+h_{22}y_j^k+h_{23}-h_{31}x_j^ky_i^k-h_{32}y_j^ky_i^k=y_i^k\\

こちらの２つの式を下記のように行列の形にします。

\begin{pmatrix}
x_j^k & y_j^k & 1 & 0 & 0 & 0 & -x_j^kx_i^k & -y_j^kx_i^k \\
0 & 0 & 0 & x_j^k & y_j^k & 1 & -x_j^ky_i^k & -y_j^ky_i^k 
\end{pmatrix} 
\begin{pmatrix}
h_{11} \\ h_{12} \\ h_{13} \\
h_{21} \\ h_{22} \\ h_{23} \\
h_{31} \\ h_{32}
\end{pmatrix}=
\begin{pmatrix}
x_i^k \\ y_i^k 
\end{pmatrix}

このように一組のマッチングペアから２つの式が得られました。ここで、マッチングペアの中から適当に４組選びます。

選ばれた特徴点のペアを $(u_i^1, u_j^1), (u_i^2, u_j^2), (u_i^3, u_j^3), (u_i^4, u_j^4)$ とすると、下記のように８つの式が導けて行列の形で表すことが出来ます。

\begin{pmatrix}
x_j^1 & y_j^1 & 1 & 0 & 0 & 0 & -x_j^1x_i^1 & -y_j^1x_i^1 \\
0 & 0 & 0 & x_j^1 & y_j^1 & 1 & -x_j^1y_i^1 & -y_j^1y_i^1 \\
x_j^2 & y_j^2 & 1 & 0 & 0 & 0 & -x_j^2x_i^2 & -y_j^2x_i^2 \\
0 & 0 & 0 & x_j^2 & y_j^2 & 1 & -x_j^2y_i^2 & -y_j^2y_i^2 \\
x_j^3 & y_j^3 & 1 & 0 & 0 & 0 & -x_j^3x_i^3 & -y_j^3x_i^3 \\
0 & 0 & 0 & x_j^3 & y_j^3 & 1 & -x_j^3y_i^3 & -y_j^3y_i^3 \\
x_j^4 & y_j^4 & 1 & 0 & 0 & 0 & -x_j^4x_i^4 & -y_j^4x_i^4 \\
0 & 0 & 0 & x_j^4 & y_j^4 & 1 & -x_j^4y_i^4 & -y_j^4y_i^4
\end{pmatrix} 
\begin{pmatrix}
h_{11} \\ h_{12} \\ h_{13} \\
h_{21} \\ h_{22} \\ h_{23} \\
h_{31} \\ h_{32}
\end{pmatrix}=
\begin{pmatrix}
x_i^1 \\ y_i^1 \\
x_i^2 \\ y_i^2 \\
x_i^3 \\ y_i^3 \\
x_i^4 \\ y_i^4 
\end{pmatrix}

未知数が８つ $( h_{11} \sim h_{32} )$ の連立方程式の形に持ってこれたので、こちらを解けばホモグラフィー行列の各要素が求まります。

というわけで４組のマッチングペアからホモグラフィー行列が計算できました。実際にはそれ以上の数のマッチングペアを使って、再投影誤差を最小化するという小難しい方法がとられたりするのですが、話が長くなるので割愛します。

RANSAC によるアウトライヤの除去

ホモグラフィー行列はマッチした特徴点ペアを何組か選んで計算することになりますが、この計算に間違ったマッチングペアを選んでしまうとホモグラフィー行列が上手く計算できず、真の値とは違う結果になってしまいます。

そこで一回で計算して決めるのではなく、何回もランダムでマッチングペアを選んでホモグラフィー行列を計算し、より誤差の少ない計算結果を採用することにします。変なペアを選んで計算されたホモグラフィー行列は、各特徴点ペアに当てはめるてみると、てんでバラバラな動きになるため、計算結果がおかしなことになっているのはすぐにわかるのです。ざっくりですが一般的な手順は次のようになります。

手順	説明
1	特徴点ペアの中から H 行列の計算に必要な数だけペアをランダムに抽出。
2	1 で選んだペアだけで H 行列を計算。
3	2 で計算した H 行列を全てのペアに当てはめる。
4	3 の結果、アウトライアっぽい動きをするペアがそれほど多くなければ、正解候補とする。
5	4 で得た正解候補の H 行列をインライアのみを使って再評価する。
6	1 ～ 5 を繰り返して最も性能の良い H 行列を採用する。

このような手法は RANSAC (Random Sample Consensus) と呼ばれており、特に統計の分野ではよく使われているそうです。最終的なマッチング処理の結果は下図のようになり、アウトライアが取り除かれていることが確認できます。

コードリーディング

今回解説した場所はソースコードでいうと Stitcher::matchImages() の後半部分です。その中でも下記の BestOf2NearestMatcher::operator () でやっている処理が、この記事で解説するメインの部分です。見かけはたった１行の処理ですが、この中ではまあまあ複雑なことをやっています。

(*features_matcher_)(features_, pairwise_matches_, matching_mask_);

そもそもパッと見ただけでは何をやっているかわかりにくいかもしれませんが、実はこれ、演算子をオーバーロードしていて、FeaturesMatcher で定義されている operator () を実行しているのです。パノラマモードの場合、メンバ変数である features_matcher_ には BestOf2NearestMatcher クラスがセットされていますが、こちらには operator () は定義されていません。というわけでその親クラスである FeaturesMatcher を見てみると、そちらの方に定義がありました。が、FeaturesMatcher には operator () が２つも定義されています。

CV_WRAP_AS(apply) void operator () (
    const ImageFeatures &features1,  /* 入力の特徴点１ */
    const ImageFeatures &features2,  /* 入力の特徴点２ */
    CV_OUT MatchesInfo& matches_info /* マッチング情報がここに出力される */ ) 
{ 
    match(features1, features2, matches_info); 
}

CV_WRAP_AS(apply2) void operator () (
    const std::vector<ImageFeatures> &features,         /* 入力の特徴点の配列 */ 
    CV_OUT std::vector<MatchesInfo>  &pairwise_matches, /* マッチング情報の配列 */ 
    const cv::UMat                   &mask = cv::UMat() /* マスク */ );

入力される引数によって処理を呼び分けているわけですが、まずは apply2 の方の operator () が呼ばれます。この中で MatchPairsBody というまたややこしいクラスが作られ、オーバーロードされた operator () が実行されて、その後にもうひとつの apply の方の operator () が呼ばれて、、、みたいな感じでとにかくソースコード上をたらい回しにされます。というわけでざっくりとフローを下記にまとめました。

FeaturesMatcher::operator () ( apply2 の方 )
 | 
 |- MatchPairsBody::operator ()
     | 
     |- for loop 全画像総当たり --------------------------
     |   |
     |   |- FeaturesMatcher::operator () ( apply の方 )
     |       |
     |       |- BestOf2NearestMatcher::match() 
     |  
     |----------------------------------------------------

とりあえず分かって欲しい事は単純で、最終的に BestOf2NearestMatcher::match() という関数が、入力した全画像に対して総当たりで実行されるということです。こちらの中身はざっくり下記のような流れとなっていますが、細かなノウハウが詰まっており、大変重要な読み応えのある部分になります。

BestOf2NearestMatcher::match() 
 |
 |- CpuMatcher::operator() 
 |    |
 |    |- CpuMatcher::match()
 |         |
 |         |- FlannBasedMatcher::knnMatch()
 |
 |- （ マッチした点を中心がゼロとなる座標に変換 ）
 | 
 |- findHomography() 
 |  
 |- （ 結果からチェックしてダメ押しでもう一回計算すべきか判断 ）
 | 
 |- findHomography() ( インライアだけで実行 )

特徴点のマッチング

BestOf2NearestMatcher::match() の最初で次のコードが実行されています。またしても演算子をオーバーロードする形です。この impl_ はコンストラクタで設定されており、 GPU 処理 (CUDA) に対応してなければ CpuMatcher がセットされています。

 (*impl_)(features1, features2, matches_info);

というわけでその中身がこちらです。

void CpuMatcher::match(
    const ImageFeatures &features1,   /* 入力の特徴点１ */ 
    const ImageFeatures &features2,   /* 入力の特徴点２ */
    MatchesInfo&        matches_info  /* マッチング情報出力 */ )
{
    matches_info.matches.clear();

    Ptr<cv::DescriptorMatcher> matcher;
    Ptr<flann::IndexParams> indexParams = makePtr<flann::KDTreeIndexParams>();
    Ptr<flann::SearchParams> searchParams = makePtr<flann::SearchParams>();
    if (features2.descriptors.depth() == CV_8U)
    {
        indexParams->setAlgorithm(cvflann::FLANN_INDEX_LSH);
        searchParams->setAlgorithm(cvflann::FLANN_INDEX_LSH);
    }
    matcher = makePtr<FlannBasedMatcher>(indexParams, searchParams);

    std::vector< std::vector<DMatch> > pair_matches;
    MatchesSet matches;

    // Find 1->2 matches
    matcher->knnMatch(features1.descriptors, features2.descriptors, pair_matches, 2);
    for (size_t i = 0; i < pair_matches.size(); ++i)
    {
        if (pair_matches[i].size() < 2)
            continue;
        const DMatch& m0 = pair_matches[i][0];
        const DMatch& m1 = pair_matches[i][1];
        if (m0.distance < (1.f - match_conf_) * m1.distance)
        {
            matches_info.matches.push_back(m0);
            matches.insert(std::make_pair(m0.queryIdx, m0.trainIdx));
        }
    }

    // Find 2->1 matches
    pair_matches.clear();
    matcher->knnMatch(features2.descriptors, features1.descriptors, pair_matches, 2);
    for (size_t i = 0; i < pair_matches.size(); ++i)
    {
        if (pair_matches[i].size() < 2)
            continue;
        const DMatch& m0 = pair_matches[i][0];
        const DMatch& m1 = pair_matches[i][1];
        if (m0.distance < (1.f - match_conf_) * m1.distance)
            if (matches.find(std::make_pair(m0.trainIdx, m0.queryIdx)) == matches.end())
                matches_info.matches.push_back(DMatch(m0.trainIdx, m0.queryIdx, m0.distance));
    }
}

knnMatch() でマッチング処理を実行します。第３引数に各特徴点のマッチングの結果が格納されますが、この時に第４引数で指定した数だけマッチング候補が得られます。ここでは候補を２つに絞っています。

void knnMatch( 
    InputArray queryDescriptors, 
    InputArray trainDescriptors,
    CV_OUT std::vector<std::vector<DMatch> >& matches, 
    int k,
    InputArray mask=noArray(), 
    bool compactResult=false ) ;

knnMatch() を実行した後の下記のコードがマッチングペアの選別を行っているところです。distance には特徴量同士の距離が入っていて、この距離が近いほうがよくマッチしていることになります。単にベクトル間の距離を見ているのではなく、第１候補との距離が第２候補との距離よりもある程度小さくないとマッチングペアとして選ばれません。このように、より信頼できるペアだけを選別して間違ってマッチングしないようにしているわけです。

const DMatch& m0 = pair_matches[i][0];
const DMatch& m1 = pair_matches[i][1];
if (m0.distance < (1.f - match_conf_) * m1.distance)
{
    matches_info.matches.push_back(m0);
    matches.insert(std::make_pair(m0.queryIdx, m0.trainIdx));
}

ホモグラフィー行列の計算と外れ値の除去

マッチングしたペアの点群をもとにホモグラフィーを計算します。が、その前に中心の座標がゼロになるように、各点の座標から画像サイズの半分の値を引き算しています。その後に findHomography() という関数でホモグラフィー行列を計算しますが、なんとこの時にオプションで RANSAC を指定することにより、外れ値まで一挙に除外してくれるのです。すごく便利。

// マッチングポイント
Mat src_points(1, static_cast<int>(matches_info.matches.size()), CV_32FC2);
Mat dst_points(1, static_cast<int>(matches_info.matches.size()), CV_32FC2);

// 画像の中心がゼロになるように座標を変換。
for (size_t i = 0; i < matches_info.matches.size(); ++i)
{
    const DMatch& m = matches_info.matches[i];
    Point2f p = features1.keypoints[m.queryIdx].pt;
    p.x -= features1.img_size.width * 0.5f;
    p.y -= features1.img_size.height * 0.5f;
    src_points.at<Point2f>(0, static_cast<int>(i)) = p;
    p = features2.keypoints[m.trainIdx].pt;
    p.x -= features2.img_size.width * 0.5f;
    p.y -= features2.img_size.height * 0.5f;
    dst_points.at<Point2f>(0, static_cast<int>(i)) = p;
}

// ホモグラフィー計算。
matches_info.H = findHomography(src_points, dst_points, matches_info.inliers_mask, RANSAC);

そして、ここからがちょっとしたノウハウ部分です。

// Find number of inliers
matches_info.num_inliers = 0;
for (size_t i = 0; i < matches_info.inliers_mask.size(); ++i)
    if (matches_info.inliers_mask[i])
        matches_info.num_inliers++;

// These coeffs are from paper M. Brown and D. Lowe. "Automatic Panoramic Image Stitching
// using Invariant Features"
matches_info.confidence = matches_info.num_inliers / (8 + 0.3 * matches_info.matches.size());

// Set zero confidence to remove matches between too close images, as they don't provide
// additional information anyway. The threshold was set experimentally.
matches_info.confidence = matches_info.confidence > 3. ? 0. : matches_info.confidence;

// Check if we should try to refine motion
if (matches_info.num_inliers < num_matches_thresh2_)
    return;

//        //
//  中略  // 
//        //

// Rerun motion estimation on inliers only
matches_info.H = findHomography(src_points, dst_points, RANSAC);

まずマッチングしたペアの数を数えて、それをもとに confidence なる値を計算しています。confidence とは「自信」とか「信用」とかいう意味です。扱いにくい変な画像が混じっているとこの先の処理での精度が落ちてしまうので、そういった画像はこの時点で間引くために confidence 値を低めに設定してあるのです。また高ければ良いというわけでもありません。特に confidence が 3 を超える場合はとても似ている画像（ほとんどカメラを回転させずに撮影した画像）とみなすことができ、追加の情報は与えてくれないので、この時点で confidence をゼロにして確実に省かれるようにしています。

続いてその先ですが、ペア数がある閾値を下回ったらそこで終了させて、もし上回っていたらインライアだけでもう一回ホモグラフィーを計算しています。このような細かい工夫でアプリとしての品質を上げているのでした。

画像の選別

BestOf2NearestMatcher::match() の最後の部分が下記のコードです。ここではうまくマッチング処理が出来なかった画像を間引いて次の処理に引き渡します。

// 確実に使える画像だけを残す。
indices_ = detail::leaveBiggestComponent(features_, pairwise_matches_, (float)conf_thresh_);

// メンバ変数を更新する。
std::vector<UMat> seam_est_imgs_subset;
std::vector<UMat> imgs_subset;
std::vector<Size> full_img_sizes_subset;
for (size_t i = 0; i < indices_.size(); ++i)
{
    imgs_subset.push_back(imgs_[indices_[i]]);
    seam_est_imgs_subset.push_back(seam_est_imgs_[indices_[i]]);
    full_img_sizes_subset.push_back(full_img_sizes_[indices_[i]]);
}
seam_est_imgs_ = seam_est_imgs_subset;
imgs_ = imgs_subset;
full_img_sizes_ = full_img_sizes_subset;

// 残っている画像が２枚以上なかったら失敗となる。
if ((int)imgs_.size() < 2)
{
    LOGLN("Need more images");
    return ERR_NEED_MORE_IMGS;
}

leaveBiggestComponent() で合成に使える画像を選別し、画像のインデックスを返しています。入力された画像の中には、他の画像とマッチングする特徴点が少なすぎて合成に使えない画像や、写っているモノがほとんど同じで合成するための追加の情報を与えてくれない画像はいらない子なのです。leaveBiggestComponent() の結果をもとにメンバ変数の seam_est_imgs_, imgs_, full_img_sizes_ を更新しています。

以上のような流れで入力された各画像の位置関係がわかるようになるのでした。

あとがき

目次
https://qiita.com/itoshogo3/items/7a3279668b24008a3761
前の記事（特徴点の話）
https://qiita.com/itoshogo3/items/0b0869ae3a119964f036
次の記事（カメラ幾何の話）
https://qiita.com/itoshogo3/items/5fd652b72b8c424f076b