はじめに

本日はふと目にしたIssueを解説します。
本記事はOpenCV アドベントカレンダー 2020、第2日目の記事です。
他の記事は目次を参照して下さい。

TL;DR

Homographyの自由度は8
対応点4組が必要

はじまりの物語

Proper extrapolation with findHomography · Issue #18254 · opencv/opencv

cv::findHomography関数に対応点3点入れてうまく行かないんだけれど！？という~~明後日の~~issue
これに対し、「そりゃそうでしょ。自由度が8なんだから」とリプライが付き、対応点が4組未満の場合に実行時エラーが投げられるPRが提出、マージされました。

calib3d: check for minimum corresponding points to calculate homography by aitikgupta · Pull Request #18444 · opencv/opencv

非常にたくさん言いたいことがあるんです。が、ここでは割愛します。
今日の主題は、「何で今まで実行時エラーにならなかったの？」という点です。

数学的側面

Homographyは、以下のように、ある平面上の点を、別の平面上の点に投影する行列のことです
一意に座標が決まるため、視点変換などでよく使われます。
以下の式で表される様に、$(x,y)$を$(x^\prime,y^\prime)$に投影します。

\begin{pmatrix}
H_{11} & H_{12} & H_{13} \\
H_{21} & H_{22} & H_{23} \\
H_{31} & H_{32} & 1
\end{pmatrix}
\begin{pmatrix}
x \\
y \\
1
\end{pmatrix}
=
s
\begin{pmatrix}
x^\prime \\
y^\prime \\
1
\end{pmatrix}

このときHomographyはスケール不変の性質があるので、最後の要素を1として固定し、変数$s$にスケールの違いを集約しました。
前述の行列の掛け算は以下のように3本の方程式に書き換えらます。

\left\{
\begin{array}{l}
H_{11}x+H_{12}y+H_{13}=sx^\prime\\
H_{21}x+H_{22}y+H_{23}=sy^\prime\\
H_{31}x+H_{32}y+1=s
\end{array}
\right.

以上の3本の方程式は$s$に関する従属性があるので、書き換えて、下記になります。

\left\{
\begin{array}{l}
H_{11}x+H_{12}y+H_{13}=(H_{31}x+H_{32}y+1)x^\prime\\
H_{21}x+H_{22}y+H_{23}=(H_{31}x+H_{32}y+1)y^\prime\\
\end{array}
\right.

一部展開して、2本の方程式を得ます。

\left\{\begin{array}{c}H_{11}x+H_{12}y+h_{13}-H_{31}xx^\prime-H_{32}yx^\prime&=x^\prime\\H_{21}x+H_{22}y+H_{23}-H_{31}xy^\prime-H_{32}yy^\prime&=y^\prime\\\end{array}\right.

この連立方程式を行列演算の形に書き換えて

\left(\begin{array}{cccccccc}x&y&1&0&0&0&-xx^\prime&-yx^\prime\\0&0&0&x&y&1&-xy^\prime&-yy^\prime\\\end{array}\right)\left(\begin{array}{c}H_{11}\\H_{12}\\H_{13}\\H_{21}\\H_{22}\\H_{23}\\H_{31}\\H_{32}\\\end{array}\right)=\left(\begin{array}{c}x^\prime\\y^\prime\\\end{array}\right)

対応点1組で、独立した方程式が2本得られます。
Homographyのパラメータは$H_{11}$から$H_{32}$まで8つあるので、$8/2=4$で4組の対応点が必要になる

\left(\begin{array}{cccccccc}
x_1&y_1&1&0&0&0&-x_1x_1^\prime&-y_1x_1^\prime\\
0&0&0&x_1&y_1&1&-x_1y_1^\prime&-y_1y_1^\prime\\
x_2&y_2&1&0&0&0&-x_2x_2^\prime&-y_2x_2^\prime\\
0&0&0&x_2&y_2&1&-x_2y_2^\prime&-y_2y_2^\prime\\
x_3&y_3&1&0&0&0&-x_3x_3^\prime&-y_3x_3^\prime\\
0&0&0&x_3&y_3&1&-x_3y_3^\prime&-y_3y_3^\prime\\
x_4&y_4&1&0&0&0&-x_4x_4^\prime&-y_4x_4^\prime\\
0&0&0&x_4&y_4&1&-x_4y_4^\prime&-y_4y_4^\prime\\
\end{array}\right)
\left(\begin{array}{c}H_{11}\\H_{12}\\H_{13}\\H_{21}\\H_{22}\\H_{23}\\H_{31}\\H_{32}\\\end{array}\right)=
\left(\begin{array}{c}
x_1^\prime\\
y_1^\prime\\
x_2^\prime\\
y_2^\prime\\
x_3^\prime\\
y_3^\prime\\
x_4^\prime\\
y_4^\prime\\
\end{array}\right)

3点しか対応点がない状態は、前述の数式の左側の行列が以下のように6行8列になり、逆行列を求めることができなくなります。

6
\begin{array}{c}
8\\
\left(\begin{array}{cccccccc}
x_1&y_1&1&0&0&0&-x_1x_1^\prime&-y_1x_1^\prime\\
0&0&0&x_1&y_1&1&-x_1y_1^\prime&-y_1y_1^\prime\\
x_2&y_2&1&0&0&0&-x_2x_2^\prime&-y_2x_2^\prime\\
0&0&0&x_2&y_2&1&-x_2y_2^\prime&-y_2y_2^\prime\\
x_3&y_3&1&0&0&0&-x_3x_3^\prime&-y_3x_3^\prime\\
0&0&0&x_3&y_3&1&-x_3y_3^\prime&-y_3y_3^\prime\\
\end{array}\right)
\end{array}

結局、3点しか無いと、大きさが6行8列である行列の逆行列を求めることになり、ランク落ちで逆行列が求められず、前述のパラメータが求められない。
なので、Homographyは4組対応点が必要。
となると、どこかで「ランク落ち」で逆行列が求められない、という実行時エラーが起きてるのではなかろうか？
その謎を探るため、探検隊はGithubの奥地へと飛んだ

実装チェック

OpenCVの実装を追ってみよう。

calib3d.hpp

CV_EXPORTS_W Mat findHomography( InputArray srcPoints, InputArray dstPoints,
                                 int method = 0, double ransacReprojThreshold = 3,
                                 OutputArray mask=noArray(), const int maxIters = 2000,
                                 const double confidence = 0.995);

引数を省略した場合、methodには0が使われる。

fundam.cpp

cv::Mat cv::findHomography( InputArray _points1, InputArray _points2,
                            int method, double ransacReprojThreshold, OutputArray _mask,
                            const int maxIters, const double confidence)
{
    CV_INSTRUMENT_REGION();

    if (method >= USAC_DEFAULT && method <= USAC_MAGSAC)
        return usac::findHomography(_points1, _points2, method, ransacReprojThreshold,
            _mask, maxIters, confidence);

ここの、USAC_DEFAULTの値は32で、全然別のアプローチでHomographyを求めているので、本記事では割愛。

calib3d.hpp

//! type of the robust estimation algorithm
enum { LMEDS  = 4,  //!< least-median of squares algorithm
       RANSAC = 8,  //!< RANSAC algorithm
       RHO    = 16, //!< RHO algorithm
       USAC_DEFAULT  = 32, //!< USAC algorithm, default settings
       USAC_PARALLEL = 33, //!< USAC, parallel version
       USAC_FM_8PTS = 34,  //!< USAC, fundamental matrix 8 points
       USAC_FAST = 35,     //!< USAC, fast settings
       USAC_ACCURATE = 36, //!< USAC, accurate settings
       USAC_PROSAC = 37,   //!< USAC, sorted points, runs PROSAC
       USAC_MAGSAC = 38    //!< USAC, runs MAGSAC++
     };

実装を更に追っていくと、今回のエラーメッセージが見つかる

fundam.cpp

    for( int i = 1; i <= 2; i++ )
    {
        Mat& p = i == 1 ? points1 : points2;
        Mat& m = i == 1 ? src : dst;
        npoints = p.checkVector(2, -1, false);
        if( npoints < 0 )
        {
            npoints = p.checkVector(3, -1, false);
            if( npoints < 0 )
                CV_Error(Error::StsBadArg, "The input arrays should be 2D or 3D point sets");
            if( npoints == 0 )
                return Mat();
            convertPointsFromHomogeneous(p, p);
        }
        // Need at least 4 point correspondences to calculate Homography
        if( npoints < 4 )
            CV_Error(Error::StsVecLengthErr , "The input arrays should have at least 4 corresponding point sets to calculate Homography");
            // ↑ここのCV_Errorで明示的に実行時エラーが投げられる
        p.reshape(2, npoints).convertTo(m, CV_32F);
    }

ここで具体的にやってるのは斉次座標だったらもとに戻し、さらに対応点の組数をカウントしている。
npointsが組数を表し、4未満だと実行時エラーになる。
この実行時エラーが無いと、どこに行き着くのか。
我々探検隊は更に奥地へと進んだ

fundam.cpp

    if( ransacReprojThreshold <= 0 )
        ransacReprojThreshold = defaultRANSACReprojThreshold;

    Ptr<PointSetRegistrator::Callback> cb = makePtr<HomographyEstimatorCallback>();

    if( method == 0 || npoints == 4 )
    {
        tempMask = Mat::ones(npoints, 1, CV_8U);
        result = cb->runKernel(src, dst, H) > 0; // 今回のissueの状況だと、ここに辿り着く
    }
    else if( method == RANSAC )
        result = createRANSACPointSetRegistrator(cb, 4, ransacReprojThreshold, confidence, maxIters)->run(src, dst, H, tempMask);
    else if( method == LMEDS )
        result = createLMeDSPointSetRegistrator(cb, 4, confidence, maxIters)->run(src, dst, H, tempMask);
    else if( method == RHO )
        result = createAndRunRHORegistrator(confidence, maxIters, ransacReprojThreshold, npoints, src, dst, H, tempMask);
    else
        CV_Error(Error::StsBadArg, "Unknown estimation method");

ここで、初めて実装らしい部分につきあたる。
methodが0なので、最初の部分が実行される。
cb->runKernelが呼ばれ、こいつが何者かと言うと、if文の直前で実体が生成された、HomographyEstimatorCallbackのrunKernelが該当する
このrunKernel、本来ならRANSACなどで使われるコールバック関数なのだが、今回は1回きり、シングルスレッドでだけ呼ばれる。
runKernel内を追っていると、最後以下の部分に到着する

fundam.cpp

        eigen( _LtL, matW, matV );
        _Htemp = _invHnorm*_H0;
        _H0 = _Htemp*_Hnorm2;
        _H0.convertTo(_model, _H0.type(), 1./_H0.at<double>(2,2) );

ここでは、擬似逆行列を求めており、これによりHomographyのパラメータを決定している。

擬似逆行列とは

今回の問題の場合みたいな6行8列といった、正方行列じゃない場合に逆行列を求める方法
今回はランク落ちしている(8個のパラメータを推定したいのに、6本しか方程式が無い)ので、擬似逆行列を使って求めた逆行列は正しい逆行列ではない。
もとのIssueに書かれていた対応点4組をもとに計算した「正しい」逆行列

\frac{1}{1000}
\left(\begin{array}{cccccccc}
-100&-100&100&100&0&100&0&-100\\
-100&0&0&0&100&0&0&0\\
1000&0&0&0&0&0&0&0\\
0&-100&0&100&0&0&0&0\\
-100&-100&100&0&100&100&-100&0\\
0&1000&0&0&0&0&0&0\\
0&-1&0&1&0&1&0&-1\\
-1&0&1&0&1&0&-1&0\\
\end{array}\right)

対応点3組から擬似逆行列を使って求めた「誤った」逆行列

\frac{1}{100010}
\left(\begin{array}{cccccccc}
-1&0&1&0&0&0\\
-10001&0&0&0&10001&0\\
100010&0&0&0&0&0\\
0&-10001&0&10001&0&0\\
0&-1&0&0&0&1\\
0&100010&0&0&0&0\\
100&0&-100&0&0&0\\
0&100&0&0&0&-100\\
\end{array}\right)

何故正しく求まらないのに擬似逆行列を使っているのか？それは対応点が4組より多く存在した場合のためである。
対応点が4組より多く存在した場合、行列が正方行列ではなくなり、単純に逆行列は求められない。
そのため、擬似逆行列を使うことになる。
擬似逆行列を使ってるため、ランク落ちしてる行列でもなんかそれっぽい値が出てきてしまうのが、今回の実装のオチ

というわけで

Homographyを3組以下の対応点で求めようとしちゃ駄目だZO！
お兄さんとの約束だ！¹
明日はhon_no_mushiさんの記事で、本記事執筆時のタイトルは「はじめてのG-API(文字描画編)」だよっ！お楽しみに！

既にお兄さんという年齢では無いというツッコミは無しの方向で ↩

対応点3点ではHomographyは求まらない

はじめに

TL;DR

はじまりの物語

数学的側面

実装チェック

擬似逆行列とは

というわけで