// ２つのグレイスケール画像の周波数成分の振幅と位相を交換した２つの画像を生成するプログラム
// プログラム名：ex_amp_phase.cpp

#include "eps-header.hh"
#include "fft-header.hh"

// 元画像１，２(src1, src2)の周波数成分の位相と振幅を交換し，処理済画像１，２(dst1, dst2)に入れる関数
void exchange_amp_phase(
	class image& src1,		// 元画像１クラス(入力)
	class image& src2,		// 元画像２クラス(入力)
	class image& dst1,		// 元画像１の振幅と元画像２の位相をもつ画像クラス(出力)
	class image& dst2)		// 元画像１の位相と元画像２の振幅をもつ画像クラス(出力)
{
	if (src1.p == NULL) {
		printf("元画像１がありません。");
		return;
	}
	if (src2.p == NULL) {
		printf("元画像２がありません。");
		return;
	}

	if (src1.get_nr_lines() != src2.get_nr_lines()) {
		printf("元画像１と２のライン数が異なります。\n");
		return;
	}
	if (src1.get_nr_pixels() != src2.get_nr_pixels()) {
		printf("元画像１と２のピクセル数が異なります。\n");
		return;
	}

	printf("\n周波数成分の振幅と位相を交換した画像を生成します。\n");

	int src_nr_lines  = src1.get_nr_lines();	// 元画像１，２のライン数   (縦の画素数)
	int src_nr_pixels = src1.get_nr_pixels();	// 元画像１，２のピクセル数 (横の画素数)

	int dst_nr_lines  = src_nr_lines;			// 処理済画像のライン数
	int dst_nr_pixels = src_nr_pixels;			// 処理済画像のピクセル数
	dst1.gmake(dst_nr_lines, dst_nr_pixels);	// 処理済画像１（元画像１の振幅と元画像２の位相をもつ画像）の配列を確保
	dst2.gmake(dst_nr_lines, dst_nr_pixels);	// 処理済画像２（元画像１の位相と元画像２の振幅をもつ画像）の配列を確保

	class fourier src1_ft(src_nr_lines, src_nr_pixels);	// 元画像１用高速フーリエ変換クラスの宣言
	class fourier src2_ft(src_nr_lines, src_nr_pixels);	// 元画像１用高速フーリエ変換クラスの宣言

	// 実数部を元画像１，２の輝度値，虚数部を0.0とする複素数を２次元配列に入れる。
	for (int iy = 0; iy < src_nr_lines; iy ++) {
		for (int ix = 0; ix < src_nr_pixels; ix ++) {
			src1_ft.p[iy][ix] = creal8((double)src1.p[iy][ix], 0.0);
			src2_ft.p[iy][ix] = creal8((double)src2.p[iy][ix], 0.0);
		}
	}
	src1_ft.Forward_FFT();		// 元画像１を２次元高速フーリエ変換
	src2_ft.Forward_FFT();		// 元画像２を２次元高速フーリエ変換

	for (int iy = 0; iy < src_nr_lines; iy ++) {
		for (int ix = 0; ix < src_nr_pixels; ix ++) {
			double re1    = real(src1_ft.p[iy][ix]);
			double im1    = imag(src1_ft.p[iy][ix]);
			double amp1   = sqrt(re1 * re1 + im1 * im1);		// 元画像１の振幅
			double phase1 = arg(src1_ft.p[iy][ix]);				// 元画像１の位相

			double re2    = real(src2_ft.p[iy][ix]);
			double im2    = imag(src2_ft.p[iy][ix]);
			double amp2   = sqrt(re2 * re2 + im2 * im2);		// 元画像２の振幅
			double phase2 = arg(src2_ft.p[iy][ix]);				// 元画像２の位相

			src1_ft.p[iy][ix] = amp1 * exp(creal8(0.0, phase2));
			src2_ft.p[iy][ix] = amp2 * exp(creal8(0.0, phase1));
		}
	}

	src1_ft.Backward_FFT();		// ２次元高速フーリエ逆変換し処理済画像１とする。
	src2_ft.Backward_FFT();		// ２次元高速フーリエ逆変換し処理済画像２とする。

	src1_ft.MakeImage(
		OUTPUT_AMP,
		FREQ_SPACE_SHUFFLING_OFF,
		OUTPUT_SCALE_ON,
		dst1.p);

	src2_ft.MakeImage(
		OUTPUT_AMP,
		FREQ_SPACE_SHUFFLING_OFF,
		OUTPUT_SCALE_ON,
		dst2.p);
}

// プログラムの開始位置
int main()
{
	class image src1;	// 元画像１クラスを宣言
	class image src2;	// 元画像２クラスを宣言
	class image dst1;	// 処理済画像１（元画像１の振幅と元画像２の位相をもつ画像）クラスを宣言
	class image dst2;	// 処理済画像２（元画像１の位相と元画像２の振幅をもつ画像）クラスを宣言

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

	// 元画像１，２(src1, src2)の周波数成分の位相と振幅を交換し，処理済画像１，２(dst1, dst2)に入れる
	exchange_amp_phase(
		src1,		// 元画像１クラス(入力)
		src2,		// 元画像２クラス(入力)
		dst1,		// 元画像１の振幅と元画像２の位相をもつ画像クラス(出力)
		dst2);		// 元画像１の位相と元画像２の振幅をもつ画像クラス(出力)

	printf("\n元画像１の振幅と元画像２の位相をもつ画像を保存します。\n");
	dst1.gsave();		// 処理済画像１をBMPファイルに保存する
						// BMPファイルを指定したい場合は，
						//   dst1.gsave("パス付きBMPファイル名");
						// とする。
	printf("\n元画像１の位相と元画像２の振幅をもつ画像を保存します。\n");
	dst2.gsave();		// 処理済画像２をBMPファイルに保存する
						// BMPファイルを指定したい場合は，
						//   dst2.gsave("パス付きBMPファイル名");
						// とする。
}