opencv

最適化手法について調べていたところ，wikipediaに粒子群最適化という方法があったので，実装してみた．

wikipediaの記述どおりにざっくりと実装．パラメタの次元数をテンプレート引数にしてある．

//[ParticleSwarmOpt.h]
//	粒子群最適化

#pragma once

#include <vector>
#include <array>

//"粒子群最適化"による目的関数の最小化
//	N_PARAM : 最適化対象パラメタの個数(次元数)
//
//[使い方]
//	コーディング：
//		これを継承して，仮想関数を実装する．
//	最適化計算：
//		1) Initialize()で初期化
//		2) Update()を任意回コールして最適化
//		3) GetCurrBestParam()で結果パラメタ値を取得
template< size_t N_PARAM >
class PSO
{
public:
	//パラメタベクトルXおよび速度Vの型
	typedef std::array<double,N_PARAM> Vec_t;
protected:
	//粒子情報(CreateInitialState()引数用)
	struct ParticleSetting
	{
		Vec_t X0;	//初期位置
		Vec_t V0;	//初期速度
		std::vector< unsigned int > GroupIndexesBelong;	//粒子が所属グループ群のindex
	};
private:
	//粒子
	struct Particle
	{
		Vec_t X;	//現在位置
		Vec_t V;	//速度
		Vec_t BestX;	//この粒子の，過去最良目的関数値になった位置
		double BestOFV;	//この粒子，過去最良目的関数値
		std::vector< struct Group* > Groups;	//この粒子が所属するグループ

		Particle(){	BestOFV = std::numeric_limits<double>::max();	}
	};

	//グループ
	struct Group
	{
		Vec_t BestX;	//このグループの，過去最良目的関数値になった位置
		double BestOFV;	//このグループの，過去最良目的関数値
		
		Group(){	BestOFV = std::numeric_limits<double>::max();	}
	};

	//-----------------------------------------------------
public:
	PSO()
	{
		SetParam( 1.0, 1.0, 1.0 );
		m_GlobalBestOFV = std::numeric_limits<double>::max();
	}

public:
	//処理パラメタ値のSet,Get
	void SetParam( double W, double C1, double C2 )
	{	m_W=W;	m_C1=C1;	m_C2=C2;	}

	double W() const {	return m_W;	}
	double C1() const {	return m_C1;	}
	double C2() const {	return m_C2;	}

	//初期状態の生成.
	//最適化処理の最初の状態を生成する．
	//Update()の繰り返しを行う前に，コールする必要がある
	bool Initialize();

	//最適化処理の1stepを実行．
	//Initialize()に成功したら，気が済むまでこれを繰り返しコールすることで最適化処理を進行させる．
	//[Ret]
	//	最良パラメタ（GetCurrBestParam()が返す値）が更新された場合はtrue．
	bool Update();

	//結果取得
	Vec_t GetCurrBestParam() const {	return m_GlobalBestX;	}
	double GetObjectiveFuncVal_at_CurrBestParam() const {	return m_GlobalBestScore;	}


	//粒子位置情報取得(状態見る用)
	void GetCurrParamsOfPartilces( std::vector<Vec_t> &rDst )
	{
		rDst.resize( m_Particles.size() );
		for( size_t i=0; i<m_Particles.size(); ++i )
		{	rDst[i] = m_Particles[i].X;	}
	}
	//粒子速度情報取得(状態見る用)
	void GetCurrVelocityOfPartilces( std::vector<Vec_t> &rDst )
	{
		rDst.resize( m_Particles.size() );
		for( size_t i=0; i<m_Particles.size(); ++i )
		{	rDst[i] = m_Particles[i].V;	}
	}

	//-----------------------------------------------------
protected:
	//乱数で 0<=ret<=1 な値を返す
	virtual double Rand0to1() = 0;

	//パラメタ値Xに対する目的関数値f(X)の計算．
	//[Args]
	//	X : 位置(パラメタ値)
	//[Ret]
	//	位置Xにおける目的関数値を返す．
	virtual double ObjectiveFunc( const Vec_t &X ) = 0;

	//初期化．
	//粒子群の初期状態と，グループの個数を引数に設定して返す．
	//[Args]
	//	rDstParticles : 粒子群初期状態受取
	//	rDst_nGroup : グループ個数受取．（正常な値は1以上）
	//[Ret]
	//	成功時はtrue，失敗時はfalse
	virtual bool CreateInitialState( std::vector<ParticleSetting> &rDstParticles, unsigned int &rDst_nGroup ) = 0;

	//-----------------------------------------------------
private:
	//粒子の移動
	void MoveParticle( Particle &P )
	{
		//粒子が属するグループ群の中で，最も良い結果を探索
		Vec_t GroupX;
		{
			double OFV = std::numeric_limits<double>::max();
			for( auto pG : P.Groups )
			{
				if( pG->BestOFV < OFV )
				{
					OFV = pG->BestOFV;
					GroupX = pG->BestX;
				}
			}
		}
		//速度の更新と移動
		for( size_t i=0; i<N_PARAM; ++i )
		{
			P.V[i] *= m_W;
			P.V[i] += m_C1 * Rand0to1() * ( P.BestX[i] - P.X[i] );
			P.V[i] += m_C2 * Rand0to1() * ( GroupX[i] - P.X[i] );

			P.X[i] += P.V[i];
		}
	}

	//目的関数値の更新
	//グローバルな最良値が更新されたらtrueを返す
	bool UpdateOFV( Particle &P )
	{
		double OFV = ObjectiveFunc( P.X );
		if( OFV >= P.BestOFV )
		{	return false;	}

		P.BestOFV = OFV;
		P.BestX = P.X;
		bool ret = false;

		for( auto pG : P.Groups )
		{
			if( OFV < pG->BestOFV )
			{
				pG->BestOFV = OFV;
				pG->BestX = P.BestX;

				if( OFV < m_GlobalBestOFV )
				{
					m_GlobalBestOFV = OFV;
					m_GlobalBestX = P.BestX;
					ret = true;
				}
			}
		}
		return ret;
	}

private:
	std::vector< Particle > m_Particles;	//粒子群
	std::vector< Group > m_Groups;	//グループ群
	Vec_t m_GlobalBestX;	//粒子群全体の，過去最良目的関数値になった位置
	double m_GlobalBestOFV;	//粒子群全体の，過去最良目的関数値

	//速度Vの更新式のパラメタ
	double m_W;	//慣性定数
	double m_C1, m_C2;	//粒子が自身のベスト位置方向にいく度合，グループのベスト位置に行く度合
};

//Initialize()
template< size_t N_PARAM >
bool PSO<N_PARAM>::Initialize()
{
	//初期情報の取得
	std::vector<ParticleSetting> PSs;
	unsigned int nGroups = 0;
	if( !CreateInitialState( PSs, nGroups ) )
	{	return false;	}

	//返ってきたデータのチェック
	if( nGroups==0 || PSs.empty() )
	{	return false;	}

	for( const auto &PS : PSs )
	{
		if( PS.GroupIndexesBelong.empty() )
		{	return false;	}

		for( unsigned int g : PS.GroupIndexesBelong )
		{
			if( g >= nGroups )
			{	return false;	}
		}
	}

	//既存データ破棄
	m_Groups.clear();
	m_Particles.clear();
	m_GlobalBestOFV = std::numeric_limits<double>::max();

	//データ生成
	m_Groups.resize( nGroups );
	m_Particles.resize( PSs.size() );

	auto iP = m_Particles.begin();
	for( const auto &PS : PSs )
	{
		iP->X = iP->BestX = PS.X0;
		iP->V = PS.V0;
		iP->Groups.reserve( PS.GroupIndexesBelong.size() );
		for( unsigned int g : PS.GroupIndexesBelong )
		{
			iP->Groups.push_back( &(m_Groups[g]) );
		}

		UpdateOFV( *iP );
		iP++;
	}

	return true;
}

//Update()
template< size_t N_PARAM >
bool PSO<N_PARAM>::Update()
{
	//粒子群の移動
	for( auto &P : m_Particles )
	{	MoveParticle( P );	}

	//目的関数計算と，｛グループ，グローバル｝情報更新
	bool bGlobalUpdated = false;
	for( auto &P : m_Particles )
	{
		bGlobalUpdated |= UpdateOFV( P );
	}
	return bGlobalUpdated;
}

				
で，２次元パラメタ空間で動作テストしてみる．つまらないテストコードは以下．
				
#include "ParticleSwarmOpt.h"

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

#include <random>
#include <iostream>

class TestOpt : public PSO<2>
{
public:
	TestOpt()
		: m_RandomEngine( std::random_device()() )
	{
	}

public:
	void Show()
	{
		cv::cvtColor( m_FuncMap, m_ShowImg, CV_GRAY2BGR );

		std::vector< Vec_t > Xs, Vs;
		GetCurrParamsOfPartilces( Xs );
		GetCurrVelocityOfPartilces( Vs );
		
		size_t n = Xs.size();
		for( size_t i=0; i<n; ++i )
		{
			cv::Point P0( cvRound(Xs[i][0]), cvRound(Xs[i][1]) );
			cv::circle( m_ShowImg, P0, 2, cv::Scalar(0,255,32), -1 );

			cv::Point P1( cvRound(Xs[i][0]+Vs[i][0]), cvRound(Xs[i][1]+Vs[i][1]) );
			cv::line( m_ShowImg, P0,P1, cv::Scalar(255,32,32) );
		}

		Vec_t BestX = GetCurrBestParam();
		cv::circle( m_ShowImg, cv::Point( cvRound(BestX[0]),cvRound(BestX[1]) ), 3, cv::Scalar(0,0,255), -1 );

		cv::imshow( "PSO<2>", m_ShowImg );
	}

protected:
	//
	virtual double ObjectiveFunc( const Vec_t &X ) override
	{
		int x = cvRound( X[0] );
		int y = cvRound( X[1] );
		if( x<0 || x>=m_FuncMap.cols || y<0 || y>=m_FuncMap.rows )
		{	return 1000;	}	//※適当にでかい値

		return m_FuncMap.at<unsigned char>( y,x );
	}

	//
	virtual double Rand0to1() override
	{	return std::uniform_real_distribution<>(0.0,1.0)( m_RandomEngine );	}

	//
	virtual bool CreateInitialState( std::vector<ParticleSetting> &rDstParticles, unsigned int &rDst_nGroup ) override
	{
		const int W = 320;
		const int H = 240;
		{//目的関数の具合をてきとーにつくる
			m_FuncMap.create( H,W, CV_8UC1 );
			for( int y=0; y<H; ++y )
			{
				unsigned char *p = m_FuncMap.ptr<unsigned char>( y );
				for( int x=0; x<W; ++x, ++p )
				{
					*p = cvRound( Rand0to1()*128 );
				}
			}

			int x = cvRound( W*0.1 + W*0.8*Rand0to1() );
			int y = cvRound( H*0.1 + H*0.8*Rand0to1() );
			cv::circle( m_FuncMap, cv::Point(x,y), 7, cv::Scalar(0), -1 );
			//m_FuncMap.at<unsigned char>(y,x) = 0;

			cv::medianBlur( m_FuncMap, m_FuncMap, 15 );

			std::cout << "正解の箇所 = (" << x << ", " << y << ")" << std::endl;
		}

		//
		const unsigned int N = 100;
		rDst_nGroup = N;
		rDstParticles.resize( N );

		const double PI = acos(-1.0);
		for( unsigned int i=0; i<N; ++i )
		{
			ParticleSetting &P = rDstParticles[i];

			P.X0[0] = Rand0to1() * (W-1);
			P.X0[1] = Rand0to1() * (H-1);
			double V = Rand0to1()*5;
			double theta = 2*PI*Rand0to1();
			P.V0[0] = V * cos( theta );
			P.V0[1] = V * sin( theta );

			P.GroupIndexesBelong.reserve( 3 );
			P.GroupIndexesBelong.push_back( i>1  ?  i-1  :  N-1 );
			P.GroupIndexesBelong.push_back( i );
			P.GroupIndexesBelong.push_back( i+1<N  ?  i+1  :  0 );
		}

		//
		m_ShowImg.create( H,W, CV_8UC3 );
		return true;
	}

private:
	std::mt19937 m_RandomEngine;
	cv::Mat m_FuncMap;
	cv::Mat m_ShowImg;
};

//
int _tmain(int argc, _TCHAR* argv[])
{
	TestOpt Test;

	Test.SetParam( 0.8, 1, 1 );
	if( !Test.Initialize() )return 0;

	std::cout << "Initial" << std::endl;
	Test.Show();
	cv::waitKey();

	unsigned int iter = 0;
	do
	{
		++iter;
		std::cout << "step " << iter;
		if( Test.Update() ){	std::cout << " *";	}
		std::cout << std::endl;
		Test.Show();
	}
	while( cv::waitKey(0)!='q' );

	return 0;
}

				
適当に乱数でグレースケール画像を作り，２次元パラメタ空間上の目的関数だということにする．黒いところほど関数値が小さい．

				・緑の●が粒子．

				・粒子から出ている青い線が粒子の速度

				・赤●はこれまでに見つけた最良のパラメタの位置．
				


				初期状態の粒子をランダムに配置した様子．

				画像左下付近のブラックホールっぽい個所が目的関数が最小になる場所であるが，初期配置では粒子がそのあたりに存在しないため，ローカルミニマムな箇所に赤●が表示されている．

				問題は粒子のグループ設定をどうするかであるが，全粒子が１つのグループというのだとなんとなくつまらないので，テストコードを見てわかるように，１粒子が３グループに属する形にしてみた．これだと，ある粒子がこれまでより良い場所を見つけても，別グループの粒子にその情報が伝達するのにある程度時間がかかるハズ．

				この状態から最適化処理開始．とりあえず1000step走らせてみた．

				↓

				

				2step目で早くも粒子の一つが偶然最小値付近に到達．

				移動速度がとても大きい粒子が多数存在しているが，元ネタの式に忠実だと多分こうなるのが正解であろう．（用途次第では最大速度制限などしてもよいか？）

				↓

				

				step6．さらに解が更新された．

				情報伝達に時間がかかるため（全粒子にこの赤●の情報がいきわたるには50stepくらい必要になる），まだ赤●の箇所に粒子が集まってくる気配はない．

				↓

				

				step108．だいぶ粒子が赤●付近に集まっているが，他のいくつかのローカルミニマムな箇所にも粒子が集まっている．

				↓

				

				step550．多くの粒子が赤●付近に吸い寄せられた．

				↓

				

				step1000．まだパラメタ空間上を飛びまわって探索している粒子が残っているのが頼もしい．
				
粒子位置の目的関数の勾配等を見ることなく粒子の移動速度が決定されるので，乱数と粒子数頼みではあるが，パラメタ空間の広範囲を探索してくれることが期待できそう．

				目的関数のJacobianやHessianが無くても動かせるので，滑降シンプレックス法と同様，「とりあえずの解」が手早く欲しいときに使えるかも．
				
ただし，終了条件には工夫が要りそうに思う．

				step数や，解の更新度合といった条件だけだと，不十分な状態で止まってしまいそうな気がするので，粒子の分布の広がりのようなものを見るなどすべきかもしれない．
				
by nakiusagi3
			

今日の話題は，OpenCVを使った機械学習．物体検出をご紹介しよう
ちなみに使ったOpenCVのバージョンは1.0．上位バージョンも大筋は変わらない．

機械学習に使うツールは、
・createsamples.exe
・haartraining.exe
なのだ．
これらはOpenCVをインストールすると、デフォルトで”C:\Program Files\OpenCV\bin”の下にできる．

１．学習データを用意する

学習に必要なデータは、学習対象の画像データ（正解画像）と、学習対象が写っていない不正解画像である。
正解画像7000、不正解画像3000位は必要と言われているらしい。
で、これらの画像から学習に使う、vecファイルなるものを作成するのだが、そこでcreatesamples.exeを使用する。
大量の正解画像を用意する方法は2つあり、まず根性で全部自分で用意する方法。そして回転等加え1つの画像から自動的に大量の画像を作成する方法である。
自力で画像を用意した場合、それらの情報を記述したファイルを用意する。
ファイルの書式は、
画像のパス　画像内の学習対象の数　学習対象領域の開始座標X Y 終了座標X Y
例えば以下の画像だと、

foloder/pogi_list_image.PNG 2 45 43 113 111 195 127 264 196
みたいな感じになる。
そんな感じですべての画像について記述したファイルを用意する。
さらに、不正解画像についてのリストファイルもつくる。こちらは仕様する全ての不正解画像のパスを記述する。

２．正解データのvecファイルを作成する

createsample.exeを使用してvecファイルを作成する。
createsamples.exe -info 正解画像リストファイル名 -vec 出力vecファイル名 -num 正解サンプル数 -w サンプル横幅 -h サンプル縦幅
最低限以上の引数があればよい。
例）createsample.exe -info positivelist.txt -vec positive.vec -num 1000 -w 24 -h 24
1つの画像から量産する場合は、画像ファイルを直接指定する。
例）createsample.exe -img positive.jpg -vec positive.vec -num 1000 -w 24 -h 24
他の引数は以下の通り

Usage: createsamples.exe
[-info (collection_file_name)]
[-img (image_file_name)]
[-vec (vec_file_name)]
[-bg (background_file_name)]：サンプル量産時の背景画像リストファイル（書式は不正解画像リストと同じ）
[-num (number_of_samples = 1000)]
[-bgcolor (background_color = 0)]
[-inv] [-randinv] [-bgthresh (background_color_threshold = 80)]
[-maxidev (max_intensity_deviation = 40)]
[-maxxangle (max_x_rotation_angle = 1.100000)]：サンプル量産時の回転範囲
[-maxyangle (max_y_rotation_angle = 1.100000)]：サンプル量産時の回転範囲
[-maxzangle (max_z_rotation_angle = 0.500000)]：サンプル量産時の回転範囲
[-show [(scale = 4.000000)]]
[-w (sample_width = 24)]：検出時の最小物体サイズになる
[-h (sample_height = 24)]：検出時の最小物体サイズになる

サンプル量産時は-bgを必ず指定したほうがいいと思う。

３．トレーニングを行う

OpenCvのbin内のhaartraining.exeを使用する。
haartraining.exe -data 出力検出器ファイル名 -vec 正解画像vecファイル名 -bg 不正解画像リストファイル名 -npos 正解サンプル数 -nneg 不正解サンプル数 -w サンプル横幅 -h サンプル縦幅
例）haartraining.exe -data out.xml -vec positive.vec -bg nagativelist.txt -npos 1000 -nneg 1000 -w 24 -h 24
他の引数は以下の通り
Usage: haartraining.exe
-data (dir_name)
-vec (vec_file_name)
-bg (background_file_name)
[-npos (number_of_positive_samples = 2000)]
[-nneg (number_of_negative_samples = 2000)]
[-nstages (number_of_stages = 14)]：多い方が精度が良いが、トレーニング時間が膨大になる
[-nsplits (number_of_splits = 1)]：多い方が精度が良いが、トレーニング時間が膨大になる
[-mem (memory_in_MB = 200)]：処理に使用可能な最大メモリ
[-sym (default)] [-nonsym]：正解画像が左右対称かどうか、非対称の場合-nonsymを指定するのだが、トレーニング時間が膨大になる
[-minhitrate (min_hit_rate = 0.995000)]
[-maxfalsealarm (max_false_alarm_rate = 0.500000)]
[-weighttrimming (weight_trimming = 0.950000)]
[-eqw]
[-mode (BASIC (default) | CORE | ALL)]
[-w (sample_width = 24)]
[-h (sample_height = 24)]
[-bt (DAB | RAB | LB | GAB (default))]
[-err (misclass (default) | gini | entropy)]
[-maxtreesplits (max_number_of_splits_in_tree_cascade = 0)]
[-minpos (min_number_of_positive_samples_per_cluster = 500)]

トレーニングが終了すると指定したファイル名のxmlファイルができる。

実験

今回は会社に何故か飾ってあるハム太郎のぬいぐるみを撮影し、学習、検出を行ってみることにした．
正解画像は1000、不正解画像は1020．
真正面から何枚か撮影し、それぞれ数百くらいずつのvecファイルを作り、複数のvecファイルを合成して1000のvecファイルを作成した。
不正解画像は，カルフォルニア工科大学様などが公開している画像データベースから拝借した様々な画像1020。
トレーニング時間は弊社のスーパーコンピュータで10時間くらい。
他のデータで試したときは3週間くらいかかっても終わらず、強制終了させたこともあった。

そして、結果発表．
(目線は検出処理してから入れましたよ)
奨励画像数よりだいぶ少なかったが、結構がんばる。正面向き、回転少なめならばかなりの精度で検出する．
検出位置が微妙に右側にずれてるのが少し気になるところ．

今日ご紹介した方法で自分の好きな物体を検出できる．
ヤクルトの小瓶にひたすら反応する「ヤクルト検出器」や岩手県の形状に反応する「岩手県検出器」，あるいは馬の顔を検出する「馬顔検出」など応用は尽きないのである．