wavelet.c

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

				}
				
				bits_written += last_refinement;
				last_refinement += next_refinement;
				
				if (bit_stat)
					bit_stat[max_bp - cur_bp] = bits_written;
			}
		}
	}
	else
	{
		bits_written += 4; // bit_write(0, 4, f); // number of bitplanes
		if (bit_stat)
			bit_stat[0] = bits_written;
	}
	
	return bits_written;
}

static int rice_encode_block(const wv_pel* src, const int num, const int skip_bp, int k, const int imax, t_bit_file* bf)
{
	int bits_written = 0;
	int max_bp, cur_bp;

	if (num == 0)
		return 0;
	
	if ((imax >> skip_bp) > 0)
	{
		max_bp = wv_log2i(imax >> skip_bp) - 1;
		bits_written += bit_write(max_bp + 1, 4, bf);
		if (max_bp >= 0)
		{
			unsigned int *refinement_cache = malloc(((num + 31) / 32) * sizeof *refinement_cache);

			bits_written += bit_write(skip_bp, 4, bf);
			for (cur_bp = wv_log2i(imax) - 1; cur_bp >= skip_bp; cur_bp--)
			{
				int	mask = 1 << cur_bp;
				int ref = (2 << cur_bp) - 1;
				int num_zeros = 0, num_refinement = 0;
				int idx, i;

				for (idx = 0; idx < num; idx++)
				{
					int cur = abs(src[idx]);

					if (cur <= ref)
					{ // not refinement
						if ((cur & mask) == 0)
							num_zeros++;
						else
						{
							i = 0;
							while (num_zeros >= (1 << k))
							{ // count the # of accumulated zeros bits
								i++;
								num_zeros -= (1 << k);
								k = min(RICE_MAX_K, k + 1); // dynamically adapt k
							}
							bits_written += bit_write(0, i, bf);
							bits_written += bit_write((1 << (k + 1)) + (num_zeros << 1) + (src[idx] < 0), 1 + k + 1, bf);
							num_zeros = 0;
							k = max(RICE_MIN_K, k - 1); // dynamically adapt k

							if (num_refinement > 0)
							{
								bits_written += num_refinement;
								for (i = 0; i < num_refinement >> 5; i++)
									bit_write(refinement_cache[i], 32, bf);
								bit_write(refinement_cache[num_refinement >> 5], num_refinement & 31, bf);
								num_refinement = 0;
							}
						}
					}
					else
					{
						refinement_cache[num_refinement >> 5] = (refinement_cache[num_refinement >> 5] << 1) + ((cur & mask) != 0);
						num_refinement++;
					}
				}
				
				i = 0;
				while (num_zeros >= (1 << k))
				{ // count the # of accumulated zeros bits
					i++;
					num_zeros -= (1 << k);
					k = min(RICE_MAX_K, k + 1); // dynamically adapt k
				}
				if (num_zeros > 0)
				{
					i++; // we have indicated too many zeros, but in the decoder we know when a bitplane is done
					k = max(RICE_MIN_K, k - 1); // dynamically adapt k
				}				
				bits_written += bit_write(0, i, bf);
				
				bits_written += num_refinement;
				for (i = 0; i < num_refinement >> 5; i++)
					bit_write(refinement_cache[i], 32, bf);
				bit_write(refinement_cache[num_refinement >> 5], num_refinement & 31, bf);
			}
			free(refinement_cache);
		}
	}
	else
		bits_written += bit_write(0, 4, bf); // number of bitplanes
	
	return bits_written;
}

static int rice_decode_block(wv_pel* dst, const int num, int k, t_bit_file* bf)
{
	if (num > 0 && dst && bf)
	{
		int cur_bp = 0;
		
		memset(dst, 0, num * sizeof *dst);

		cur_bp = bit_read(4, bf);
		if (cur_bp > 0)
		{
			int rlr_bits_left = num;
			int skip_bp;
			wv_pel sign_flag;

			skip_bp = bit_read(4, bf);
			for (cur_bp--; cur_bp >= 0; cur_bp--)
			{
				int num_zeros = 0;
				int	idx = 0;
				int refinement_left = num - rlr_bits_left;
				int cur_rlr_bits_left = rlr_bits_left;
				
				// arranged so that after all the bitplanes are read, this flag ends up in the sign-bit position (MSB)
				sign_flag = 1 << ((8 * sizeof (wv_pel)) - (1 + cur_bp));
				
				while (cur_rlr_bits_left > 0)
				{
					if (bit_read_single(bf) == 0)
					{
						num_zeros += 1 << k;
						cur_rlr_bits_left -= 1 << k;
						if (cur_rlr_bits_left < 0)
						{
							k = max(RICE_MIN_K, k - 1);
							break;
						}
						else
							k = min(RICE_MAX_K, k + 1);
					}
					else
					{
						int added_zeros;
						unsigned char encoded_sign;
						
						// read sign bit as well and increase # of 0 by 1 (2 as it's shifted >> 1), the extra zero will become the '1' bit
						added_zeros = bit_read(k + 1, bf) + 2;
						k = max(RICE_MIN_K, k - 1);
						encoded_sign = added_zeros & 1; // remember the sign
						added_zeros = min(added_zeros >> 1, cur_rlr_bits_left); // safety
						cur_rlr_bits_left -= added_zeros;
						rlr_bits_left--; // one less RLR-encoded bit
						num_zeros += added_zeros;

						for (; num_zeros > 0; idx++)
						{
							if (dst[idx] != 0)
								dst[idx] = (dst[idx] << 1) + bit_read_single(bf); // refinement
							else
								num_zeros--;
						}
						dst[idx - 1] = 0x01 + sign_flag * encoded_sign; // as we were looking for an extra zero which ended the loop
					}
				}
				// and read the last refinement bits (no need to write zeros ;))
				for (; idx < num; idx++)
					if (dst[idx] != 0)
						dst[idx] = (dst[idx] << 1) + bit_read_single(bf);
			}
			
			sign_flag--; // set all the bits below the mask
			if (skip_bp > 0)
			{ // have to dequantize as well
				int i;

				skip_bp = 1 << skip_bp;
				for (i = 0; i < num; i++)
				{
					int val = dst[i] & sign_flag;

					val = (val > 0) ? ((val + val + 1) * skip_bp) >> 1 : 0;
					if (dst[i] < 0)
						dst[i] = -val; // set sign
					else
						dst[i] = val;
				}
			}
			else
			{ // only need to "properly" set the sign
				int i;

				for (i = 0; i < num; i++)
					if (dst[i] < 0)
						dst[i] = (wv_pel)(-(dst[i] & sign_flag)); // set sign
			}
			return num;
		}
	}
	return 0;
}


// returns the next largest float
#define flarger(fl) ((fl) == 0.0f ? FLT_MIN : (fl) + (fl) * FLT_EPSILON)

//----------------------------------------------------------------------------
// wv_init_channel	initialises a channel for wavelet encoding and returns
//				a handle to it (used for finding compression parameters)
// Width		the _original_ width of the channel (also see Channel)
// Height		the _original_ height of the channel (also see Channel)
// Channel		pixel data of the channel. this _has_ to be stored with a
//				width and height that are the next-bigger power of 2 to Width
//				and Height given above. a suitable extension should be done
//				on the boundary. Unused wavelet coefficients are set to 0.
//				e.g. source channel 720x480, then Width = 720
//				and Height = 480, but the data in Channel has to be 1024x512!
// Lossless		should be set to > 0 if you only want to compress
//				losslessly (greatly speeds up computations)
//				if < 0, the error computation is approximated (but allows
//				for lossy compression [at reduced quality in some cases])
// ReorderTable	pointer to a pointer to the reorder-table. if the pointer the
//				pointers points to (huh?) is NULL, a new table is generated
//				and its address written (so we can share a single table)
// Progress		function callback used for progress display, can be NULL
// UserData		passed to the Progress function
//----------------------------------------------------------------------------
WAVELET_DLL_API t_wv_cchannel* WAVELET_DLL_CC wv_init_channel(const int Width, const int Height, const wv_pel* Channel,
	const int Lossless, int*** ReorderTable, wv_progress_function Progress, void* UserData)
{
	t_wv_cchannel* cc = NULL;

	if (Width >= 1 && Height >= 1 && Channel)
	{
		int i;
		wv_pel *channel;
		int cur, end;

		cc = malloc(sizeof *cc);

		cc->owidth = Width;
		cc->oheight = Height;
		cc->width = 1 << wv_log2i(cc->owidth - 1);
		cc->height = 1 << wv_log2i(cc->oheight - 1);
		cc->size = cc->width * cc->height;
		cc->max_levels = max(0, wv_log2i(min(cc->width, cc->height)) - 2);

		cc->reorder_table = (ReorderTable && *ReorderTable) ? *ReorderTable : init_reorder_tables(cc->width, cc->owidth, cc->height, cc->oheight, cc->max_levels, 0);
		if (ReorderTable && *ReorderTable == NULL)
			*ReorderTable = cc->reorder_table;
		
		channel = malloc(cc->size * sizeof *channel);

		cc->reordered_channel = malloc(cc->size * sizeof *cc->reordered_channel);
		memcpy(channel, Channel, cc->size * sizeof *channel); // use Channel as temp-storage
		
		wavelet_transform(channel, cc->width, cc->height, cc->max_levels);

		// reorder now sets the irrelevant wavelet coefficients to 0
		reorder(cc->reordered_channel, channel, cc->reorder_table, cc->width, cc->height, cc->max_levels);

		cc->num_blocks = cc->max_levels + 1;
		cc->size_min_block = (cc->width >> (cc->num_blocks - 1)) * (cc->height >> (cc->num_blocks - 1));

		cc->final_mse = 65536.0f;

		cc->blocks = malloc(cc->num_blocks * sizeof *cc->blocks);
		end = 0;
		for (i = 0; i < cc->num_blocks; i++)
		{ // do this before, so we know how much work we have to do
			t_wv_cblock* bl = &cc->blocks[i];
			int j;
			
			bl->offset = (i == 0) ? 0 : (cc->size_min_block << (2 * (i - 1)));
//			bl->size = (i == 0) ? cc->size_min_block : ((cc->size_min_block << (2 * i)) - bl->offset); // original
			bl->size = (i == 0) ? cc->size_min_block : 3 * cc->reorder_table[i - 1][0]; // do not write zeroes
			bl->max_value = 0;
			for (j = bl->offset; j < bl->offset + bl->size; j++)
				bl->max_value = max(bl->max_value, abs(cc->reordered_channel[j]));

			bl->max_bpp = wv_log2i(bl->max_value);
			bl->max_quant = ((Lossless > 0) || i == 0) ? 0 : bl->max_bpp;
			end += bl->max_quant + 1;
		}

		// build dictionary
		end--;
		cur = 0;
		for (i = 0; i < cc->num_blocks; i++)
		{
			t_wv_cblock* bl = &cc->blocks[i];
			int* bit_stat;
			int j;

			bit_stat = malloc(bl->max_bpp * sizeof *bit_stat); // we don't need more

			if (bl->max_bpp > 0)
				rice_encode_size(cc->reordered_channel + bl->offset, bl->size, 0, bl->max_value, bit_stat);
			
			bl->quant = malloc((bl->max_quant + 1) * sizeof *bl->quant);
			for (j = 0; j <= bl->max_quant; j++)
			{
				t_wv_quantised_block* ebl = &bl->quant[j];

				if (j == bl->max_bpp)
					ebl->compressed_size = 4;
				else
					ebl->compressed_size = bit_stat[bl->max_bpp - j - 1];
				
				ebl->mse = estimate_error((Lossless < 0) ? NULL : channel, cc->width, cc->owidth, cc->height, cc->oheight, cc->num_blocks - i, 1 << j);
				if (j > 0)
					ebl->mse = max(ebl->mse, flarger(cc->blocks[i].quant[j - 1].mse));
				
//				printf("block: %i q: %i mse: %f bits: %i\n", i, j, ebl->mse, ebl->compressed_size);
				if (Progress)
					Progress(cur, end, UserData);
				cur++;
			}
			
			free(bit_stat);
		}
		free(channel);
	}
	return cc;
}


//----------------------------------------------------------------------------
// wv_done_channel	frees a cchannel and all its associated data
// FreeReorderTable	!= 0, if you want to free the ReorderTable associated with
//					the cchannel (if a shared table is used, do this _once_!)
//----------------------------------------------------------------------------
WAVELET_DLL_API void WAVELET_DLL_CC wv_done_channel(t_wv_cchannel* CC, const int FreeReorderTable)
{
	if (CC)
	{
		int i;

		for (i = 0; i < CC->num_blocks; i++)
			free(CC->blocks[i].quant);
		free(CC->blocks);

		if (FreeReorderTable)
			free_reorder_tables(CC->reorder_table, CC->max_levels);
		free(CC->reordered_channel);

		free(CC);
	}
}


// --------------- QUANTIZER ROUTINES ---------------

static float recursive_selector_bits(t_wv_cchannel* cc, const int block_idx, const int bits_left, const float mse, unsigned char *went_down)
{
	int i, need_left;
	t_wv_cblock* bl;
	float best_lmse;
	float min_mse, max_mse;

	if (went_down)
		*went_down = 0;

	if (block_idx >= cc->num_blocks)
		return mse;

	min_mse = max_mse = mse;
	for (i = block_idx; i < cc->num_blocks; i++)
	{
		min_mse += cc->blocks[i].quant[0].mse;
		max_mse += cc->blocks[i].quant[cc->blocks[i].max_quant].mse;
	}
	
	bl = &cc->blocks[block_idx];	
	
	if (min_mse > bl->best_mse)
		return max_mse + 65536.0f; // only try this if it can possibly be better (with infite bits)
	
	best_lmse = max_mse;
	need_left = ((cc->num_blocks - 1) - block_idx) * 4; // need 4 bits for each block to indicate num_bp