// グレイスケール画像を幾何変換するプログラム
// 内挿方法：最近隣内挿法
// プログラム名：affine_nn.cpp

#include "eps-header.hh"

enum enum_coord {
	COORD_ORIGIN,		// 処理済画像の座標系は，元画像の座標系と同じ(処理済画像は，元画像と同じ大きさ)
	COORD_MINMAX};		// 処理済画像の位置にあわせて，はみ出ないように座標系を決める

// 　　元画像の座標 (sx, sy)
// 処理済画像の座標 (dx, dy)
//
//  幾何変換行列　（同次座標表記）
//   ┌       ┐
//   │a  b  c│
//   │d  e  f│
//   │0  0  1│
//   └       ┘
//
//（１）元画像の座標を処理済画像の座標に幾何変換
//   ┌  ┐   ┌    ┐┌  ┐   ┌ ┐
//   │dx│   │a  b││sx│   │c│
//   │  │ = │    ││  │ + │ │
//   │dy│   │d  e││sy│   │f│
//   └  ┘   └    ┘└  ┘   └ ┘
//
//（２）処理済画像の座標を元画像の座標に幾何変換
//   ┌  ┐              ┌     ┐┌      ┐
//   │sx│        1     │ e -b││dx - c│
//   │  │ = -----------│     ││      │
//   │sy│   (a*e - b*d)│-d  a││dy - f│
//   └  ┘              └     ┘└      ┘

// ３×３行列と同次座標を乗算する関数
void mul_coord(
	double mat[3][3],			// ３×３行列(入力)
	double   sx, double   sy,	// 同次座標(入力)
	double * dx, double * dy)	// 乗算した座標(出力)
{
	* dx = mat[0][0] * sx + mat[0][1] * sy + mat[0][2];
	* dy = mat[1][0] * sx + mat[1][1] * sy + mat[1][2];
}

// 元画像(src)を幾何変換して，処理済画像(dst)に入れる関数
// 内挿方法：最近隣内挿法
void affine_nn(
	class image& src,			// 元画像クラス(入力)
	class image& dst,			// 処理済画像クラス(出力)
	double affine_mat[3][3],	// 幾何変換行列(入力)
	enum  enum_coord coord)		// 座標系の指定(入力)
{
	if (src.p == NULL) {
		printf("元画像がありません。");
		return;
	}

	printf("\nグレイスケール画像を幾何変換します。\n");
	printf("内挿方法：最近隣内挿法\n");

	int nr_src_lines  = src.get_nr_lines();		// 元画像のライン数   (縦の画素数)
	int nr_src_pixels = src.get_nr_pixels();	// 元画像のピクセル数 (横の画素数)

	// 元画像の四隅の座標を幾何変換し，その最大値と最小値を求める。
	double dx, dx_min, dx_max;
	double dy, dy_min, dy_max;

	mul_coord(affine_mat, 0.0, 0.0, &dx_min, &dy_min);
	dx_max = dx_min;
	dy_max = dy_min;

	mul_coord(affine_mat, 0.0, (double)(nr_src_lines), &dx, &dy);
	if (dx < dx_min) dx_min = dx;	if (dx > dx_max) dx_max = dx;
	if (dy < dy_min) dy_min = dy;	if (dy > dy_max) dy_max = dy;

	mul_coord(affine_mat, (double)(nr_src_pixels), 0.0, &dx, &dy);
	if (dx < dx_min) dx_min = dx;	if (dx > dx_max) dx_max = dx;
	if (dy < dy_min) dy_min = dy;	if (dy > dy_max) dy_max = dy;

	mul_coord(affine_mat, (double)(nr_src_pixels), (double)(nr_src_lines), &dx, &dy);
	if (dx < dx_min) dx_min = dx;	if (dx > dx_max) dx_max = dx;
	if (dy < dy_min) dy_min = dy;	if (dy > dy_max) dy_max = dy;

	printf("幾何変換行列：\n");
	printf("+-                            -+\n");
	printf("|  %8.4lf %8.4lf %8.4lf  |\n", affine_mat[0][0], affine_mat[0][1], affine_mat[0][2]);
	printf("|  %8.4lf %8.4lf %8.4lf  |\n", affine_mat[1][0], affine_mat[1][1], affine_mat[1][2]);
	printf("|  %8.4lf %8.4lf %8.4lf  |\n", affine_mat[2][0], affine_mat[2][1], affine_mat[2][2]);
	printf("+-                            -+\n");
	printf("幾何変換前の画像のピクセル座標の範囲：[%8.2lf, %8.2lf]\n", 0.0, (double)(nr_src_pixels));
	printf("　　　　　　　　　　ライン座標の範囲：[%8.2lf, %8.2lf]\n", 0.0, (double)(nr_src_lines));
	printf("幾何変換後の画像のピクセル座標の範囲：[%8.2lf, %8.2lf]\n", dx_min, dx_max);
	printf("　　　　　　　　　　ライン座標の範囲：[%8.2lf, %8.2lf]\n", dy_min, dy_max);

	int nr_dst_lines;
	int nr_dst_pixels;
	if (coord == COORD_ORIGIN) {
		printf("\n元画像と同じ座標系で幾何変換画像を作成します。\n");
		nr_dst_pixels = nr_src_pixels;
		nr_dst_lines  = nr_src_lines;
		dx_min = 0.0;
		dy_min = 0.0;

	} else if (coord == COORD_MINMAX) {
		printf("\n幾何変換後すべての画素がはみ出ないように幾何変換画像を作成します。\n");
		nr_dst_lines  = (int)(dy_max - dy_min + 0.5);
		if (nr_dst_lines == 0) nr_dst_lines = 1;
		nr_dst_pixels = (int)(dx_max - dx_min + 0.5);
		if (nr_dst_pixels == 0) nr_dst_pixels = 1;
	}
	printf("作成する画像の大きさは，横 %d ピクセル，縦 %d ラインになります。\n", nr_dst_pixels, nr_dst_lines);

	dst.gmake(nr_dst_lines, nr_dst_pixels);			// 処理済画像の配列を確保

	double det = affine_mat[0][0] * affine_mat[1][1] - affine_mat[0][1] * affine_mat[1][0];
	if (fabs(det) < 1.0e-20) {
		printf("幾何変換行列の逆行列が存在しないため，幾何変換できませんでした。\n");
		return;
	}

	// 最近隣内挿法によって幾何変換する。
	for (int iy = 0; iy < nr_dst_lines; iy ++) {
		for (int ix = 0; ix < nr_dst_pixels; ix ++) {
			double dxx = (ix + dx_min) - affine_mat[0][2];
			double dyy = (iy + dy_min) - affine_mat[1][2];
			double sx = (  affine_mat[1][1] * dxx - affine_mat[0][1] * dyy) / det;
			double sy = (- affine_mat[1][0] * dxx + affine_mat[0][0] * dyy) / det;
			int isx = (int)(sx + 0.5);
			int isy = (int)(sy + 0.5);
			if (isx < 0 || isx >= nr_src_pixels ||
			    isy < 0 || isy >= nr_src_lines) {
				dst.p[iy][ix] = MIN_GRAY_LEVEL;
			} else {
				// 最近隣内挿法
				dst.p[iy][ix] = src.p[isy][isx];
			}
		}
	}
}

// プログラムの開始位置
int main()
{
	class image src;	// 元画像クラスを宣言
	class image dst;	// 処理済画像クラスを宣言

	// 幾何変換行列：平行移動(tx, ty)
	double tx = 10.0;
	double ty = 50.0;
	double mat_move[3][3] = {
		{1.0, 0.0,  tx},
		{0.0, 1.0,  ty},
		{0.0, 0.0, 1.0}
	};

	// 幾何変換行列：拡大・縮小率(sx, sy)
	double sx = 4.0;
	double sy = 4.0;
	double mat_scale[3][3] = {
		{ sx, 0.0, 0.0},
		{0.0,  sy, 0.0},
		{0.0, 0.0, 1.0}
	};

	// 幾何変換行列：回転角(th) [rad]
	double th = 45.0 * PI / 180.0;
	double mat_rotate[3][3] = {
		{ cos(th), sin(th), 0.0},
		{-sin(th), cos(th), 0.0},
		{     0.0,     0.0, 1.0}
	};

	// 幾何変換行列：反転
	double ax = 0.0;	// 	反転の基準点座標
	double ay = 0.0;	//
	double mat_trans[3][3] = {
		{-1.0,  0.0, 2.0 * ax},
		{ 0.0, -1.0, 2.0 * ay},
		{ 0.0,  0.0,      1.0}
	};

	// 幾何変換行列：スキュー
	double th_h = 15.0 * PI / 180.0;	// 水平方向のスキュー角 [rad]
	double th_v = 35.0 * PI / 180.0;	// 垂直方向のスキュー角 [rad]
	double mat_skew[3][3] = {
		{       1.0, tan(th_h), 0.0},
		{ tan(th_v),       1.0, 0.0},
		{       0.0,       0.0, 1.0}
	};

	src.gload();		// 元画像をBMPファイルから読み込む
						// BMPファイルを指定したい場合は，
						//   src.gload("パス付きBMPファイル名");
						// とする。

	// 元画像(src)を幾何変換して，処理済画像(dst)に入れる。
	// 内挿方法：最近隣内挿法
	affine_nn(
		src,			// 元画像クラス(入力)
		dst,			// 処理済画像クラス(出力)
		mat_rotate,		// 幾何変換行列(入力)
		COORD_ORIGIN);	// 処理済画像の座標系の指定(入力)	元画像と同じにするとき       COORD_ORIGIN
						//									処理済画像の大きさにするとき COORD_MINMAX

	dst.gsave();		// 処理済画像をBMPファイルに保存する
						// BMPファイルを指定したい場合は，
						//   dst.gsave("パス付きBMPファイル名");
						// とする。
}