Go to Qiita Advent Calendar 2024 Top
LoginSignup
23
12

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

OpenCVAdvent Calendar 2020

Day 2
@tomoaki_teshima(Aki Teshima)

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

Last updated at Posted at 2020-12-01

はじめに

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の奥地へと飛んだ

実装チェック

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);
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のパラメータを決定している。

擬似逆行列とは

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

\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！
  • お兄さんとの約束だ！1
  • 明日はhon_no_mushiさんの記事で、本記事執筆時のタイトルは「はじめてのG-API(文字描画編)」だよっ！お楽しみに！

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

23
12
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
23
12

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?