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


/*首先必須注意的是，在小波變換的時(shí)候采用的是串行工作方式
（即，小波系數(shù)存放在原圖象的存儲(chǔ)區(qū)內(nèi)），因此不得不對(duì)小波系數(shù)重新排序。
這樣一來，對(duì)于高速緩存的利用是十分低效率的，卻可以節(jié)省不少的空間。*/
 
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; //高頻細(xì)節(jié)空間 
	
	*det = ptr[1] - (ptr[2] + ptr[0]) / 2; //細(xì)節(jié)預(yù)處理 
	*avg = ptr[0] + (det[0] + det[0]) / 4;  
	ptr += 2; det++; avg++;

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


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
		// 將結(jié)果存放在一行之下，并且工作方式是從底端往上 
		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]; //矩陣轉(zhuǎn)置 
			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]; //再轉(zhuǎn)置 
		}
	}
	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];  //存儲(chǔ)最后一列 
		}
		
		// 行 
		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開始小波系數(shù)恢復(fù) 
		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));
}

/*標(biāo)題5版式*****************************************************************
函數(shù)輸入：	int max_val：	對(duì)數(shù)函數(shù)的自變量（比如log2i(0) = 0）
功能描述：	對(duì)給定的一個(gè)整數(shù)作底為2的對(duì)數(shù)運(yùn)算
函數(shù)輸出：	函數(shù)返回對(duì)數(shù)運(yùn)算的取整結(jié)果
***************************************************************************/
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;
}


 
//本函數(shù)利用比特交錯(cuò)將重排序索引轉(zhuǎn)換稱(x,y)形式的坐標(biāo)
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++)
		{
						//確定一塊小波系數(shù)在坐標(biāo)值的大小
			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;
			}
			// 現(xiàn)在求和所有未用的1（原本應(yīng)該全部為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);
	}
}

//本函數(shù)用來重排序小波變換系數(shù)，輸出是自定義的表
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--)
		{
			//確定一塊小波系數(shù)在坐標(biāo)值的大小
			int block_size_x = width >> l; 
			int block_size_y = height >> l; 

			//交替輸出LH端/HL端的小波系數(shù)
			idx = tables[l - 1];
			for (i = 0; i < block_size_y * block_size_x; i++)
			{
				*dst++ = src[*idx + block_size_x]; // LH端的小波系數(shù)以逐行順序存放 
				*dst++ = src[*idx + block_size_y * width]; // HL小波系數(shù)也按逐行順序存放 
				idx++;
			}
			//現(xiàn)在是右邊下端的塊 
			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);
}

//本函數(shù)用于將重排序函數(shù)所自定義的交替順序恢復(fù)為直接可用的數(shù)據(jù)
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--)
		{
			//確定一塊小波系數(shù)在坐標(biāo)值的大小
			int block_size_x = width >> l;  
			int block_size_y = height >> l;  

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