2016-12

この2年ほど忙しすぎて会社のHPをほったらかしにしていたのですが、WEBの更新がない会社って大丈夫だろうか？っていう時代だし、何か書かないといけないねという強迫観念に駆られまして、約2年ぶりに記事を書いております。
久しぶりの記事の内容は何にするかというのは、これもまた悩ましいところですが、世間の関心のある先端技術はよそにまかせて、自分が先日困ったことを話題に取り上げて久しぶりの記事投稿再開の1号としたいと思います。

WindowsのアプリをC++でゴリゴリ書くのは、もう時代遅れかもしれないのですが、画像を扱うアプリを書くときには、今もMFCは大事な選択しだと思うのです。
2016年の年末になった今日も、16年前とほとんど変わらないMFCを使ってWindowsアプリを書いているというのは、考えてみるとなかなかMFCすごいんじゃないかって思うんですが同意してくれる人があんまりいませんね。
日進月歩のソフトウエア技術のなかにあっても、未だに約20年前とほぼ同じ仕組みをつかって最新のWindowsでも動作するハイパフォーマンスソフトを書けるのですから私のような中年のソフト屋さんにはMFCは手放せません。

MFCを使って画像処理アプリを書いていると画像を扱うIplImageとかcv::MatとかをCArchive.Serializeしちゃいたい時がよくある訳なんですが、IplImage->imageDataをそのままWriteしちゃったりすると、ものすごいファイルが大きくなるわけです。だから多少の画像の劣化は妥協して、JPEG圧縮したうえで, CArchiveにSerializeしたいなと思うことになるんですが、いざやってみようかとおもうと、libjpeg使わないとダメかなあとか、結構めんどくさそうな感じになってくるんですね。
手軽にだれかが作ってくれたIplImageをCArchiveしてくれるコードがないかなあと探してみたりしたのですが、今の時代にMFCを使ってる人ってあんまりいないみたいで、MFC and IplImage and CArchiveなどという情報は、Google先生もズバリな答えを教えてくれません。仕方なしに、ほかの断片的な情報をつなぎ合わせて、Gdiplusの機能を使って IplImageをCArchiveに保管するコードを書いてみましたのでご紹介します。

IplImageをCArchiveに保存する関数がwrite_iplimage_with_jpeg()です。
Jpeg圧縮後のデータを保存するためのメモリーをGlobalAllocで確保していますが、生画像サイズよりは小さくなるだろうということで、生画像サイズ丸ごと確保しています。
その後にGdiplusに圧縮してもらうため、メモリーブロックからIStremハンドルを取得しています。GetEncoderClsid()は、MSDNのサンプルから頂きました。
IplImage2GdiPlusBitmap()は、IplImageをGdiplusのBitmapに変換する自作関数です。”image/jpeg”を、”image/png”にしたらPNG形式でSerializeできます。
IplImage2GdiPlusBitmap()も記事の後半に貼り付けておきます。

	//IplImageをJPEG圧縮してCArchiveに保存する
	void write_iplimage_with_jpeg(IplImage* ipl, CArchive& ar, int quality)
	{
		int w = ipl->width, h = ipl->height;
		HGLOBAL hJpegBuffWrite = ::GlobalAlloc(GMEM_MOVEABLE, w * h * 3);
		BYTE * pJpegBuffWrite = (BYTE*)::GlobalLock(hJpegBuffWrite);
		IStream*isImageWrite;
		CreateStreamOnHGlobal(hJpegBuffWrite, TRUE, &isImageWrite);/* メモリブロックからIStreamを作成 */
		CLSID jpgClsid;
		GetEncoderClsid(L"image/jpeg", &jpgClsid);

		ULARGE_INTEGER pos;
		LARGE_INTEGER dis;
		dis.QuadPart = 0;

		Bitmap* pBmp = IplImage2GdiPlusBitmap(ipl);

		Gdiplus::EncoderParameters p;
		p.Count = 1;
		p.Parameter[0].Guid = EncoderQuality;
		p.Parameter[0].Type = EncoderParameterValueTypeLong;
		p.Parameter[0].NumberOfValues = 1;
		p.Parameter[0].Value = &quality;

		pBmp->Save(isImageWrite, &jpgClsid, &p);
		delete pBmp;

		isImageWrite->Seek(dis, STREAM_SEEK_CUR, &pos);
		TRACE("after seek=%d\n", pos.QuadPart);

		BYTE* pJpegBuff = (BYTE*)::GlobalLock(hJpegBuffWrite);
		SIZE_T jpeg_length = (SIZE_T)pos.QuadPart;

		ar << jpeg_length;
		ar.Write(pJpegBuff, jpeg_length);

		isImageWrite->Release();

		::GlobalUnlock(hJpegBuffWrite);
		::GlobalFree(hJpegBuffWrite);
	}

				
続いて、CArchiveからIplImageを読み出すコードです。
				
	//CArchiveに保存されたJPEG圧縮画像をIplImageに読み出す
	IplImage* read_iplimage_with_jpeg(CArchive& ar)
	{
		//ar >> 
		SIZE_T jpeg_length;
		ar >> jpeg_length;
		HGLOBAL hJpegBuff = ::GlobalAlloc(GMEM_MOVEABLE, jpeg_length);
		BYTE * pJpegBuff = (BYTE*)::GlobalLock(hJpegBuff);
		ar.Read(pJpegBuff, jpeg_length);//JPEGデータの読み込み
		::GlobalUnlock(hJpegBuff);
		IStream* isImage;
		CreateStreamOnHGlobal(hJpegBuff, TRUE, &isImage);/* メモリブロックからIStreamを作成 */
		Gdiplus::Bitmap bmp(isImage);

		IplImage* new_image = GdiPlusBitmap2IplImage(&bmp);
		TRACE("(%d,%d)\n", bmp.GetWidth(), bmp.GetHeight());

		::GlobalFree(hJpegBuff);
		return new_image;
	}

				
IplImageからGdiplusのBitmapに変換する関数とその逆の関数です。GdiPlusBitmap2IplImage()の戻り値のIplImage*は、使用後にReleaseするのを忘れてはいけません。
				
	Gdiplus::Bitmap*  IplImage2GdiPlusBitmap(IplImage* pIplImg)
	{
		if (!pIplImg)
			return NULL;

		Gdiplus::Bitmap *pBitmap = new Gdiplus::Bitmap(pIplImg->width, pIplImg->height, PixelFormat24bppRGB);
		if (!pBitmap)
			return NULL;

		Gdiplus::BitmapData bmpData;
		Gdiplus::Rect rect(0, 0, pIplImg->width, pIplImg->height);
		pBitmap->LockBits(&rect, Gdiplus::ImageLockModeWrite, PixelFormat24bppRGB, &bmpData);
		BYTE *pByte = (BYTE*)bmpData.Scan0;

		if (pIplImg->widthStep == bmpData.Stride) //likely  
			memcpy(bmpData.Scan0, pIplImg->imageDataOrigin, pIplImg->imageSize);

		pBitmap->UnlockBits(&bmpData);
		return pBitmap;
	}

	IplImage* GdiPlusBitmap2IplImage(Gdiplus::Bitmap* bmp)
	{
		auto format = bmp->GetPixelFormat();
		if (format != PixelFormat24bppRGB)
			return nullptr;

		Gdiplus::Rect rcLock(0, 0, bmp->GetWidth(), bmp->GetHeight());
		Gdiplus::BitmapData bmpData;

		bmp->LockBits(&rcLock, Gdiplus::ImageLockModeRead, format, &bmpData);

		int buffSz = bmpData.Stride * bmpData.Height;
		int depth = 8, channel = 3;
		IplImage* cvImage = cvCreateImage(cvSize(rcLock.Width, rcLock.Height), depth, channel);
		const unsigned char* src = static_cast<unsigned char*>(bmpData.Scan0);
		std::copy(src, src + buffSz, cvImage->imageData);

		bmp->UnlockBits(&bmpData);
		if (bmpData.Stride<0)
			cvFlip(cvImage, cvImage);

		return cvImage;
	}


				
MSDNのサンプルから引っ張ってきた関数です。
				
	static int GetEncoderClsid(const WCHAR* format, CLSID* pClsid)
	{
		UINT  num = 0;          // number of image encoders
		UINT  size = 0;         // size of the image encoder array in bytes

		ImageCodecInfo* pImageCodecInfo = NULL;

		GetImageEncodersSize(&num, &size);
		if (size == 0)
			return -1;  // Failure

		pImageCodecInfo = (ImageCodecInfo*)(malloc(size));
		if (pImageCodecInfo == NULL)
			return -1;  // Failure

		GetImageEncoders(num, size, pImageCodecInfo);

		for (UINT j = 0; j < num; ++j)
		{
			if (wcscmp(pImageCodecInfo[j].MimeType, format) == 0)
			{
				*pClsid = pImageCodecInfo[j].Clsid;
				free(pImageCodecInfo);
				return j;  // Success
			}
		}

		free(pImageCodecInfo);
		return -1;  // Failure
	}

				
T.Kamata
			

最適化手法について調べていたところ，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