wavelet.c

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

					errs[VDI][RI][RI] = ERR(VD(ofs));
					errs[DDI][RI][RI] = ERR(DD(ofs));
					

					// x even, y even
					diff = (errs[DDI][LI][LI] + errs[DDI][LI][CI] + errs[DDI][CI][LI] + errs[DDI][CI][CI]
						- 4 * errs[HDI][CI][LI] - 4 * errs[HDI][CI][CI]
						- 4 * errs[VDI][LI][CI] - 4 * errs[VDI][CI][CI]) / 16.0;
					sum2 += diff * diff;
					
					// x odd, y even
					diff = (errs[DDI][LI][LI] + errs[DDI][LI][RI] - 6 * errs[DDI][LI][CI] + errs[DDI][CI][LI] + errs[DDI][CI][RI] - 6 * errs[DDI][CI][CI]
						- 4 * errs[HDI][CI][LI] - 4 * errs[HDI][CI][RI] + 24 * errs[HDI][CI][CI]
						- 4 * errs[VDI][LI][RI] - 4 * errs[VDI][LI][CI] - 4 * errs[VDI][CI][RI] - 4 * errs[VDI][CI][CI]) / 32.0;
					sum2 += diff * diff;

					// x even, y odd
					diff = (errs[DDI][LI][LI] + errs[DDI][LI][CI] + errs[DDI][RI][LI] + errs[DDI][RI][CI] - 6 * errs[DDI][CI][LI] - 6 * errs[DDI][CI][CI]
						- 4 * errs[HDI][RI][LI] - 4 * errs[HDI][RI][CI] - 4 * errs[HDI][CI][LI] - 4 * errs[HDI][CI][CI]
						- 4 * errs[VDI][LI][CI] - 4 * errs[VDI][RI][CI] + 24 * errs[VDI][CI][CI]) / 32.0;
					sum2 += diff * diff;
					
					// x odd, y odd
					diff = (errs[DDI][LI][LI] + errs[DDI][LI][RI] - 6 * errs[DDI][LI][CI] + errs[DDI][RI][LI] + errs[DDI][RI][RI] + - 6 * errs[DDI][RI][CI] - 6 * errs[DDI][CI][LI] - 6 * errs[DDI][CI][RI] + 36 * errs[DDI][CI][CI]
						- 4 * errs[HDI][RI][LI] - 4 * errs[HDI][RI][RI] + 24 * errs[HDI][RI][CI] - 4 * errs[HDI][CI][LI] - 4 * errs[HDI][CI][RI] + 24 * errs[HDI][CI][CI]
						- 4 * errs[VDI][LI][RI] - 4 * errs[VDI][LI][CI] - 4 * errs[VDI][RI][RI] - 4 * errs[VDI][LI][CI] + 24 * errs[VDI][CI][RI] + 24 * errs[VDI][CI][CI]) / 64.0;
					sum2 += diff * diff;
				}
				#undef HDI
				#undef VDI
				#undef DDI
				#undef LI
				#undef CI
				#undef RI
			}
			#undef OFS
			#undef HD
			#undef VD
			#undef DD
			#undef ERR
			return (float)(sum2 / (float)((block_size_x - unused_size_x) * (block_size_y - unused_size_y) * 4));
		}
		else // estimate errors
		{
			int numq = quant / 2 - 1;
			double est_error;

			est_error = (double)quant / 4.0; // 1 / q * (1 / 2)^2

			// Sum[i^2, {i, n}] = 1/6 n (1 + n) (1 + 2n)
			est_error += (double)(numq * (1 + numq) * (1 + 2 * numq)) / (3.0 * (double)quant);

			return (float)est_error;
		}
	}
	return 0.0f;
}



//----------------------------------------------------------------------------
// wv_log2i		returns the integer logarithm to the base 2 of a given integer
//				(i.e. the number of bits needed to encode that number)
// max_val		number whose log is wanted (log2i(0) = 0)
//----------------------------------------------------------------------------
WAVELET_DLL_API int WAVELET_DLL_CC wv_log2i(int max_val)
{
	int i;

	for (i = 0; max_val > 0; max_val >>= 1, i++) ;
	return i;
}


// this converts a reorder-index into an (x,y) coordinate to fetch from
// using bit-interleaving (fractal z-curve)
static void get_zig_zag_idx(const int idx, const int maxx, const int maxy, int* xx, int* yy)
{
	int xbits = wv_log2i(maxx) - 1;
	int ybits = wv_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--;
			}
		}
	}
//	printf("x: %i y: %i -> idx: %i x: %i, y: %i\n", x, y, i, *xx, *yy);
}

static int* make_single_table(const int width, const int owidth, const int height, const int oheight, const int l, const int r)
{
	int block_size_x = width >> (l + r); // how big is a block ("level") of coefficients
	int block_size_y = height >> (l + r); // in x and y?
	int *tab = malloc((1 + block_size_x * block_size_y) * sizeof *tab);
	int *cur_table = tab + 1;
	int unused_size_x = (width - owidth) >> (l + r), unused_size_y = (height - oheight) >> (l + r);
	int i, j, pitch;
			
	pitch = 1 << l;
	tab[0] = 0; // number of non-zeroed entries
	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 * pitch * (width >> r) + j * pitch;
			tab[0]++; // increase number of non-zeroed entries
		}
	}
	// now add the unused ones (which should be all zero, anyway)
	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 * pitch * (width >> r) + j * pitch;
	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 * pitch * (width >> r) + j * pitch;
	return tab;
}

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

	if (width > 2 && height > 2)
	{
		int l;
		
		tables = malloc(levels * sizeof *tables);
		for (l = 1; l <= levels; l++)
			tables[levels - l] = make_single_table(width, owidth, height, oheight, l, reduce);
	}
	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, num;
		const int *idx;
		
		// output base average
		num = tables[0][0];
		idx = tables[0] + 1;
		for (i = 0; i < num; i++)
			*dst++ = src[*idx++];
		memset(dst, 0, ((width >> levels) * (height >> levels) - num) * sizeof *dst);
		dst += (width >> levels) * (height >> levels) - num;
		
		for (l = levels; l >= 1; l--)
		{
			int block_size_x = width >> l; // how big is a block ("level") of coefficients
			int block_size_y = height >> l; // in x and y?
			int pitch = 1 << (l - 1);

			// output HD coefficients 
			num = tables[levels - l][0];
			idx = tables[levels - l] + 1;
			for (i = 0; i < num; i++)
			{
				*dst++ = src[*idx + pitch]; // HD coeff in row by row order
				*dst++ = src[*idx + pitch * width]; // VD coeff in row by row order
				*dst++ = src[*idx + pitch * width + pitch]; // DD coeff in row by row order
				idx++;
			}
			memset(dst, 0, 3 * (block_size_y * block_size_x - num) * sizeof *dst);
			dst += 3 * (block_size_y * block_size_x - num);
		}
	}
	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;

		// output base average
		idx = tables[0] + 1;
		for (i = 0; i < (width >> levels) * (height >> levels); i++)
			dst[*idx++] = *src++;
		for (l = levels; l >= 1; l--)
		{
			int block_size_x = width >> l; // how big is a block ("level") of coefficients
			int block_size_y = height >> l; // in x and y?
			int pitch = 1 << (l - 1);

			// output HD coefficients 
			idx = tables[levels - l] + 1;
			for (i = 0; i < block_size_y * block_size_x; i++)
			{
				dst[*idx + pitch] = *src++; // HD coeff in "z" order
				dst[*idx + pitch * width] = *src++; // VD coeff in in "z" order
				dst[*idx + pitch * width + pitch] = *src++; // DD coeff in in "z" order
				idx++;
			}
		}
	}
	else
		memcpy(dst, src, width * height * sizeof *dst);
}


//----------------------------------------------------------------------------
// wv_mse_to_psnr	converts mse (mean square error) into psnr (peak signal to
//				noise ratio (in dB))
//				returns 0.0 if mse <= 0.0f (instead of "infinity")
// mse			mean square error to be expressed as psnr
//----------------------------------------------------------------------------
WAVELET_DLL_API float WAVELET_DLL_CC wv_mse_to_psnr(const float mse)
{
	if (mse > 0.0f)
		return (float)(20.0 * log10(255.0 / sqrt(mse)));
	return 0.0f;
}


//----------------------------------------------------------------------------
// wv_psnr_to_mse	converts psnr (peak signal to noise ratio (in dB)) into a
//				mse (mean square error)
// psnr			psnr to be expressed as mean square error
//----------------------------------------------------------------------------
WAVELET_DLL_API float WAVELET_DLL_CC wv_psnr_to_mse(const float psnr)
{
	if (psnr > 0.0f)
		return (float)pow(255.0 / pow(10.0, psnr / 20.0), 2.0);
	return 65025.0f;
}


//----------------------------------------------------------------------------
// wv_calc_psnr	calculates the difference between the two given channels
//				and returns its psnr
// a			channel a
// b			channel b
// width		width of the rectangle to be compared
// height		height of the rectangle to be compared
// pitch		pitch of the channels in wv_pel units
//				mse (mean square error)
// pmse			address of a float where the mean-square error is written into
//				(if != NULL)
//----------------------------------------------------------------------------
WAVELET_DLL_API float WAVELET_DLL_CC wv_calc_psnr(const wv_pel *a, const wv_pel *b, const int width, const int height, const int pitch, float *pmse)
{
	double sum2 = 0.0f;
	int i, j;

	for (i = 0; i < height; i++)
		for (j = 0; j < width; j++)
		{
			double diff = (double)(a[i * pitch + j] - b[i * pitch + j]);

			sum2 += diff * diff;
		}

	if (sum2 > 0.0)
	{
		float mse = (float)(sum2 / (width * height)); // mean squared error
		
		if (pmse)
			*pmse = mse;
		return wv_mse_to_psnr(mse);
	}
	else
		*pmse = 0.0f;
	return 0.0f;
}


// --------------- RICE ENCODING ROUTINES ---------------

#define RICE_MIN_K 0
#define RICE_MAX_K 30
// 2 * 4 bits for num_bp and skip_bp, and 3 bits for a 1 + k + sign
#define RICE_MIN_SPACE (2 * 4 + 3)

static int rice_encode_size(const wv_pel* src, int num, int k, const int imax, int* bit_stat)
{
	int bits_written = 0;
	int cur_bp, max_bp;

	if (num == 0)
		return 0;
	
	if (imax > 0)
	{
		max_bp = cur_bp = wv_log2i(imax) - 1;
		bits_written += 4; // bit_write(num_bp, 4, f);
		if (max_bp >= 0)
		{
			int last_refinement, next_refinement;

			bits_written += 4; // bit_write(skip_bp, 4, f);
			last_refinement = 0;

			for (; cur_bp >= 0; cur_bp--)
			{
				int	mask = 1 << cur_bp;
				int ref = (2 << cur_bp) - 1;
				int num_zeros = 0;
				int idx;

				next_refinement = 0;
				for (idx = 0; idx < num; idx++)
				{
					int cur = abs(src[idx]);
					
					if (cur <= ref)
					{
						if ((cur & mask) == 0)
							num_zeros++;
						else
						{
							while (num_zeros >= (1 << k))
							{
								bits_written++; // += bit_write_single(0, f); // emit '0'
								num_zeros -= (1 << k);
								k = min(RICE_MAX_K, k + 1); // dynamically adapt k
							}
							bits_written += 1 + k + 1;
							num_zeros = 0;
							k = max(RICE_MIN_K, k - 1); // dynamically adapt k
							next_refinement++;
						}
					}
				}
				
				while (num_zeros >= (1 << k))
				{
					bits_written++; // += bit_write_single(0, f); // emit '0'
					num_zeros -= (1 << k);
					k = min(RICE_MAX_K, k + 1); // dynamically adapt k
				}
				if (num_zeros > 0)
				{
					bits_written++; // += bit_write_single(0, f); // emit '0'
					// we have indicated too many zeros, but in the decoder we know when a bitplane is done
					num_zeros = 0;
					k = max(RICE_MIN_K, k - 1); // dynamically adapt k