/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Port to C from C++
 *
 * Mow-Song, Ng 2/9/2002
 * msng@mmu.edu.my
 * http://www.pesona.mmu.edu.my/~msng
 *
 * I do not claim copyright to the code, but if you use them or modify them,
 * please drop me a mail.
 *
 */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Origianl copyright info */
/*---------------------------------------------------------------------------*/
// Baseline Wavelet Transform Coder Construction Kit
//
// Geoff Davis
// gdavis@cs.dartmouth.edu
// http://www.cs.dartmouth.edu/~gdavis
//
// Copyright 1996 Geoff Davis 9/11/96
//
// Permission is granted to use this software for research purposes as
// long as this notice stays attached to this software.
//
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

#include "wtransform.h"


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
WTRANSFORM *WaveletTransformAlloc(WAVELET *wavelet, FIMAGE *image, int nsteps, 
											 int symmetric)
{
	WTRANSFORM *WTransform;

	if ((WTransform=(WTRANSFORM *)malloc(sizeof(WTRANSFORM)))==NULL){
		return NULL;
	}
	
	WTransform->wavelet=wavelet;
	WTransform->nsteps=nsteps;
	WTransform->symmetric=symmetric;

	WTransform->value=NULL;

	if (image!=NULL){
		WTransform->hsize=image->xsize;
		WTransform->vsize=image->ysize;
		WaveletTransformForward(WTransform, image, wavelet, nsteps, symmetric);
	}
	else{
		WTransform->hsize=WTransform->vsize=0;
	}
	return WTransform;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
WTRANSFORM *WaveletTransformAllocBlank(WAVELET *wavelet, 
								int hsize, int vsize, int nsteps, int symmetric)
{	
	//int i;
	WTRANSFORM *WTransform;

	if ((WTransform=(WTRANSFORM *)malloc(sizeof(WTRANSFORM)))==NULL){
		return NULL;
	}
	
	WTransform->hsize=hsize;
	WTransform->vsize=vsize;
	WTransform->wavelet=wavelet;
	WTransform->nsteps=nsteps;
	WTransform->symmetric=symmetric;

	//- Strange, now we set them back to 0 & -1
	//WTransform->nsteps=0;
	//WTransform->symmetric=-1;
	WaveletTransformInit(WTransform);
	
	//for (i=0; i<hsize*vsize; i++){
	//	WTransform->value[i]=0;
	//}

	return WTransform;
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
WTRANSFORM *WaveletTransformAllocCopy(WTRANSFORM *WTransformSrc)
{
	int i;
	WTRANSFORM *WTransform;

	if ((WTransform=(WTRANSFORM *)malloc(sizeof(WTRANSFORM)))==NULL){
		return NULL;
	}
	
	WTransform->wavelet = WTransformSrc->wavelet;
	WTransform->hsize = WTransformSrc->hsize;
	WTransform->vsize = WTransformSrc->vsize;
	WTransform->nsteps = WTransformSrc->nsteps;
	WTransform->symmetric = WTransformSrc->symmetric;
	
	if (WTransformSrc->value == NULL) {
		WTransform->value = NULL;
	} 
	else {
		WaveletTransformInit(WTransform);
		for (i = 0; i < WTransform->hsize*WTransform->vsize; i++)
			WTransform->value[i] = WTransformSrc->value[i];
	}
	
	return WTransform;

}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletTransformDealloc(WTRANSFORM *WTransform)
{
	WaveletTransformFreeAll(WTransform);
	free(WTransform);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
Real WaveletTransformGetValue(WTRANSFORM *WaveletTransform, int scale, 
										int orientation, int x, int y)
{
	int idx;
	if (scale==0){
		/* special case */
		idx = 0;
	}
	else{
		idx = 3*scale - 2 + orientation;
	}
	
	return (WaveletTransform->subbandPtr[idx]
			[y*WaveletTransform->subbandHSize[idx] + x]);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletTransformSetValue(WTRANSFORM *WaveletTransform, int scale, 
										int orientation, int x, int y, Real val)
{
	int idx;
	if (scale==0){
		/* special case */
		idx = 0;
	}
	else{
		idx = 3*scale - 2 + orientation;
	}
	
	WaveletTransform->subbandPtr[idx]
			[y*WaveletTransform->subbandHSize[idx] + x] = val;
	
	return;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletTransformForward(WTRANSFORM *WTransform, FIMAGE *image, 
									  WAVELET *wavelet, int nsteps, int symmetric)
{
	Real *temp;

	// clear out old info and set up subband pointers
	WaveletTransformFreeAll(WTransform);
	WTransform->hsize = image->xsize;
	WTransform->vsize = image->ysize;
	WTransform->wavelet = wavelet;
	WTransform->nsteps = nsteps;
	WTransform->symmetric = symmetric;
	WaveletTransformInit (WTransform);
	
	temp = (Real *)calloc(WTransform->hsize*WTransform->vsize, sizeof(Real));

	WaveletTransform2D(WTransform->wavelet, image->pixelLinear, 
		temp, WTransform->hsize, WTransform->vsize,
		WTransform->nsteps, WTransform->symmetric);

	// linearize data
	MallatToLinear(WTransform, temp);
	free(temp);

}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletTransformInvert(WTRANSFORM *WTransform, FIMAGE *invertedImage)
{
	Real *temp;
	
	temp = (Real *)calloc(WTransform->hsize*WTransform->vsize, sizeof(Real));
	
	/* put data in Mallat format*/
	/* note that the transform->value is always in linear format */
	LinearToMallat (WTransform, temp);
	
	WaveletInvert2D(WTransform->wavelet, temp, invertedImage->pixelLinear,
		WTransform->hsize, WTransform->vsize, WTransform->nsteps, WTransform->symmetric);
	free(temp);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
FIMAGE *WaveletTransformInvertThis(WTRANSFORM *WTransform)
{
	FIMAGE *image;
	
	image = FImageAlloc(WTransform->hsize, WTransform->vsize);

   WaveletTransformInvert(WTransform, image);

   return image;
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void MallatToLinear(WTRANSFORM *WTransform, Real *mallat)
{
	int i, j, k;
	int *lowHSize, *lowVSize;
	
   /* arrays of top left corner coordinates for subbands */
	lowHSize = (int *) calloc(WTransform->nsteps, sizeof(int));
	lowVSize = (int *) calloc(WTransform->nsteps, sizeof(int));
	
	/* highsed scale size */
	lowHSize[WTransform->nsteps-1] = (WTransform->hsize+1)/2;
	lowVSize[WTransform->nsteps-1] = (WTransform->vsize+1)/2;
	
	for (i = WTransform->nsteps-2; i >= 0; i--) {
		lowHSize[i] = (lowHSize[i+1]+1)/2;
		lowVSize[i] = (lowVSize[i+1]+1)/2;
	}
	
	// move transformed image (in Mallat order) into linear array structure
	// special case for LL subband
	for (j = 0; j < WTransform->subbandVSize[0]; j++){
		for (i = 0; i < WTransform->subbandHSize[0]; i++){
			WTransform->subbandPtr[0][j*WTransform->subbandHSize[0]+i] 
            = mallat[j*WTransform->hsize+i];
		}
	}
		
	for (k = 0; k < WTransform->nsteps; k++) {
		for (j = 0; j < WTransform->subbandVSize[k*3+1]; j++){
			for (i = 0; i < WTransform->subbandHSize[k*3+1]; i++){
					WTransform->subbandPtr[k*3+1][j*WTransform->subbandHSize[k*3+1]+i] = 
					mallat[j*WTransform->hsize+(lowHSize[k]+i)];
			}
		}
		
		for (j = 0; j < WTransform->subbandVSize[k*3+2]; j++){
			for (i = 0; i < WTransform->subbandHSize[k*3+2]; i++){
				WTransform->subbandPtr[k*3+2][j*WTransform->subbandHSize[k*3+2]+i] = 
						mallat[(lowVSize[k]+j)*WTransform->hsize+i];
			}
		}
		
		for (j = 0; j < WTransform->subbandVSize[k*3+3]; j++){
			for (i = 0; i < WTransform->subbandHSize[k*3+3]; i++){
				WTransform->subbandPtr[k*3+3][j*WTransform->subbandHSize[k*3+3]+i] = 
					mallat[(lowVSize[k]+j)*WTransform->hsize+(lowHSize[k]+i)];
			}