// グレイスケール画像の特定周波数成分を調整するプログラム
// プログラム名：fft_adjust.cpp

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

// 特定周波数成分の調整用構造体
struct st_freq_element {
	int		wave_no_line;		// ライン方向の波数
	int		wave_no_pixel;		// ピクセル方向の波数
	double	amp_factor;			// 元の振幅にかける係数
	double	amp_bias;			// 元の振幅に加える数値
};

// 元画像(src)の特定周波数領域成分を調整する関数
void img_adj_freq(
	class image& src,					// 元画像クラス(入力)
	class image& dst,					// 調整した画像クラス(出力)
	int	nr_freq,						// 特定周波数成分の数
	struct st_freq_element * adj_freq)	// 特定周波数成分の配列へのポインタ
{
	if (src.p == NULL) {
		printf("元画像がありません。");
		return;
	}

	printf("\n指定された周波数成分のノイズを画像に加えます。\n");
	int nr_lines  = src.get_nr_lines();		// 元画像のライン数   (縦の画素数)
	int nr_pixels = src.get_nr_pixels();	// 元画像のピクセル数 (横の画素数)
	dst.gmake(nr_lines, nr_pixels);			// 処理済画像の配列を確保

	class fourier data(nr_lines, nr_pixels);	// 高速フーリエ変換クラスの宣言

	// 実数部を元画像の輝度値，虚数部を0.0とする複素数を２次元配列に入れる。
	for (int iy = 0; iy < nr_lines; iy ++) {
		for (int ix = 0; ix < nr_pixels; ix ++) {
			data.p[iy][ix] = creal8((double)src.p[iy][ix], 0.0);
		}
	}

	data.Forward_FFT();		// ２次元高速フーリエ変換

	for (int n = 0; n < nr_freq; n ++) {
		printf("\n%4d番目の特定周波数成分(ライン方向波数:%4d, ピクセル方向波数：%4d, 元の振幅にかける係数:%lf，加える数値:%lf)を加えています。\n", 
			n,
			adj_freq[n].wave_no_line,
			adj_freq[n].wave_no_pixel,
			adj_freq[n].amp_factor,
			adj_freq[n].amp_bias);
		int    iy     = adj_freq[n].wave_no_line;
		int    ix     = adj_freq[n].wave_no_pixel;
		double factor = adj_freq[n].amp_factor;
		double bias   = adj_freq[n].amp_bias;
		if (iy >= 0) {
			iy = (iy % nr_lines);
		} else {
			iy = -((-iy) % nr_lines);
			if (iy !=0) iy += nr_lines;
		}

		if (ix >= 0) {
			ix = (ix % nr_pixels);
		} else {
			ix = -((-ix) % nr_pixels);
			if (ix != 0) ix += nr_pixels;
		}

		double amp_org = sqrt(norm(data.p[iy][ix]));
		double amp_dst = amp_org * factor + bias;
		printf("２次元周波数空間の配列（ライン方向:%4d, ピクセル方向:%4d）の振幅%lfを%lfに調整します。\n", iy, ix, amp_org, amp_dst);
		if (amp_org >= 1.0e-20) {
			data.p[iy][ix] *= (amp_dst / amp_org);
		} else {
			data.p[iy][ix] = creal8(amp_dst, 0.0);
		}
	}

	data.Backward_FFT();	// ２次元高速フーリエ逆変換
	for (int iy = 0; iy < nr_lines; iy ++) {
		for (int ix = 0; ix < nr_pixels; ix ++) {
			data.p[iy][ix] /= (nr_lines * nr_pixels);
		}
	}

	for (int iy = 0; iy < nr_lines; iy ++) {
		for (int ix = 0; ix < nr_pixels; ix ++) {
			double pp = sqrt(norm(data.p[iy][ix]));
			if (pp >= (double)MAX_GRAY_LEVEL) pp = (double)MAX_GRAY_LEVEL;
			if (pp <= (double)MIN_GRAY_LEVEL) pp = (double)MIN_GRAY_LEVEL;
			dst.p[iy][ix] = (int)(pp + 0.5);
		}
	}
}

// プログラムの開始位置
int main()
{
		// 実数部を対象とする場合，ライン方向の波数とピクセル方向の波数を，それぞれ符号を逆にした周波数成分を組みにする。
	struct st_freq_element adj_freq[] = {
		// ライン方向の波数，ピクセル方向の波数，振幅にかける係数，振幅に加える数値
		{                -2,                 +3,             10.0,             0.0},
		{                +2,                 -3,             10.0,             0.0},

		{                +4,                +10,              1.0,       1000000.0},
		{                -4,                -10,              1.0,       1000000.0}
	};
	int nr_freq = sizeof(adj_freq) / sizeof(struct st_freq_element);	// 周波数成分の数

	class image src;	// 元画像クラスを宣言
	class image dst;	// 処理済画像クラスを宣言

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

	// 元画像(src)の特定周波数領域成分を調整する。
	img_adj_freq(
		src,			// 元画像クラス(入力)
		dst,			// 調整した画像クラス(出力)
		nr_freq,		// 特定周波数成分の数
		adj_freq);		// 特定周波数成分の配列へのポインタ

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