image.c

基于小波的图像压缩

			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;p<im->planes;p++) {
		fmr = im->data[p];
		tor = subb->data[p];
		for(y=0;y<im->height;y++) {
			a = y>>levels;
			b = y - (a<<levels);
			/** y = (a<<levels) + b **/
			toy = a + ysize*b;
			for(x=0;x<im->width;x++) {
				a = x>>levels;
				b = x - (a<<levels);
				tox = a + xsize*b;
				tor[toy][tox] = fmr[y][x];
			}
		}
	}
}
void deshuffleImageBad(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;p<im->planes;p++) {
		fmr = im->data[p];
		tor = subb->data[p];
		for(y=0;y<im->height;y++) {
			a = y>>levels;
			b = y - (a<<levels);
			/** y = (a<<levels) + b **/
			toy = a + ysize*b;
			for(x=0;x<im->width;x++) {
				a = x>>levels;
				b = x - (a<<levels);
				tox = a + xsize*b;
				fmr[y][x] = tor[toy][tox];
			}
		}
	}
}

void shuffleImage(image *im,image *subb,int levels)
{
int p,bx,by,l,tobx,toby,w,h;
int sw,sh,offx,offy,dsize,i,j,block;
int **dctr,**subr;

	w = im->width;
	h = im->height;
	block = 1<<levels;

	for(p=0;p<im->planes;p++) {
		dctr = im->data[p];
		subr = subb->data[p];
		for(by=0;by<h;by+=block) {
			for(bx=0;bx<w;bx+=block) {
				//l == -1 case :
				subr[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<<l;	//length of the dct subband
					for(i=0;i<dsize;i++) {
						for(j=0;j<dsize;j++) {
							subr[0 + offy + j][sw + offx + i] = dctr[ 0   + by + j][dsize + bx + i];
							subr[sh+ offy + j][sw + offx + i] = dctr[dsize+ by + j][dsize + bx + i];
							subr[sh+ offy + j][ 0 + offx + i] = dctr[dsize+ by + j][0     + bx + i];
						}	
					}
				}
			}
		}
	}
}

void deshuffleImage(image *im,image *subb,int levels)
{
int p,bx,by,l,tobx,toby,w,h,block;
int sw,sh,offx,offy,dsize,i,j;
int **dctr,**subr;

	w = im->width;
	h = im->height;
	block = 1<<levels;

	for(p=0;p<im->planes;p++) {
		dctr = im->data[p];
		subr = subb->data[p];
		for(by=0;by<h;by+=block) {
			for(bx=0;bx<w;bx+=block) {
				//l == -1 case :
				dctr[by][bx] = subr[by>>levels][bx>>levels];
				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<<l;	//length of the dct subband
					for(i=0;i<dsize;i++) {
						for(j=0;j<dsize;j++) {
							dctr[ 0   + by + j][dsize + bx + i] = subr[0 + offy + j][sw + offx + i];
							dctr[dsize+ by + j][dsize + bx + i] = subr[sh+ offy + j][sw + offx + i];
							dctr[dsize+ by + j][0     + bx + i] = subr[sh+ offy + j][ 0 + offx + i];
						}	
					}
				}
			}
		}
	}
}

void shuffleImageFloat(imageFloat **imptr)
{
int p,bx,by,l,tobx,toby,w,h;
int sw,sh,offx,offy,dsize,i,j;
imageFloat *dctim,*subb;
double **dctr,**subr;

	dctim = *imptr;
	if ( (subb = newImageFloatFromFloat(dctim)) == NULL ) return;

	w = dctim->width;
	h = dctim->height;

	for(p=0;p<dctim->planes;p++) {
		dctr = dctim->data[p];
		subr = subb->data[p];
		for(by=0;by<h;by+=8) {
			for(bx=0;bx<w;bx+=8) {
				//l == -1 case :
				subr[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<<l;	//length of the dct subband
					for(i=0;i<dsize;i++) {
						for(j=0;j<dsize;j++) {
							subr[0 + offy + j][sw + offx + i] = dctr[ 0   + by + j][dsize + bx + i];
							subr[sh+ offy + j][sw + offx + i] = dctr[dsize+ by + j][dsize + bx + i];
							subr[sh+ offy + j][ 0 + offx + i] = dctr[dsize+ by + j][0     + bx + i];
						}	
					}
				}
			}
		}
	}

	freeImageFloat(dctim);
	*imptr = subb;
}

void deshuffleImageFloat(imageFloat **imptr)
{
int p,bx,by,l,tobx,toby,w,h;
int sw,sh,offx,offy,dsize,i,j;
imageFloat *dctim,*subb;
double **dctr,**subr;

	subb = *imptr;
	if ( (dctim = newImageFloatFromFloat(subb)) == NULL ) return;

	w = dctim->width;
	h = dctim->height;

	for(p=0;p<dctim->planes;p++) {
		dctr = dctim->data[p];
		subr = subb->data[p];
		for(by=0;by<h;by+=8) {
			for(bx=0;bx<w;bx+=8) {
				//l == -1 case :
				dctr[by][bx] = subr[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<<l;	//length of the dct subband
					for(i=0;i<dsize;i++) {
						for(j=0;j<dsize;j++) {
							dctr[ 0   + by + j][dsize + bx + i] = subr[0 + offy + j][sw + offx + i];
							dctr[dsize+ by + j][dsize + bx + i] = subr[sh+ offy + j][sw + offx + i];
							dctr[dsize+ by + j][0     + bx + i] = subr[sh+ offy + j][ 0 + offx + i];
						}	
					}
				}
			}
		}
	}

	freeImageFloat(subb);
	*imptr = dctim;
}

#define MAX_DIFF 256
#define PSNR_MAX 	(48.165)	//(10*log10(256^2))

double imagePSNR(image *src,image *comp)
{
long diffs[MAX_DIFF+1];
int diff,i,tot,totsq,pnum;
double mse,psnr;

	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] ++;
		}
	}

	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;p<im->planes;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;p<im->planes;p++) {
				ch = im->data[p][0][r] + databyte_mid;
				putminmax(ch,databyte_min,databyte_max);
				fputub(ch,outFile);
			}
		}
	} else {

		for(p=0;p<im->planes;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;p<im->planes;p++) {
				ch = im->data[p][0][r] + databyte_mid;
				putminmax(ch,databyte_min,databyte_max);
				fputuw(ch,outFile);
			}
		}
	} else {
		int *plane;

		for(p=0;p<im->planes;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++;
	}
}