/**
 *
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include "dct.h"

#include "qcif.h"
#include "quantizer.h"

/****************************************************************************/

/* Pi */
#define PI 3.14159265

/****************************************************************************/

#ifndef FALSE
#define FALSE 0
#endif

#ifndef TRUE
#define TRUE 1
#endif

/****************************************************************************/

/* cosine coefficient values */
#define C1 (cos (PI       / 16))
#define C2 (cos ((2 * PI) / 16))
#define C3 (cos ((3 * PI) / 16))
#define C4 (cos ((4 * PI) / 16))
#define C5 (cos ((5 * PI) / 16))
#define C6 (cos ((6 * PI) / 16))
#define C7 (cos ((7 * PI) / 16))

/* DCT addition nodes, form of N<col><row> */
#define N11 (f[0] + f[7])
#define N12 (f[1] + f[6])
#define N13 (f[2] + f[5])
#define N14 (f[3] + f[4])
#define N15 (f[3] - f[4])
#define N16 (f[2] - f[5])
#define N17 (f[1] - f[6])
#define N18 (f[0] - f[7])

#define N21 (N11 + N14)
#define N22 (N12 + N13)
#define N23 (N12 - N13)
#define N24 (N11 - N14)
#define N25 (N15)
#define N26 ((N17 - N16) * C4)
#define N27 ((N17 + N16) * C4)
#define N28 (N18)

#define N31 ((N21 + N22) * C4)
#define N32 ((N21 - N22) * C4)
#define N33 (N23 * C6 + N24 * C2)
#define N34 (N24 * C6 - N23 * C2)
#define N35 (N25 + N26)
#define N36 (N25 - N26)
#define N37 (N28 - N27)
#define N38 (N28 + N27)

/* IDCT addition nodes, form of IN<col><row> */
#define IN31 (F[0])
#define IN32 (F[4])
#define IN33 (F[2])
#define IN34 (F[6])
#define IN35 (F[1] * C7 - F[7] * C1)
#define IN36 (F[5] * C3 - F[3] * C5)
#define IN37 (F[5] * C5 + F[3] * C3)
#define IN38 (F[1] * C1 + F[7] * C7)

#define IN21 ((IN31 + IN32) * C4)
#define IN22 ((IN31 - IN32) * C4)
#define IN23 (IN33 * C6 - IN34 * C2)
#define IN24 (IN34 * C6 + IN33 * C2)
#define IN25 (IN35 + IN36)
#define IN26 (IN35 - IN36)
#define IN27 (IN38 - IN37)
#define IN28 (IN38 + IN37)

#define IN11 (IN21 + IN24)
#define IN12 (IN22 + IN23)
#define IN13 (IN22 - IN23)
#define IN14 (IN21 - IN24)
#define IN15 (IN25)
#define IN16 ((IN27 - IN26) * C4)
#define IN17 ((IN27 + IN26) * C4)
#define IN18 (IN28)

/****************************************************************************/

/* 1-D DCT across a row or column, horizontal should be TRUE or FALSE */
void compute_1d_dct (DctBlock *block, int horizontal, int row_col)
{
	double f[8];
	double F2[8];
	int    i;
	
	for (i = 0; i < 8; i++)
	{
		if (horizontal)
			f[i] = block->f[row_col][i];
		else
			f[i] = block->f[i][row_col];
	}
	
	/* calculate transforms using common addition nodes */
	F2[0] = N31;
	F2[4] = N32;
	F2[2] = N33;
	F2[6] = N34;
	F2[1] = N35 * C7 + N38 * C1;
	F2[5] = N36 * C3 + N37 * C5;
	F2[3] = N37 * C3 - N36 * C5;
	F2[7] = N38 * C7 - N35 * C1;
	
	for (i = 0; i < 8; i++)
	{
		if (horizontal)
			block->f[row_col][i] = F2[i] / 2.0;
		else
			block->f[i][row_col] = F2[i] / 2.0;
	}
}

/* 1-D IDCT across a row or column, horizontal should be TRUE or FALSE */
void compute_1d_idct (DctBlock *block, int horizontal, int row_col)
{
	double f[8];
	double F[8];
	int    i;
	
	for (i = 0; i < 8; i++)
	{
		if (horizontal)
			F[i] = block->f[row_col][i] * 2.0;
		else
			F[i] = block->f[i][row_col] * 2.0;
	}
	
	/* calculate transforms using common addition nodes */
	f[0] = IN11 + IN18;
	f[1] = IN12 + IN17;
	f[2] = IN13 + IN16;
	f[3] = IN14 + IN15;
	f[4] = IN14 - IN15;
	f[5] = IN13 - IN16;
	f[6] = IN12 - IN17;
	f[7] = IN11 - IN18;
	
	/* note that when converting back, we scale by a factor of 2/N = 2/8 */
	for (i = 0; i < 8; i++)
	{
		if (horizontal)
			block->f[row_col][i] = f[i] / 4.0;
		else
			block->f[i][row_col] = f[i] / 4.0;
	}
}

/****************************************************************************/

void dct_block_free (DctBlock *block)
{
	free (block);
}

DctBlock *dct_block_init ()
{
	DctBlock *block   = (DctBlock *)malloc (sizeof (DctBlock));
	int       i, j;
	
	for (i = 0; i < 8; i++)
	{
		for (j = 0; j < 8; j++)
			block->f[i][j] = 0.0;
	}
	
	return block;
}

DctBlock *dct_block_new_char (char *row0, char *row1, char *row2, char *row3, 
                              char *row4, char *row5, char *row6, char *row7, 
                              int offset)
{
	char   *crow[8];
	double  drow[8][8];
	int     i, j;
	
	/* for looping convenience */
	crow[0] = row0;
	crow[1] = row1;
	crow[2] = row2;
	crow[3] = row3;
	crow[4] = row4;
	crow[5] = row5;
	crow[6] = row6;
	crow[7] = row7;
	
	for (i = 0; i < 8; i++)
	{
		for (j = 0; j < 8; j++)
			drow[i][j] = (double)crow[i][j];
	}
	
	return dct_block_new_double (drow[0], drow[1], drow[2], drow[3], 
	                             drow[4], drow[5], drow[6], drow[7], 
	                             offset);
}

DctBlock *dct_block_new_double (double *row0, double *row1, double *row2, 
                                double *row3, double *row4, double *row5, 
                                double *row6, double *row7, int offset)
{
	DctBlock *block   = dct_block_init ();
	double   *drow[8];
	int       i, j;
	
	/* for looping convenience */
	drow[0] = row0;
	drow[1] = row1;
	drow[2] = row2;
	drow[3] = row3;
	drow[4] = row4;
	drow[5] = row5;
	drow[6] = row6;
	drow[7] = row7;
	
	for (i = 0; i < 8; i++)
	{
		for (j = 0; j < 8; j++)
			block->f[i][j] = drow[i][j + offset];
	}
	
	return block;
}

void dct_block_print (DctBlock *block)
{
	int i;
	
	printf ("DCT Block:\n");
	
	for (i = 0; i < 8; i++)
		printf ("[ %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f ]\n", 
		        block->f[i][0], block->f[i][1], block->f[i][2], 
		        block->f[i][3], block->f[i][4], block->f[i][5], 
				block->f[i][6], block->f[i][7]);
}

void dct_block_quantize (DctBlock *block)
{
	int i, j;
	
	for (i = 0; i < 8; i++)
	{
		for (j = 0; j < 8; j++)
			block->f[i][j] = quantize (block->f[i][j]);
	}
}

/****************************************************************************/

DctFrame *dct_frame_init ()
{
	return (DctFrame *)malloc (sizeof (DctFrame));
}

DctFrame *dct_frame_new (double **components)
{
	
}

/****************************************************************************/

void compute_dct (DctBlock *block)
{
	int i;
	
	/* computer 1-D horizontal DCTs */
	for (i = 0; i < 8; i++)
		compute_1d_dct (block, TRUE, i);
	
	/* computer 1-D vertical DCTs */
	for (i = 0; i < 8; i++)
		compute_1d_dct (block, FALSE, i);
}

void compute_idct (DctBlock *block)
{
	int i;

	/* compute 1-D vertical IDCTs */
	for (i = 0; i < 8; i++)
		compute_1d_idct (block, FALSE, i);

	/* compute 1-D horizontal IDCTs */
	for (i = 0; i < 8; i++)
		compute_1d_idct (block, TRUE, i);
}

/****************************************************************************/