coder_frac.c

好东西呢

#define FIX_ERRORS
#define DO_LOG

#define MAX_SCALE_0
//#define MAX_SCALE_1	
	/** allowing scales hurts about 0.1 bpp,
	*		fixes and escapes are both reduce by about 200 (1400 -> 1200)
	*		but coding the scale is about 1 bit per match even though it's
	*		peaked around scale == 0
	***/

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

a Fractal-Inspired vector quantizer.

	we beat VQ by about 0.05 when quant < 5
	VQ stomps us by about 0.5 (!!!) when quant >= 10
		presumably this is because when quant is large, compressiblity is high,
		and we are sending about 3 bits for the index regardless.

	Many todos are needed before this becomes a competitive coder.

todo :

	1. perhaps our biggest problem is that in the very bottom HH levels, we
		use parents for contexts which are much noisier than we are.  Allowing
		scale down isn't the solution.  Perhaps the best solution is #3
		Another way would be to use the neighboring bands of the same size as
		context, instead of the parent.

	2.	A. skip over redundant vectors in the idx_to_vec search
		B. skip over vectors that differ by a scale factor (if scaling)
			(2B is probably not worth the effort)

	3. search also already-sent vectors in current plane (frac/vq hybrid)
			Perhaps do context-inspired switching between Frac and VQ !?
			(switch to frac after success in VQ, switch to VQ after failure in Frac)

	4. instead of just taking the first vec under error, do a full R-D search (?)

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

#include <stdio.h>
#include <stdlib.h>
#include <crblib/inc.h>
#include <crblib/arithc.h>
#include <crblib/scontext.h>
#include <crblib/intmath.h>

typedef struct {
	int A,B,C,D;
} vector;

extern int tune_param;

#ifdef DO_LOG
static int matches=0,escapes=0,fixes=0,nscale0=0,nscale1=0;
#endif

#ifdef DO_LOG
#define LOG(x)	if(0) ; else { x; }
#else
#define LOG(x)
#endif

#define VEC_ESCAPE			8		/** best at == 0 !! just doing order-0 ! **/

#define VEC_CNTXMAX			3		/** helps **/
#define VEC_CONTEXTS		(VEC_CNTXMAX+1)
#define VEC_CONTEXT(cntx)	min(VEC_CNTXMAX,intlog2(cntx))

#define VEC_TOTMAX		5000
#define VEC_INC			15
#define VEC_ALPHABET	(VEC_ESCAPE+1)

#define LIT_ESCAPE		30
#define LIT_TOTMAX		5000
#define LIT_INC			5
#define LIT_ALPHABET	(LIT_ESCAPE+1)

#include "coder.h"

void coderFrac_encodeBand(coder *me,int *band,int w,int h,int fullw,int *parent);
void coderFrac_decodeBand(coder *me,int *band,int w,int h,int fullw,int *parent);

typedef struct {
	scontext **vec_o1,*lit_o0,*delta_o0,*scale_o0;
	arithInfo * ari;
	int max_err;
} myInfo;

void coderFrac_init(coder *c)
{
myInfo *d;
int i;

	if ( (d = new(myInfo)) == NULL )
		errexit("alloc failed");

	c->data = d;

	d->ari = c->arith;

	d->max_err = 2;
#ifdef FIX_ERRORS
	if ( d->max_err != 2 ) errexit("max err must be 2 to fix error");
#endif

	if ( (d->vec_o1 = newarray(scontext *,VEC_CONTEXTS)) == NULL )
		errexit("alloc failed");

	for(i=0;i<VEC_CONTEXTS;i++) {
		if ( (d->vec_o1[i] = scontextCreate(c->arith,VEC_ALPHABET,0,VEC_TOTMAX,VEC_INC,true)) == NULL )
			errexit("ozero init failed");
	}

	if ( (d->lit_o0 = scontextCreate(c->arith,LIT_ALPHABET,0,LIT_TOTMAX,LIT_INC,true)) == NULL )
		errexit("ozero init failed");

	if ( (d->delta_o0 = scontextCreate(c->arith,3,0,10000,30,true)) == NULL )
		errexit("ozero init failed");

	if ( (d->scale_o0 = scontextCreate(c->arith,3,0,10000,30,true)) == NULL )
		errexit("ozero init failed");
}

void coderFrac_free(coder *c)
{

#ifdef DO_LOG
if ( matches+escapes > 0 ) {
	printf("matches = %d, fixes = %d, escapes = %d\n",matches,fixes,escapes);
	printf("nscale0 = %d,nscale>0 = %d \n",nscale0,nscale1);
	matches = fixes = escapes = 0; nscale0=nscale1=0;
}
#endif

	if ( c->data ) {
		myInfo *d;
		d = c->data;
		if ( d->vec_o1 ) { int i;
			for(i=0;i<VEC_CONTEXTS;i++)
				if ( d->vec_o1[i] ) scontextFree(d->vec_o1[i]);
			free(d->vec_o1);
		}
		if ( d->lit_o0 ) scontextFree(d->lit_o0);
		if ( d->delta_o0 ) scontextFree(d->delta_o0);
		if ( d->scale_o0 ) scontextFree(d->scale_o0);
		free(d);
		c->data = NULL;
	}
}

coder coderFrac = {
		"fractal",
		coderFrac_init,
		coderFrac_free,
		coderFrac_encodeBand,
		coderFrac_decodeBand
	};


static void idx_to_vec_init(int *zplane,int zwidth,int zheight,int zfullw);
static void idx_to_vec(vector *into,int idx,int x,int y);

static int diffVecScaled(vector *x,vector *y,int *scaleptr)	/** x can be scaled down **/
{
int err0,err1,err2;
	err0 =	(x->A - y->A)*(x->A - y->A) + (x->B - y->B)*(x->B - y->B) +
			(x->C - y->C)*(x->C - y->C) + (x->D - y->D)*(x->D - y->D);
#ifdef MAX_SCALE_0
	*scaleptr = 0; return err0;
#endif
	err1 =	((x->A >>1) - y->A)*((x->A >>1) - y->A) + ((x->B >>1) - y->B)*((x->B >>1) - y->B) +
			((x->C >>1) - y->C)*((x->C >>1) - y->C) + ((x->D >>1) - y->D)*((x->D >>1) - y->D);
#ifdef MAX_SCALE_1
	if ( err0 <= err1 ) {
		*scaleptr = 0; return err0;
	} else {
		*scaleptr = 1; return err1;
	}
#endif
	err2 =	((x->A >>2) - y->A)*((x->A >>2) - y->A) + ((x->B >>2) - y->B)*((x->B >>2) - y->B) +
			((x->C >>2) - y->C)*((x->C >>2) - y->C) + ((x->D >>2) - y->D)*((x->D >>2) - y->D);
	if ( err0 <= err1 && err0 <= err2 ) {
		*scaleptr = 0; return err0;
	} else if ( err1 <= err2 ) {
		*scaleptr = 1; return err1;
	} else {
		*scaleptr = 2; return err2;
	}
}

static void encode_esc(scontext *o0,arithInfo *ari,int val)
{
if ( val == 0 ) { scontextEncode(o0,0); return; }
else { 
	int v = abs(val);
	if ( v < LIT_ESCAPE ) scontextEncode(o0,v);
	else {
		scontextEncode(o0,LIT_ESCAPE);
		encode_m1(ari,v - LIT_ESCAPE);
	}
	if ( isneg(val) ) arithBit(ari,1);
	else arithBit(ari,0);
	}
}

static int decode_esc(scontext *o0,arithInfo *ari)
{
int val;
	val = scontextDecode(o0);
	if ( val == 0 ) return 0;
	else if ( val == LIT_ESCAPE ) {
		val += decode_m1(ari);
	}
	if ( arithGetBit(ari) ) val = -val;
return val;
}

static void codeDelta(myInfo *mi,vector *delta,int err)
{
int pos,sign,sign2;

	scontextEncode(mi->delta_o0,err);

	switch(err) {
		case 0: 
			break;
		case 1:
			// 8 values : 2 bit pos, 1 bit sign

#if 0
pos = abs(delta->A) + abs(delta->B) + abs(delta->C) + abs(delta->D);
if ( pos > 1 ) errexit("err > 1");
#endif

			if ( delta->A ) { pos=0; sign = delta->A; }
			if ( delta->B ) { pos=1; sign = delta->B; }
			if ( delta->C ) { pos=2; sign = delta->C; }
			if ( delta->D ) { pos=3; sign = delta->D; }

			if ( sign == -1 ) pos += 4;
			arithEncode(mi->ari,pos,pos+1,8);
			break;
		case 2:
			// 6 positions, 4 sign values

			if ( delta->A ) { sign = delta->A; 
				if ( delta->B ) { pos=0; sign2 = delta->B; }
				else if ( delta->C ) { pos=1; sign2 = delta->C; }
				else if ( delta->D ) { pos=2; sign2 = delta->D; }
				else errexit("should not get here");
			} else if ( delta->B ) { sign = delta->B; 
				if ( delta->C ) { pos=3; sign2 = delta->C; }
				else if ( delta->D ) { pos=4; sign2 = delta->D; }
				else errexit("should not get here");
			} else { pos =5; sign = delta->C; sign2 = delta->D; }

			if ( sign == - 1) pos += 6;
			if ( sign2 == - 1) pos += 12;

			arithEncode(mi->ari,pos,pos+1,24);
			break;
	}
}

static void decodeDelta(myInfo * mi,vector *vec) /** add the delta onto vec **/
{
int err,pos,sign,sign2;

err = scontextDecode(mi->delta_o0);

	switch(err) {
		case 0: 
			break;
		case 1:
			pos = arithGet(mi->ari,8); arithDecode(mi->ari,pos,pos+1,8);
			if ( pos&4 ) { pos -=4; sign = -1; }
			else sign = 1;
			switch(pos) {
				case 0: vec->A += sign; break;
				case 1: vec->B += sign; break;
				case 2: vec->C += sign; break;
				case 3: vec->D += sign; break;
			}
			break;
		case 2:
			pos = arithGet(mi->ari,24); arithDecode(mi->ari,pos,pos+1,24);
			if ( pos >= 12) { pos -= 12; sign2 = -1; } else sign2 = 1;
			if ( pos >= 6) { pos -= 6; sign = -1; } else sign = 1;

			switch(pos) {
				case 0: vec->A += sign; vec->B += sign2; break;
				case 1: vec->A += sign; vec->C += sign2; break;
				case 2: vec->A += sign; vec->D += sign2; break;
				case 3: vec->B += sign; vec->C += sign2; break;
				case 4: vec->B += sign; vec->D += sign2; break;
				case 5: vec->C += sign; vec->D += sign2; break;
				default: errexit("should not get here");
			}
			break;
	}
}

void coderFrac_encodeBand(coder *me,int *band,int width,int height,int fullw,int *parent)
{
int x,y,cntx,idx,err,scale;
int best_idx,best_err,best_scale;
int *dp,*pp,*dpn;
vector vec,vs;
myInfo *mi = (myInfo *)(me->data); 
arithInfo *ari = mi->ari;

	idx_to_vec_init(parent,width>>1,height>>1,fullw);

	dp = band;	pp = parent;
	for(y=0;y<height;y+=2) {
		dpn = dp + fullw;
		if ( coder_timetostop(me) ) { coder_didstop(me,y); return; }
		for(x=0;x<width;x+=2) {
			cntx = abs(pp[x>>1]);
			cntx = VEC_CONTEXT(cntx);

			vec.A = dp[x];	vec.B = dp[x+1];
			vec.C = dpn[x];	vec.D = dpn[x+1];

			best_idx = VEC_ESCAPE; best_err = mi->max_err + 1;
			for(idx=0;idx<VEC_ESCAPE;idx++) {
				idx_to_vec(&vs,idx,x>>1,y>>1);
				err = diffVecScaled(&vs,&vec,&scale);
				if ( err < best_err ) {
					best_err = err; best_idx = idx; best_scale = scale;
					if ( err == 0 ) break;
				}
			}
				
			if ( best_err <= mi->max_err ) {
				idx_to_vec(&vs,best_idx,x>>1,y>>1);
				LOG(matches++);
				LOG( if (best_scale ==0) nscale0++; else nscale1++; );
				scontextEncode(mi->vec_o1[cntx],best_idx);
				scontextEncode(mi->scale_o0,best_scale);
				vs.A >>= best_scale; vs.B >>= best_scale;
				vs.C >>= best_scale; vs.D >>= best_scale;
#ifdef FIX_ERRORS
				if ( best_err > 0 ) LOG(fixes++);
				vec.A -= vs.A;	vec.B -= vs.B;
				vec.C -= vs.C;	vec.D -= vs.D;
				codeDelta(mi,&vec,best_err);
#else // FIX_ERRORS
				if ( best_err != 0 ) {	/** must put in the used values for context coding **/
					dp[x] =  vs.A;	dp[x+1]= vs.B;
					dpn[x] = vs.C;	dpn[x+1]=vs.D;
				}
#endif
			} else {
				LOG(escapes++);
				scontextEncode(mi->vec_o1[cntx],VEC_ESCAPE);
				encode_esc(mi->lit_o0,ari,vec.A); encode_esc(mi->lit_o0,ari,vec.B);
				encode_esc(mi->lit_o0,ari,vec.C); encode_esc(mi->lit_o0,ari,vec.D);
			}
		}
		pp += fullw;
		dp += fullw + fullw;
	}
}

void coderFrac_decodeBand(coder *me,int *band,int width,int height,int fullw,int *parent)
{
int x,y,cntx,idx,scale;
int *dp,*pp,*dpn;
vector vec;
myInfo *mi = (myInfo *)(me->data); 
arithInfo *ari = mi->ari;

	idx_to_vec_init(parent,width>>1,height>>1,fullw);

	dp = band;	pp = parent;
	for(y=0;y<height;y+=2) {
		dpn = dp + fullw;
		if ( coder_timetostopd(me,y) ) return;
		for(x=0;x<width;x+=2) {	/** x & y are the parent's location *2 **/
			cntx = abs(pp[x>>1]);
			cntx = VEC_CONTEXT(cntx);

			idx = scontextDecode(mi->vec_o1[cntx]);

			if ( idx == VEC_ESCAPE ) {
				dp[x]	= decode_esc(mi->lit_o0,ari);
				dp[x+1] = decode_esc(mi->lit_o0,ari);
				dpn[x] 	= decode_esc(mi->lit_o0,ari);
				dpn[x+1]= decode_esc(mi->lit_o0,ari);
			} else {
				idx_to_vec(&vec,idx,x>>1,y>>1);
				scale = scontextDecode(mi->scale_o0);
				vec.A >>= scale; vec.B >>= scale;
				vec.C >>= scale; vec.D >>= scale;
#ifdef FIX_ERRORS
				decodeDelta(mi,&vec);
#endif
				dp[x]  = vec.A;	dp[x+1]  = vec.B;
				dpn[x] = vec.C;	dpn[x+1] = vec.D;
			}
		}
		pp += fullw;
		dp += fullw + fullw;
	}
}

/******** idx_to_vec section : *********/

static int *plane,width,height,fullw;

static void idx_to_vec_init(int *zplane,int zwidth,int zheight,int zfullw)
{
plane = zplane;
width = zwidth;
height = zheight;
fullw = zfullw;
}

static void idx_to_vec(vector *into,int idx,int x,int y)
{
int offx,offy,size;
int *row;

#define STEP 2		/** <> 1 or 2 ? **/

	/** spiral out **/
	size = STEP;
	offx = 0; offy = 0;
	for(;;) {
		while(offx<size) {
			if ( idx == 0 ) goto got_off;
			idx--; offx += STEP;
		}
		while(offy<size) {
			if ( idx == 0 ) goto got_off;
			idx--; offy += STEP;
		}
		while(offx> -size) {
			if ( idx == 0 ) goto got_off;
			idx--; offx -= STEP;
		}
		while(offy> -size) {
			if ( idx == 0 ) goto got_off;
			idx--; offy -= STEP;
		}
		size += STEP;
	}

	got_off:

	while ( offx > width ) offx -= width;
	while ( offx < 0 ) offx += width;
	while ( offy > height ) offy -= height;
	while ( offy < 0 ) offy += height;

x += offx; y += offy;
row = plane + y*fullw;
into->A = row[x]; into->B = row[x+1];
row += fullw;
into->C = row[x]; into->D = row[x+1];

}