/***
 * transform.c
 ***
 * soucast knihovny libopen
 * 2007-04-29
 * xbarin02@stud.fit.vutbr.cz  
 ***
 * soubor s fcemi, kt. provadeji transformace (2D (I)DWT a 2D (I)DCT)
 */

// pomocne deklarace
void fwt(double *const_input, double *output);
void fwt_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x);

/***
 * dopredna 2D DWT
 */ 
void fwt2D(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH])
{
	for(unsigned y=0; y<LENGTH; y++)
		fwt(const_input[y],output[y]);
	for(unsigned x=0; x<LENGTH; x++)
		fwt_col(output,output,x);
}

/***
 * dopredna 1D/sloupcova DWT (vyuziva radkovou)
 */ 
void fwt_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x)
{
	static double var_input[LENGTH];
	for(unsigned y=0; y<LENGTH; y++)
		var_input[y] = const_input[y][x];
	static double var_output[LENGTH];
	fwt(var_input,var_output);
	for(unsigned y=0; y<LENGTH; y++)
		output[y][x] = var_output[y];
}

/***
 * dopredna 1D/radkova DWT (haarova)
 */
void fwt(double *const_input, double *output)
{
	static double input[LENGTH];
	memcpy(input,const_input,sizeof(double)*LENGTH);
	for(int length = LENGTH >> 1; ; length >>= 1)
	{
		for(int i = 0; i < length; i++)
		{
			double sum = (input[i*2]+input[i*2+1])/2;
			double difference = (input[i*2]-input[i*2+1])/2;
			output[i] = sum;
			output[length+i] = difference;
		}
		if(length == 1) 
			return;
		memcpy(input,output,sizeof(double)*length);
	}
}

// pomocne deklarace
void iwt(double *const_input, double *output);
void iwt_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x);

/***
 * inverzni 2D DWT
 */ 
void iwt2D(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH])
{
	for(unsigned x=0; x<LENGTH; x++)
		iwt_col(const_input,output,x);
	for(unsigned y=0; y<LENGTH; y++)
		iwt(output[y],output[y]);
}

/***
 * inverzni 1D/sloupcova DWT (vyuziva radkovou)
 */ 
void iwt_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x)
{
	static double var_input[LENGTH];
	for(unsigned y=0; y<LENGTH; y++)
		var_input[y] = const_input[y][x];
	static double var_output[LENGTH];
	iwt(var_input,var_output);
	for(unsigned y=0; y<LENGTH; y++)
		output[y][x] = var_output[y];
}

/***
 * inverzni 1D/radkova DWT (haarova)
 */ 
void iwt(double *const_input, double *output)
{
	static double input[LENGTH];
	memcpy(input,const_input,sizeof(double)*LENGTH);
	for(int length = 1; ; length <<= 1)
	{
		for(int i = 0; i < length; i++)
		{
			double a = (input[i] - input[length+i]);
			double b = (input[i] + input[length+i]);
			output[i*2] = b;
			output[i*2+1] = a;
		}
		if(length == LENGTH >> 1) 
			return;
		memcpy(input,output,sizeof(double)*length*2);
	}
}

// ----------------------------------------------------------------------------

// pomocne deklarace
void fct(double *input, double *output);
void fct_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x);

/***
 * dopredna 2D DCT (optimalizovany posun)
 */ 
void fct2D(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH])
{
	for(unsigned y=0; y<LENGTH; y++)
		fct(const_input[y],output[y]);
	for(unsigned x=0; x<LENGTH; x++)
		fct_col(output,output,x);
	output[0][0] -= 8192.0;
}

/***
 * dopredna 1D DCT/sloupcova (vyuziva radkovou)
 */ 
void fct_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x)
{
	static double var_input[LENGTH];
	for(unsigned y=0; y<LENGTH; y++)
		var_input[y] = const_input[y][x];
	static double var_output[LENGTH];
	fct(var_input,var_output);
	for(unsigned y=0; y<LENGTH; y++)
		output[y][x] = var_output[y];
}

// pomocna tabulka pro optimalizaci (predpocita pomocne koeficienty)
double fct_tbl[LENGTH][LENGTH];

/***
 * dopredna 1D DCT/radkova (s optimalizaci, bez posunu)
 */ 
void fct(double *input, double *output)
{
	static int first = 0;
	if(0 == first)
	{
		for(int h=0; h<LENGTH; h++)
			for(int j=0; j<LENGTH; j++)
				fct_tbl[h][j] = cos( (M_PI/LENGTH)*h*(j+0.5) );
		first++;
	}

	for(int h=0; h<LENGTH; h++)
	{
		double sum = 0;
		for(int j=0; j<LENGTH; j++)
			sum += input[j]*fct_tbl[h][j];
		output[h] = sum;
	}
}

// pomocne deklarace
void ict(double *input, double *output);
void ict_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x);

/***
 * inverzni 2D DCT (opt. posun)
 */ 
void ict2D(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH])
{
	const_input[0][0] += 8192.0;
	double tmp_output[LENGTH][LENGTH];
	for(unsigned x=0; x<LENGTH; x++)
		ict_col(const_input,tmp_output,x);
	for(unsigned y=0; y<LENGTH; y++)
		ict(tmp_output[y],output[y]);
}

/***
 * inverzni 1D/sloupcova DCT
 */ 
void ict_col(double const_input[LENGTH][LENGTH], double output[LENGTH][LENGTH], int x)
{
	static double var_input[LENGTH];
	for(unsigned y=0; y<LENGTH; y++)
		var_input[y] = const_input[y][x];
	static double var_output[LENGTH];
	ict(var_input,var_output);
	for(unsigned y=0; y<LENGTH; y++)
		output[y][x] = var_output[y];
}

// pro predpocitane pomocne koeficienty
double ict_tbl[LENGTH][LENGTH];

/***
 * inverzni 1D DCT (opt., bez posunu)
 */ 
void ict(double *input, double *output)
{
	static int first = 0;
	if(0 == first)
	{
		for(int h=0; h<LENGTH; h++)
			for(int j=1; j<LENGTH; j++)
				ict_tbl[h][j] = cos( (M_PI/LENGTH)*j*(h+0.5) );
		first++;
	}

	for(int h=0; h<LENGTH; h++)
	{
		double sum = 0;
		for(int j=1; j<LENGTH; j++)
			sum += input[j]*ict_tbl[h][j];
		sum += 0.5*input[0];
		sum *= 2/(double)LENGTH;
		output[h] = sum;
	}
}

/***
 * ladici fce pro hromadny posun matice
 */ 
void shift2D(double input2D[LENGTH][LENGTH], double output2D[LENGTH][LENGTH], double shift)
{
	for(unsigned y=0; y<LENGTH; y++)
	{
		for(unsigned x=0; x<LENGTH; x++)
			output2D[y][x] = input2D[y][x] + shift;
	}
}

/***
 * ladici fce pro vypis matice (printf)
 */ 
void print(double a[LENGTH][LENGTH])
{
	for(unsigned y=0; y<LENGTH; y++)
	{
		for(unsigned x=0; x<LENGTH; x++)
			printf("%8.3f ",a[y][x]);
		printf("\n");
	}
	printf("\n");
}