dpcm.c

语音图像压缩dpcm编码源程序

#define LOG(x)	if (0) ; else x;
//#define LOG(x)	if (1) ; else x;

#define NEW_PREDICTOR
//#define NEW_SHAPES	

/*******

NEW_SHAPES hurts or helps about 0.01 depending on image!!! : this indicates our
choice of contexts & escaping is still far from perfect

4-1-98

beart.256  :   65536 ->   29762 =  3.633 bpb =  2.202 to 1 ,  57996 byps
checa.256  :   65536 ->   25013 =  3.053 bpb =  2.620 to 1 ,  61248 byps
adjani.256 :   65536 ->   34964 =  4.268 bpb =  1.874 to 1 ,  56013 byps
fenn.256   :   65536 ->   27378 =  3.342 bpb =  2.393 to 1 ,  57487 byps
jen.256    :   65536 ->   22434 =  2.738 bpb =  2.921 to 1 ,  60124 byps

*******/

/*******

did , do better? :

		 use a high order PPM first, and then drop to the current system
			only when a match isn't made

			load PPM contexts like

								6
							 425
							31%

			to make an order-6 and then escape down.  I don't want to directly
			use old PPM code because we only need to code from a context if it
			is good (e.g. deterministic), and we almost surely should never use
			order-1 or order-2 PPM, I believe this guy here will be better.

didn't work/help ? :

	1. 15-bit shape context escape down to 8-bit.  Why not?  The larger seems
			to fill up and get used very quickly (90%!!) 	need to log stats somehow

	2.	keep a min & max delta seen in each shape context, as an
				extra measure of the confidence.
			(it seems we should be able to pick up deterministic contexts this way!)

	3. 45 & 135 edge predictors from ISO proposal

-----------

I get about the same compression as CALIC, but I'm throwing a lot more muscle
at the problem : things are imperfect here

we lose about 0.05 bpc from the poor startup (first 2 lines)

*******/

#include <stdlib.h>
#include <stdio.h>
#include <crbinc/inc.h>
#include <crbinc/fileutil.h>
#include <crbinc/cindcatr.h>
#include <crbinc/timer.h>
#include <crbinc/bbitio.h>
#include <crbinc/arithc.h>
#include <crbinc/context.h>
#include <crbinc/o0coder.h>

#include "imppm.h"

#define cleanup(str) if ( 0 ) ; else { errputs(str); exit(10); }

#define cap_pred(pred)	if ( pred > 0xFF ) pred = 0xFF; else if ( pred < 0 ) pred = 0; else ;

int log_shift_helped = 0;
int log_num_smooth = 0;
int log_num_weak_grad = 0;
int log_used_big_shape = 0;
int log_used_ppm = 0;
int log_small_delta=0;

int tune_param = 0;

#define DELTA_SMALL_FLAG 20

		// gradient cuttoffs
#define GRAD1_DIFF	50	// <- doesn't matter much or help much
#define GRAD1_LO	20
#define GRAD2_DIFF	10	// <- hard to tune
#define GRAD2_LO	50

#define SHAPE_CONTOUR1 5	// abs > contour makes a new context
#define SHAPE_CONTOUR2 14	// !!! these makes a strong difference
#define SHAPE_CONTOUR3 80

#define DIFFA1	5
#define DIFFA2	15
#define DIFFA3	25
#define DIFFA4	42
#define DIFFA5	60
#define DIFFA6	85
#define DIFFA7	140

#define DIFFB1	5
#define DIFFB2	10
#define DIFFB3	30

#define CONTEXT_SCALE			3000
#define CONTEXT_SCALE_O0		1024

#define SHAPES_PER_BUCKET		100

#define SHAPE_CONTEXT_SHIFT			8
#define SHAPE_CONTEXTS				(1<<SHAPE_CONTEXT_SHIFT)

#define SHAPE_CONTEXT_BIG_SHIFT		16
#define SHAPE_CONTEXTS_BIG			(1<<SHAPE_CONTEXT_BIG_SHIFT)

#define CODE_CONTEXTS_A_SHIFT		3
#define CODE_CONTEXTS_A	 			(1<<CODE_CONTEXTS_A_SHIFT)
#define CODE_CONTEXTS_B_SHIFT		2
#define CODE_CONTEXTS_B	 			(1<<CODE_CONTEXTS_B_SHIFT)

#define CODE_CONTEXTS (1<<(CODE_CONTEXTS_A_SHIFT + CODE_CONTEXTS_B_SHIFT + 1))

ipi * imppminfo = NULL;
ubyte **raw = NULL;
ubyte *comp = NULL;		//packed deltas
ubyte *blankline = NULL;
FILE *rawF = NULL,*compF = NULL;
struct FAI * FAI = NULL;
struct BitIOInfo * BII = NULL;
bool coding;
ozero *order0 = NULL;
context **order1 = NULL;
int num_planes,rawpad;

int shape_delta[SHAPE_CONTEXTS];
int shape_seen[SHAPE_CONTEXTS];
int shape_min[SHAPE_CONTEXTS];
int shape_max[SHAPE_CONTEXTS];
int shape_big_delta[SHAPE_CONTEXTS_BIG];
int shape_big_seen[SHAPE_CONTEXTS_BIG];
int shape_big_min[SHAPE_CONTEXTS_BIG];
int shape_big_max[SHAPE_CONTEXTS_BIG];

void ExitFunc(void);
void encode(long sym,long context);
long decode(long context);

int main(int argc,char *argv[])
{
char *t1,*t2;
int width,height;
int rawsize,complen,i;
clock_t startclock,stopclock;

	if ( argc < 6 ) {
		errputs("dcpm <e/d> <raw file> <comp file> <width> <height> [planes] [tune param]");
		exit(10);
	}

	if ( atexit(ExitFunc) )
	  cleanup("Couldn't install exit trap!");

	if ( toupper(argv[1][0]) == 'E' || toupper(argv[1][0]) == 'C' )
		coding = TRUE;
	if ( toupper(argv[1][0]) == 'D' || toupper(argv[1][0]) == 'U' )
		coding = FALSE;

	if ( coding ) { t1 = "rb"; t2 = "wb"; }
	else					{ t1 = "wb"; t2 = "rb"; }

	if ( ( rawF = fopen(argv[2],t1) ) == NULL )
		cleanup("open raw failed");

	if ( ( compF = fopen(argv[3],t2) ) == NULL )
		cleanup("open comp failed");

	width = atol(argv[4]);
	height = atol(argv[5]);


	if ( argc > 6 )
		num_planes = atol(argv[6]);
	else
		num_planes = 3;

	if ( argc > 7 )
		tune_param = atol(argv[7]);

	rawsize = width * height;
	rawpad = width * 3;

	if ( coding )
		complen = ((rawsize * num_planes)*9)>>3;
	else
		if ( (complen = FileLengthofFH(compF)) == FileLengthofFH_Error )
			cleanup("FileLen error!");

	if ( (raw = malloc(sizeofpointer*num_planes)) == NULL )
		cleanup("malloc failed");

	for(i=0;i<num_planes;i++) {
		if ( (raw[i] = malloc(rawsize + rawpad)) == NULL )
			cleanup("malloc failed");
		MemClear(raw[i],rawpad);
		raw[i] += rawpad;
	}

	if ( (comp = malloc(complen + 1024)) == NULL )
		cleanup("malloc failed");

	if ( (blankline = malloc(width)) == NULL )
		cleanup("malloc failed");
	MemClear(blankline,width);

	if ( (BII = BitIO_Init(comp)) == NULL )
		cleanup("BitIOInit failed!");

	if ( (FAI = FastArithInit(BII)) == NULL )
		cleanup("FastArithInit failed!");

	for(i=0;i<SHAPE_CONTEXTS;i++) {
		shape_delta[i] = 0;
		shape_min[i] = 999; shape_max[i] = -999;
		shape_seen[i] = 1;
	}

	for(i=0;i<SHAPE_CONTEXTS_BIG;i++) {
		shape_big_delta[i] = 0;
		shape_big_min[i] = 999; shape_big_max[i] = -999;
		shape_big_seen[i] = 1;
	}

	if ( coding ) {

		for(i=0;i<num_planes;i++) {
			if ( !FReadOk(rawF,raw[i],rawsize) )
				errputs("fread short! continuing..");
		}

		FastArithEncodeCInit(FAI);
	} else {
		if ( ! FRead(compF,comp,complen) )
			errputs("fread short! continuing..");

		BitIO_InitRead(BII);
		FastArithDecodeCInit(FAI);
	}

	if ( (imppminfo = iOpen(FAI)) == NULL )
		cleanup("Image_PPM_Init failed!");

	if ( (order0 = ozeroCreateMax(FAI,512,CONTEXT_SCALE_O0)) == NULL )
		cleanup("Order0_Init failed!");

	if ( (order1 = AllocMem(CODE_CONTEXTS*sizeofpointer,MEMF_CLEAR)) == NULL )
		cleanup("Order1_Init failed!");

	startclock = clock();

	/*** precondition the order0 ***/
	for(i=0;i<12;i++) {
		ozeroAdd(order0,256 + i);
		ozeroAdd(order0,256 - i);
	}
	for(i=0;i<4;i++) {
		ozeroAdd(order0,256 + i);
		ozeroAdd(order0,256 - i);
	}

	/*#*/ {
	ubyte *plane,*line,*prevline,*prevline2;
	int pnum,x,y;
	int delta,pred,pred_diff,value;
	int code_context,shape_context,shape_context_big,shape_shift;
	int w,ww,nw,n,ne,nn,nne;
	int dx,dy,predx,predy,d45,d135;
	ulong imppm_cntx1,imppm_cntx2;

	/** do first two lines of all planes **/
	for(pnum=0;pnum<num_planes;pnum++) {
		plane = raw[pnum];

		// y = 0;
		line = plane;

		// code x=0,1
		if ( coding ) {
			ozeroEncode(order0, line[0] + 128);
			ozeroEncode(order0, line[1] - line[0] + 256);
		} else {
			line[0] = ozeroDecode(order0) - 128;
			line[1] = ozeroDecode(order0) + line[0] - 256;
		}

		for(x=2;x<width;x++) {
			w = line[x-1]; ww = line[x-2];
			pred = w + w - ww;
			cap_pred(pred);
			if ( coding ) {
				ozeroEncode(order0,line[x] - pred + 256);
			} else {
				line[x] = ozeroDecode(order0) - 256 + pred;
			}
		}

		// y = 1;
		line = plane + width;
		prevline = plane;

		// code x=0,1
		if ( coding ) {
			ozeroEncode(order0, line[0] - prevline[0] + 256);
			pred = line[0] + prevline[1] - prevline[0];
			cap_pred(pred);
			ozeroEncode(order0, line[1] - pred + 256);
		} else {
			line[0] = ozeroDecode(order0) + prevline[0] - 256;
			pred = line[0] + prevline[1] - prevline[0];
			cap_pred(pred);
			line[1] = ozeroDecode(order0) + pred - 256;
		}

		for(x=2;x<width;x++) {
			w = line[x-1]; ww = line[x-2];
			n = prevline[x]; nw = prevline[x-1];
			if ( x == (width-1) ) ne = n;
			else ne = prevline[x+1];

			pred = (3*(w + ((w - ww)>>2) + ((ne - nw)>>1)) + n)>>2;	// vector from x
			cap_pred(pred);

			if ( coding ) {
				delta = line[x] - pred;
				ozeroEncode(order0,delta + 256);
			} else {
				delta = ozeroDecode(order0) - 256;
				line[x] = pred + delta;
			}
		}
	}

	for(pnum=0;pnum<num_planes;pnum++) {
		plane = raw[pnum];
		for(y=2;y<height;y++) {
			fprintf(stderr,"%d/%d\r",y,height);
			line = plane + y*width;
			prevline = line - width;
			prevline2 = prevline - width;
			for(x=0;x<width;x++) {

				/*{**** load in the neighbors ******/

				n = prevline[x];
				nn = prevline2[x];

				if ( x == 0 ) {
					w = ww = nw = prevline[0];
				} else if ( x == 1 ) {
					w = ww = line[0];
					nw = prevline[0];
				} else {
					w = line[x-1]; ww = line[x-2];
					nw = prevline[x-1];
				}

				if ( x != (width-1) ) {
					nne = prevline2[x+1];
					ne = prevline[x+1];
				} else {
					nne = nn; ne = n;
				}

				/*}{**** make predictions *******/

				/**** constants in here were tuned on my .256 set **/

				predx = w + ((w - ww)>>2) + ((ne - nw)>>1);	// vector from x