coder.c

小波变换算法

/**

<> todo : if (!coder->encodeBand)
	provide a generic stub which calls encodeBandBP several times
	(same for BandZT)
	this is a pain in the ass, and perhaps obsolete; the calls to encodeBand
		should probably all go through encodeImage, which is the master.

***/

#define CHECK_EOF

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <crblib/inc.h>
#include <crblib/arithc.h>
#include <crblib/codeutil.h>

#include "coder.h"
#include "image.h"
#include "subbands.h"

int tune_param = 0;

extern coder coderNone, coderOzero, coderOone, coderSigMap,coderSM2,coderSMo0 ,
	coderBitPlane,coderBP,coderBPsorted,coderBPbin,coderBPB2,coderBPBF, coderVQ, coderFrac, coderO1_1,coderO2, coderO1_SB,
	coderSH_O1,coderSH_O1SB,coderZF,coderBPZT,coderBPZT2,coderBPBFZT,coderNOP;

const num_coders = 24;	// (sizeof(coder_list)/sizeofpointer)-1
const coder * coder_list[] = {
		&coderBP,
		&coderNOP,
		&coderBPsorted,
		&coderBitPlane,
		&coderBPZT,
		&coderBPZT2,
		&coderBPBF,
		&coderBPBFZT,
		&coderBPbin,
		&coderBPB2,
		&coderZF,
		&coderSigMap,
		&coderSM2,
		&coderSMo0,
		&coderO1_1,
		&coderO1_SB,
		&coderO2,
		&coderOone,
		&coderOzero,
		&coderVQ,
		&coderFrac,
		&coderNone,
		&coderSH_O1,
		&coderSH_O1SB,
		NULL
};

#define SAFE_PAD (1<<16)

extern wavelet * newWavelet(const image *im,int levels)
{
wavelet *w;
	if ( (w = new(wavelet)) == NULL ) return NULL;
	w->width = im->width;
	w->height = im->height;
	w->planes = im->planes;
	w->levels = levels;
	w->complen = im->tot_bytes;
	if ( (w->comp = malloc(w->complen)) == NULL ) {
		freeWavelet(w); return NULL;
	}

	w->stopline = -1;
	w->im = im;

return w;
}
extern void freeWavelet(wavelet *w)
{
	if ( w ) {
		freeSubbands(w->subband_root);
		smartfree(w->qi);
		smartfree(w->comp);
		free(w);
	}
}


coder * coder_create_read(wavelet *w)
{
coder *ret;

	if ( (ret = new(coder)) == NULL ) return NULL;
	if ( w->coder_template ) memcpy(ret,w->coder_template,sizeof(coder));

	ret->w = w;

	w->stoplen = min(w->complen,w->stoplen);
	
	if ( (ret->arith = arithInit()) == NULL ){
		coder_destroy(ret); return NULL;
	}
	arithDecodeInit(ret->arith,ret->w->comp);

return ret;
}

coder * coder_create_write(const coder *template,wavelet *w,int stoplen)
{
coder *ret;

	if ( (ret = new(coder)) == NULL ) return NULL;
	if ( template ) memcpy(ret,template,sizeof(coder));

	w->coder_template = template;
 
	ret->w = w;

	w->stoplen = stoplen;

	if ( (ret->arith = arithInit()) == NULL ){
		coder_destroy(ret); return NULL;
	}
	arithEncodeInit(ret->arith,ret->w->comp);

return ret;
}

void coder_flush_read(coder *c)
{
#ifdef CHECK_EOF
int got;
got = arithGet(c->arith,79);
if ( got == 43 ) arithDecode(c->arith,43,44,79);
else errputs("warning : didn't get EOF");
#endif // EOF

arithDecodeDone(c->arith);
}
void coder_flush_write(coder *c)
{
#ifdef CHECK_EOF
arithEncode(c->arith,43,44,79);
#endif //EOF

c->w->complen = arithEncodeDone(c->arith);
}

void coder_destroy(coder *c)
{
	if ( c ) {
		if ( c->arith ) arithFree(c->arith);
		free(c);
	}
}

/************* routines for DPCM coding the top plane 

as of v1.5 this coder is quite bad-ass , even beats things like CALIC on your teeny
images (because context coders take too long to rev up) :

we do something cheezy like all the standard image coders.  Since the LL band
is so small, we don't have time to adapt a context coder.  Also, the different LL
bands of images have widely varying statistics, so we can't just use the plain
order -1 coder.

instead we try to make a local estimate of the mean prediction error (MPE).  We then
code using a static coder :

	0  		+ (1 < MPE)
	01 		+ (MPE < 3*MPE )
	001 	+ (3*MPE < 7*MPE )
	0001	+ (7*MPE < 15 * MPE)
	...

it would make more sense to use a Guassian probability distribution
sent to the arithcoder, with MPE as the standard dev, but the cumprobs
for a gaussian are the error function.  you would think we could
handle that by tabulation, but there are intricacies with coding
(sym/mpe) as a real number vs. an integer. (eg. hard to define the
"next" and "prev" symbols to get the neighboring cumprobs)

------

our sign coder is even cheezier.  We build the N+W sign neighbors context.
if N == W we have confidence, and use a binary adaptive coder on (sign == neighbor)
else we just send sign raw

***************/

#define EST_ERROR_SHIFT			3
#define EST_ERROR_PAD			0
#define EST_ERROR(grad,prev)	(((grad + grad + prev)>>EST_ERROR_SHIFT) + 2 + EST_ERROR_PAD)
#define ERROR_CODER_MULT		2		// tunes to 1 !!! almost an 0.2 b gain (on the LL)

#define INC 5

#define bitModel(bit,P0,PT)		do { PT += INC; if (!(bit)) P0 += INC;  if ( PT > 1000 ) { PT >>= 1; P0 >>= 1; P0++; PT += 2; } } while(0)
#define bitEnc(bit,ari,P0,PT)	do { arithEncBit(ari,P0,PT,bit);	bitModel(bit,P0,PT); } while(0)
#define bitDec(bit,ari,P0,PT)	do { bit = arithDecBit(ari,P0,PT);	bitModel(bit,P0,PT); } while(0)

static int signs_p0,signs_pt;

void coder_encodesign(arithInfo *ari,bool sign,bool W,bool N)
{
	if ( N == W ) {
		if ( ! N ) sign = !sign;
		bitEnc(sign,ari,signs_p0,signs_pt);
	} else {
		arithEncBitRaw(ari,sign);
	}
}
bool coder_decodesign(arithInfo *ari,bool W,bool N)
{
	if ( N == W ) {
		bool sign;
		bitDec(sign,ari,signs_p0,signs_pt);
		return sign ? N : !N;
	} else {
		return arithDecBitRaw(ari);
	}
}

void coder_encodeDPCM(coder *c,int *plane,int width,int height,int rowpad)
{
int x,y,pred,grad,val,sign,prev_err,est,fullw;
arithInfo * ari = c->arith;
int *ptr,*pline;

	fullw = width + rowpad;
	signs_p0 = 10; signs_pt = 20;
	ptr = plane; prev_err = 99;
	for(y=0;y<height;y++) {
		for(x=0;x<width;x++) {
			pline = ptr - fullw;
			if ( y == 0 ) { 
				if ( x == 0 ) { pred = 0; grad = 99; }
				else if ( x == 1 ) { pred = ptr[-1]; grad = 99; }
				else { pred = (ptr[-1] + ptr[-2])>>1; grad = abs(ptr[-1] - ptr[-2]); }
			} else if ( x == 0 ) { pred = (pline[0] + pline[1])>>1; grad = abs(pline[0] - pline[1]); 
			} else if ( x == width-1 ) {
				pred = (ptr[-1] + pline[0])>>1;
				grad = abs(ptr[-1] - pline[0]);
			} else {
				pred = (ptr[-1]*3 + pline[0]*3 + ptr[-1] + pline[1])>>3;
				grad = max( abs(ptr[-1] - ptr[-1]) , abs( pline[0] - pline[1]) );
			}

			val = (*ptr) - pred;

			if ( val < 0 ) { sign = 1; val = -val; }
			else sign = 0;
	
			est = EST_ERROR(grad,prev_err);
			cu_putMulting_ari(val,ari,est,ERROR_CODER_MULT);

			if ( val > 0 ) {
				if ( x == 0 || y == 0 ) {
					arithEncBitRaw(ari,sign);
				} else {
					coder_encodesign(ari,sign,isneg(ptr[-1]),isneg(pline[0]));
				}
			}

			ptr++;
			prev_err = val;
		}
		ptr += rowpad;
	}
}

void coder_decodeDPCM(coder *c,int *plane,int width,int height,int rowpad)
{
int x,y,pred,grad,val,prev_err,est,fullw;
arithInfo * ari = c->arith;
int *ptr,*pline;

	fullw = width + rowpad;
	signs_p0 = 10; signs_pt = 20;
	ptr = plane; prev_err = 99;
	for(y=0;y<height;y++) {
		for(x=0;x<width;x++) {
			pline = ptr - fullw;
			if ( y == 0 ) { 
				if ( x == 0 ) { pred = 0; grad = 99; }
				else if ( x == 1 ) { pred = ptr[-1]; grad = 99; }
				else { pred = (ptr[-1] + ptr[-2])>>1; grad = abs(ptr[-1] - ptr[-2]); }
			} else if ( x == 0 ) { pred = (pline[0] + pline[1])>>1; grad = abs(pline[0] - pline[1]); 
			} else if ( x == width-1 ) {
				pred = (ptr[-1] + pline[0])>>1;
				grad = abs(ptr[-1] - pline[0]);
			} else {
				pred = (ptr[-1]*3 + pline[0]*3 + ptr[-1] + pline[1])>>3;
				grad = max( abs(ptr[-1] - ptr[-1]) , abs( pline[0] - pline[1]) );
			}

			est = EST_ERROR(grad,prev_err);
			val = cu_getMulting_ari(ari,est,ERROR_CODER_MULT);

			prev_err = val;

			if ( val > 0 ) {
				if ( x == 0 || y == 0 ) {
					if ( arithDecBitRaw(ari) ) val = -val;
				} else {
					if ( coder_decodesign(ari,isneg(ptr[-1]),isneg(pline[0])) )
						val = -val;
				}
			}

			*ptr++ = val + pred;
		}
		ptr += rowpad;
	}
}

/*****

the LL DPCM no longer uses the _m1 routines , but some of
the coders use them to do order (-1) coding.

******/

#define M1_THRESH_1 3
#define M1_THRESH_2 8
#define M1_THRESH_3 30
#define M1_THRESH_4 128

void encode_m1(arithInfo * ari,int sym)
{
	if ( sym < M1_THRESH_1 ) {
		arithBit(ari,0);
		arithEncode(ari,sym,sym+1,M1_THRESH_1);
	} else {
		arithBit(ari,1); sym -= M1_THRESH_1;
		if ( sym < M1_THRESH_2 ) {
			arithBit(ari,0);
			arithEncode(ari,sym,sym+1,M1_THRESH_2);
		} else {
			arithBit(ari,1); sym -= M1_THRESH_2;
			if ( sym < M1_THRESH_3 ) {
				arithBit(ari,0);
				arithEncode(ari,sym,sym+1,M1_THRESH_3);
			} else {
				int escape = M1_THRESH_4;

				arithBit(ari,1); sym -= M1_THRESH_3;

				while( sym >= escape ) {
					arithEncode(ari,escape,escape+1,escape+1);
					sym -= escape;	escape += escape;
					if ( escape > ari->safeProbMax ) escape = ari->safeProbMax;
				}
				arithEncode(ari,sym,sym+1,escape+1);
			}
		}
	}
}

int decode_m1(arithInfo *ari)
{
int sym,cur;
sym=0;
	if ( arithGetBit(ari) == 0 ) {
		cur = arithGet(ari,M1_THRESH_1);
		arithDecode(ari,cur,cur+1,M1_THRESH_1);
		sym+=cur;
	} else {
		sym += M1_THRESH_1;
		if ( arithGetBit(ari) == 0 ) {
			cur = arithGet(ari,M1_THRESH_2);
			arithDecode(ari,cur,cur+1,M1_THRESH_2);
			sym+=cur;
		} else {
			sym += M1_THRESH_2;
			if ( arithGetBit(ari) == 0 ) {
				cur = arithGet(ari,M1_THRESH_3);
				arithDecode(ari,cur,cur+1,M1_THRESH_3);
				sym+=cur;
			} else {
				int escape = M1_THRESH_4>>1;

				sym += M1_THRESH_3;
				do {
					escape += escape;
					if ( escape > ari->safeProbMax ) escape = ari->safeProbMax;
					cur = arithGet(ari,escape+1);
					arithDecode(ari,cur,cur+1,escape+1);
					sym += cur;
				} while( cur == escape );
			}
		}
	}
return sym;
}