wavelet.c

c++ Builder中网络图象渐进传输实例

// --------------- SINGLE-CHANNEL ROUTINES ---------------
// wavelet image compression
// encodes channels (which necessarily have to have
// power of 2 dimensions) using a wavelet-transform and
// bit-plane run-length coding into a progressive bit-stream
// (c) 2001-2003 Daniel Vollmer ([e-mail])

// 2.7	-	removed MMX optimisations for wavelet transforms and made code even faster ;)
//		-	removed unused MaxBits parameter from wv_init_channel
//		-	changed bitstream format (order in which bits are written) and removed writing
//			unneccessary zeros at the end of each block
//		-	changed yuv transform slightly (Cr / CB are now centred around 0, not 128), as we're writing
//			the sign in any case
//		-	changed colorspace conversion to be in-place
//		-	fixed bug in raw_load if file was too small
//		-	misc optimisations to bit-files
//		-	added wv_ prefix to log2i, mse_to_psnr, psnr_to_mse
//		-	changed the # of iterations for the multi-channel size selector (now bails out earlier)
//		-	new function to return the header of an image (wv_read_header), changed layout of t_wv_dchannels
//		-	changed decompression to work for (hopefully) all invalid data w/o overwriting
//			anything in memory
//		-	wv_init_decode_channels now accepts an extra reduction parameter (return a scaled down version of image)
//			(see -dr parameter in wako.exe)
// 2.6	-	reordering now clears unused coefficients as well
//		-	init_channel now accepts a NEGATIVE Lossless parameter to indicate that the
//			(faster) approximate error estimation is to be used
//		-	optimised the wavelet-transforms (including MMX versions)
//		-	wrote seperate size-counter for rice-encoding
// 2.5	-	changed reordering again slightly, now blocks (VD, HD, DD) are transmitted
//			completely interleaved, reordering lowest level unnecessarily fixed
//		-	fixed bug with blocks consisting of all 0s
// 2.1	-	file slightly smaller -> reordered '0' coefficients to the end, incompatible

/* 	This program is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA	02111-1307	USA
*/

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

#define wv_MIN_HEADER_VERSION 1
#define wv_CUR_HEADER_VERSION 1
#define wv_MIN_BITSTREAM_VERSION 1
#define wv_CUR_BITSTREAM_VERSION 1

#define DETAIL(ptr, pitch) ((ptr)[0] -= ((ptr)[-(pitch)] + (ptr)[(pitch)]) / 2)
#define DETAIL_LAST(ptr, pitch) ((ptr)[0] -= (ptr)[-(pitch)])
#define AVERAGE(ptr, pitch) ((ptr)[-(pitch)] += ((ptr)[0] + (ptr)[-2 * (pitch)]) / 4)
#define AVERAGE_FIRST(ptr, pitch) ((ptr)[-(pitch)] += (ptr)[0] / 2)

static void wavelet_row(wv_pel *dst, const int pitch, const int count)
{ // biorthogonal (2,2) Cohen-Daubechies-Feauveau (lifting)
	int i;
	wv_pel *ptr;
	
	ptr = dst + pitch; // odd
	
	DETAIL(ptr, pitch);
	AVERAGE_FIRST(ptr, pitch);
	for (i = count - 2; i > 0; i--)
	{
		ptr += 2 * pitch;
		DETAIL(ptr, pitch);
		AVERAGE(ptr, pitch);
	}
	ptr += 2 * pitch;
	DETAIL_LAST(ptr, pitch);
	AVERAGE(ptr, pitch);
}

static void wavelet_even_no_row_first(wv_pel *dst, const int pitch, const int rowpitch, const int count)
{ // biorthogonal (2,2) Cohen-Daubechies-Feauveau (lifting)
	int i;
	wv_pel *rowptr;
	
	rowptr = dst - rowpitch;
	
	for (i = 2 * count; i > 0; i--)
	{
		// column (first)
		DETAIL(rowptr, rowpitch);
		AVERAGE_FIRST(rowptr, rowpitch);
		rowptr += pitch;
	}
}

static void wavelet_even_no_row(wv_pel *dst, const int pitch, const int rowpitch, const int count)
{ // biorthogonal (2,2) Cohen-Daubechies-Feauveau (lifting)
	int i;
	wv_pel *rowptr;
	
	rowptr = dst - rowpitch;
	
	for (i = 2 * count; i > 0; i--)
	{
		// column (generic)
		DETAIL(rowptr, rowpitch);
		AVERAGE(rowptr, rowpitch);
		rowptr += pitch;
	}
}

static void wavelet_no_row_last(wv_pel *dst, const int pitch, const int rowpitch, const int count)
{ // biorthogonal (2,2) Cohen-Daubechies-Feauveau (lifting)
	int i;
	wv_pel *rowptr;
	
	rowptr = dst;
	
	for (i = 2 * count; i > 0; i--)
	{
		// column (last)
		DETAIL_LAST(rowptr, rowpitch);
		AVERAGE(rowptr, rowpitch);
		rowptr += pitch;
	}
}

// dst must contain the original (untransformed) image
// also width & height must be divisible by 2^(levels+1)
static void wavelet_transform(wv_pel *dst, const int width, const int height, const int levels)
{
	int l;

	for (l = 0; l < levels; l++)
	{
		int wmid = width >> (l + 1), hmid = height >> (l + 1);
		int pitch = 1 << l, rowpitch = width << l;
		int i;

		// precalculate the first three rows
		wavelet_row(dst, pitch, wmid);
		wavelet_row(dst + pitch * width, pitch, wmid);
		wavelet_row(dst + 2 * pitch * width, pitch, wmid);
		// do the third row (and the 1st two columns)
		wavelet_even_no_row_first(dst + 2 * pitch * width, pitch, rowpitch, wmid);
		for (i = 3 * pitch; i < height - pitch; i += 2 * pitch)
		{
			wavelet_row(dst + i * width, pitch, wmid); // odd row
			wavelet_row(dst + (i + pitch) * width, pitch, wmid); // even row
			wavelet_even_no_row(dst + (i + pitch) * width, pitch, rowpitch, wmid); // even row + previous two columns
		}
		wavelet_row(dst + i * width, pitch, wmid); // last (odd) row
		wavelet_no_row_last(dst + i * width, pitch, rowpitch, wmid); // and the last two columns
	}
}


#define INV_AVERAGE_FIRST(ptr, pitch) ((ptr)[0] -= (ptr)[(pitch)] / 2)
#define INV_AVERAGE(ptr, pitch) ((ptr)[0] -= ((ptr)[-(pitch)] + (ptr)[(pitch)]) / 4)
#define INV_DETAIL(ptr, pitch) ((ptr)[-(pitch)] += ((ptr)[-2 * (pitch)] + (ptr)[0]) / 2)
#define INV_DETAIL_LAST(ptr, pitch) ((ptr)[(pitch)] += (ptr)[0])

static void inverse_wavelet_row(wv_pel *dst, const int pitch, const int count)
{
	int i;
	wv_pel *ptr;
	
	ptr = dst; // even (average)
	
	INV_AVERAGE_FIRST(ptr, pitch);
	for (i = count - 1; i > 0; i--)
	{
		ptr += 2 * pitch;
		INV_AVERAGE(ptr, pitch);
		INV_DETAIL(ptr, pitch);
	}
	INV_DETAIL_LAST(ptr, pitch);
}

static void inverse_wavelet_no_row_first(wv_pel *dst, const int pitch, const int rowpitch, const int count)
{
	int i;
	wv_pel *rowptr;
	
	rowptr = dst;
	
	for (i = 2 * count; i > 0; i--)
	{
		// column (first)
		INV_AVERAGE_FIRST(rowptr, rowpitch);
		rowptr += pitch;
	}
}

static void inverse_wavelet_no_row_evenm1(wv_pel *dst, const int pitch, const int rowpitch, const int count)
{
	int i;
	wv_pel *rowptr;
	
	// column
	rowptr = dst;
	
	for (i = 2 * count; i > 0; i--)
	{
		// column
		INV_AVERAGE(rowptr, rowpitch);
		INV_DETAIL(rowptr, rowpitch);
		rowptr += pitch;
	}
}

static void inverse_wavelet_no_row_last(wv_pel *dst, const int pitch, const int rowpitch, const int count)
{
	int i;
	wv_pel *rowptr;
	
	// column
	rowptr = dst;
	
	for (i = 2 * count; i > 0; i--)
	{
		// column
		INV_DETAIL_LAST(rowptr, rowpitch);
		rowptr += pitch;
	}
}


// dst must contain the wavelet transform of the image
// also width & height must be divisible by 2^(levels+1)
static void inverse_wavelet_transform(wv_pel *dst, const int width, const int height, const int levels)
{
	int l;

	for (l = levels - 1; l >= 0; l--)
	{
		int wmid = width >> (l + 1), hmid = height >> (l + 1);
		int pitch = 1 << l, rowpitch = width << l;
		int i;

		// preculate first row (w/ vertical transform)
		inverse_wavelet_no_row_first(dst, pitch, rowpitch, wmid);
		for (i = 2 * pitch; i < height; i += 2 * pitch)
		{
			inverse_wavelet_no_row_evenm1(dst + i * width, pitch, rowpitch, wmid);
			inverse_wavelet_row(dst + (i - 2 * pitch) * width, pitch, wmid);
			inverse_wavelet_row(dst + (i - pitch) * width, pitch, wmid);
		}
		inverse_wavelet_no_row_last(dst + (i - 2 * pitch) * width, pitch, rowpitch, wmid);
		inverse_wavelet_row(dst + (i - 2 * pitch) * width, pitch, wmid);
		inverse_wavelet_row(dst + (i - pitch) * width, pitch, wmid);
	}
}


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

		for (i = 0; i < num; i++)
			*dst++ /= quant;
	}
}*/


/*int dead_zone_dequant(const int val, const int fac)
{ // has been multiplied by 10 to avoid floating point calculations entirely
	if (val > 0)
		return ((val * 10 + 5) * fac) / 10; // (val + 0.5) * fac
	else if (val < 0)
		return ((val * 10 - 5) * fac) / 10; // (val - 0.5) * fac
	else 
		return 0;
}

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

		// start dequantizing w/ T1
		for (i = 0; i < num; i++)
			*dst++ = (wv_pel)dead_zone_dequant(*dst, quant);
	}
}*/

#define dead_zone_dequant(v, f) (((v) > 0) ? ((((v) << 1) + 1) * (f)) >> 1 : ((v) < 0) ? -((((-(v) << 1) + 1) * (f)) >> 1) : 0)

// for the why and how of the function see "wavelet_analysis.nb". why it isn't exact (which i think it should be), ... i don't know!
static float estimate_error(const wv_pel *a, const int width, const int owidth, const int height, const int oheight, const int level, const int quant)
{
	if (quant > 1)
	{
		if (a)
		{
			int block_size_x = width >> level; // how big is a block ("level") of coefficients
			int block_size_y = height >> level; // in x and y?
			int unused_size_x = (width - owidth) >> level;
			int unused_size_y = (height - oheight) >> level;
			int pitch = 1 << level;
			int HDofs = (pitch >> 1), VDofs = (pitch >> 1) * width;
			int y, x;
			double sum2 = 0.0;

			#define OFS(x, y) ((y) * pitch * width + (x) * pitch)
			#define HD(ofs) a[(ofs) + HDofs]
			#define VD(ofs) a[(ofs) + VDofs]
			#define DD(ofs) a[(ofs) + HDofs + VDofs]
			#define ERR(v) ((v) - dead_zone_dequant((v) / quant, quant))
			
			for (y = 0; y < block_size_y - unused_size_y; y++)
			{
				#define HDI 0
				#define VDI 1
				#define DDI 2
				#define LI 0
				#define CI 1
				#define RI 2
				
				int yli = max(0, y - 1), yri = min(block_size_y - 1, y + 1); // y index left, right
				int errs[DDI + 1][RI + 1][RI + 1];
				int ofs;

				// we copy the error in the coefficients into this array and then only update the ones we need (after shifting them)
				// this way then inner loop only references fixed offsets into errs
				ofs = OFS(0, yli);
				errs[HDI][LI][CI] = errs[HDI][LI][RI] = ERR(HD(ofs));
				errs[VDI][LI][CI] = errs[VDI][LI][RI] = ERR(VD(ofs));
				errs[DDI][LI][CI] = errs[DDI][LI][RI] = ERR(DD(ofs));
				ofs = OFS(0, y);
				errs[HDI][CI][CI] = errs[HDI][CI][RI] = ERR(HD(ofs));
				errs[VDI][CI][CI] = errs[VDI][CI][RI] = ERR(VD(ofs));
				errs[DDI][CI][CI] = errs[DDI][CI][RI] = ERR(DD(ofs));
				ofs = OFS(0, yri);
				errs[HDI][RI][CI] = errs[HDI][RI][RI] = ERR(HD(ofs));
				errs[VDI][RI][CI] = errs[VDI][RI][RI] = ERR(VD(ofs));
				errs[DDI][RI][CI] = errs[DDI][RI][RI] = ERR(DD(ofs));
						
				for (x = 0; x < block_size_x - unused_size_x; x++)
				{
					int xri = min(block_size_x - 1, x + 1); // right
					double diff;

					// shift
					errs[HDI][LI][LI] = errs[HDI][LI][CI];
					errs[HDI][CI][LI] = errs[HDI][CI][CI];
					errs[HDI][RI][LI] = errs[HDI][RI][CI];

					errs[HDI][LI][CI] = errs[HDI][LI][RI];
					errs[HDI][CI][CI] = errs[HDI][CI][RI];
					errs[HDI][RI][CI] = errs[HDI][RI][RI];
					// shift
					errs[VDI][LI][LI] = errs[VDI][LI][CI];
					errs[VDI][CI][LI] = errs[VDI][CI][CI];
					errs[VDI][RI][LI] = errs[VDI][RI][CI];

					errs[VDI][LI][CI] = errs[VDI][LI][RI];
					errs[VDI][CI][CI] = errs[VDI][CI][RI];
					errs[VDI][RI][CI] = errs[VDI][RI][RI];
					// shift
					errs[DDI][LI][LI] = errs[DDI][LI][CI];
					errs[DDI][CI][LI] = errs[DDI][CI][CI];
					errs[DDI][RI][LI] = errs[DDI][RI][CI];

					errs[DDI][LI][CI] = errs[DDI][LI][RI];
					errs[DDI][CI][CI] = errs[DDI][CI][RI];
					errs[DDI][RI][CI] = errs[DDI][RI][RI];
					
					// fill
					ofs = OFS(xri, yli);
					errs[HDI][LI][RI] = ERR(HD(ofs));
					errs[VDI][LI][RI] = ERR(VD(ofs));
					errs[DDI][LI][RI] = ERR(DD(ofs));
					ofs = OFS(xri, y);
					errs[HDI][CI][RI] = ERR(HD(ofs));
					errs[VDI][CI][RI] = ERR(VD(ofs));
					errs[DDI][CI][RI] = ERR(DD(ofs));
					ofs = OFS(xri, yri);
					errs[HDI][RI][RI] = ERR(HD(ofs));