twitter アイコンから or OpenCVで Hybrid images 生成

twitter アイコンで Hybrid images が作れるページを公開しました。

http://luigitefu.xrsp.net/cgi-bin/hybridtwicon/

・Hybrid images とは

Hybrid images とは人間の視覚特性を利用し、(近視 or 縮小 or 遠くから見る)場合と、(not近視 or 拡大 or 近くから見る)場合で違って見える画像のことです。
MITで考案され、以下のページと論文に多数の具体例や原理が紹介されています。
　　　Hybrid Images @ MIT
　　　Hybrid Images - A. Oliva, A. Torralba, P.G. Schyns, ACM Transactions on Graphics, vol. 25-3, pages 527-530, 2006.

今回は2つの画像を入力とすることで、縮小したときと、拡大したときに、それぞれの画像が強く見えるような画像を生成しています。

それぞれの画像の低周波成分（ぼかした画像）と、高周波成分（エッジ部分）だけを組み合わせるとHybrid images が生成されます。
具体的な処理としては、周波数領域において、以下のαブレンドのような計算を行っています。ただし、ここでGはガウシアンカーネルです。
Output = Input1・G + Input2・(1-G)　
このガウシアンカーネルのパラメータ（標準偏差）を調節することで、各画像の強さや、どの程度の距離の人にどういった画像を見せたいか、といった調節が可能です。画像のエッジの強さ等によっても適切なパラメータが変わってくると思います。

twitterアイコンだけでなく、任意の画像で Hybrid imagesを生成したい場合は、以下のページが参考になると思います（MATLAB か Octave が必要になりますが）。
　　　hybrid imagesの作り方 - デー
ちなみに、ここで紹介されているインデックスとイカ娘で作られたHybrid imageは必見じゃなイカ？

・OpenCV で Hybrid images

ついでに、C++ / OpenCV による Hybrid images を生成するプログラムのソースコードを以下に晒しておきます。
かなり古臭い書き方をしているので、多分OpenCV1.0以上であれば動くと思います。

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// HybridImages.cpp  by @Luigitefu
//
// 2つの入力画像からHybrid imagesを生成します．
//
// Hybrid imagesの詳細については以下のWebページと論文を参照してください．
//    http://cvcl.mit.edu/hybridimage.htm
//    Hybrid Images - A. Oliva, A. Torralba, P.G. Schyns, ACM Transactions on Graphics, vol. 25-3, pages 527-530, 2006. 
//
// コンパイルには，OpenCVが必要です．
// OpenCVによるDCTに関しては，以下のWebページを参考にしています．
//    http://opencv.jp/sample/discrete_transforms.html 
//    http://opencv.jp/opencv-1.0.0/document/opencvref_cxcore_discrete.html
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>
#include <cv.h>
#include <highgui.h>

// OpenCV ライブラリ読み込み
#pragma comment(lib,"cv.lib")
#pragma comment(lib,"cxcore.lib")
#pragma comment(lib,"highgui.lib")


/////////////////////////
////  パラメータ類   
/////////////////////////
#define SIGMA 7.5            // ガウシアンカーネルの標準偏差パラメータ．これが大きいとINPUT1がより強く見える
#define SIDE_MAX 400        // INPUT1のどちらか1辺のサイズがこれより大きければリサイズ
#define INPUT1 "input1.png" // 縮小したときに強く見える画像
#define INPUT2 "input2.png" // 拡大したときに強く見える画像
#define OUTPUT "output.png" // 出力されるHybrid images


////////////////////////////////////
//// ガウシアンカーネル作成
////////////////////////////////////
void myGauss(
	CvMat* gau1,
	CvMat* gau2,
	const float sig
	){
		float gaui,gaud;
		int xs,ys,x,y;
		float xm,ym;
		int cx,cy;
		CvMat *tmp = 0;
		CvMat q1stub, q2stub;
		CvMat q3stub, q4stub;
		CvMat d1stub, d2stub;
		CvMat d3stub, d4stub;
		CvMat *q1, *q2, *q3, *q4;
		CvMat *d1, *d2, *d3, *d4;
		float ss2=sig*sig*2;

		xs=gau1->width;
		ys=gau1->height;
		xm=(float)((float)xs/2.0);
		ym=(float)((float)ys/2.0);
		for(x=0;x<xs;x++){
			for(y=0;y<ys;y++){
				gaud = (xm-x)*(xm-x)+(ym-y)*(ym-y);
				gaui = exp(-gaud/ss2);
				gau1->data.fl[(gau1->width*y+x)*2] = gaui;
				gau1->data.fl[(gau1->width*y+x)*2+1] = gaui;
				gau2->data.fl[(gau2->width*y+x)*2] = 1-gaui;
				gau2->data.fl[(gau2->width*y+x)*2+1] = 1-gaui;
			}
		}

		////////////////////////////////////////////////////////////
		// カーネルの象限を入れ替える処理
		///////////////////////////////////////////////////////////
		tmp = cvCreateMat (ys / 2, xs / 2,CV_32FC2  );
		cx=(int)xm;
		cy=(int)ym;

		q1 = cvGetSubRect (gau1, &q1stub, cvRect (0, 0, cx, cy));
		q2 = cvGetSubRect (gau1, &q2stub, cvRect (cx, 0, cx, cy));
		q3 = cvGetSubRect (gau1, &q3stub, cvRect (cx, cy, cx, cy));
		q4 = cvGetSubRect (gau1, &q4stub, cvRect (0, cy, cx, cy));
		d1 = cvGetSubRect (gau1, &d1stub, cvRect (0, 0, cx, cy));
		d2 = cvGetSubRect (gau1, &d2stub, cvRect (cx, 0, cx, cy));
		d3 = cvGetSubRect (gau1, &d3stub, cvRect (cx, cy, cx, cy));
		d4 = cvGetSubRect (gau1, &d4stub, cvRect (0, cy, cx, cy));
		cvCopy (q3, tmp, 0);
		cvCopy (q1, q3, 0);
		cvCopy (tmp, q1, 0);
		cvCopy (q4, tmp, 0);
		cvCopy (q2, q4, 0);
		cvCopy (tmp, q2, 0);
		q1 = cvGetSubRect (gau2, &q1stub, cvRect (0, 0, cx, cy));
		q2 = cvGetSubRect (gau2, &q2stub, cvRect (cx, 0, cx, cy));
		q3 = cvGetSubRect (gau2, &q3stub, cvRect (cx, cy, cx, cy));
		q4 = cvGetSubRect (gau2, &q4stub, cvRect (0, cy, cx, cy));
		d1 = cvGetSubRect (gau2, &d1stub, cvRect (0, 0, cx, cy));
		d2 = cvGetSubRect (gau2, &d2stub, cvRect (cx, 0, cx, cy));
		d3 = cvGetSubRect (gau2, &d3stub, cvRect (cx, cy, cx, cy));
		d4 = cvGetSubRect (gau2, &d4stub, cvRect (0, cy, cx, cy));
		cvCopy (q3, tmp, 0);
		cvCopy (q1, q3, 0);
		cvCopy (tmp, q1, 0);
		cvCopy (q4, tmp, 0);
		cvCopy (q2, q4, 0);
		cvCopy (tmp, q2, 0);
}


//////////////////////////////////////////
//// 1チャンネルでの Hybrid images 生成
//////////////////////////////////////////
void HybridImages1ch(
	CvMat* src1,
	CvMat* src2, 
	CvMat* dft_A,
	CvMat* dft_B, 
	const CvMat* gau1, 
	const CvMat* gau2, 
	const CvMat* srcIm )
{
	// 入力画像と虚数配列をマージして複素数平面を構成
	CvMat* complex1 = cvCreateMat( src1->rows, src1->cols, CV_32FC2 ); 
	CvMat* complex2 = cvCreateMat( src1->rows, src1->cols, CV_32FC2 );
	cvMerge(src1, srcIm, NULL, NULL,complex1);
	cvMerge(src2, srcIm, NULL, NULL,complex2);

	CvMat tmp;
	// 複素数平面をdft_A, dft_Bにコピーし，残りの行列右側部分を0で埋めた後，離散フーリエ変換を行う
	cvGetSubRect( dft_A, &tmp, cvRect(0,0,complex1->cols,complex1->rows));
	cvCopy( complex1, &tmp );
	if(dft_A->cols > complex1->cols){
		cvGetSubRect( dft_A, &tmp, cvRect(complex1->cols,0,dft_A->cols - complex1->cols,complex1->rows));
		cvZero( &tmp );
	}
	cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complex1->rows );

	
	cvGetSubRect( dft_B, &tmp, cvRect(0,0,complex1->cols,complex1->rows));
	cvCopy( complex2, &tmp );
	if(dft_A->cols > complex1->cols){
		cvGetSubRect( dft_B, &tmp, cvRect(complex1->cols,0,dft_B->cols - complex1->cols,complex1->rows));
		cvZero( &tmp );
	}
	cvDFT( dft_B, dft_B, CV_DXT_FORWARD, complex1->rows );

	// 周波数領域において，ガウシアンカーネルGを用いて以下の処理を行う
	// Output = Input1・G + Input2・(1-G)
	cvMul(dft_A,gau1,dft_A);
	cvMul(dft_B,gau2,dft_B);
	cvAdd(dft_A,dft_B,dft_A);

	// 離散フーリエ逆変換し，実数成分をsrc1に，虚数成分をsrc2に格納
	cvDFT( dft_A, dft_A, CV_DXT_INV_SCALE, src1->rows ); // 上部のみを計算する
	cvGetSubRect( dft_A, &tmp, cvRect(0,0,src1->cols,src1->rows) );
	cvCopy( &tmp, complex1 );
	cvSplit(complex1,src1,src2,0,0);
	cvReleaseMat(&complex1);
	cvReleaseMat(&complex2);
}


//////////////////////////////
////  Hybrid images 生成
//////////////////////////////
void HybridImages(
	CvMat* src1,
	const CvMat* src2, 
	const float sig )
{
	//行列確保
	int dft_M = cvGetOptimalDFTSize( src1->height ); 
	int dft_N = cvGetOptimalDFTSize( src1->width  ); 
	CvMat* src1b = cvCreateMat( src1->rows, src1->cols, CV_32FC1 ); 
	CvMat* src1g = cvCreateMat( src1->rows, src1->cols, CV_32FC1 ); 
	CvMat* src1r = cvCreateMat( src1->rows, src1->cols, CV_32FC1 ); 
	CvMat* src2b = cvCreateMat(src2->rows, src2->cols, CV_32FC1 ); 
	CvMat* src2g = cvCreateMat(src2->rows, src2->cols, CV_32FC1 ); 
	CvMat* src2r = cvCreateMat(src2->rows, src2->cols, CV_32FC1 ); 
	CvMat* srcIm = cvCreateMat(src2->rows, src2->cols, CV_32FC1 ); 
	CvMat* dft_A = cvCreateMat( dft_M, dft_N, CV_32FC2 ); 
	CvMat* dft_B = cvCreateMat( dft_M, dft_N, CV_32FC2 ); 
	CvMat* dft_tmp = cvCreateMat( dft_M, dft_N, CV_32FC1 ); 
	CvMat* gau1 = cvCreateMat( dft_M, dft_N, CV_32FC2 );
	CvMat* gau2 = cvCreateMat( dft_M, dft_N, CV_32FC2 );
	myGauss(gau1,gau2,sig);
	cvZero(srcIm);

	//各チャンネルごとに処理
	cvSplit(src1,src1b,src1g,src1r,0);
	cvSplit(src2,src2b,src2g,src2r,0);
	HybridImages1ch(src1b, src2b, dft_A, dft_B, gau1, gau2,srcIm);
	HybridImages1ch(src1g, src2g, dft_A, dft_B, gau1, gau2,srcIm);
	HybridImages1ch(src1r, src2r, dft_A, dft_B, gau1, gau2,srcIm);
	cvMerge(src1b, src1g, src1r, NULL,src1);

	//メモリの解放
	cvReleaseMat(&src1b);
	cvReleaseMat(&src1g);
	cvReleaseMat(&src1r);
	cvReleaseMat(&src2b);
	cvReleaseMat(&src2g);
	cvReleaseMat(&src2r);
	cvReleaseMat(&srcIm);
	cvReleaseMat(&dft_A);
	cvReleaseMat(&dft_B);
	cvReleaseMat(&gau1);
	cvReleaseMat(&gau2);
}

/////////////////////////////////////////
//// 画像を読み込んでHybrid imagesを保存
/////////////////////////////////////////
int HybridImages_main(
	const char* fname1,
	const char* fname2,
	const char* outname,
	float sig )
{
	IplImage *img1,*img2,*result,*dst1,*dst2, *res,img_hdr;
	int width,height;
	CvMat istub1, *src1,istub2,*src2;

	//INPUT1
	if((img1=cvLoadImage(fname1 ,CV_LOAD_IMAGE_COLOR ))==0){
		return -1;
	}
	if(img1->width > img1->height){
		if(img1->width > SIDE_MAX){
			width=SIDE_MAX;
			height = SIDE_MAX* img1->height / img1->width;
		}else{
			width =img1->width;
			height = img1->height;		
		}	
	}else{
		if(img1->height > SIDE_MAX){
			height=SIDE_MAX;
			width = SIDE_MAX* img1->width / img1->height;
		}else{
			width =img1->width;
			height = img1->height;		
		}
	}
	dst1 =cvCreateImage(cvSize(width,height),IPL_DEPTH_8U, 3);
	cvResize(img1,dst1,CV_INTER_CUBIC);
	cvReleaseImage(&img1);
	img1 = cvCreateImage(cvSize(width,height), IPL_DEPTH_32F, 3);
	cvConvertScale(dst1, img1);
	src1 = cvGetMat(img1, &istub1);
	cvReleaseImage(&dst1);
	
	//INPUT2
	if( (img2 = cvLoadImage(fname2, CV_LOAD_IMAGE_COLOR ) )==0){
		return -1;
	}
	dst2 =cvCreateImage(cvSize(width,height),IPL_DEPTH_8U, 3);
	cvResize(img2,dst2,CV_INTER_CUBIC); 
	cvReleaseImage(&img2);
	img2 = cvCreateImage(cvSize(width,height), IPL_DEPTH_32F, 3);
	cvConvertScale(dst2, img2);
	src2 = cvGetMat(img2, &istub2);
	cvReleaseImage(&dst2);

	//Hybrid images 生成処理を呼び出し，結果画像を保存
	HybridImages(src1, src2, sig);
	result = cvCreateImage(cvSize(width,height), IPL_DEPTH_8U, 3);
	res = cvGetImage(src1, &img_hdr);
	cvConvertScaleAbs(res, result);
	if( (cvSaveImage(outname,result))==0){
		return -1;
	}
	return 1;
}


int main(int argc, char* argv[]){
	HybridImages_main(INPUT1,INPUT2,OUTPUT,SIGMA);
	return 1;
}