#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 <stdio.h>
#include <math.h>
#include <crblib/inc.h>
#include <crblib/fileutil.h>
#include <crblib/imutil.h>

#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;p<planes;p++) {

		if ( (im->data[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;p<planes;p++) {

		if ( (im->data[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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<fm->planes;p++) {
		fmrows = fm->data[p];
		torows = to->data[p];
		for(y=0;y<h;y++) {
			torow = (torows[y+toy])+tox;
			fmrow = (fmrows[y+fmy])+fmx;
			for(x=w;x--;) {
				*torow++ = *fmrow++;
			}
		}
	}
}

void patchImageFromFloat(imageFloat *fm,image *to)
{
int p,r;
int *rptr; double *fptr;

	for(p=0;p<to->planes;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;p<to->planes;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;p<im->planes;p++) {
		for(l=1;l<=levels;l++) {
			w = (im->width) >> l;
			h = (im->height) >> l;
			rows = &(im->data[p][h]);
			for(y=0;y<h;y++) {
				for(x=y+1;x<w;x++) {
					z			= rows[x][y];
					rows[x][y]	= rows[y][x];
					rows[y][x]	= z;
				}
			}
		}
	}
}

void dwtImage(image *im,imageFloat *imf,int levels,bool inverse)
{
int p;

if ( im->width != imf->width || im->height != imf->height ) return;

for(p=0;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;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;p<im->planes;p++) {
			for(y=0;y<im->height;y+=8) {
				for(x=0;x<im->width;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;p<im->planes;p++) {
			for(y=0;y<im->height;y+=8) {
				for(x=0;x<im->width;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;p<im->planes;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 :