dpcm.c

傅立叶变换和小波变换是图像压缩的重要工具。该代大戏是利用小波变换进行图像压缩。

/**********

The WaveVideo DPCM coder:

1. code near-lossless like any good dpcm, with a quantizer in the error domain
	(not pixel domain)
2. since the dpcm is very small, smoothness doesn't really work; CALIC and such hurts
3. keep copies of the previous two frames of the subband
		don't ever do deltas, but use the previous frame as part of the predictor
		if (and only if) it's stable over the last two frames
4. coding is done like JPEG-LS : you track an average value and then code golomb using
		a rung for the unary and then some raw bits
5. NOTEZ: the dpcm band is the only one that is quantized to ints without
		the 'packed sign' format; we just use plain signed ints!

todos :

1. should do two different models based on perr
	(a model here is a unary rung & a tot & as count)

**********/

#include "dpcm.h"
#include "mathutil.h"
#include "wavelet.h"
#include <crblib/rungae.h>

//these constants aren't tuned, but seem reasonable:
#define USE_PERR_TOLERANCE				(4)
#define PERR_DONT_USE_DPCM_TOLERANCE	(8)

//#define DO_LOGS

void dpcmCode(SubBand *sb,int levels,int loss,arithInfo *ari,bool forward)
{
int x,y,quantizer,frame;
bool doStore;
int *bptr;
uint tot,count;
int *lptr1,*lptr2;
rung_t rung;
#ifdef DO_LOGS
uint startCompLen;
uint used_pred_dpcm=0,
	 used_pred_last=0,
	 used_pred_both=0;
#endif

	assert( loss > 0 );

	#if 0 //{
	loss = (loss+1)>>1;
	#endif //}

	doStore = ! sb->wavelet->speculative;
	quantizer = (loss*2 + 1);

	rungModelInit(&rung);

	tot = (32 << levels) / quantizer;
	count = 1;

#ifdef DO_LOGS
	startCompLen = arithTellEncPos(ari);
#endif

	frame = sb->wavelet->uses;

	bptr = sb->band;
	SubBand_MakeLastBands(sb);

	for(y = 0;y<sb->height;y++)
	{
		bptr  = ((int *)sb->band) + y * sb->stride;
		lptr1 = ((int *)sb->lastBand1) + y * sb->width;
		lptr2 = ((int *)sb->lastBand2) + y * sb->width;

		for(x=0;x<sb->width;x++)
		{
		int pred,err,val;
		int pred_dpcm,perr_dpcm;
		uint unary,k;

			if ( y == 0 )
			{
				if ( x == 0 )
				{
					if ( frame == 0 )
					{
						if ( forward )
						{
							val = *bptr;
							val = packsigns( val );
							assert( val < 65536 );
							arithEncBits(ari, (val>>8)&0xFF ,8);
							arithEncBits(ari, (val)&0xFF ,8);
							val = *bptr;
						}
						else
						{
							val = arithDecBits(ari,8);
							val = (val<<8) + arithDecBits(ari,8);
							val = unpacksigns( val );
							*bptr = val;
						}
						goto didcode;
					}

					pred_dpcm = lptr1[0];
					perr_dpcm = 64;
				}
				else
				{
					pred_dpcm = bptr[-1];
					perr_dpcm = 32;
				}
			}
			else
			{
			int N;

				N = bptr[- sb->stride];

				if ( x == 0 ) 
				{
					pred_dpcm = N;
					if ( y > 1 )	
					{
					int NN;
						NN = bptr[- 2*sb->stride];
						perr_dpcm =	ABS(N-NN);
					}
					else
					{
					int NE;
						NE = bptr[1 - sb->stride];
						perr_dpcm =	ABS(N-NE) + 1;
					}
				}
				else
				{
				int W;
					W = bptr[-1];
					pred_dpcm = (N+W+1)>>1;
					perr_dpcm = ABS(N-W);
				}
				
			}
		
			if ( frame < 2 )
			{
				pred = pred_dpcm;
			}
			else
			{
			int pred_last,perr_last;

				pred_last = lptr1[x];
				perr_last = ABS( lptr1[x] - lptr2[x] );

				if ( perr_dpcm < USE_PERR_TOLERANCE || perr_dpcm < perr_last )
				{
					pred = pred_dpcm;
					
				#ifdef DO_LOGS
					used_pred_dpcm++;
				#endif
				}
				else if ( perr_last < USE_PERR_TOLERANCE && perr_dpcm > PERR_DONT_USE_DPCM_TOLERANCE )
				{
					pred = pred_last;
				#ifdef DO_LOGS
					used_pred_last++;
				#endif
				}
				else
				{
					pred = (pred_dpcm + pred_last + 1)>>1;
				#ifdef DO_LOGS
					used_pred_both++;
				#endif
				}
			}

			//if ( pred < 0 ) pred = 0; // never predict negative

			// @@ should do two different models based on perr
			//	(a model here is a unary rung & a tot & as count)

			k = ilog2ceil(tot/count);
			if ( k > 14 ) k = 14; // arith max
			
			if ( forward )
			{
				val = *bptr;

				err = val - pred;

				if ( err < 0 )
				{
					err = (err - loss) / quantizer;
				}
				else
				{
					err = (err + loss) / quantizer;
				}

				// get the restored val
				val = pred + err*quantizer;
				
				assert( ABS(val - (*bptr)) <= loss );

				*bptr = val;

				if ( err < 0 )
				{
					arithEncBitRaw(ari, 1);
					err = (-err) - 1;
				}
				else
				{
					arithEncBitRaw(ari, 0);
				}

				assert( err >= 0 && err < 65536 );
			
				tot += err;

				unary = err>>k;

				rungEncodeUnary(ari,unary,&rung);
				arithEncBits(ari, err & ((1UL<<k)-1) ,k);
			}
			else
			{
			bool isNeg;
				isNeg = arithDecBitRaw(ari);

				unary = rungDecodeUnary(ari,&rung);

				err = (unary << k) + arithDecBits(ari,k);

				assert( err >= 0 );

				tot += err;

				if ( isNeg )
				{
					err = - (err + 1);
				}

				val = pred + err*quantizer;
				
				*bptr = val;
			}

			if ( ++count == 32 )
			{
				tot = (tot + 1)>>1;
				count >>= 1;
			}

			didcode :

			if ( doStore )
			{
				lptr2[x] = val;
			}

			bptr++;
		}
	}

	if ( doStore )
	{
		swap( sb->lastBand1, sb->lastBand2 );
	}

#ifdef DO_LOGS //{
	if ( forward && doStore )
	{
	uint compLen;
	float frac;
		frac = 1.0f / (sb->width * sb->height);
		compLen = arithTellEncPos(ari) - startCompLen;
		printf("dpcm coded LL in %d bytes = %1.2f bpp, (%1.1f,%1.1f,%1.1f)\n",
			compLen,(compLen*8.0f*frac),
			used_pred_dpcm*100.0f*frac,
			used_pred_both*100.0f*frac,
			used_pred_last*100.0f*frac);
	}
#endif //}

}