#include "MRegressionLinear.h"

//-------------------------------------------------------------------------------
// コンストラクタ　RGB同時計算向けに係数次元を3倍する
CMRegressionLinear::CMRegressionLinear(void) : CMRegression(LINEARIZE_DIM * 3)
{
}

//-------------------------------------------------------------------------------
// デストラクタ
CMRegressionLinear::~CMRegressionLinear(void)
{
}

//-------------------------------------------------------------------------------
// 推定行列の算出に使う値を返す
double CMRegressionLinear::GetX(const CvMat *mat, const int sample, const int index)
{
	CvScalar v = cvGet2D(mat, sample, index / LINEARIZE_DIM);
	return pow(v.val[0], (double)(index % LINEARIZE_DIM));
}

//-------------------------------------------------------------------------------
// 係数の重回帰推定
bool CMRegressionLinear::CalcCoef(const CvMat *macbethX)
{
	// 無彩色パッチを抜き出す(19〜24番)
	CvMat rgbX = cvMat(NUM_PATCH_LINEARIZE, macbethX->cols, CV_64FC1, NULL);
	cvGetRows(macbethX, &rgbX, 18, 24);

	double reflectance[] = REFLECTANCE_LIST;
	CvMat y = cvMat(NUM_PATCH_LINEARIZE, 1, CV_64FC1, reflectance);
	CvMat x = cvMat(NUM_PATCH_LINEARIZE, 1, CV_64FC1, NULL);
	CvMat *coef = cvCreateMat(m_Dim, rgbX.cols, CV_64FC1);
	cvSetZero(coef);

	// 色別に処理
	m_Dim = LINEARIZE_DIM;	// 一時的に係数次元を1色分に落とす
	for (int col = 0; col < rgbX.cols; col ++)
	{
		// １色取り出す
		cvGetCol(&rgbX, &x, col);

		// 色別に係数の算出
		CALL(CMRegression::CalcCoef(&x, &y));

		// 算出した係数,誤差をコピーする
		for(int i = 0; i < m_Coef->rows; i ++) 
			cvmSet(coef, i + col * LINEARIZE_DIM, col, GcvmGet(m_Coef, i, 0));

		cvReleaseData(&x);
	}

	// RGB同時計算向けに係数次元を3倍に戻す
	m_Dim = LINEARIZE_DIM * 3;

	// 係数行列をメンバー変数に保存する
	SAFE_RELEASEMAT(m_Coef);
	m_Coef = (CvMat*)cvClone(coef);
#ifdef SHOW_REGRESSION_COEF
	GShowMat(m_Coef, "CoefLinearize", "%9.6f");

	// 誤差の計算
	double reflectance3[] = REFLECTANCE_LIST3;
	CvMat y3 = cvMat(NUM_PATCH_LINEARIZE, 3, CV_64FC1, reflectance3);
	CALL(CMRegression::CalcError(&rgbX, &y3));
	GShowMat(m_Error, "Estimation error", "%9.6f");
#endif // SHOW_REGRESSION_COEF

	// 行列の解放
	SAFE_RELEASEMAT(coef);
	cvReleaseData(&rgbX);

	return true;
}

//-------------------------------------------------------------------------------
// フィッティングカーブの算出
bool CMRegressionLinear::CalcCurve()
{
	// 入力値行列作成
	const int Steps = 256;
	CvMat *fitData = cvCreateMat(Steps, m_Dim, CV_64FC1);
	for (int sample = 0; sample < Steps; sample ++)
	{
		for (int i = 0; i < m_Dim; i++)
		{
			cvmSet(fitData, sample, i, 
				pow((double)sample, (double)(i % LINEARIZE_DIM)));
		}
	}

	// フィッティングカーブの算出
	CvMat *fitCurve = cvCreateMat(Steps, m_Coef->cols, CV_64FC1);
	cvMatMul(fitData, m_Coef, fitCurve);
	GShowMat(fitCurve, "Fitting Curve", "%7.4f");

	// 行列の解放
	SAFE_RELEASEMAT(fitData);
	SAFE_RELEASEMAT(fitCurve);

	return true;
}

//-------------------------------------------------------------------------------
// 変換した行列データを生成する
CvMat *CMRegressionLinear::GenConvert(const CvMat *data)
{
	// GetXを使わないで計算する
	return GenConvertNoGetX(data);
	// 従来の処理へ
//	return CMRegression::GenConvert(data);
	// RGBチャネル毎に計算する（遅い）
//	return GenConvertEachCh(data);
}

//-------------------------------------------------------------------------------
// 変換した行列データを生成する（RGB別に処理するバージョン：遅かった）
CvMat *CMRegressionLinear::GenConvertEachCh(const CvMat *data)
{
	// m_Dim を一時的に変更
	m_Dim = LINEARIZE_DIM;

	// 行列の確保
	CvMat *convert = cvCreateMat(data->rows, data->cols, CV_64FC1);
	CvMat *sCoef = cvCreateMat(LINEARIZE_DIM, 1, CV_64F);
	CvMat *sConvert = cvCreateMat(data->rows, 1, CV_64F);

	// 各色の処理
	for (int col = 0; col < COLOR; col ++)
	{
		// 係数抜き出し
		for (int i = 0; i < LINEARIZE_DIM; i ++)
			cvmSet(sCoef, i, 0, cvmGet(m_Coef, col * LINEARIZE_DIM + i, col));

		// データ抜き出し
		CvMat sData = cvMat(data->rows, 1, data->type);
		cvGetCol(data, &sData, col);

		// データ行列作成
		CvMat *matData = GenDataMat(&sData);

		// 推定値の算出
		cvMatMul(matData, sCoef, sConvert);

		// 推定値を出力行列に書き込む
		for (int i = 0; i < data->rows; i ++)
			cvmSet(convert, i, col, cvmGet(sConvert, i, 0));

		SAFE_RELEASEMAT(matData);
	}

	// m_Dim を戻す
	m_Dim = LINEARIZE_DIM * 3;

	// 計算に使用した行列の解放
	SAFE_RELEASEMAT(sCoef);
	SAFE_RELEASEMAT(sConvert);

	return convert;
}

//-------------------------------------------------------------------------------
// 変換した行列データを生成する（GetXを使わないバージョン）
CvMat *CMRegressionLinear::GenConvertNoGetX(const CvMat *data)
{
	assert(LINEARIZE_DIM == 4);

	// 行列の確保
	CvMat *convert = cvCreateMat(data->rows, data->cols, CV_64FC1);

	// 係数の取得
	double b0[COLOR], b1[COLOR], b2[COLOR], b3[COLOR];
	for (int col = 0; col < data->cols; col ++)
	{
		b0[col] = cvmGet(m_Coef, col * LINEARIZE_DIM    , col);
		b1[col] = cvmGet(m_Coef, col * LINEARIZE_DIM + 1, col);
		b2[col] = cvmGet(m_Coef, col * LINEARIZE_DIM + 2, col);
		b3[col] = cvmGet(m_Coef, col * LINEARIZE_DIM + 3, col);
	}

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
 	for (int sample = 0; sample < data->rows; sample ++)
	{
		for (int col = 0; col < data->cols; col ++)
		{
			CvScalar v = cvGet2D(data, sample, col);
//			double x = cvmGet(data, sample, col);
			double y = b0[col] + b1[col] * v.val[0] 
				+ b2[col] * v.val[0] * v.val[0]
				+ b3[col] * v.val[0] * v.val[0] * v.val[0];
			cvmSet(convert, sample, col, y);
		}
	}

	return convert;
}

//-------------------------------------------------------------------------------
// 画像を変換しデータを生成する（係数行列が計算済みであることが前提）
IplImage *CMRegressionLinear::GenConvert(const IplImage *dataImg)
{
	return CMRegression::GenConvert(dataImg);
}