4点同士の対応点情報からホモグラフィを線型に求める

#contents

*Mat findHomography(InputArray src, InputArray dst, int method, double ransacReprojThreshold, OutputArray mask, const int maxIters, const double confidence) [#ff050f2a]

-methodで指定した方法で、srcとdstの点同士を変換するホモグラフィを求める

**引数 [#t640f091]
-src:Mat型の変換前の座標
-dst:Mat型の変換後の座標
-method:int型の内部計算方法。省略するとデフォルト0
-ransacReprojThreshold:double型のしきい値。省略するとデフォルト3
-mask:Mat型の出力。対応点として使われたかどうか。
-maxIters:int型の入力。最大試行回数。省略するとデフォルト2000
-confidence:double型の入力。省略するとデフォルト0.995

**返り値 [#q451fed8]
-Mat型のホモグラフィ

*Mat findHomography(InputArray src, InputArray dst, OutputArray mask, int method, double ransacReprojThreshold); [#b7d69a14]
-findHomographyのオーバーロード
-内部的に前者のfindHomographyをcallする

**引数 [#rbf1bd92]
-オーバーロードなので、前述のfindHomographyに準ずる

**返り値 [#q66d07e7]
-Mat型のホモグラフィ

*解説 [#h57bf74d]
-srcはn*2もしくはn*3, dstはn*2もしくはn*3, homographyは3*3のサイズの行列
-n*2はn行2列(n rows, 2 columns), n*3はn行3列(n rows, 3 columns)，幅が2か3で高さがn
-nは対応点の数
-&mimetex(\left(\begin{array}{cc} x_{1} & y_{1} \\ x_{2} & y_{2} \\ \vdots & \vdots \\ x_{n} & y_{n} \end{array}\right));か&mimetex(\left(\begin{array}{ccc} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ \vdots & \vdots & \vdots \\ x_{n} & y_{n} & 1 \end{array}\right));の形で渡すとOK
-srcを&mimetex(\bf{x});, dstを&mimetex(\bf{x}^\prime);, homographyを&mimetex(\rm{H});とすると
-&mimetex(\bf{x}^\prime \simeq \rm{H}\bf{x});
-となるhomographyを返す
-methodは以下の4つが引数として使用できる
--0:全点を使って計算する。ノイズが乗った場合は自乗誤差が最小になるような行列を求める。指定を省略した場合はこの手法が使われる
--LMEDS(4):LMEDS法を使って誤った対応を除いてホモグラフィを求める
--RANSAC(8):RANSAC法を使って誤った対応を除いてホモグラフィを求める
--RHO(16):RHO法を使って誤った対応を除いてホモグラフィを求める(要OpenCV3.0以降)
--詳細は後述
-maskに点数と同じ数だけ要素を持った行列を渡すと，推定後に除去された点の情報が返ってくる
--1か0が入って返ってくるが，1が入ってる点だけが計算に利用された
--出力としてだけ使用できる。渡した値はOpenCVによって読まれることは無い。
-maxItersを指定すると、収束するまでの最大試行回数を設定できる。
--コメント文とドキュメントには2000回が最大試行回数と書いてあるが、ソースコードを読む限り別段2000回以上の試行もできる(意味があるかはわからぬ)
-confidenceには0-1の範囲のdoubleを指定する
--各対応点のconfidenceでのしきい値
--省略すると0.995
--この値より信頼性の低い対応点は棄却される


**methodについて [#p357c302]
-対応点にはだいたい誤対応が含まれる
-通常の行列演算に誤対応が含まれていると、行列全体にノイズの影響が出る
-誤対応（ドキュメントでoutlier、外れ値）を除去するために、LMeDS、RANSAC、RHOなどの手法を用いる
--LMeDS:再投影誤差の中間値が最も小さくなるような対応を初期値として選択する
--RANSAC:対応関係の点数が最も多くなるような対応を初期値として選択する
--RHO:RHO法を使って推定した対応を初期値として選択する。要OpenCV3.0以降
-RANSACとRHOは誤対応だらけでも、対応関係を推定できる。そのためにはどれくらい近ければ対応していると認めるかしきい値を定める必要がある
--パラメータのransacReprojThresholdがしきい値にあたる
-LMeDSはしきい値をあたえないでも対応関係を推定する
--ただし、誤対応が50%を超えている場合は推定がうまくいかない
-RHO手法以外を指定した場合、それぞれの推定で初期推定を求めたうえで、最適化でさらなるフィッティングが行われる
--OpenCV3.0以降ではLM法を使ったフィッティングが行われる
--0を指定した場合も4点より多い(5点以上)の対応点があればフィッティングが行われる
--何故か、RHOを指定した場合だけLM法を使わない。なにゆえだ。。
-OpenCV 2.4系列までは、対応点は、2*N、N*2、3*N、N*3のシングルチャンネル、もしくは1*Nの2チャンネルor3チャンネルの行列でも解釈してくれる
-OpenCV 3.0系列ではcolumn、つまり列数が2、もしくは3である必要がある
--つまり、ラスタスキャンした時に、x0、y0、x1、y1の順もしくはx0、y0、1、x1、y1、1の順にアクセスできる必要がある。
--バグな気がするんだけれどなぁ。

**導出 [#n9a33393]
-以下の導出は一切知らなくてもこの関数は使える
-変換前を&mimetex(\bf{x});, 変換後を&mimetex(\bf{x}^\prime);，ホモグラフィを&mimetex(\rm{H});として，次式で表せる
-&mimetex(\rm{H}\bf{x}=s\bf{x}^\prime);
-ここで，ホモグラフィ行列&mimetex(\rm{H});も斉次座標&mimetex(\bf{x});，&mimetex(\bf{x}^\prime);も，スケーリングの不定性がある
-言い換えれば，右辺だけ2倍しても等式が成り立ってしまう((厳密には等式(=)でなく，相似&mimetex(\(\simeq\));が成り立つ．解説の節では関係性を=ではなく&mimetex(\simeq);で表している．等式に書き換えるためにはスケールファクタが必要))
-その不定性を表すために，右辺にスケールファクターとして&mimetex(s);を導入した
-各ベクトル，行列を要素ごとに書くと，次式のようになる

#mimetex(\left(\begin{array}{ccc}h_{11}&h_{12}&h_{13}\\h_{21}&h_{22}&h_{23}\\h_{31}&h_{32}&1\\\end{array}\right)\left(\begin{array}{c}x\\y\\1\\\end{array}\right)=\left(\begin{array}{c}sx^\prime\\sy^\prime\\s\\\end{array}\right))

-ホモグラフィのスケーリングを固定するために，&mimetex(h_{33}=1);とする
-この行列式を分解すると，下記の3本の方程式が導出される

#mimetex(\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})

-ただ，3本の方程式それぞれはsに関して従属性があるので，独立な方程式は2本得られる
-3本の式から&mimetex(s);を消去して次の2式を得る

#mimetex(\left\{\begin{array}{c}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})

-&mimetex(\rm{H});のパラメータを含む項だけを左辺に移して整理して，次の2式を得る

#mimetex(\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})

-これを行列の形で表すと，次式になる

#mimetex(\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))

-この連立方程式を解けばホモグラフィの8つのパラメータ&mimetex(h_{11}-h_{32});が求まる
-1つの対応点から得られる方程式は2本
-よって8÷2で4組の対応点が必要となる
-今&mimetex(\bf{x}_i=\left(\begin{array}{c}x_i & y_i\\\end{array}\right));という点を4組導入すると，次式のように書き換えられる

#mimetex(\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))

-両辺それぞれ行列&mimetex(\rm{A});，&mimetex(\rm{X});，&mimetex(\rm{B});の行列で書き換えると次式になる

#mimetex(\rm{A}\rm{X}=\rm{B})

-ここで&mimetex(\rm{A});，&mimetex(\rm{B});は既知，&mimetex(\rm{X});が未知である
-&mimetex(\rm{A});が正方行列の場合
--[[&mimetex(\rm{A}^{-1}\rm{B}=\rm{X});>逆行列の計算]]となるので，パラメータ(=ホモグラフィ)が求められる
-&mimetex(\rm{A});が非正方行列の場合
--4点以上からホモグラフィを推定する場合は，&mimetex(\rm{A});は非正方行列になる
--その場合は[[擬似逆行列>逆行列の計算]]を利用して求めれば良い
--必要数以上ある対応点の組から，誤対応をLMeDs推定やRANSAC推定を利用して除去する

*サンプルコード [#q568ce4e]
-findHomography と [[warpPerspective>ホモグラフィによる画像の変換]] を使ったサンプルコード
-対応点からホモグラフィを推定し、本の表紙を、もとある本の上に重畳する
-ロバスト推定のための各手法の効果を示すために、ノイズを加えた場合も表示する
#geshi(c++,number){{
#include <opencv2/calib3d/calib3d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgcodecs/imgcodecs.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>

const char inputForegroundImageFilename[] = "file/path/to/opencv/samples/data/box.png";
const char inputBackgroundImageFilename[] = "file/path/to/opencv/samples/data/box_in_scene.png";
const unsigned char kHigh = 255;
const double kSigma = 1.0;

// in order of {x0,y0,1,   x1,y1,1,   x2,y2,1,   x3,y3,1,}
// corresponding points in box.png
// feature points extracted using SURF
const double srcPoints[] = {
	224, 53,1,145,161,1,240,185,1,204,135,1,159,174,1,200,63, 1,236,154,1,229,163,1,170,170,1,
	218,155,1,237,152,1, 83, 54,1,188,166,1, 95,145,1,201,193,1,132,140,1,180,159,1,150, 46,1,
	180,183,1,180,132,1,155,152,1,145, 93,1, 63,132,1, 89,101,1,148, 95,1,228,130,1,210,147,1,
	133,175,1,154, 44,1, 70, 89,1,167,146,1,158,167,1,166,150,1,159, 85,1, 53, 73,1, 37, 91,1,
	 90,117,1,107,163,1,
};

// corresponding points in box_in_scene.png
// feature points extracted using SURF
const double dstPoints[] = {
	243, 53,1,173,249,1,222,269,1,207,239,1,179,257,1,211,199,1,223,250,1,218,256,1,186,256,1,
	213,251,1,223,250,1,152,191,1,195,254,1,147,237,1,338,324,1,168,234,1,192,252,1,188,190,1,
	190,264,1,194,235,1,187,298,1,181,213,1,133,229,1,149,214,1,349,320,1,218,256,1,209,247,1,
	163,255,1,188,190,1,176,239,1,185,245,1,278,156,1,185,245,1,370,288,1,136,198,1,149,214,1,
	149,230,1,123,264,1,
};
const unsigned int cPoints = 38;
const cv::Scalar lineColor = cv::Scalar(kHigh,kHigh,kHigh,0);

// draw lines between corresponding points
// concatenate them horizontally to show it in one image
void drawCorrespondingLine(const cv::Mat& foreground, const cv::Mat& background, cv::Mat& result, const cv::Mat& src, const cv::Mat& dst, const cv::Mat Mask)
{
	using namespace cv;
	int heightDiff = background.rows - foreground.rows;
	if(heightDiff < 0)
	{
		// to connect in one image, foreground image must be equall or shorter than background in height
		result = background.clone();
		return;
	}

	Mat object;
	Mat space;
	if(heightDiff > 0)
	{
		// add a blank space below the foreground
		space = Mat(heightDiff, foreground.cols, CV_8UC1);
		vconcat(foreground, space, object);
	}
	else
	{
		// foreground and background are same height
		// just copy the foreground
		object = foreground.clone();
	}
	// concatenate them horizontally
	hconcat(object, background, result);

	// draw lines
	for(int iColumn = 0;iColumn < dst.rows;iColumn++)
	{
		if(Mask.data[iColumn] != 0)
		{
			Point objPoint  = Point((int)src.at<double>(iColumn, 0),                 (int)src.at<double>(iColumn, 1));
			Point backPoint = Point((int)dst.at<double>(iColumn, 0)+foreground.cols, (int)dst.at<double>(iColumn, 1));
			line(result, objPoint, backPoint, lineColor);
		}
	}
}

// estimate the homography and show the estimated result
void estimateHomography(cv::Mat& src, cv::Mat& dst, int method, const char* windowName)
{
	using namespace cv;

	Mat pointMask  = src.clone();
	Mat homography = Mat(3, 3, CV_64FC1);
	// estimate the homography
	homography = findHomography(src, dst, method, 3.0f, pointMask);

	Mat foreground = imread(inputForegroundImageFilename, IMREAD_GRAYSCALE);
	Mat background = imread(inputBackgroundImageFilename, IMREAD_GRAYSCALE);
	Mat result     = background.clone();
	Mat whiteImage = Mat().zeros(foreground.size(), CV_8UC1);
	Mat maskImage  = Mat().zeros(background.size(), CV_8UC1);
	Mat correspondence;

	// draw lines between corresponding lines
	drawCorrespondingLine(foreground, background, correspondence, src, dst, pointMask);

	// warp the image based on homography
	warpPerspective(foreground, result, homography, result.size(),INTER_LINEAR, BORDER_TRANSPARENT);

	// resize the mask to show which points are used
	resize(pointMask, pointMask, Size(20, background.rows), 0, 0, INTER_NEAREST);
	pointMask.convertTo(pointMask, CV_8U);
	pointMask *= kHigh;

	// concatenate them horizontally
	std::vector<Mat> inputLists;
	inputLists.push_back(correspondence);
	inputLists.push_back(result.clone());
	inputLists.push_back(pointMask);
	hconcat(inputLists, result);

	// show the image
	namedWindow(windowName);
	imshow(windowName, result);
	waitKey(0);

	destroyAllWindows();

};

int main(int argc, char **argv){

	using namespace cv;
	Mat src        = Mat(cPoints, 3, CV_64FC1, (void*)srcPoints);
	Mat dst        = Mat(cPoints, 3, CV_64FC1, (void*)dstPoints);

	// estimate homography using default method
	estimateHomography(src, dst, 0,      "default");
	// estimate homography using RANSAC estimation method
	estimateHomography(src, dst, RANSAC, "RANSAC");
	// estimate homography using LMeDS estimation method
	estimateHomography(src, dst, LMEDS,  "LMeDS");
	// estimate homography using RHO estimation method
	estimateHomography(src, dst, RHO,    "RHO");

	// estimate homography with some noise on the coordinates
	// in order to show the tolerance against the noise
	Mat srcNoiseAdded = src.clone();
	Mat dstNoiseAdded = dst.clone();
	RNG noise;

	// add gaussian noise
	for(int i = 0;i < src.rows;i++)
	{
		srcNoiseAdded.at<double>(i,0) += noise.gaussian(kSigma);
		srcNoiseAdded.at<double>(i,1) += noise.gaussian(kSigma);
		dstNoiseAdded.at<double>(i,0) += noise.gaussian(kSigma);
		dstNoiseAdded.at<double>(i,1) += noise.gaussian(kSigma);
	}

	// estimate the homography using each method
	estimateHomography(srcNoiseAdded, dstNoiseAdded, RANSAC, "RANSAC");
	estimateHomography(srcNoiseAdded, dstNoiseAdded, LMEDS,  "LMeDS");
	estimateHomography(srcNoiseAdded, dstNoiseAdded, RHO,    "RHO");


	// estimate homography with wrong correspondence
	// add ranndom correspondence coordinates and see if estimation can reject them
	srcNoiseAdded = src.clone();
	dstNoiseAdded = dst.clone();

	// fill the noisePart with random coordinates
	// concatenate them using push_back()
	Mat noisePart = Mat(20, 3, CV_64FC1);
	noise.fill(noisePart, RNG::UNIFORM, Mat::zeros(1, 1, CV_64FC1), Mat::ones(1, 1, CV_64FC1) * 200);
	srcNoiseAdded.push_back(noisePart);
	noise.fill(noisePart, RNG::UNIFORM, Mat::zeros(1, 1, CV_64FC1), Mat::ones(1, 1, CV_64FC1) * 200);
	dstNoiseAdded.push_back(noisePart);

	// estimate the homography using each method
	estimateHomography(srcNoiseAdded, dstNoiseAdded, RANSAC, "RANSAC");
	estimateHomography(srcNoiseAdded, dstNoiseAdded, LMEDS,  "LMeDS");
	estimateHomography(srcNoiseAdded, dstNoiseAdded, RHO,    "RHO");

	return 0;
}
}}
**実行結果 [#qd9aead3]
-誤対応を含んだ推定結果
-ハードコードされている対応点の座標には、誤対応点が38点中6点含まれている
-一番右側に表示されている白い帯は、ホモグラフィ推定に使われた対応点を表す(白=1=使用された、黒=0=誤対応)
--通常の推定結果
---たった6点の誤対応のために推定が上手く行かない
---誤対応を推定しないので、Maskはすべての要素が1になる
#ref(result000-noEstimation.png)
--RANSAC手法を使用した推定結果
---誤対応を取り除いて、正しくホモグラフィが推定されている
#ref(result001-RANSAC.png)
--LMeDs手法を使用した推定結果
---誤対応を取り除いて、正しくホモグラフィが推定されている
---RANSAC、LMeDs、RHOの3つの中で、LMeDsが一番数多くの正しい対応点を使用している
#ref(result002-LMeD.png)
--RHO手法を使用した推定結果
---誤対応を取り除いて、正しくホモグラフィが推定されている
#ref(result003-RHO.png)
-ノイズを加えた場合
-各座標に&mimetex(\sigma=1.0);のノイズを加えてみた
--RANSAC手法を使用した推定結果
#ref(result004-RansacWithNoise.png)
--LMeDs手法を使用した推定結果
#ref(result005-LmedWithNoise.png)
--RHO手法を使用した推定結果
#ref(result006-RhoWithNoise.png)
--一応RANSACが良い推定結果を出しているが、たまたまであり、これはロバスト推定ではなく、LM法による最適化のおかげ
--ノイズの乗り方では、下図の様に、RANSACを使用してもホモグラフィの推定に失敗することもある
#ref(result004.png)
-誤対応を大量に含む場合
-乱数で発生させた適当な座標を対応点に追加してみた
--RANSAC手法を使用した推定結果
---誤対応を取り除いて、正しくホモグラフィが推定されている
---右側の白い帯が特徴点が使われたか否かを表しているが、帯の下半分が黒くなっている
---つまり誤対応として除去されたことを意味する
#ref(result007-RansacWithWrongCorrespondence.png)
--LMeDs手法を使用した推定結果
---もともとあった6点の誤対応に加え、20点の誤対応を追加したため、ホモグラフィの推定に失敗している
---右側の帯を見ると、誤対応関係を多数含んで無理やりホモグラフィを推定している
#ref(result008-LmedWithWrongCorrespondence.png)
--RHO手法を使用した推定結果
---誤対応を取り除いて、正しくホモグラフィが推定されている
---RANSAC同様、右側の帯下半分が黒くなっている
#ref(result009-RhoWithWrongCorrespondence.png)
-RHOとRANSACの違いは検証できず
-対応点関係が十分にあるのならば、RANSACやRHOで弾いてしまうのが既知
-ただし、誤対応の数が十分に少ないとわかっているならば、LMeDs手法の方も使える。

*実体ファイル [#me219ce5]
**OpenCV 3.0系列 [#b1e4507f]
-findHomography
--modules/calib3d/include/opencv2/calib3d/calib3d.hpp
--modules/calib3d/src/fundam.cpp
**OpenCV 2.4系列 [#u29d903a]
--modules/calib3d/include/opencv2/calib3d/calib3d.hpp
--modules/calib3d/src/fundam.cpp

*その他 [#c794c260]
-2.4系列だと、findHomographyはCインタフェースのラッパー

ジャンル[[:OpenCV]][[:OpenCV 2.4]][[:OpenCV 3.0]][[:OpenCV 3.1]]準拠
トップ   編集 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS