ezdct.c

easy discrete time transform implemented in C

	free_coders();

	for(pnum=0;pnum<num_planes;pnum++) {
		MemClear(transout[pnum],rawsize*sizeof(uword));
	}

	init_coders_read();

	/* decode */
	/*#*/ {
	uword *transline,*transplane;
	int x,y,b;
	int context,block,parent,parent_sign;
	int bitpn,bitmask,donemask,nextmask,A,B,C,D;
	int trans_top,top_bitpn;

	top_bitpn = arithGet(FAI,TRANS_MAX_BPN);
	arithDecode(FAI,top_bitpn,top_bitpn+1,TRANS_MAX_BPN);
	trans_top = 1<<top_bitpn;

	nextmask=0; bitpn = top_bitpn;
	for(bitmask = trans_top;bitmask>=TRANS_LOWBIT;bitmask>>=1) {
		donemask = nextmask; bitpn--;
		nextmask = donemask + bitmask;
		for(b=0;b<16;b++) { /** 16 4x4 blocks **/
			for(pnum=0;pnum<num_planes;pnum++) {
				transplane = transout[pnum];
				for(y=0;y<height;y += DCTLINE) {
					transline = transplane + y*width;
					for(x=0;x<width;x += DCTLINE) {

						A = ZAG_INDEX(transline,x, (b<<2)+0 );
						B = ZAG_INDEX(transline,x, (b<<2)+1 );
						C = ZAG_INDEX(transline,x, (b<<2)+2 );
						D = ZAG_INDEX(transline,x, (b<<2)+3 );

						if ( b == 0 ) {
							parent = TRANS_MAX;
						} else {
							parent = ZAG_INDEX(transline,x,b);
						}

						MAKE_CONTEXT(context,(parent&nextmask),((A & donemask)?1:0),((B & donemask)?1:0),((C & donemask)?1:0),((D & donemask)?1:0));
						
						block = decode(context);

						 // only code from parent's sign if it has been sent
						if ( parent & nextmask )	parent_sign = parent & 1;
						else						parent_sign = 2;

						if ( block & 1 ) { A += bitmask; if ( !(A & donemask) ) A += decodeSign(parent_sign); }
						if ( block & 2 ) { B += bitmask; if ( !(B & donemask) ) B += decodeSign(parent_sign); }
						if ( block & 4 ) { C += bitmask; if ( !(C & donemask) ) C += decodeSign(parent_sign); }
						if ( block & 8 ) { D += bitmask; if ( !(D & donemask) ) D += decodeSign(parent_sign); }

						ZAG_INDEX(transline,x, (b<<2)+0 ) = A;
						ZAG_INDEX(transline,x, (b<<2)+1 ) = B;
						ZAG_INDEX(transline,x, (b<<2)+2 ) = C;
						ZAG_INDEX(transline,x, (b<<2)+3 ) = D;
					}

				if ( (BitIO_GetPos(BII)) > req_len )
					goto break_decoding; /** don't waste our time any more **/
				}
			}
		}
#ifdef BITPLANE_ADAPT
		coder_adapt();
#endif
	}

	break_decoding:
		/** push values up to midpoints **/

#ifdef DECODE_MIDPOINTS
		bitmask >>= 1; bitmask--;
		if ( bitmask >= TRANS_LOWBIT ) {
			for(pnum=0;pnum<num_planes;pnum++) { transline = trans[pnum];
				for(x=rawsize;x--;) {
					if ( *transline > 0 ) *transline += bitmask;
					transline++;
				}
			}
		}
#endif // DECODE_MIDPOINTS

#ifdef DPCM_TRANSFORM /** undo it **/
	/*#*/ {
	int val,pred; uword *prevline;
	for(pnum=0;pnum<num_planes;pnum++) {
		transplane = transout[pnum];
		for(y=0;y<height;y += DCTLINE) {
			transline = transplane + y*width;
			prevline  = transline - DCTLINE*width;
			for(x=0;x<width;x += DCTLINE) {
				val = trans_data(transline[x]);

				if ( x == 0 && y == 0 ) {	pred = 0;
				} else if ( y == 0 ) {		pred = trans_data(transline[x-DCTLINE]);
				} else if ( x == 0 ) {		pred = trans_data(prevline[x]);
				} else { int a,b;
					a = trans_data(transline[x-DCTLINE]);
					b = trans_data(prevline[x]);
					pred = (a + b)>>1;
				}

				val += pred;

				transline[x] = data_trans(val);
			}
		}
	}
	/*#*/ }
#endif //DPCM_TRANSFORM

	/*#*/ }

	free_coders();

#endif // NO_CODING /*}*/

	idct_image(rawout,transout,width,height,num_planes,qtable);

return coded_len;
}

void coder_adapt(void)
{
int context;
	for(context=0;context<CODE_CONTEXTS;context++) {
		scontextHalve(order1[context]);
		scontextHalve(order1[context]);
	}
}

void encode(int sym,int context)
{
scontextEncode(order1[context],sym);
}

void encodeSign(bool sign,int parent)
{
#ifdef CODE_SIGNS
	scontextEncode(signorder0[parent],sign);
#else
	arithBit(FAI,sign);
#endif
}

int decode(int context)
{
return scontextDecode(order1[context]);
}
bool decodeSign(int prev)
{
#ifdef CODE_SIGNS
	return scontextDecode(signorder0[prev]);
#else
	return arithGetBit(FAI);
#endif
}

void ExitFunc(void)
{
int i;

free_coders();

smartfree(comp);

if ( raw ) {
	for(i=0;i<num_planes;i++) {
		smartfree(raw[i]);
	}
	free(raw);
}
if ( rawout ) {
	for(i=0;i<num_planes;i++) {
		smartfree(rawout[i]);
	}
	free(rawout);
}
if ( trans ) {
	for(i=0;i<num_planes;i++) {
		smartfree(trans[i]);
	}
	free(trans);
}
if ( transout ) {
	for(i=0;i<num_planes;i++) {
		smartfree(transout[i]);
	}
	free(transout);
}


if ( rawF ) fclose(rawF);

}

int mse_image(ubyte **original,ubyte **vs,int width,int height,int num_planes)
{
int diffs[256];
int diff,i,totsq,pnum;
int rawsize;

	rawsize = width*height;

	MemClear(diffs,256*sizeof(int));

	for(pnum=0;pnum<num_planes;pnum++) {
		ubyte *rptr,*vptr;
		rptr = original[pnum]; vptr = vs[pnum];
		for(i=rawsize;i--;) {
			diff = *rptr++ - *vptr++;
			if ( diff < 0 ) diff = -diff;
			diffs[diff] ++;
		}
	}

	totsq = 0;
	for(i=1;i<256;i++) {
		if ( diffs[i] > 0 ) {
			totsq += i*i*diffs[i];
		}
	}

return totsq;
}

void TheImageAnalyzer(ubyte **original,ubyte **im2,
					int num_planes,int width,int height,
					float ratio,FILE *sio)
{
int diffs[256];
int diff,i,tot,totsq,max,pnum,j;
int rawsize,totsize;
float mse,me,rmse,mse_percep;

rawsize = width*height;
totsize = width*height*num_planes;

	MemClear(diffs,256*sizeof(int));

	for(pnum=0;pnum<num_planes;pnum++) {
		ubyte *rptr,*vptr;
		rptr = original[pnum]; vptr = im2[pnum];
		for(i=rawsize;i--;) {
			diff = *rptr++ - *vptr++;
			if ( diff < 0 ) diff = -diff;
			diffs[diff] ++;
		}
	}

	tot = totsq = max = 0;
	for(i=1;i<256;i++) {
		if ( diffs[i] > 0 ) {
			max = i;
			tot += i * diffs[i];
			totsq += i*i * diffs[i];
		}
	}

	me = (float)tot/totsize;
	mse = (float)totsq/totsize;
	rmse = sqrt(mse);

	for(pnum=0;pnum<num_planes;pnum++) {
		int x,y,av1,av2;
		ulong totds;
		ubyte *line1,*pline1,*nline1;
		ubyte *line2,*pline2,*nline2;
		line1 = original[pnum]; nline1 = line1+width;
		line2 = im2[pnum]; nline2 = line2+width;
		totds = 0;
		for(y=1;y<(width-1);y++) {
			pline1 = line1; line1 = nline1; nline1 += width;
			pline2 = line2; line2 = nline2; nline2 += width;
			for(x=1;x<(width-1);x++) {
				av1 = line1[x] + line1[x-1] + line1[x+1] + pline1[x] + nline1[x];
				av2 = line2[x] + line2[x-1] + line2[x+1] + pline2[x] + nline2[x];
				diff = (av1 - av2); diff *= diff;
				if ( abs( line1[x-1] - line1[x+1] ) < 5 &&
 					 abs( pline1[x] - nline1[x] ) < 5 ) totds += diff + diff;
				else totds += diff;
			}
		}

		mse_percep = (float)totds/(25.0*totsize);
	}

#ifdef VERBOSE
	fprintf(sio,"error: av= %.2f,max= %d,mse= %.3f,rmse= %.2f,psnr= %.2f,perc= %.2f\n",me,max,mse,rmse,(PSNR_MAX - 10*log10(mse)),mse_percep);
	fprintf(sio,"performance: MSE = %f , RMSE = %f , percep = %f\n",(ratio/mse),(ratio/rmse),(ratio/mse_percep));
		/** use MSE in high error regime, RMS in low error **/
#else
	fprintf(sio,"RMSE=%f,PSNR=%f\n",rmse,(PSNR_MAX - 10*log10(mse)));;
#endif

}

void dct_image(ubyte **raw,uword **trans,int width,int height,int num_planes,
		int * quant_table)
{
ubyte *plane,*line,*lptr;
uword *transline,*transplane,*tptr;
RAWDATA rawblock[DCTBLOCK],*rptr;
DCTDATA dctblock[DCTBLOCK],*bptr;
int x,y,pnum,i;

	dct_init(quant_table);

	for(pnum=0;pnum<num_planes;pnum++) {
		plane = raw[pnum];
		transplane = trans[pnum];
		for(y=0;y<height;y += DCTLINE) {
			line = plane + y*width;
			transline = transplane + y*width;
			for(x=0;x<width;x += DCTLINE) {

				rptr = rawblock;
				for(i=0;i<DCTLINE;i++) {
					lptr = line + x + i*width;
					unroll_dctline(*rptr++ = *lptr++);
				}
	
				dct(rawblock,dctblock);

				bptr = dctblock;
				for(i=0;i<DCTLINE;i++) {
					tptr = transline + x + i*width;
					unroll_dctline(*tptr++ = data_trans(*bptr); bptr++);
				}
			}
		}
	}


}


void idct_image(ubyte **rawout,uword **trans,int width,int height,int num_planes,
		int * quant_table)
{
ubyte *plane,*line,*lptr;
uword *transline,*transplane,*tptr;
RAWDATA rawblock[DCTBLOCK],*rptr;
DCTDATA dctblock[DCTBLOCK],*bptr;
int x,y,pnum,i;

	idct_init(quant_table);

	for(pnum=0;pnum<num_planes;pnum++) {
		plane = rawout[pnum];
		transplane = trans[pnum];
		for(y=0;y<height;y += DCTLINE) {
			line = plane + y*width;
			transline = transplane + y*width;
			for(x=0;x<width;x += DCTLINE) {

				bptr = dctblock;
				for(i=0;i<DCTLINE;i++) {
					tptr = transline + x + i*width;
					unroll_dctline(*bptr++ = trans_data(*tptr); tptr++;);
				}

				idct(dctblock,rawblock);

				rptr = rawblock;
				for(i=0;i<DCTLINE;i++) {
					lptr = line + x + i*width;
					unroll_dctline(*lptr++ = *rptr++);
				}
	
			}
		}
	}

}

void init_allocs(void)
{
int i;
	if ( (raw = malloc(sizeofpointer*num_planes)) == NULL )
		cleanup("malloc failed");
	if ( (rawout = malloc(sizeofpointer*num_planes)) == NULL )
		cleanup("malloc failed");
	if ( (trans = malloc(sizeofpointer*num_planes)) == NULL )
		cleanup("malloc failed");
	if ( (transout = malloc(sizeofpointer*num_planes)) == NULL )
		cleanup("malloc failed");

	for(i=0;i<num_planes;i++) {
		if ( (raw[i] = malloc(rawsize)) == NULL )
			cleanup("malloc failed");
		if ( (rawout[i] = malloc(rawsize)) == NULL )
			cleanup("malloc failed");
		if ( (trans[i] = newarray(uword,rawsize)) == NULL )
			cleanup("malloc failed");
		if ( (transout[i] = newarray(uword,rawsize)) == NULL )
			cleanup("malloc failed");
	}

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

void init_coders_write(void)
{
	MemClear(comp,complen);

	init_coders();

	FastArithEncodeCInit(FAI);
}

void init_coders_read(void)
{
	init_coders();

	BitIO_InitRead(BII);

	FastArithDecodeCInit(FAI);
}

void init_coders(void)
{
int i;
	if ( (BII = BitIO_Init(comp)) == NULL )
		cleanup("BitIOInit failed!");
	if ( (FAI = FastArithInit(BII)) == NULL )
		cleanup("FastArithInit failed!");

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

	for(i=0;i<CODE_CONTEXTS;i++) {
		if ( (order1[i] = scontextCreate(FAI,CODE_ALPHABET,0,
				ORDER1_TOTMAX,ORDER1_INC,true)) == NULL )
			cleanup("context creation failed!");
	}

#ifdef CODE_SIGNS
	if ( (signorder0 = AllocMem(SIGN_CONTEXTS*sizeofpointer,MEMF_CLEAR)) == NULL )
		cleanup("Order1_Init failed!");
	for(i=0;i<SIGN_CONTEXTS;i++) {
		if ( (signorder0[i] = scontextCreate(FAI,2,
				0,SIGNORDER0_TOTMAX,SIGNORDER0_INC, true)) == NULL )
			cleanup("Order0_Init failed!");
	}
#endif
}

int done_coders(void)
{
int complen;
	FastArithEncodeCDone(FAI);
	complen = 0;
	complen += BitIO_FlushWrite(BII) - 2;
return complen;
}

void free_coders(void)
{
int i;

if ( order1 ) {
	for(i=0;i<CODE_CONTEXTS;i++)
		if ( order1[i] ) scontextFree(order1[i]);
	free(order1);
	order1 = NULL;
}
#ifdef CODE_SIGNS
if ( signorder0 ) {
	for(i=0;i<SIGN_CONTEXTS;i++)
		if ( signorder0[i] ) scontextFree(signorder0[i]);
	free(signorder0);
	signorder0 = NULL;
}
#endif

if ( FAI ) FastArithCleanUp(FAI); FAI = NULL;
if ( BII ) BitIO_CleanUp(BII); BII = NULL;

}