dct.c

关于数字图象处理的ＤＣＴ变换的源代码

#ifdef DEBUG
if ( sym >= ORDER0_ALPHABET ) dbf();
if ( context >= CODE_CONTEXTS ) dbf();
#endif

	if ( order1[context] == NULL )
		if ( (order1[context] = scontextCreate(FAI,ORDER1_ALPHABET,0,
				ORDER1_TOTMAX,ORDER1_INC,true)) == NULL )
			cleanup("context creation failed!");

#ifdef USE_EOBS
	if ( sym == VAL_EOB ) {
		scontextEncode(order1[context],ORDER1_EOB);
	} else 
#endif
		if ( sym < ORDER1_ESCAPE ) {
		scontextEncode(order1[context],sym);
	} else {
		scontextEncode(order1[context],ORDER1_ESCAPE);
		sym -= ORDER1_ESCAPE;

		if ( ! contextEncode(order0,sym) ) {
#ifdef FLAT_ORDER_M1
			encode_m1_flat(sym);
#else
#ifdef JPEG_ORDER_M1
			encode_m1_jpeg(sym);
#else
			encode_m1_custom(sym);
#endif
#endif
			LOG(log_nm1++);
		}
	}

}

void encode_m1_flat(int sym) // order -1
{
	while( sym >= FLAT_ORDER_M1_ESCAPE ) {
		arithEncode(rawFAI,FLAT_ORDER_M1_ESCAPE,FLAT_ORDER_M1_ESCAPE+1,FLAT_ORDER_M1_ESCAPE+1);
		sym -= FLAT_ORDER_M1_ESCAPE;
	}
arithEncode(rawFAI,sym,sym+1,(FLAT_ORDER_M1_ESCAPE+1));
}

void encode_m1_jpeg(int sym) // order -1
{
int topbit,top;

	topbit = 0;
	top = sym;
	while(top)
		{ topbit++; top>>=1; }

	if ( topbit == 0 ) {
		BitIO_WriteZeroBit(rawBII);
		BitIO_WriteZeroBit(rawBII);
	} else {
		top = 1<<(topbit-1);
		sym -= top;
		if ( topbit < 3 ) {
			BitIO_WriteZeroBit(rawBII);
			BitIO_WriteBit(rawBII,1);
			BitIO_WriteBit(rawBII,topbit&1);
		} else {
			BitIO_WriteBit(rawBII,1);
			topbit -= 3;
			if ( topbit >= 7 ) {
				BitIO_WriteBit(rawBII,1);
				BitIO_WriteBit(rawBII,1);
				BitIO_WriteBit(rawBII,1);
				topbit -= 7;
				while(topbit) {	BitIO_WriteBit(rawBII,1); topbit--; }
				BitIO_WriteZeroBit(rawBII);
			} else {
				BitIO_WriteBit(rawBII,(topbit&4));
				BitIO_WriteBit(rawBII,(topbit&2));
				BitIO_WriteBit(rawBII,(topbit&1));
			}
		}
#ifdef DEBUG
		if ( sym >= top ) cleanup(" sym >= top , you blew it ");
#endif
		arithEncode(rawFAI,sym,sym+1,top);
	}
}


void encode_m1_custom(int sym) // order -1
{
	if ( sym < M1_THRESH_1 ) {
		BitIO_WriteZeroBit(rawBII);
		arithEncode(rawFAI,sym,sym+1,M1_THRESH_1);
	} else {
		BitIO_WriteBit(rawBII,1); sym -= M1_THRESH_1;
		if ( sym < M1_THRESH_2 ) {
			BitIO_WriteZeroBit(rawBII);
			arithEncode(rawFAI,sym,sym+1,M1_THRESH_2);
		} else {
			BitIO_WriteBit(rawBII,1); sym -= M1_THRESH_2;
			if ( sym < M1_THRESH_3 ) {
				BitIO_WriteZeroBit(rawBII);
				arithEncode(rawFAI,sym,sym+1,M1_THRESH_3);
			} else {
				BitIO_WriteBit(rawBII,1); sym -= M1_THRESH_3;
				while( sym >= (M1_THRESH_4) ) {
					arithEncode(rawFAI,M1_THRESH_4,M1_THRESH_4+1,M1_THRESH_4+1);
					sym -= (M1_THRESH_4);
				}
				arithEncode(rawFAI,sym,sym+1,M1_THRESH_4+1);
			}
		}
	}
}

void encodeSign(bool sign,int parent)
{
LOG( log_nsign++; );
#ifdef CODE_SIGNS
LOG( if ( parent == 2 ) log_nsign_noc++; else if ( parent == sign ) log_nsign_same++; else log_nsign_diff++; );
	
	scontextEncode(signorder0[parent],sign);
#else
	BitIO_WriteBit(signBII,sign);
#endif
}

void ExitFunc(void)
{
int i;

free_coders();

smartfree(comp);
smartfree(signcomp);
smartfree(rawcomp);

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 ( rawF ) fclose(rawF);

}


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, performanse(rmse) = %f\n",rmse,(PSNR_MAX - 10*log10(mse)),(ratio/rmse));;
#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(int num_planes,int rawsize,int complen)
{
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");

	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] = malloc(sizeof(uword)*rawsize)) == NULL )
			cleanup("malloc failed");
	}

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

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

	if ( (signBII = BitIO_Init(signcomp)) == NULL )
		cleanup("BitIOInit failed!");
	if ( (signFAI = FastArithInit(signBII)) == NULL )
		cleanup("FastArithInit failed!");

	if ( (rawBII = BitIO_Init(rawcomp)) == NULL )
		cleanup("BitIOInit failed!");
	if ( (rawFAI = FastArithInit(rawBII)) == NULL )
		cleanup("FastArithInit failed!");

	FastArithEncodeCInit(FAI);
	FastArithEncodeCInit(signFAI);
	FastArithEncodeCInit(rawFAI);

	if ( (order0 = contextCreateMax(FAI,ORDER0_ALPHABET,ORDER0_TOTMAX)) == NULL )
		cleanup("Order0_Init failed!");

	if ( (order1 = AllocMem(CODE_CONTEXTS*sizeofpointer,MEMF_CLEAR)) == NULL )
		cleanup("Order1_Init 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(signFAI,2,
				0,SIGNORDER0_TOTMAX,SIGNORDER0_INC, true)) == NULL )
			cleanup("Order0_Init failed!");
	}
#endif
}

int done_coders(void)
{
int complen;
	FastArithEncodeCDone(FAI);
	FastArithEncodeCDone(signFAI);
	FastArithEncodeCDone(rawFAI);
	complen = 0;
	complen += BitIO_FlushWrite(BII) - 2;
	complen += BitIO_FlushWrite(signBII) - 2;
	complen += BitIO_FlushWrite(rawBII) - 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 ( order0 ) contextFree(order0); order0 = NULL;

if ( FAI ) FastArithCleanUp(FAI); FAI = NULL;
if ( BII ) BitIO_CleanUp(BII); BII = NULL;
if ( signFAI ) FastArithCleanUp(signFAI); signFAI = NULL;
if ( signBII ) BitIO_CleanUp(signBII); signBII = NULL;
if ( rawFAI ) FastArithCleanUp(rawFAI); rawFAI = NULL;
if ( rawBII ) BitIO_CleanUp(rawBII); rawBII = NULL;

}