2010-06-30

ウォーリーの途中ですが、次へ進む前に、ひとつだけツールが必要だということに気がつきました。なので、ソースも含めて、公開します。

興味のある方は、こちらをご覧ください。
http://denshika.cc/open/20100630/

（免責）
使って何か起きても、責任とれませんので、あしからず。

（使い方）

SSIM.zipをダウンロードして、解凍して、Cドライブの直下などにフォルダごと入れます。
Cドライブ直下に入れたとして、A.tifとB.tifのSSIMを計算したい場合、こんな感じです。
3番目のMAP.tifについては、次回、ウォーリーの続きで説明します。
今回のプログラムでは、TIFF、JPEGなど一般的なフォーマットが対象です。

/* 
 *
 * http://mehdi.rabah.free.fr/SSIM/
 * を参考に、作りましたよー
 * でも、プログラムはアマチュアレベルなので、
 * 細かい部分は、直してくださーい。
 * http://denshikA.cc
 * 
 * 動かすためには、OpenCVが必要
 * http://opencv.jp/
 * 
 */

/* オリジナルメッセージ
 *
 * The equivalent of Zhou Wang's SSIM matlab code using OpenCV.
 * from http://www.cns.nyu.edu/~zwang/files/research/ssim/index.html
 * The measure is described in :
 * "Image quality assessment: From error measurement to structural similarity"
 * C++ code by Rabah Mehdi. http://mehdi.rabah.free.fr/SSIM
 *
 * This implementation is under the public domain.
 * @see http://creativecommons.org/licenses/publicdomain/
 * The original work may be under copyrights. 
 *
 */

#include <iostream>
//std::maxとstd::minを使えるようにするため、こうすると良いと書いてあった
//http://d.hatena.ne.jp/pyopyopyo/20100329/p1
#define NOMINMAX 
//
#include "opencv\cv.h"	
#include "opencv\highgui.h"

/* 
 * SSIM_Cal [a.tif] [b.tif] [map.tif]
 * TIFFを２つ読み込んで、SSIM_MAPを保存して、計算結果を返す
 * 読み込んだ画像がカラーの場合、グレースケール変換
 * 変換のウェイト：Y <- 0.299*R + 0.587*G + 0.114*B
 */
double SSIM_Cal(char* argv1, char* argv2, char* argv3)
{
	// デフォルトさんですね
	const double K1 = 0.01;
	const double K2 = 0.03;
	const int L = 255;
	const double C1 = (K1 * L) * (K1 * L); //6.5025 C1 = (K(1)*L)^2;
	const double C2 = (K2 * L) * (K2 * L); //58.5225 C2 = (K(2)*L)^2;

	//　画像の宣言
	IplImage
		*img1=NULL, *img2=NULL, *img1_img2=NULL,
		*img1_temp=NULL, *img2_temp=NULL,
		*img1_L=NULL, *img2_L=NULL,
		*img1_sq=NULL, *img2_sq=NULL,
		*mu1=NULL, *mu2=NULL,
		*mu1_sq=NULL, *mu2_sq=NULL, *mu1_mu2=NULL,
		*sigma1_sq=NULL, *sigma2_sq=NULL, *sigma12=NULL,
		*ssim_map=NULL, *ssim_map_for_display=NULL, 
		*temp1=NULL, *temp2=NULL, *temp3=NULL;
	

	//　引数にある二つの画像を読み込み
	img1_temp = cvLoadImage(argv1,CV_LOAD_IMAGE_GRAYSCALE);
	img2_temp = cvLoadImage(argv2,CV_LOAD_IMAGE_GRAYSCALE);
	if(img1_temp==NULL || img2_temp==NULL)
		return -1;
	//　縦横のサイズ取得
	int x=img1_temp->width, y=img1_temp->height;

	//　ディスタンス（ダウンサンプリング）の設定、短い方の辺を256にする
	int rate_downsampling = std::max(1, int((std::min(x,y) / 256) + 0.5));

	//　チャンネルは1で固定、DEPTHは途中計算のため32F
	int nChan=1, d=IPL_DEPTH_32F;

	//　ダウンサンプリング前のサイズ
	CvSize size_L = cvSize(x, y);
	//　ダウンサンプリング後のサイズ
	CvSize size = cvSize(x / rate_downsampling, y / rate_downsampling);

	//　画像の格納先確保
	img1_L = cvCreateImage( size_L, d, nChan);
	img2_L = cvCreateImage( size_L, d, nChan);
	img1 = cvCreateImage( size, d, nChan);
	img2 = cvCreateImage( size, d, nChan);

	//　8bit整数で読み込み→32bit浮動小数点数にスケール変換
	cvConvert(img1_temp, img1_L);
	cvConvert(img2_temp, img2_L);
	//　tempの解放
	cvReleaseImage(&img1_temp);
	cvReleaseImage(&img2_temp);

	//　ダウンサンプリング
	cvResize(img1_L, img1);
	cvResize(img2_L, img2);
	//　Lの解放
	cvReleaseImage(&img1_L);
	cvReleaseImage(&img2_L);

	//　2乗の格納先確保	
	img1_sq = cvCreateImage( size, d, nChan);
	img2_sq = cvCreateImage( size, d, nChan);
	img1_img2 = cvCreateImage( size, d, nChan);
	
	//　2乗の計算
	cvPow( img1, img1_sq, 2 );
	cvPow( img2, img2_sq, 2 );
	cvMul( img1, img2, img1_img2, 1 );

	//　μの格納先確保
	mu1 = cvCreateImage( size, d, nChan);
	mu2 = cvCreateImage( size, d, nChan);
	mu1_sq = cvCreateImage( size, d, nChan);
	mu2_sq = cvCreateImage( size, d, nChan);
	mu1_mu2 = cvCreateImage( size, d, nChan);

	//　σの格納先確保
	sigma1_sq = cvCreateImage( size, d, nChan);
	sigma2_sq = cvCreateImage( size, d, nChan);
	sigma12 = cvCreateImage( size, d, nChan);

	//　tempの格納先確保
	temp1 = cvCreateImage( size, d, nChan);
	temp2 = cvCreateImage( size, d, nChan);
	temp3 = cvCreateImage( size, d, nChan);

	//ssimマップの格納先確保
	ssim_map = cvCreateImage( size, d, nChan);


	//////////////////////////////////////////////////////////////////////////
	// 予備的な計算ですね

	//　μの計算時に、ガウシアンフィルターを適用
	cvSmooth( img1, mu1, CV_GAUSSIAN, 11, 11, 1.5 );
	cvSmooth( img2, mu2, CV_GAUSSIAN, 11, 11, 1.5 );
	
	//　μを2乗
	cvPow( mu1, mu1_sq, 2 );
	cvPow( mu2, mu2_sq, 2 );
	cvMul( mu1, mu2, mu1_mu2, 1 );

	//　σの計算時に、ガウシアンフィルターを適用
	cvSmooth( img1_sq, sigma1_sq, CV_GAUSSIAN, 11, 11, 1.5 );
	cvAddWeighted( sigma1_sq, 1, mu1_sq, -1, 0, sigma1_sq );
	
	cvSmooth( img2_sq, sigma2_sq, CV_GAUSSIAN, 11, 11, 1.5 );
	cvAddWeighted( sigma2_sq, 1, mu2_sq, -1, 0, sigma2_sq );

	cvSmooth( img1_img2, sigma12, CV_GAUSSIAN, 11, 11, 1.5 );
	cvAddWeighted( sigma12, 1, mu1_mu2, -1, 0, sigma12 );
	

	//////////////////////////////////////////////////////////////////////////
	// 計算しましょう

	// (2*mu1_mu2 + C1)
	cvScale( mu1_mu2, temp1, 2 );
	cvAddS( temp1, cvScalarAll(C1), temp1 );

	// (2*sigma12 + C2)
	cvScale( sigma12, temp2, 2 );
	cvAddS( temp2, cvScalarAll(C2), temp2 );

	// ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
	cvMul( temp1, temp2, temp3, 1 );

	// (mu1_sq + mu2_sq + C1)
	cvAdd( mu1_sq, mu2_sq, temp1 );
	cvAddS( temp1, cvScalarAll(C1), temp1 );

	// (sigma1_sq + sigma2_sq + C2)
	cvAdd( sigma1_sq, sigma2_sq, temp2 );
	cvAddS( temp2, cvScalarAll(C2), temp2 );

	// ((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2))
	cvMul( temp1, temp2, temp1, 1 );

	// ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2))
	cvDiv( temp3, temp1, ssim_map, 1 );

	//　32bit浮動小数点数→8bit整数にスケール変換
	d=IPL_DEPTH_8U;
	ssim_map_for_display = cvCreateImage( size, d, nChan);
	cvConvertScale(ssim_map, ssim_map_for_display, 255.0, 0.0);
	//　ssimマップを保存
	cvSaveImage(argv3, ssim_map_for_display);
	//　画像の解放
	cvReleaseImage(&ssim_map_for_display);
	
	//　ssimマップをスカラーに変換
	CvScalar index_scalar = cvAvg( ssim_map );

	//　画像の解放
	cvReleaseImage(&img1);
	cvReleaseImage(&img2);
	cvReleaseImage(&ssim_map);

	return index_scalar.val[0];
}

int main(int argc, char** argv)
{
	if(argc!=4) {
		std::cout << "こんな感じでお願いよ　SSIM.exe a.tif b.tif map.tif" << std::endl;
		return -1;
	}
	double SSIM = SSIM_Cal(argv[1],argv[2],argv[3]);
	// 結果の出力
	std::cout << "SSIMは" << SSIM << "でした。" << std::endl;

	return 0;
}

電子化

SSIMを計算するプログラムのソース