dpcm.c

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

				predy = n + ((n - nn) + (w - nw) + (ne - nne))/3;	// vector from y

				dx = abs(w - ww) + abs(ne - n) + abs(n - nw);
				dy = abs(n - nn) + abs(w - nw) + abs(ne - nne);
				
#ifdef NEW_PREDICTOR	/*{*/
				/**** the idea of this tree is to not use the vector from a direction
						which has a large gradient.  that is, e.g. for a hard black line on
						a white background, we try to pick up the line ;
					the smooth version of this is actually worse:
						pred = (dy* predx + dx* predy) / (dx + dy);
				 *****/

					   if ( (dy - dx) > GRAD1_DIFF && dx < GRAD1_LO ) {	// horizontal line
					/** almost zero pixels use the high-gradient predictor **/
					pred = predx;
				} else if ( (dx - dy) > GRAD1_DIFF && dy < GRAD1_LO ) { // vertical line
					pred = predy;
				} else if ( (dy - dx) > GRAD2_DIFF && dx < GRAD2_LO ) { // weak horizontal
					pred = (3*predx + predy)>>2;
					LOG(log_num_weak_grad++);
				} else if ( (dx - dy) > GRAD2_DIFF && dy < GRAD2_LO ) { // weak vertical
					pred = (3*predy + predx)>>2;
					LOG(log_num_weak_grad++);
				} else {						// no strong lines (or both strong lines)
					pred = (predx + predy)>>1;
					LOG(log_num_smooth++);
				}

#else	/*}{*/

				/**** ISO CALIC predictors ****/

				if ( dx + dy > 32 ) 
					pred = (dx*predx + dy*predy)/(dx + dy);
				else if ( dy - dx > 12 )
					pred = (2*predx + predy)/3;
				else if ( dx - dy > 12 )
					pred = (2*predy + predx)/3;
				else				
					pred = (predx + predy)>>1;

			/****
				if ( dx + dy > 32 ) 
					pred = (dx*w + dy*n)/(dx + dy);
				else if ( dy - dx > 12 )
					pred = (2*w + n)/3;
				else if ( dx - dy > 12 )
					pred = (2*n + w)/3;
				else				
					pred = ((n + w)>>1);
			****/

#endif	/** gradient adjustor type ****}*/

				if ( y > 2 ) {
					/*** now do corrections for 45/135 gradients: ***/
					// uses nww and nnw	(not loaded)
						
					d45  = abs(w - prevline[x-2]) + abs(ne - nn) + abs(n - prevline2[x-1]);
					d135 = abs(w - n) + abs(nw - nn) + abs(n - nne);

					if ( d45 - d135 > 32 )
						pred += ((ne - nw)>>3);
					else if ( d45 - d135 > 16 )
						pred += ((ne - nw)>>4);
					else if ( d45 - d135 < 32 )
						pred -= ((ne - nw)>>3);
					else if ( d45 - d135 < 16 )
						pred -= ((ne - nw)>>4);
				}
				
				cap_pred(pred);

				/*}{***** make the context ******/

				/* shape context picks up deltas to detect patterns, plus the top bits of
						value to see large differences in patterns correlated to pixel value ****/

#ifdef NEW_SHAPES
				shape_context = 0;
				if ( abs(pred - w) <= 1 ) {
					shape_context <<= 2;
					if ( abs(pred - ww) > 1 ) shape_context ++;
				} else {
					shape_context <<= 1;
					shape_context ++; shape_context <<= 1;
					if ( pred < w ) shape_context ++;
				}
				if ( abs(pred - n) <= 1 ) {
					shape_context <<= 2;
					if ( abs(pred - nn) > 1 )
						shape_context ++;
				} else {
					shape_context <<= 1;
					shape_context ++; shape_context <<= 1;
					if ( pred < n ) shape_context ++;
				}
				if ( abs(pred - nw) <= 1 ) {
					shape_context <<= 2;
					if ( abs(pred - ne) > 1 )
						shape_context ++;
				} else {
					shape_context <<= 1;
					shape_context ++; 
					shape_context <<= 1;
					if ( pred < nw ) shape_context ++;
				}
				shape_context <<= 2;
				shape_context += (pred>>6);

				shape_context_big = shape_context;

				shape_context_big <<= 2;
				pred_diff = abs(w - pred);
					 if ( pred_diff > SHAPE_CONTOUR3 ) shape_context_big += 1;
				else if ( pred_diff > SHAPE_CONTOUR2 ) shape_context_big += 2;
				else if ( pred_diff > SHAPE_CONTOUR1 ) shape_context_big += 3;
				shape_context_big <<= 2;
				pred_diff = abs(n - pred);
					 if ( pred_diff > SHAPE_CONTOUR3 ) shape_context_big += 1;
				else if ( pred_diff > SHAPE_CONTOUR2 ) shape_context_big += 2;
				else if ( pred_diff > SHAPE_CONTOUR1 ) shape_context_big += 3;
				shape_context_big <<= 1;
				if ( predx > pred ) shape_context_big ++;
				shape_context_big <<= 1;				
				if ( predy >= pred ) shape_context_big ++;

#else /** OLD_SHAPES **/

				shape_context_big = 0;
				if ( w  > pred ) shape_context_big += 1<<0;
				if ( n  > pred ) shape_context_big += 1<<1;
				if ( ww > pred ) shape_context_big += 1<<2;
				if ( nn > pred ) shape_context_big += 1<<3;
				if ( nw > pred ) shape_context_big += 1<<4;
				if ( ne > pred ) shape_context_big += 1<<5;
				pred_diff = abs(w - pred);
					 if ( pred_diff > SHAPE_CONTOUR3 ) shape_context_big += 1<<6;
				else if ( pred_diff > SHAPE_CONTOUR2 ) shape_context_big += 2<<6;
				else if ( pred_diff > SHAPE_CONTOUR1 ) shape_context_big += 3<<6;
				pred_diff = abs(n - pred);
					 if ( pred_diff > SHAPE_CONTOUR3 ) shape_context_big += 1<<8;
				else if ( pred_diff > SHAPE_CONTOUR2 ) shape_context_big += 2<<8;
				else if ( pred_diff > SHAPE_CONTOUR1 ) shape_context_big += 3<<8;
				if ( predx > pred ) shape_context_big += 1<<10;
				if ( predy >= pred ) shape_context_big += 1<<11;
				shape_context_big <<= 2;
				shape_context_big += (pred>>6);

				shape_context = 0;
				if ( w  > pred ) shape_context += 1<<0;
				if ( n  > pred ) shape_context += 1<<1;
				if ( ww > pred ) shape_context += 1<<2;
				if ( nn > pred ) shape_context += 1<<3;
				if ( nw > pred ) shape_context += 1<<4;
				if ( ne > pred ) shape_context += 1<<5;
				shape_context <<= 2;
				shape_context += (pred>>6);

#endif /* SHAPES */

				if ( shape_context >= SHAPE_CONTEXTS )
					errputs("shape context too big");
				if ( shape_context_big >= SHAPE_CONTEXTS_BIG )
					errputs("shape context big too big");

				/** code context weighs our confidence in the prediction, so we use a
						statistical spread to fit this confidence **/

				code_context = 0; 

			/*** <> we can do a better job on making code_context *******/

				// part A : 3 bits

				pred_diff = dx + dy;

				     if ( pred_diff < DIFFA1 ) code_context += 0;
				else if ( pred_diff < DIFFA2 ) code_context += 1;
				else if ( pred_diff < DIFFA3 ) code_context += 2;
				else if ( pred_diff < DIFFA4 ) code_context += 3;
				else if ( pred_diff < DIFFA5 ) code_context += 4;
				else if ( pred_diff < DIFFA6 ) code_context += 5;
				else if ( pred_diff < DIFFA7 ) code_context += 6;
				else						   code_context += 7;

				code_context <<= CODE_CONTEXTS_B_SHIFT;

				// part B : 2 bits
				// large deltas previously imply large deltas next
					 if ( delta < DIFFB1 ) code_context += 0;
				else if ( delta < DIFFB2 ) code_context += 1;
				else if ( delta < DIFFB3 ) code_context += 2;
				else					   code_context += 3;

				/*** done with code_context ; now make the "shape_shift" (see below) ***/

				if ( shape_big_seen[shape_context_big] < 4 && ( shape_seen[shape_context] > shape_big_seen[shape_context_big] ) ) {
					shape_shift = shape_delta[shape_context]/shape_seen[shape_context];
					delta = shape_max[shape_context] - shape_min[shape_context];
					if ( delta < 0 ) delta = DELTA_SMALL_FLAG+1;
				} else {
					shape_shift = shape_big_delta[shape_context_big]/
											shape_big_seen[ shape_context_big];
					delta = shape_big_max[shape_context_big] - shape_big_min[shape_context_big];
					if ( delta < 0 ) delta = DELTA_SMALL_FLAG+1;
					LOG(log_used_big_shape++);
				}

				/** flag code_context if this shape_context is a "good" predictor (small range) **/

				code_context <<= 1;
				if ( delta < DELTA_SMALL_FLAG ) {
					code_context++;
					LOG(log_small_delta++);
				}

				/*}{***** (de)code it ***********/

				/*** we shift the delta by shape_delta[] so try to center the deltas on zero,
						so that all source distributions coded with the same code_context are
						at least centered on the same spot ;	(shape_shift is loaded above)
						we also do code_delta *= signof(shape_delta) , as a semi-heuristic fix
				 *****/

				/*** first try the PPM ***/

				imppm_cntx1 = w + (n<<8) + (nw<<16) + (ne<<24);
				imppm_cntx2 = ww + (nn<<8);

				if ( coding ) {
					delta = line[x] - pred;

					if ( ! iEncode(imppminfo,imppm_cntx1,imppm_cntx2,line[x],pred_diff) ) {
						int code_delta;
						code_delta = delta - shape_shift;
						if ( shape_shift < 0 ) code_delta *= -1;

						LOG( if ( abs(code_delta) < abs(delta) ) log_shift_helped++ );

						encode(code_delta + 256,code_context);
					} else {
						LOG(log_used_ppm++);
					}
				} else {
					if ( (value = iDecode(imppminfo,imppm_cntx1,imppm_cntx2,pred_diff)) == -1 ) {
						int code_delta;
						code_delta = decode(code_context) - 256;
						if ( shape_shift < 0 ) code_delta *= -1;
						delta = code_delta + shape_shift;
						line[x] = delta + pred;
				
						iDecodeGot(imppminfo,imppm_cntx1,imppm_cntx2,line[x]);
					} else {
						line[x] = value;
						delta = value - pred;
					}
				}

				shape_delta[shape_context] += delta;
				if ( ++shape_seen[shape_context] == SHAPES_PER_BUCKET ) {
					shape_delta[shape_context] >>= 1;
					shape_seen[ shape_context] >>= 1;
				}
				if ( delta < shape_min[shape_context] ) shape_min[shape_context] = delta;
				if ( delta > shape_max[shape_context] ) shape_max[shape_context] = delta;

				shape_big_delta[shape_context_big] += delta;
				if ( ++shape_big_seen[shape_context_big] == SHAPES_PER_BUCKET ) {
					shape_big_delta[shape_context_big] >>= 1;
					shape_big_seen[ shape_context_big] >>= 1;
				}
				if ( delta < shape_big_min[shape_context_big] ) shape_big_min[shape_context_big] = delta;
				if ( delta > shape_big_max[shape_context_big] ) shape_big_max[shape_context_big] = delta;

				/**}**/
			}
		}
	}
	errputs("");

	/*#*/ }

	stopclock = clock();

	if ( coding ) {
		FastArithEncodeCDone(FAI);
		complen = BitIO_FlushWrite(BII);
		if ( ! FWriteOk(compF,comp,complen) )
			errputs("fwrote short! continuing..");
	} else {
		FastArithDecodeCDone(FAI);

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

	printf("%-12s :",argv[2]);
	TheCompressionIndicator(rawsize*num_planes,complen,stdout);
	printf(", %6ld byps\n",NumPerSec(rawsize*num_planes,stopclock-startclock));

	if ( coding) {
		int totsize;
		totsize = rawsize*num_planes;
		LOG(printf("log_used_big_shape = %f %%\n",(float)log_used_big_shape*100/(float)totsize));
		LOG(printf("log_num_smooth = %f %%\n",(float)log_num_smooth*100/(float)totsize));
		LOG(printf("log_num_weak_grad = %f %%\n",(float)log_num_weak_grad*100/(float)totsize));
		LOG(printf("log_used_ppm = %f %%\n",(float)log_used_ppm*100/(float)totsize));
		LOG(printf("log_small_delta = %f %%\n",(float)log_small_delta*100/(float)totsize));
		LOG(printf("log_shift_helped = %f %%\n",(float)log_shift_helped*100/(float)totsize));
	}

return 0;
}

void encode(long sym,long context)
{

	if ( order1[context] == NULL )
		if ( (order1[context] = contextCreateMax(FAI,512,CONTEXT_SCALE)) == NULL )
			cleanup("context creation failed!");

	if ( ! contextEncode(order1[context],sym) ) {
		ozeroEncode(order0,sym);
	}

return;
}

long decode(long context)
{
long sym;

	if ( order1[context] == NULL )
		if ( (order1[context] = contextCreateMax(FAI,512,CONTEXT_SCALE)) == NULL )
			cleanup("context creation failed!");

	if ( ( sym = contextDecode(order1[context]) ) == -1 ) {
		sym = ozeroDecode(order0);
		contextAdd(order1[context],sym);
	}

return(sym);
}

void ExitFunc(void)
{
int i;

if ( imppminfo ) iFree(imppminfo);

if ( order1 ) {
	for(i=0;i<CODE_CONTEXTS;i++)
		if ( order1[i] ) contextFree(order1[i]);
	free(order1);
}
if ( order0 ) ozeroFree(order0);

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

smartfree(blankline);

if ( raw ) {
	for(i=0;i<num_planes;i++) {
		if ( raw[i] ) {
			raw[i] -= rawpad;
			free(raw[i]);
		}
	}
	free(raw);
}

smartfree(comp);

if ( rawF ) fclose(rawF);
if ( compF) fclose(compF);

}