www.pudn.com > wavecode.rar > image.c

 
#define AVE_DEL_YUV	// turn U&V to (U+V) , (U-V) 
/** works dandy; it visually works at compactifying magnitudes : U & V are very similar 
* also compactifies noise into the third band 
* <> needs more testing to make sure this isn't image-anomalous 
**/ 
 
#include  
#include  
#include  
#include  
#include  
 
#include "image.h" 
#include "quantim.h" 
 
#include "haar.h" 
#include "dwt.h" 
#include "dct.h" 
#include "spt.h" 
#include "cdf22.h" 
#include "cdf24.h" 
#include "bcw3.h" 
#include "d4.h" 
#include "f97.h" 
#include "b97.h" 
#include "l97.h" 
 
#if 1	// words are Intel 
 
#define fgetuw fgetuwi 
#define fputuw fputuwi 
 
#endif	// intel words 
 
image * newImage(int width, int height,int planes) 
{ 
int p,y; 
image * im; 
int **rows; 
 
	if ( (im = new(image)) == NULL ) return NULL; 
 
	im->width = width; 
	im->height = height; 
	im->planes = planes; 
	im->plane_size = width*height; 
	im->plane_bytes = (im->plane_size)*sizeof(int); 
	im->tot_bytes = im->plane_bytes * planes; 
	im->tot_size = im->plane_size * planes; 
 
	if ( (im->data = newarray(int **,planes)) == NULL ) { 
		freeImage(im); return NULL; 
	} 
 
	for(p=0;pdata[p] = newarray(int *,height+1)) == NULL ) { 
			freeImage(im); return NULL; 
		} 
		rows = im->data[p]; 
 
		if ( (rows[0] = malloc(im->plane_bytes + height)) == NULL ) { 
			freeImage(im); return NULL; 
		} 
 
		for (y = 1; y <= height; y++) 
		    rows[y] = rows[y-1] + width; 
 
  	} 
 
return im; 
} 
 
imageFloat * newImageFloat(int width, int height,int planes) 
{ 
int p,y; 
imageFloat * im; 
double **rows; 
 
	if ( (im = new(image)) == NULL ) return NULL; 
 
	im->width = width; 
	im->height = height; 
	im->planes = planes; 
	im->plane_size = width*height; 
	im->plane_bytes = (im->plane_size)*sizeof(double); 
	im->tot_bytes = im->plane_bytes * planes; 
	im->tot_size = im->plane_size * planes; 
 
	if ( (im->data = newarray(double **,planes)) == NULL ) { 
		freeImageFloat(im); return NULL; 
	} 
 
	for(p=0;pdata[p] = newarray(double *,height+1)) == NULL ) { 
			freeImageFloat(im); return NULL; 
		} 
		rows = im->data[p]; 
 
		if ( (rows[0] = malloc(im->plane_bytes + height)) == NULL ) { 
			freeImageFloat(im); return NULL; 
		} 
 
		for (y = 1; y <= height; y++) 
		    rows[y] = rows[y-1] + width; 
 
  	} 
 
return im; 
} 
 
image * copyImage(image *im) 
{ 
image *new; 
int p; 
 
	if ( (new = newImage(im->width,im->height,im->planes)) == NULL ) 
		return NULL; 
 
	for(p=0;pplanes;p++) { 
		memcpy(new->data[p][0],im->data[p][0],im->plane_bytes); 
	} 
 
return new; 
} 
 
image * newImageFromFloat(imageFloat *im) 
{ 
image *new; 
 
	if ( (new = newImage(im->width,im->height,im->planes)) == NULL ) 
		return NULL; 
 
return new; 
} 
 
imageFloat * newImageFloatFromFloat(imageFloat *im) 
{ 
imageFloat *new; 
 
	if ( (new = newImageFloat(im->width,im->height,im->planes)) == NULL ) 
		return NULL; 
 
return new; 
} 
 
image * newImageFrom(image *im) 
{ 
image *new; 
 
	if ( (new = newImage(im->width,im->height,im->planes)) == NULL ) 
		return NULL; 
 
return new; 
} 
imageFloat * newImageFloatFrom(image *im) 
{ 
imageFloat *new; 
 
	if ( (new = newImageFloat(im->width,im->height,im->planes)) == NULL ) 
		return NULL; 
 
return new; 
} 
 
 
void freeImage(image *im)  
{ 
	if (im ) { 
		if ( im->data ) { 
			int p; 
			for(p=0;pplanes;p++) { 
				if ( im->data[p] ) { 
					if (im->data[p][0] ) 
						free(im->data[p][0]); 
					free(im->data[p]); 
				} 
			} 
			free(im->data); 
		} 
		free(im); 
	} 
} 
 
void freeImageFloat(imageFloat *im)  
{ 
	if (im ) { 
		if ( im->data ) { 
			int p; 
			for(p=0;pplanes;p++) { 
				if ( im->data[p] ) { 
					if (im->data[p][0] ) 
						free(im->data[p][0]); 
					free(im->data[p]); 
				} 
			} 
			free(im->data); 
		} 
		free(im); 
	} 
} 
 
void capBandByte(int *band,int width,int height,int fullw,int cap) 
{ 
int x,y,v,minv,maxv; 
int *dp; 
int rowPad = fullw - width; 
double scale; 
 
	minv = 999; maxv = -999; 
	dp = band; 
	for(y=height;y--;) { 
		for(x=width;x--;) { 
			v = *dp++; 
			if ( v < minv ) minv = v; 
			if ( v > maxv ) maxv = v; 
		} 
		dp += rowPad; 
	} 
 
	maxv = abs(maxv); 
	if ( abs(minv) > abs(maxv) ) maxv = abs(minv); 
	if ( maxv < cap ) return; 
 
	scale = (double)cap/maxv; 
	dp = band; 
	for(y=height;y--;) { 
		for(x=width;x--;) { 
			*dp = (*dp)*scale; dp++; 
		} 
		dp += rowPad; 
	} 
} 
 
void capImageBands(image *im,int levels,int cap) 
{ 
int p,l,sizeX,sizeY; 
int **rows; 
 
	for(p=0;pplanes;p++) { 
		rows = im->data[p]; 
		for (l = levels; l > 0; l--) { 
			sizeX = im->width >> l; 
			sizeY = im->height >> l; 
 
			if (l == levels)  
				capBandByte(rows[0], sizeX, sizeY, im->width,cap); 
 
			capBandByte(rows[0] + sizeX, sizeX, sizeY, im->width,cap); 
			capBandByte(rows[sizeY], sizeX, sizeY, im->width,cap); 
			capBandByte(rows[sizeY]+sizeX, sizeX, sizeY, im->width,cap); 
		} 
	} 
} 
 
void subtractImage(image *im,image *by) 
{ 
int p,r; 
int *ptr,*bptr; 
	for(p=0;pplanes;p++) { 
		ptr = im->data[p][0]; 
		bptr =by->data[p][0]; 
		for(r=im->plane_size;r--;) { 
			*ptr++ -= *bptr++; 
		} 
	} 
} 
 
void addImage(image *im,image *by) 
{ 
int p,r; 
int *ptr,*bptr; 
	for(p=0;pplanes;p++) { 
		ptr = im->data[p][0]; 
		bptr =by->data[p][0]; 
		for(r=im->plane_size;r--;) { 
			*ptr++ += *bptr++; 
		} 
	} 
} 
 
void zeroImage(image *im) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		memclear(im->data[p][0],im->plane_bytes); 
} 
 
void patchImage(image *fm,image *to,int fmx,int fmy,int w,int h,int tox,int toy) 
{ 
int x,y,p; 
int **fmrows,**torows; 
int *fmrow,*torow; 
	for(p=0;pplanes;p++) { 
		fmrows = fm->data[p]; 
		torows = to->data[p]; 
		for(y=0;yplanes;p++) { 
		fptr = fm->data[p][0]; rptr = to->data[p][0]; 
		for(r=to->plane_size;r--;) *rptr++ = *fptr++; 
	} 
} 
 
void patchImageFloatFrom(image *fm,imageFloat *to) 
{ 
int p,r; 
int *rptr; double *fptr; 
 
	for(p=0;pplanes;p++) { 
		fptr = to->data[p][0]; rptr = fm->data[p][0]; 
		for(r=to->plane_size;r--;) *fptr++ = *rptr++; 
	} 
} 
 
void transposeLHs(image *im,int levels) 
{ 
int p,x,y,l,w,h; 
int **rows,z; 
 
	/*** 
	* 
	*	This is a beautiful free compression win! 
	* 
	*	<> the SPT and DWT do a transpose right before this 
	*		instead of two transposes, do zero 
	*	this should be absorbed into each transform in any case 
	* 
	***/ 
 
	assert( im->width == im->height ); 
 
	for(p=0;pplanes;p++) { 
		for(l=1;l<=levels;l++) { 
			w = (im->width) >> l; 
			h = (im->height) >> l; 
			rows = &(im->data[p][h]); 
			for(y=0;ywidth != imf->width || im->height != imf->height ) return; 
 
for(p=0;pplanes;p++) 
	waveletTransform2D(im->data[p],imf->data[p],im->width,im->height,levels,inverse); 
 
} 
 
void cdf22ImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		cdf22_2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void cdf22Image(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	cdf22ImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
 
 
void f97ImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		f97_2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void f97Image(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	f97ImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
void b97ImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		b97_2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void b97Image(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	b97ImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
void l97ImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		l97_2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void l97Image(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	l97ImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
 
void cdf24ImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		cdf24_2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void cdf24Image(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	cdf24ImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
 
void bcw3ImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		bcw3_2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void bcw3Image(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	bcw3ImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
 
void d4ImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
	for(p=0;pplanes;p++) 
		d4_2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void d4Image(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	d4ImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
 
void sptImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
 
/** special-case for the situation where quantizer == 1 
	instead of	:	spt -> float -> quantized ints 
	we do		: 	spt -> ints 
***/ 
 
	for(p=0;pplanes;p++) 
		sp_Transform2D(im->data[p],im->width,im->height,levels,inverse); 
} 
 
void sptImage(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
/** it's unnatural to copy the results of the spt into 
*	a float buffer, but we do it for similarity to the 
*	other transforms 
**/ 
 
	if ( inverse ) patchImageFromFloat(imF,im); 
 
	sptImageInt(im,levels,inverse); 
 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
 
void dctImage(image *im,imageFloat * imf,bool inverse) 
{ 
int p,x,y,i; 
int block[64],*bptr,*line; 
double blockf[64],*fptr,*fline; 
 
	if ( inverse ) { 
		idct_init();	 
		/** imf -> im  **/ 
 
		for(p=0;pplanes;p++) { 
			for(y=0;yheight;y+=8) { 
				for(x=0;xwidth;x+=8) { 
					fptr = blockf; 
					for(i=0;i<8;i++) { 
						fline = imf->data[p][y+i] + x; 
						*fptr++ = *fline++; *fptr++ = *fline++; *fptr++ = *fline++; *fptr++ = *fline++; 
						*fptr++ = *fline++; *fptr++ = *fline++; *fptr++ = *fline++; *fptr++ = *fline++; 
					} 
 
					idct(blockf,block); 
				 
					bptr = block; 
					for(i=0;i<8;i++) { 
						line = im->data[p][y+i] + x; 
						*line++ = *bptr++; *line++ = *bptr++; *line++ = *bptr++; *line++ = *bptr++; 
						*line++ = *bptr++; *line++ = *bptr++; *line++ = *bptr++; *line++ = *bptr++; 
					} 
				} 
			} 
		} 
 
	} else { 
		dct_init();			 
		/** im  -> imf **/ 
 
		for(p=0;pplanes;p++) { 
			for(y=0;yheight;y+=8) { 
				for(x=0;xwidth;x+=8) { 
					bptr = block; 
					for(i=0;i<8;i++) { 
						line = im->data[p][y+i] + x; 
						*bptr++ = *line++; *bptr++ = *line++; *bptr++ = *line++; *bptr++ = *line++; 
						*bptr++ = *line++; *bptr++ = *line++; *bptr++ = *line++; *bptr++ = *line++; 
					} 
 
					dct(block,blockf); 
 
					fptr = blockf; 
					for(i=0;i<8;i++) { 
						fline = imf->data[p][y+i] + x; 
						*fline++ = *fptr++; *fline++ = *fptr++; *fline++ = *fptr++; *fline++ = *fptr++; 
						*fline++ = *fptr++; *fline++ = *fptr++; *fline++ = *fptr++; *fline++ = *fptr++; 
					}				 
				} 
			} 
		} 
	} 
 
} 
 
void haarImageInt(image *im,int levels,bool inverse) 
{ 
int p; 
 
	/** <> 
	*	by far the slowest part of the Haar transform is 
	*	the copy & shuffle.  This could be eliminated entirely 
	*	if we adapted our coders to work on the tree-type data 
	**/ 
 
	if ( inverse ) { 
		image * temp; 
		temp = copyImage(im); 
		deshuffleImage(im,temp,levels); 
		freeImage(temp); 
	} 
 
	for(p=0;pplanes;p++) { 
		haar_Transform2D(im->data[p],im->width,im->height,levels,inverse); 
	} 
 
	if ( !inverse ) { 
		image * temp; 
		temp = copyImage(im); 
		shuffleImage(temp,im,levels); 
		freeImage(temp); 
	} 
 
} 
void haarImage(image *im,imageFloat *imF,int levels,bool inverse) 
{ 
	if ( inverse ) patchImageFromFloat(imF,im); 
	haarImageInt(im,levels,inverse); 
	if (!inverse ) patchImageFloatFrom(im,imF); 
} 
 
void transformImage(image *raw,imageFloat *trans,int levels,bool inverse,int transformN) 
{ 
	switch(transformN) { 
		case TRANS_SPT : 
			sptImage(raw,trans,levels,inverse); 
			break; 
		case TRANS_HAAR : 
			haarImage(raw,trans,levels,inverse); 
			break; 
		case TRANS_DWT : 
			dwtImage(raw,trans,levels,inverse); 
			break; 
		case TRANS_DCT : 
			assert( levels == 3 ); 
			dctImage(raw,trans,inverse); 
			break; 
		case TRANS_CDF22 : 
			cdf22Image(raw,trans,levels,inverse); 
			break; 
		case TRANS_CDF24 : 
			cdf24Image(raw,trans,levels,inverse); 
			break; 
		case TRANS_BCW3 : 
			bcw3Image(raw,trans,levels,inverse); 
			break; 
		case TRANS_D4 : 
			d4Image(raw,trans,levels,inverse); 
			break; 
		case TRANS_F97 : 
			f97Image(raw,trans,levels,inverse); 
			break; 
		case TRANS_B97 : 
			b97Image(raw,trans,levels,inverse); 
			break; 
		case TRANS_L97 : 
			l97Image(raw,trans,levels,inverse); 
			break; 
		default: 
			errexit("tried to do invalid transform number"); 
			break; 
	} 
} 
 
void transformImageInt(image *im,int levels,bool inverse,int transformN) 
{ 
	switch(transformN) { 
		case TRANS_SPT : 
			sptImageInt(im,levels,inverse); 
			break; 
		case TRANS_HAAR : 
			haarImageInt(im,levels,inverse); 
			break; 
		case TRANS_CDF22 : 
			cdf22ImageInt(im,levels,inverse); 
			break; 
		case TRANS_CDF24: 
			cdf24ImageInt(im,levels,inverse); 
			break; 
		case TRANS_BCW3: 
			bcw3ImageInt(im,levels,inverse); 
			break; 
		case TRANS_D4: 
			d4ImageInt(im,levels,inverse); 
			break; 
		case TRANS_F97: 
			f97ImageInt(im,levels,inverse); 
			break; 
		case TRANS_B97: 
			b97ImageInt(im,levels,inverse); 
			break; 
		case TRANS_L97: 
			l97ImageInt(im,levels,inverse); 
			break; 
		default: 
			errexit("tried to do invalid transformInt number"); 
			break; 
	} 
} 
 
void transQuad(int *band,int w,int h,int fullw,bool inverse,int transformN) 
{	 
	switch(transformN) { 
		case TRANS_SPT : 
			spQuad(band,w,h,fullw,inverse); 
			break; 
		case TRANS_CDF22 : 
			cdfQuad(band,w,h,fullw,inverse); 
			break; 
		case TRANS_BCW3 : 
			bcwQuad(band,w,h,fullw,inverse); 
			break; 
		case TRANS_D4 : 
			d4Quad(band,w,h,fullw,inverse); 
			break; 
		case TRANS_F97 : 
			f97Quad(band,w,h,fullw,inverse); 
			break; 
		case TRANS_B97 : 
			b97Quad(band,w,h,fullw,inverse); 
			break; 
		case TRANS_L97 : 
			l97Quad(band,w,h,fullw,inverse); 
			break; 
		default: 
			errexit("tried to do invalid transQuad number"); 
			break; 
	} 
} 
 
/** 
* shuffle takes dct-style locally treed data 
*	turns into subbands 
* this "shuffle" basically scales up each 8x8 block to 
*		take up the full image, then bumps them so they 
*		all lock together. 
*	this is *not* quite the right shuffle to tree-structure dct. 
*		we need to keep plaquettes within a dct-block together spatially  
*	(use the -w option to see what I mean) 
**/ 
 
void shuffleImageBad(image *im,image *subb,int levels) 
{ 
int p,x,y,tox,toy,a,b,xsize,ysize; 
int **fmr,**tor; 
 
	xsize = im->width>>levels; 
	ysize = im->height>>levels; 
	for(p=0;pplanes;p++) { 
		fmr = im->data[p]; 
		tor = subb->data[p]; 
		for(y=0;yheight;y++) { 
			a = y>>levels; 
			b = y - (a<width;x++) { 
				a = x>>levels; 
				b = x - (a<width>>levels; 
	ysize = im->height>>levels; 
	for(p=0;pplanes;p++) { 
		fmr = im->data[p]; 
		tor = subb->data[p]; 
		for(y=0;yheight;y++) { 
			a = y>>levels; 
			b = y - (a<width;x++) { 
				a = x>>levels; 
				b = x - (a<width; 
	h = im->height; 
	block = 1<planes;p++) { 
		dctr = im->data[p]; 
		subr = subb->data[p]; 
		for(by=0;by>levels][bx>>levels] = dctr[by][bx]; 
				for(l=0;l< levels;l++) { 
					sw = w>>(levels-l); 
					sh = h>>(levels-l);	//subband sizes 
					offx = bx>>(levels-l); 
					offy = by>>(levels-l); //offsets in the subband 
					dsize = 1<width; 
	h = im->height; 
	block = 1<planes;p++) { 
		dctr = im->data[p]; 
		subr = subb->data[p]; 
		for(by=0;by>levels][bx>>levels]; 
				for(l=0;l>(levels-l); 
					sh = h>>(levels-l);	//subband sizes 
					offx = bx>>(levels-l); 
					offy = by>>(levels-l); //offsets in the subband 
					dsize = 1<width; 
	h = dctim->height; 
 
	for(p=0;pplanes;p++) { 
		dctr = dctim->data[p]; 
		subr = subb->data[p]; 
		for(by=0;by>3][bx>>3] = dctr[by][bx]; 
				for(l=0;l<=2;l++) { 
					sw = w>>(3-l); 
					sh = h>>(3-l);	//subband sizes 
					offx = bx>>(3-l); 
					offy = by>>(3-l); //offsets in the subband 
					dsize = 1<width; 
	h = dctim->height; 
 
	for(p=0;pplanes;p++) { 
		dctr = dctim->data[p]; 
		subr = subb->data[p]; 
		for(by=0;by>3][bx>>3]; 
				for(l=0;l<=2;l++) { 
					sw = w>>(3-l); 
					sh = h>>(3-l);	//subband sizes 
					offx = bx>>(3-l); 
					offy = by>>(3-l); //offsets in the subband 
					dsize = 1<planes);pnum++) { 
		int *rptr,*vptr; 
		rptr = src->data[pnum][0];  
		vptr = comp->data[pnum][0]; 
		for(i=(src->plane_size);i--;) { 
			diff = *rptr++ - *vptr++; 
			if ( diff < 0 ) diff = -diff; 
			if ( diff > MAX_DIFF ) diff = MAX_DIFF; 
			diffs[diff] ++; 
		} 
	} 
 
	totsq = 0; 
	for(i=1;i<=MAX_DIFF;i++) { 
		if ( diffs[i] > 0 ) { 
			totsq += i*i * diffs[i]; 
		} 
	} 
 
	mse = (float)totsq/(src->tot_size); 
 
	if ( mse < 0.0001 ) psnr = 999.0; 
	else psnr = (PSNR_MAX - 10*log10(mse)); 
 
return psnr; 
} 
 
void imageCompare(image *src,image *comp,FILE *out) 
{ 
long diffs[MAX_DIFF+1]; 
int diff,i,tot,totsq,max,pnum; 
float mse,me,rmse,psnr; 
 
	if ( (src->tot_bytes != comp->tot_bytes) ) { 
		fprintf(out,"Images are not of same size! cannot compare.\n"); 
		return; 
	} 
 
	MemClearLongs(diffs,MAX_DIFF+1); 
 
	for(pnum=0;pnum<(src->planes);pnum++) { 
		int *rptr,*vptr; 
		rptr = src->data[pnum][0];  
		vptr = comp->data[pnum][0]; 
		for(i=(src->plane_size);i--;) { 
			diff = *rptr++ - *vptr++; 
			if ( diff < 0 ) diff = -diff; 
			if ( diff > MAX_DIFF ) diff = MAX_DIFF; 
			diffs[diff] ++; 
		} 
	} 
 
	tot = totsq = max = 0; 
	for(i=1;i<=MAX_DIFF;i++) { 
		if ( diffs[i] > 0 ) { 
			max = i; 
			tot += i * diffs[i]; 
			totsq += i*i * diffs[i]; 
		} 
	} 
 
	if ( tot == 0 ) return; 
 
	me = (float)tot/(src->tot_size); 
	mse = (float)totsq/(src->tot_size); 
	rmse = sqrt(mse); 
 
	if ( mse < 0.0001 ) psnr = 999.0; 
	else psnr = (PSNR_MAX - 10*log10(mse)); 
 
	fprintf(out,"error: av= %.2f,max= %d,mse= %.3f,rmse= %.2f,psnr= %.2f\n",me,max,mse,rmse,psnr); 
} 
 
int imageAverage(image *im) 
{ 
int p,r; 
int *dp; 
int tot=0; 
 
	for(p=0;pplanes;p++) { 
		dp = im->data[p][0]; 
		for(r=im->plane_size;r--;) { 
			tot += *dp++; 
		} 
	} 
 
return ( tot / im->tot_size); 
} 
 
#define databyte_max	0xFF 
#define databyte_min	0 
#define databyte_mid	0x80	// pointless? 
//#define databyte_mid	0 
 
#define dataword_max	0xFFFF 
#define dataword_min	0 
#define dataword_mid	0x8080	// pointless? 
 
bool readImageFileBytes(char *name,image *im,bool interleaved) 
{ 
int p, r; 
FILE *inFile; 
 
	if ((inFile = fopen(name, "rb")) == NULL) 
		return false; 
 
	if ( interleaved && im->planes > 1 ) { 
 
		for(r=0;r<(im->plane_size);r++) { 
			for(p=0;p<(im->planes);p++) { 
				im->data[p][0][r] = (int)(fgetub(inFile) - databyte_mid); 
			} 
		} 
	} else { 
		int *plane; 
 
		for(p=0;p<(im->planes);p++) { 
			plane = im->data[p][0]; 
			for(r=(im->plane_size);r--;) { 
				*plane++ = (int)(fgetub(inFile) - databyte_mid); 
			} 
		} 
	} 
 
	fclose(inFile); 
 
return true; 
} 
 
bool writeImageFileBytes(char *name,image *im,bool interleaved) 
{ 
int p, r, ch; 
int *plane; 
FILE *outFile; 
 
	if ((outFile = fopen(name, "wb")) == NULL) 
		return false; 
 
	if ( interleaved && im->planes > 1 ) { 
		for(r=0;r<(im->plane_size);r++) { 
			for(p=0;pplanes;p++) { 
				ch = im->data[p][0][r] + databyte_mid; 
				putminmax(ch,databyte_min,databyte_max); 
				fputub(ch,outFile); 
			} 
		} 
	} else { 
 
		for(p=0;pplanes;p++) { 
			plane = im->data[p][0]; 
			for(r=im->plane_size;r--;) { 
				ch = *plane++ + databyte_mid; 
				putminmax(ch,databyte_min,databyte_max); 
				fputub(ch,outFile); 
			} 
		} 
	} 
 
	fclose(outFile); 
 
return true; 
} 
 
bool readImageFileWords(char *name,image *im,bool interleaved) 
{ 
int p, r, ch; 
int *plane; 
FILE *inFile; 
 
	if ((inFile = fopen(name, "rb")) == NULL) 
		return false; 
 
	if ( interleaved && im->planes > 1 ) { 
 
		for(r=0;r<(im->plane_size);r++) { 
			for(p=0;p<(im->planes);p++) { 
				im->data[p][0][r] = (int)(fgetuw(inFile) - databyte_mid); 
			} 
		} 
	} else { 
		int *plane; 
 
		for(p=0;p<(im->planes);p++) { 
			plane = im->data[p][0]; 
			for(r=(im->plane_size);r--;) { 
				*plane++ = (int)(fgetuw(inFile) - databyte_mid); 
			} 
		} 
	} 
 
	fclose(inFile); 
 
return true; 
} 
 
bool writeImageFileWords(char *name,image *im,bool interleaved) 
{ 
int p, r, ch; 
FILE *outFile; 
 
	if ((outFile = fopen(name, "wb")) == NULL) 
		return false; 
 
	if ( interleaved && im->planes > 1 ) { 
		for(r=0;r<(im->plane_size);r++) { 
			for(p=0;pplanes;p++) { 
				ch = im->data[p][0][r] + databyte_mid; 
				putminmax(ch,databyte_min,databyte_max); 
				fputuw(ch,outFile); 
			} 
		} 
	} else { 
		int *plane; 
 
		for(p=0;pplanes;p++) { 
			plane = im->data[p][0]; 
			for(r=im->plane_size;r--;) { 
				ch = *plane++ + databyte_mid; 
				putminmax(ch,databyte_min,databyte_max); 
				fputuw(ch,outFile); 
			} 
		} 
	} 
 
	fclose(outFile); 
 
return true; 
} 
 
void RGBtoYUV(image *im) 
{ 
int i; 
int *p0,*p1,*p2; 
int A; 
 
	if ( im->planes < 3 ) return; 
 
	p0 = im->data[0][0]; 
	p1 = im->data[1][0]; 
	p2 = im->data[2][0]; 
	for(i=im->plane_size;i--;) { 
		RGB_to_YUV(*p0,*p1,*p2,p0,p1,p2); 
#ifdef AVE_DEL_YUV 
		A = ((*p1 + *p2)>>1); 
		*p2 -= *p1; 
		*p1 = A; 
#endif 
		p0++; p1++; p2++; 
	} 
} 
 
void YUVtoRGB(image *im) 
{ 
int i; 
int *p0,*p1,*p2; 
 
	if ( im->planes < 3 ) return; 
 
	p0 = im->data[0][0]; 
	p1 = im->data[1][0]; 
	p2 = im->data[2][0]; 
	for(i=im->plane_size;i--;) { 
#ifdef AVE_DEL_YUV 
		*p1 = *p1 - ((*p2 +1)>>1); 
		*p2 += *p1; 
#endif 
		YUV_to_RGB(*p0,*p1,*p2,p0,p1,p2); 
		p0++; p1++; p2++; 
	} 
}