OpenCV Advent Calendar 2015 26日目

はじめに†

OpenCV Advent Calendar 2015は昨日無事終了しました&note{fukushima1981-advent-calendar-25:OpenCV Advent Calendar 2015まとめ - Qiita, 25日目, 2015-12-25公開, 2015-12-25閲覧};が、なんとなくあたらしくマイナーなネタを見つけたので新しく書いておきます。
OpenCVに入ってるRNG(Random Number Generator)クラスで乱数を生成する話です

↑

乱数をOpenCV内で使う場面†

OpenCV内では、関数が想定通りの挙動を示すかどうかテストをします。
その際、大量にデータを生成するのに使われるのがRNGクラスです。

↑

RNG の実装を見てみよう†

RNGクラスはcore.hpp と rand.cpp内に記述があります

Random number generator. It encapsulates the state (currently, a 64-bit
integer) and has methods to return scalar random values and to fill
arrays with random values. Currently it supports uniform and Gaussian
(normal) distributions. The generator uses Multiply-With-Carry
algorithm, introduced by G. Marsaglia (
<http://en.wikipedia.org/wiki/Multiply-with-carry> ).
Gaussian-distribution random numbers are generated using the Ziggurat
algorithm ( <http://en.wikipedia.org/wiki/Ziggurat_algorithm> ),
introduced by G. Marsaglia and W. W. Tsang.

というふうにコメント内に説明があります

日本語でおｋということなので、ざっくり和訳してみます。

乱数生成器。内部に状態を保持し（状態は64bit整数）、
乱数スカラ値を返し、もしくは行列を乱数で埋める。
現在サポートされてるのは一様分布とガウシアン（正規）分布。
生成器は内部でMultiply-With-Carryアルゴリズムを使う。
G. Marsaglia ( <http://en.wikipedia.org/wiki/Multiply-with-carry> ) 導出。
ガウシアン分布の乱数はZigguratのアルゴリズム
( <http://en.wikipedia.org/wiki/Ziggurat_algorithm> ) を実装している。

ざっくりまとめると、一様分布と正規分布の乱数を生成できるとのことです。
ちなみに、一様分布の生成アルゴリズムMultiply with carryアルゴリズム&note{wikipedia-multiply-with-carry:Multiply-with-carry - Wikipedia, the free encyclopedia, 2015-12-17 14:30版, 2015-12-23閲覧};とZigguratアルゴリズム&note{wikipedia-ziggura-algorithm:Ziggurat algorithm - Wikipedia, the free encyclopedia, 2015-11-04 02:57版, 2015-12-23閲覧};に関してはWikipediaをご参照下さい
それにしてもジグラットって。バーチャロンを彷彿させる。
早速、使ってみましょう。

↑

一様分布†

一様分布は、seedでインスタンスを生成した後、a.next()で新しい乱数を生成生成します。
#geshi(C++){{

#include <iostream>

#include <iomanip>

#include <opencv2/core/core.hpp>

const uint64 initState = 0x1234354435;

int main() {

	using namespace cv;
	using namespace std;
	RNG a(initState);
	for(int i = 0;i < 10;i++)
	{
		unsigned int randomNumber = a.next();
		cout << "0x" << std::hex << std::setw(8) << std::setfill('0') << randomNumber << endl;
	}
	return 0;

} }}

出力結果

0xc4802924
0x6f670ec9
0x4ec21e44
0xe40ce671
0x5fcd40f9
0x3090cb5e
0x5d55db31
0x7145ce59
0xc5efc30e
0x7215b472

と、こんな感じに32bitの乱数が生成されます。
charやshortが必要だったらマスキングしましょう。
この一様分布のスゴいところは、実装が掛け算とマスキングだけで行われていることです。
なので、OpenCVを使いながらも、libファイルとリンクしないでも、ヘッダだけで実装が完結しています。(スゴい！)(小並感)
#geshi(C++){{ inline RNG::RNG(uint64 _state) { state = _state ? _state : 0xffffffff; } inline unsigned RNG::next() {
```
   state = (uint64)(unsigned)state* /*CV_RNG_COEFF*/ 4164903690U + (unsigned)(state >> 32);
   return (unsigned)state;
```
} }}

なんと以上。
内部の状態変数と掛け算と足し算を行うことで新しい数を生成します。
実装が簡単なだけあって、速度もそこそこ出ます。

試しに1000000回生成してみました。

#geshi(c++){{

	int64 start = cv::getTickCount();
	for(unsigned int i = 0;i < iteration;i++)
	{
		a.next();
	}
	int64 finish = getTickCount();
	std::cout << "0x" << std::hex << a.next() << std::endl;
	std::cout << (finish - start) * second2ms / getTickFrequency() << "ms for generating " << std::dec << cIteration << " random #" << std::endl;
	std::cout << ((finish - start) * second2ms / getTickFrequency()  / cIteration) << "ms / 1 random #" << std::endl;

}}

実行結果

0xd497e65c
1.5647ms for generating 1000000 random #
1.5647e-006ms / 1 random #

用途によりますが、これぐらい速いんだったら様々な場面で活躍するのではないかと思います。
参考までに、ヒストグラムを作ってみました。

↑

正規分布†

正規分布は、gaussianメソッドを呼ぶことでdouble型の値が返って来ます
このgaussianメソッドにはsigmaを渡すことができ、分散がsigmaになるような正規分布に則った値が生成されます。
実際に使ってみましょう。
#geshi(c++){{ const unsigned int cIteration = 8000000; const double second2ms = 1000.0f; const double kSigma = 3.0f;

void checkEachBitVarianceGaussian(uint64 init, unsigned int iteration) {

	using namespace cv;
	RNG a(init);
	double foo = a.gaussian(kSigma);
	int64 start = cv::getTickCount();
	for(unsigned int i = 0;i < iteration;i++)
	{
		a.gaussian(kSigma);
	}
	int64 finish = getTickCount();
	foo = a.gaussian(kSigma);
	std::cout << (finish - start) * second2ms / getTickFrequency() << "ms for generating " << cIteration << " gaussian distributed random #" << std::endl;
	std::cout << ((finish - start) * second2ms / getTickFrequency()  / cIteration) << "ms / 1 random #" << std::endl;
	std::cout << foo << std::endl;

} }}

実行結果

24.3736ms for generating 1000000 gaussian distributed random #
2.43736e-005ms / 1 random #
-0.965316

一様分布の場合は1回あたりの生成が平均1.5647e-006ms でしたが、ガウシアン分布の場合は2.43736e-005msです。10-20倍ぐらいガウシアンの方が遅いです。（一様分布が速いというべきか？）

ついでに正規分布の様子をヒストグラムで見てみましょう。

#geshi(c++){{

	using namespace cv;
	RNG a(init);
	unsigned int count[60]; //ヒストグラムのカウント
	memset(count, 0, sizeof(unsigned int)*60);
	// ヒストグラム作成
	for(unsigned int i = 0;i < iteration;i++)
	{
		double foo = a.gaussian(kSigma);
		int index = (int)(foo * 3 + 30);
		index = index <  0  ? 0  : index;
		index = index >= 60 ? 59 : index;
		count[index]++;
	}
	// ヒストグラム出力
	for(int i = 0;i < 60;i++)
	{
		std::cout << std::setw(2) << i << ' ' << ((double)(i-30)/3.0f) << '\t' << std::setw(7) << count[i] << std::endl;
	}

}}

でグラフにしてみました。
なかなか綺麗な分布ですね。

↑

OpenCV内での本来の使い方†

前述の通り、本来は大量のテストデータを生成するために使われています。
#geshi(c++){{ Mat randomMat(RNG& rng, Size size, int type, double minVal, double maxVal, bool useRoi) {
```
   Size size0 = size;
   if( useRoi )
   {
       size0.width += std::max(rng.uniform(0, 10) - 5, 0);
       size0.height += std::max(rng.uniform(0, 10) - 5, 0);
   }
```

   Mat m(size0, type);

   rng.fill(m, RNG::UNIFORM, minVal, maxVal); // ここで配列mを乱数で埋める
   if( size0 == size )
       return m;
   return m(Rect((size0.width-size.width)/2, (size0.height-size.height)/2, size.width, size.height));

} }}