www.pudn.com > wavecode.rar > 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
#include
#include
#include
#include
#include
#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>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>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;
}