#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include "wavelet.h"


/*首先必須注意的是，在小波變換的時候采用的是串行工作方式
（即，小波系數存放在原圖象的存儲區內），因此不得不對小波系數重新排序。
這樣一來，對于高速緩存的利用是十分低效率的，卻可以節省不少的空間。*/
 
static void wavelet_row(wv_pel *dst, const wv_pel *src, const int len)
{ 
	//利用的是雙正交小波提升
	int i, mid;
	const wv_pel *ptr;
	wv_pel *avg, *det;

	mid = len / 2;

	ptr = src;
	avg = dst; //直流平均值 
	det = dst + mid; //高頻細節空間 
	
	*det = ptr[1] - (ptr[2] + ptr[0]) / 2; //細節預處理 
	*avg = ptr[0] + (det[0] + det[0]) / 4;  
	ptr += 2; det++; avg++;

	for (i = 1; i < mid - 1; i++)
	{
	//左邊不存在一個奇鄰居，右邊提升兩次
		*det = ptr[1] - (ptr[2] + ptr[0]) / 2;
		*avg = ptr[0] + (det[0] + det[-1]) / 4;
		ptr += 2; det++; avg++;
	}
	//右邊不存在一個偶鄰居，左邊提升兩次 
	*det = ptr[1] - (ptr[0] + ptr[0]) / 2; 
	*avg = ptr[0] + (det[0] + det[-1]) / 4; //最近平均值計算 
}


static void inverse_wavelet_row(wv_pel *dst, const wv_pel *src, const int len)
{
	int i, mid;
	wv_pel *ptr;
	const wv_pel *avg, *det;

	mid = len / 2;

	ptr = dst;
	avg = src;  
	det = src + mid; 
	//邊界情況
	*ptr = avg[0] - det[0] / 2; 
	ptr += 2; det++; avg++;
	
	for (i = 0; i < mid - 1; i++)
	{
		ptr[0] = avg[0] - (det[0] + det[-1]) / 4;
		ptr[-1] = det[-1] + (ptr[0] + ptr[-2]) / 2;
		ptr += 2; det++; avg++;
	}
	ptr[-1] = det[-1] + ptr[-2]; //邊界 
}

 
static void wavelet_transform(wv_pel *dst, const int width, const int height, const int levels)
{
	int l;
	wv_pel *tmp_lastrow = malloc((width + height + height) * sizeof *tmp_lastrow); // 寬 
	wv_pel *tmp_column_in = tmp_lastrow + width; // 高 
	wv_pel *tmp_column_out = tmp_lastrow + width + height;  //高 
	
	for (l = 0; l < levels; l++)
	{
		int w = width >> l, h = height >> l;
		int i;
		
		// 行 
		wavelet_row(tmp_lastrow, dst + (h - 1) * width, w); // put last row in tmp_row
		// 將結果存放在一行之下，并且工作方式是從底端往上 
		for (i = (h - 2) * width; i >= 0; i -= width)
			wavelet_row(dst + i + width, dst + i, w);

		// 列 
		for (i = 0; i < w; i++)
		{
			int j;

			for (j = 1; j < h; j++)
				tmp_column_in[j - 1] = dst[j * width + i]; //矩陣轉置 
			tmp_column_in[h - 1] = tmp_lastrow[i]; //放入最近一行 
			wavelet_row(tmp_column_out, tmp_column_in, h);
			for (j = 0; j < h; j++)
				dst[j * width + i] = tmp_column_out[j]; //再轉置 
		}
	}
	free(tmp_lastrow);
}

 
static void inverse_wavelet_transform(wv_pel *dst, const int width, const int height, const int levels)
{
	int l;
	wv_pel *tmp_lastrow = malloc((width + height + height) * sizeof *tmp_lastrow);  
	wv_pel *tmp_column_in = tmp_lastrow + width;  
	wv_pel *tmp_column_out = tmp_lastrow + width + height; 
	
	for (l = levels - 1; l >= 0; l--)
	{ 
		int w = width >> l, h = height >> l;
		int i;

		//列 
		for (i = 0; i < w; i++)
		{
			int j;

			for (j = 0; j < h; j++)
				tmp_column_in[j] = dst[j * width + i];  
			inverse_wavelet_row(tmp_column_out, tmp_column_in, h);
			for (j = 0; j < h - 1; j++)
				dst[(j + 1) * width + i] = tmp_column_out[j];  
			tmp_lastrow[i] = tmp_column_out[h - 1];  //存儲最后一列 
		}
		
		// 行 
		for (i = 0; i < (h - 1) * width; i += width)
			inverse_wavelet_row(dst + i, dst + i + width, w);
		inverse_wavelet_row(dst + (h - 1) * width, tmp_lastrow, w);
	}
	free(tmp_lastrow);
}


static void quantize(wv_pel* dst, const int num, const int T)
{
	if (T > 1)
	{
		int i;

		for (i = 0; i < num; i++)
			*dst++ /= T;
	}
}
 
#define dead_zone_dequant(v, f) (((v) > 0) ? (((v) * 10 + 5) * (f)) / 10 : ((v) < 0) ? (((v) * 10 - 5) * (f)) / 10 : 0)

static void dequantize(wv_pel* dst, const int num, const int T)
{
	if (T > 1)
	{
		int i;

		//以閾值T1開始小波系數恢復 
		for (i = 0; i < num; i++)
			*dst++ = dead_zone_dequant(*dst, T);
	}
}


static float log2(const double max_val)
{
	return (float)(log10(max_val) / log10(2.0));
}

/*標題5版式*****************************************************************
函數輸入：	int max_val：	對數函數的自變量（比如log2i(0) = 0）
功能描述：	對給定的一個整數作底為2的對數運算
函數輸出：	函數返回對數運算的取整結果
***************************************************************************/
WAVELET_DLL_API int WAVELET_DLL_CC log2i(int max_val)
{
	int i;
	for (i = 0; max_val > 0; max_val >>= 1, i++) ;
	return i;
}


 
//本函數利用比特交錯將重排序索引轉換稱(x,y)形式的坐標
static void get_zig_zag_idx(const int idx, const int maxx, const int maxy, int* xx, int* yy)
{
	int xbits = log2i(maxx) - 1;
	int ybits = log2i(maxy) - 1;
	int mask;

	mask = 1 << (xbits + ybits - 1);
	*xx = *yy = 0;
	if (ybits >= xbits)
	{
		while (ybits > 0)
		{
			*yy <<= 1;
			*yy |= (idx & mask) > 0;
			mask >>= 1;
			ybits--;
			if (xbits > 0)
			{
				*xx <<= 1;
				*xx |= (idx & mask) > 0;
				mask >>= 1;
				xbits--;
			}
		}
	}
	else
	{
		while (xbits > 0)
		{
			*xx <<= 1;
			*xx |= (idx & mask) > 0;
			mask >>= 1;
			xbits--;
			if (ybits > 0)
			{
				*yy <<= 1;
				*yy |= (idx & mask) > 0;
				mask >>= 1;
				ybits--;
			}
		}
	}
 
}

static int** init_reorder_tables(const int width, const int owidth, const int height, const int oheight, const int levels)
{
	int** tables = NULL;

	if (width > 2 && height > 2)
	{
		int l;
		
		tables = malloc((levels + 1) * sizeof *tables);
		for (l = 1; l <= levels + 1; l++)
		{
						//確定一塊小波系數在坐標值的大小
			int block_size_x = width >> l; 
			int block_size_y = height >> l; 
			int unused_size_x, unused_size_y;
			int *cur_table;
			int i, j;

			tables[l - 1] = cur_table = malloc(block_size_x * block_size_y * sizeof *tables[l - 1]);
						
			unused_size_x = (width - owidth) >> l;
			unused_size_y = (height - oheight) >> l;
			for (i = 0; i < block_size_x * block_size_y; i++)
			{
				int k;

				get_zig_zag_idx(i, block_size_x, block_size_y, &j, &k);
				if (j <= (block_size_x - unused_size_x) && k <= (block_size_y - unused_size_y))
					*cur_table++ = k * width + j;
			}
			// 現在求和所有未用的1（原本應該全部為0）  
			for (i = 0; i < min(block_size_y, block_size_y - unused_size_y + 1); i++)
				for (j = block_size_x - unused_size_x + 1; j < block_size_x; j++)
					*cur_table++ = i * width + j;
			for (i = block_size_y - unused_size_y + 1; i < block_size_y; i++)
				for (j = 0; j < block_size_x; j++)
					*cur_table++ = i * width + j;
		}
	}
	return tables;
}


static void free_reorder_tables(int** tables, const int levels)
{
	if (tables)
	{
		int i;

		for (i = 0; i <= levels; i++)
			free(tables[i]);
		free(tables);
	}
}

//本函數用來重排序小波變換系數，輸出是自定義的表
static void reorder(wv_pel *dst, const wv_pel* src, const int** tables, const int width, const int height, const int levels)
{
	if (tables)
	{
		int i, l;
		const int *idx;
		
		//輸出基頻均值 
		idx = tables[levels];
		for (i = 0; i < (width >> (levels + 1)) * (height >> (levels + 1)); i++)
			*dst++ = src[*idx++];
		for (l = levels + 1; l >= 1; l--)
		{
			//確定一塊小波系數在坐標值的大小
			int block_size_x = width >> l; 
			int block_size_y = height >> l; 

			//交替輸出LH端/HL端的小波系數
			idx = tables[l - 1];
			for (i = 0; i < block_size_y * block_size_x; i++)
			{
				*dst++ = src[*idx + block_size_x]; // LH端的小波系數以逐行順序存放 
				*dst++ = src[*idx + block_size_y * width]; // HL小波系數也按逐行順序存放 
				idx++;
			}
			//現在是右邊下端的塊 
			idx = tables[l - 1];
			for (i = 0; i < block_size_x * block_size_y; i++)
				*dst++ = src[*idx++ + block_size_y * width + block_size_x];
		}
	}
	else
		memcpy(dst, src, width * height * sizeof *dst);
}

//本函數用于將重排序函數所自定義的交替順序恢復為直接可用的數據
static void unreorder(wv_pel* dst, const wv_pel* src, const int** tables, const int width, const int height, const int levels)
{
	if (tables)
	{
		int i, l;
		const int *idx;

		//輸出基頻均值 
		idx = tables[levels];
		for (i = 0; i < (width >> (levels + 1)) * (height >> (levels + 1)); i++)
			dst[*idx++] = *src++;
		for (l = levels + 1; l >= 1; l--)
		{
			//確定一塊小波系數在坐標值的大小
			int block_size_x = width >> l;  
			int block_size_y = height >> l;  

			//交替輸出LH端/HL端的小波系數 
			idx = tables[l - 1];
			for (i = 0; i < block_size_y * block_size_x; i++)
			{
				dst[*idx + block_size_x] = *src++; //LH端的小波系數按“鋸齒”順序存放 
				dst[*idx + block_size_y * width] = *src++; //HL端的系數同樣按“鋸齒”順序存放