dct.c

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

	while(nactive>0) {
		for(b=0;b<DCTBLOCK;b++) {
			if ( ! active[b] ) continue;
#ifdef VERBOSE
			fprintf(stderr,"%d\r",b); fflush(stderr);
#endif
			if ( qtable[b] == 1 )
				{ nactive--; active[b] = false; continue; } //** do this better

			inc = qtable[b]>>4; if ( inc < 1 ) inc = 1;

			r_c = code_image(raw,rawout,width,height,num_planes,qtable,1000,0, NULL,NULL,NULL,NULL);
			d_c=0; for(i=0;i<num_planes;i++) { ubyte *rptr,*vptr; rptr = raw[i]; vptr = rawout[i];
			for(j=rawsize;j--;) { d = *rptr++ - *vptr++; 
			d_c += d*d; } }
			d_c = sqrt((float)d_c);

			qtable[b] += inc;
			r_p = code_image(raw,rawout,width,height,num_planes,qtable,1000,0, NULL,NULL,NULL,NULL);
			d_p=0; for(i=0;i<num_planes;i++) { ubyte *rptr,*vptr; rptr = raw[i]; vptr = rawout[i];
			for(j=rawsize;j--;) { d = *rptr++ - *vptr++; 
			d_p += d*d; } }
			d_p = sqrt((float)d_p);

			qtable[b] -= (inc<<1);
			r_m = code_image(raw,rawout,width,height,num_planes,qtable,1000,0, NULL,NULL,NULL,NULL);
			d_m=0; for(i=0;i<num_planes;i++) { ubyte *rptr,*vptr; rptr = raw[i]; vptr = rawout[i];
			for(j=rawsize;j--;) { d = *rptr++ - *vptr++; 
			d_m += d*d; } }
			d_m = sqrt((float)d_m);

			qtable[b] += inc;

			r_p -= r_c;
			d_p -= d_c;
			r_m -= r_c;
			d_m -= d_c;

			if ( d_p <= 0 && r_p < 0  ) { // go plus
				qtable[b] += inc;
				tot_wiggle++;
				num_wiggles++;
#ifdef VERBOSE
				printf("  p: %d/%d , m: %d/%d\n",r_p,d_p,r_m,d_m);
#endif
			} else if ( d_m <= 0 && r_m < 0 ) {
				qtable[b] -= inc;
				tot_wiggle--;
				num_wiggles++;
#ifdef VERBOSE
				printf("  p: %d/%d , m: %d/%d\n",r_p,d_p,r_m,d_m);
#endif
			} else {
				active[b] = false;
				nactive--;
			}
		}
	}
#ifdef VERBOSE
		errputs("");
#endif

#ifdef VERBOSE
		printf("tot_wiggle = %d, num_wiggles = %d\n",tot_wiggle,num_wiggles);
#endif

	} while(0);
#endif // FIDDLE_QUANTIZERS

	RESET_LOG();

	startclock = clock();

	code_image(raw,rawout,width,height,num_planes,qtable,J_sought_s5,max_distortion,
		&comp_dpcm_len,&comp_ac_len,&comp_sign_len,&comp_raw_len);

	stopclock = clock();

	comp_tot_len = comp_sign_len + comp_ac_len + comp_raw_len + comp_dpcm_len;

	printf("%-13s :",FilePart(argv[1]));
	TheCompressionIndicator(totsize,comp_tot_len,stdout);
	printf(", %6ld byps\n",NumPerSec(totsize,stopclock-startclock));

#ifdef VERBOSE
	printf("%-13s :","DC dpcm");TheCompressionIndicator((totsize>>6),comp_dpcm_len,stdout); puts("");
	printf("%-13s :","AC order-0-1");  TheCompressionIndicator(totsize-(totsize>>6),comp_ac_len,stdout); puts("");
	LOG( printf("%-13s :","AC order(-1)");TheCompressionIndicator(log_nm1,comp_raw_len,stdout); puts(""); );
	LOG( printf("%-13s :","signs"); TheCompressionIndicator((log_nsign - log_nsign_dpcm)/8,comp_sign_len,stdout); puts(""); );
#endif

	TheImageAnalyzer(raw,rawout,num_planes,width,height,((float)totsize/comp_tot_len),stdout);

#ifdef DO_LOG
	printf("n0 = %f %% , n1 = %f %%\n",(float)log_n0*100/(float)(totsize),(float)log_n1*100/(float)(totsize));
	printf("zero_len = %d , av = %f , smashed = %d, av gain = %f\n",log_zero_len,(float)log_zero_len/(totsize>>6),log_zero_len_smashed,(float)log_zero_len_gained/log_zero_len_smashed);
	printf("0p0 = %f %% , 0p1 = %f %% , 0pbig = %f %%\n",(float)log_0p0*100/(float)(log_0p0+log_0p1+log_0pbig),(float)log_0p1*100/(float)(log_0p0+log_0p1+log_0pbig),(float)log_0pbig*100/(float)(log_0p0+log_0p1+log_0pbig));
#ifdef CODE_SIGNS
	printf("nsign = %.4f %% , same = %.4f %% , diff = %.4f %% , noc  = %.4f %%\n",(float)log_nsign*100/(float)(totsize),(float)log_nsign_same*100/(float)(log_nsign),(float)log_nsign_diff*100/(float)(log_nsign),(float)log_nsign_noc *100/(float)(log_nsign) );
#endif

	if ( log_zero_len_smashed > 0 )
		printf("max RoD = %.4f , min RoD = %.4f\n",log_max_RoDs5/32.0,log_min_RoDs5/32.0);

	if ( ( rawF = fopen("out.raw","wb") ) == NULL )
		cleanup("open rawout failed");
	for(i=0;i<num_planes;i++) {
		if ( !FWriteOk(rawF,rawout[i],rawsize) )
			errputs("fwrote short! continuing..");
	}
	fclose(rawF); rawF = NULL;
#endif // DO_LOG

return 0;
}

typedef struct { int RoDs5,x,y,len,zlen,D; } smash_data;
int smash_data_cmp(const void *v1,const void *v2)
{
return( (((*(smash_data **)v1))->RoDs5) - (((*(smash_data **)v2))->RoDs5) );
}

int code_image(ubyte **raw,ubyte **rawout,int width,int height,int num_planes,
	int * qtable,int J_sought_s5,int max_distortion,
	int * comp_dpcm_len_ptr,int *comp_ac_len_ptr,int *comp_sign_len_ptr,int *comp_raw_len_ptr)
{
int coded_len,comp_dpcm_len,comp_ac_len,comp_sign_len,comp_raw_len;
int pnum;
bool smash_tail;

if ( J_sought_s5 > 999 || max_distortion < 10 ) smash_tail = false;
else smash_tail = true;

	init_coders();

	dct_image(raw,trans,width,height,num_planes,qtable);

#define TRANS_ZAG(z) (transline[zag_x[z] + x + zag_y[z]])
#define TRANS_ZAG_PREVB(z) (transline[zag_x[z] + x-DCTLINE + zag_y[z]])

#ifdef USE_EOBS
	/* look for zero-tail ; this is a transform, not coding */
if ( smash_tail ) {

	LOG(log_max_RoDs5=0; log_min_RoDs5=99999;);

	for(pnum=0;pnum<num_planes;pnum++) {
		uword *transline,*transplane;
		int x,y,z,val,i;
		int totR,totD;
		int nblocks,b;
		smash_data * sd;
		smash_data ** sdp;

		nblocks = (width*height)/DCTBLOCK;

		if ( (sd = malloc(nblocks*sizeof(smash_data))) == NULL )
			cleanup("malloc smash_data failed!");
		if ( (sdp = malloc(nblocks*sizeofpointer)) == NULL ) {
			free(sd); cleanup("malloc smash_data failed!");
		}
		for(b=0;b<nblocks;b++) {
			sdp[b] = &sd[b];
		}

		transplane = trans[pnum];

		b= -1;

		/** first count R/D for all blocks **/
		for(y=0;y<height;y += DCTLINE) {
			transline = transplane + y*width;
			for(x=0;x<width;x += DCTLINE) {
				b++;

				sd[b].x = x;
				sd[b].y = y;

				z = 63;
				while( z > 0 && TRANS_ZAG(z) == 0 ) {
					z--;
				} //zag(z) is not zero
				sd[b].zlen = (63-z);
				if ( z == 0 ) {
					sd[b].RoDs5 = 0;
					sd[b].len = 63;
					sd[b].D = 0;
					continue;
				}

				totD = 0; totR = 0;
				do {
					int gain,s,R,D;
					s = z;
					val = TRANS_ZAG(s);
					D = (val*val*qtable[zigzag[z]]*qtable[zigzag[z]])>>2;
					s--;
					while( s > 0 && TRANS_ZAG(s) == 0 ) s--;
					gain = (z-s);

					if ( val == 1 )	R = gain + 2;
					else R = gain + intlog2(val+1) + intlog2(intlog2(val+1)) + 1;

					/* J = R/D , J > J_quant ? */

					if ( R > ((J_sought_s5 * intsqrt(D))>>5) ) {
						z = s; // points at nonzero one
						totR += R;
						totD += D;
					} else {
						break;
					}
				} while( z > 0 );
				sd[b].len = (63-z);
				if ( sd[b].len == 0 || totD == 0 ) { 
					sd[b].D = 0; sd[b].RoDs5 = 0;
				} else {
					sd[b].D = totD;
					sd[b].RoDs5 = (totR<<5)/intsqrt(totD);
				}
			}
		}

		/** sort blocks Id's by R/D **/

		qsort(sdp,nblocks,sizeofpointer,smash_data_cmp);

		/** do smashing in order from top R/D **/
		totD = 0;
		for(b=nblocks-1;b>=0;b--) {
			if ( sdp[b]->len == 0 ) continue;

			y = sdp[b]->y; x = sdp[b]->x;
			transline = transplane + y*width;

			if ( totD < max_distortion && sdp[b]->D > 0 ) {
				LOG( if ( sdp[b]->RoDs5 > log_max_RoDs5 ) log_max_RoDs5 = sdp[b]->RoDs5; );
				LOG( if ( sdp[b]->RoDs5 < log_min_RoDs5 ) log_min_RoDs5 = sdp[b]->RoDs5; );

				z = 63;
				for(i=(sdp[b]->len);i--;) {
					LOG( if ( TRANS_ZAG(z) > 0 ) log_zero_len_smashed ++; );
					TRANS_ZAG(z) = TRANS_DONT_CODE;
					z--;
				}
#ifdef DEBUG
				if ( z < 0 || z > 62 ) cleanup("z out of range, you doof");
#endif
				TRANS_ZAG(z+1) = TRANS_EOB;
				totD += sdp[b]->D;
				LOG( log_zero_len += (sdp[b]->len) );
				LOG( log_zero_len_gained += (sdp[b]->len - sdp[b]->zlen); );
			} else if ( (i = (sdp[b]->zlen)) > 0 ) {
				for(z=63;i--;) {
					TRANS_ZAG(z) = TRANS_DONT_CODE;
					z--;
				}
				TRANS_ZAG(z+1) = TRANS_EOB;
				LOG(log_zero_len += (sdp[b]->zlen) );
			}
		}

#ifdef VERBOSE
		printf("max_distortion = %d , totd = %d\n",max_distortion,totD);
#endif

		free(sd);
		free(sdp);
	}
} else { /** don't smash tails **/
	uword *transline,*transplane;
	int x,y,z,val;

	for(pnum=0;pnum<num_planes;pnum++) {
		transplane = trans[pnum];
		for(y=0;y<height;y += DCTLINE) {
			transline = transplane + y*width;
			for(x=0;x<width;x += DCTLINE) {
				if ( TRANS_ZAG(63) > 0 ) continue;
				TRANS_ZAG(63) = TRANS_DONT_CODE;
				z = 62;
				while( z > 0 && TRANS_ZAG(z) == 0 ) {
					TRANS_ZAG(z) = TRANS_DONT_CODE;
					z--;
				}
				TRANS_ZAG(z+1) = TRANS_EOB;

				LOG(log_zero_len += (63-z) );
			}
		}
	}
}
#endif // USE_EOBS

#ifndef NO_CODING // *{*

	arithEncode(FAI,0,1,100); /** send the "quality" factor **/

	/* DPCM the 0ths */
	/*#*/ {
	uword *transline,*transplane,*prevline,sign,context;
	int x,y,pred,prev_sign,val;

	prev_sign = 2;

	for(pnum=0;pnum<num_planes;pnum++) {
		transplane = trans[pnum];
		for(y=0;y<height;y += DCTLINE) {
			prevline  = transline;
			transline = transplane + y*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;

				if ( val == 0 ) sign = 2;
				else if ( val < 0 ) { sign = 1; val = -val; }
				else sign = 0;

				/** warning : val can now be up to 2* order0_alphabet ! ;
					encode_m1 is the safest way to go. **/

#ifdef DPCM_ORDER_M1
				encode_m1_custom(val);
#else
#ifdef DPCM_ORDER_0
				if ( ! contextEncode(order0,val) )
					encode_m1_custom(val);
#else
				encode(val,(CODE_CONTEXTS-1));
					// using context=0 here hurts a huge amount!
#endif
#endif //orders

				if ( val > 0 ) encodeSign(sign,prev_sign);

#ifdef CODE_SIGNS
				prev_sign = sign;
#endif
			}
		}
	}
	/*#*/ }

	comp_dpcm_len = BitIO_GetPos(BII) + BitIO_GetPos(rawBII) + BitIO_GetPos(signBII);
	BitIO_ResetArray(BII,comp); BitIO_ResetArray(rawBII,rawcomp); BitIO_ResetArray(signBII,signcomp);
	LOG( log_nsign_dpcm = log_nsign;);

	/* encode the frequencies (AC) */
	/*#*/ {
	uword *transline,*transplane,val,sign,context;
	int prev_val,prev_sign,prev;
	int x,y,pnum,i,t;

	for(i=1;i<DCTBLOCK;i++) {
		for(pnum=0;pnum<num_planes;pnum++) {
			transplane = trans[pnum];
			for(y=0;y<height;y += DCTLINE) {
				transline = transplane + y*width;
				for(x=0;x<width;x += DCTLINE) {

					/** we scan through all spacial locations of
						a certain frequency (so arithcoders adapt)
						but use the previous zigzag as context **/

					if ( (t = TRANS_ZAG(i)) == TRANS_DONT_CODE )
						continue;

					prev = TRANS_ZAG(i-1);

					sign = t & 1; val = (t>>1) + sign;
					if ( val == 0 || t == TRANS_EOB ) sign = 2;

					prev_sign = prev & 1; prev_val = (prev>>1) + prev_sign;
					if ( prev_val == 0 ) prev_sign = 2;

				LOG( if ( val == 0 ) log_n0 ++; else if ( val == 1 ) log_n1 ++; );
				LOG( if ( prev_val == 0 ) { if (val == 0) log_0p0 ++; else if (val==1) log_0p1++; else log_0pbig++; } );

					for(context=0;;context++) //** could use direct lookup table
						if ( prev_val <= context_buckets[context] )	
							break;

					encode(val,context);
					if ( sign != 2 ) encodeSign(sign,prev_sign);
				}
			}
		}
	}
	/*#*/ }

#endif // *}* ALL CODING

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

	coded_len = done_coders() + comp_dpcm_len;

	comp_ac_len = BitIO_GetPos(BII) - 2;
	comp_sign_len = BitIO_GetPos(signBII) - 2;
	comp_raw_len = BitIO_GetPos(rawBII) - 2;

	if ( comp_dpcm_len_ptr ) *comp_dpcm_len_ptr = comp_dpcm_len;
	if ( comp_ac_len_ptr ) *comp_ac_len_ptr = comp_ac_len;
	if ( comp_sign_len_ptr ) *comp_sign_len_ptr = comp_sign_len;
	if ( comp_raw_len_ptr ) *comp_raw_len_ptr = comp_raw_len;

	free_coders();

return coded_len;
}

void encode(int sym,int context)
{