dpcm.sav

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

/**********

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
	// @@ 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;
uint *bptr,*lptr1,*lptr2,tot,count;
rung_t rung;
rung_t signrung;
#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);
	rungModelInit(&signrung);

	for(x=64;x--;) // precondition towards positive; negative is very rare
	{
		rungModelBit(&signrung,0);
	}

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

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

	// there's funny >>1 stuff all over cuz the lowest bit is the sign

	frame = sb->wavelet->uses;

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

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

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

			if ( y == 0 )
			{
				if ( x == 0 )
				{
					if ( frame == 0 )
					{
						if ( forward )
						{
							val = *bptr;
							assert( val < 65536 );
							arithEncBits(ari, (val>>8)&0xFF ,8);
							arithEncBits(ari, (val)&0xFF ,8);
						}
						else
						{
							val = arithDecBits(ari,8);
							val = (val<<8) + arithDecBits(ari,8);
							*bptr = val;
						}
						tot += ((val>>1) + loss) / quantizer;
						count ++;
						goto didcode;
					}

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

				N = bptr[- sb->stride]>>1;

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

				pred_last = lptr1[x];
				perr_last = ABS( (int)lptr1[x] - (int)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
				}
			}

			// @@ 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;
				isNeg = (val&1);
				val >>= 1;

				err = (int)val - (int)pred;

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

				sval = (int)pred + err*quantizer;
				if ( sval < 0 )
				{
					assert( sval > -quantizer );
					sval = 0;
				}
				if ( sval == 0 )
					isNeg = 0;
				val = (sval<<1) + isNeg;
				
				assert( ! sval || (val&1) == (*bptr&1) );
				assert( ABS((int)(val>>1) - (int)(*bptr>>1)) <= 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);
				
				if ( *bptr )
				{
					// the sign is *almost* always positive..
					rungModelEncBit(ari,isNeg,&signrung);
				}
			}
			else
			{
				isNeg = arithDecBitRaw(ari);

				unary = rungDecodeUnary(ari,&rung);

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

				assert( err >= 0 );

				tot += err;

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

				sval = (int)pred + err*quantizer;
				if ( sval < 0 )
				{
					assert( sval > -quantizer );
					sval = 0;
				}
				val = (sval<<1);
				
				if ( val )
				{
					val += rungModelDecBit(ari,&signrung);
				}

				*bptr = val;
			}

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

			didcode :

			if ( doStore )
			{
				lptr2[x] = *bptr >> 1;
			}

			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 //}

}