wavelet.c

Jepg2000的压缩解压

 
#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端的系数同样按“锯齿”顺序存放