Home > Tags > 技術
技術
粒子群最適化
最適化手法について調べていたところ,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; }で,2次元パラメタ空間で動作テストしてみる.つまらないテストコードは以下.
#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; }適当に乱数でグレースケール画像を作り,2次元パラメタ空間上の目的関数だということにする.黒いところほど関数値が小さい.
・緑の●が粒子.
・粒子から出ている青い線が粒子の速度
・赤●はこれまでに見つけた最良のパラメタの位置.
初期状態の粒子をランダムに配置した様子.
画像左下付近のブラックホールっぽい個所が目的関数が最小になる場所であるが,初期配置では粒子がそのあたりに存在しないため,ローカルミニマムな箇所に赤●が表示されている.
問題は粒子のグループ設定をどうするかであるが,全粒子が1つのグループというのだとなんとなくつまらないので,テストコードを見てわかるように,1粒子が3グループに属する形にしてみた.これだと,ある粒子がこれまでより良い場所を見つけても,別グループの粒子にその情報が伝達するのにある程度時間がかかるハズ.
この状態から最適化処理開始.とりあえず1000step走らせてみた.
↓
2step目で早くも粒子の一つが偶然最小値付近に到達.
移動速度がとても大きい粒子が多数存在しているが,元ネタの式に忠実だと多分こうなるのが正解であろう.(用途次第では最大速度制限などしてもよいか?)
↓
step6.さらに解が更新された.
情報伝達に時間がかかるため(全粒子にこの赤●の情報がいきわたるには50stepくらい必要になる),まだ赤●の箇所に粒子が集まってくる気配はない.
↓
step108.だいぶ粒子が赤●付近に集まっているが,他のいくつかのローカルミニマムな箇所にも粒子が集まっている.
↓
step550.多くの粒子が赤●付近に吸い寄せられた.
↓
step1000.まだパラメタ空間上を飛びまわって探索している粒子が残っているのが頼もしい.粒子位置の目的関数の勾配等を見ることなく粒子の移動速度が決定されるので,乱数と粒子数頼みではあるが,パラメタ空間の広範囲を探索してくれることが期待できそう.
目的関数のJacobianやHessianが無くても動かせるので,滑降シンプレックス法と同様,「とりあえずの解」が手早く欲しいときに使えるかも.ただし,終了条件には工夫が要りそうに思う.
step数や,解の更新度合といった条件だけだと,不十分な状態で止まってしまいそうな気がするので,粒子の分布の広がりのようなものを見るなどすべきかもしれない.by nakiusagi3
- Comments (Close): 0
- Trackbacks (Close): 0
画像から対象を検出して数える検討
上司「同じ物が複数写っている画像に対して,ユーザがその中の1個を『これ』って指定したら残りのやつを全部見つけるソフトの検討しろ」
私「条件が自由すぎて難し…
上司「いいからやれ.これ画像例な.」
という話になっており,泣きそうです.
例として渡されたが画像はコレ↓
どうやら目的は,こういった画像上で対称物体の個数をカウントしたいという感じのようです.
(撮影対象は南瓜の種の表面に謎の粉が付着したMADE IN CHINAな菓子(?)ですが,面倒なので本稿ではこれを単に「豆」と呼称することにします.)
与えられる情報は正例1個だけなので機械学習みたいな方法はとれません,というかそもそも特定の対象(例えば「人」とか「車」とか)を検出するという話ではないために,対象物に対する前提知識が全く使えません.
指定された領域をテンプレートとしたテンプレートマッチングくらいしか方法が思い浮かばないので,そういった方向で検討開始.
【検討内容】
画像上を走査し,各位置でのZernikeモーメント特徴と正解領域の同特徴とを比較してそれっぽいところを見つける.
なんでZernikeモーメント?→特徴量が回転不変とのことなので,豆がいろんな方向を向いて写っていることに対処してくれるかな?ということで採用.
豆のサイズにはいくらか変動があるので,走査窓のスケールを変えたりする必要がありそうですが,その辺のことは後回しにして,まずは単一のスケールだけで走査してみました.
特徴量比較結果がどの程度であれば「それっぽい」のか?という判断基準が必要ですが,適当に結果値を閾値で切るのではあまりうまくいかなかったので,正解画像のモーメントを元にしてそこから「明らかに不正解っぽい」モーメントを複数作って(=でっちあげて),「それらよりも正解の方に近いこと」という条件で判定してみました.
また,比較に際しては,n=m=0の要素は用いないことにしました.
太い赤丸が指定された正解領域で,細いのが検出結果です.
この例では,半径方向の次数(? 正式には何と呼ぶのか不明)nは n=0~6まで(かつ半径方向の次数(?)m≧0 を,n+m≦6 に制限)のZernike多項式を用いて特徴ベクトルを計算しています.
赤丸が豆の中央をきちんと捉えていないものが多く,「見つけた」と言っていいのかどうか判断しかねる結果.
(複数の走査スケールで処理すれば改善するだろうか?)
あと,Zernikeモーメントからの特徴量は画像輝度を反転しても変わらないようで,最初何も考えないでやったら「中央が背景で周辺部が豆」な個所が大量に出てきてしまい途方にくれましたが,そのような個所についてはn=2,m=0のモーメントの符号で判断して除外することにしました.
ちなみに,正解をもう一つ指定して,2回の処理の結果を統合してみると↓のようになります.
2回目の指定領域が緑太丸で,検出結果が緑細丸です.
ただし,2回の処理結果の中で位置が重複したものについては,より良さそうな方を残すことにし,除外されたものは青で描画してあります.
正解画像のモーメントから画像を再構成してみるとこんな感じ↓(画像は2倍拡大してある.1回目のが左,2回目のが右)
かろうじて正円ではなくて楕円な雰囲気が感じ取れますが,豆なのか何なのかわからない様相となっています.
より次数の高い項までを用いればそれだけ高周波成分の情報が増えるので豆な形に近づくのですが,やりすぎると今度は豆のテクスチャ(模様)の差が大きく影響するようになってくるので,増やせばいいというわけでもなさそうで,この辺の決め方をどうするのかについては難しい問題に思います.
他の画像(ネットで拾ってきた猫画像)での例
やはりなんとも言い難い結果です.
あと,どうでもいい話ですが,この豆画像,グレースケールで見るとMy Favorite菓子である みながわ製菓のとうがらしの種 にそっくり.
(売っている店が近場にほとんど無くて困っています! オンラインショップで買えるけど,価格に占める送料の割合が大きすぎるのがどうにも…)
by nakiusagi3
- Comments (Close): 0
- Trackbacks (Close): 0
Home > Tags > 技術