/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* 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.
 *
 */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Original 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 "wavelet.h"

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
WAVELET *WaveletAlloc(FILTERSET *filterset)
{
	WAVELET *wavelet;

	if((wavelet=(WAVELET *)malloc(sizeof(WAVELET)))==NULL){
		return NULL;
	}

	wavelet->analysisLow=filterset->analysisLow;
	wavelet->analysisHigh=filterset->analysisHigh;
	wavelet->synthesisLow=filterset->synthesisLow;
	wavelet->synthesisHigh=filterset->synthesisHigh;
	wavelet->symmetric=filterset->symmetric;

	// amount of space to leave for padding vectors for symmetric extensions
	wavelet->npad = max(wavelet->analysisLow->size, wavelet->analysisHigh->size);

	return wavelet;
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletDealloc(WAVELET *wavelet)
{
	free(wavelet);
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletTransform1D(WAVELET *wavelet, Real *input, Real *output, int size,
								int nsteps, int symExt)
{
	int i, currentIndex=0;
	Real *data[2];
	int lowSize=size, highSize;
	
	// If form of extension unspecified, default to symmetric
   // extensions for symmetrical filters and periodic extensions for
   // asymmetrical filters
	
	//- Explicitly set symExt = -1
   if (symExt == -1){
		symExt = wavelet->symmetric;
	}
	
	// data[0] and data[1] are padded with npad entries on each end
   data [0] = (Real *)calloc(2*wavelet->npad+size, sizeof(Real));
   data [1] = (Real *)calloc(2*wavelet->npad+size, sizeof(Real));
	
	
   for (i = 0; i < size; i++){
		data[currentIndex][wavelet->npad+i] = input[i];
	}
	
	while(nsteps--){
		if (lowSize<=2 && wavelet->symmetric == 1){
			WaveletWarning("Reduce # of transform steps or increase signal size.\n");
			WaveletWarning("or switch to periodic extension.\n");
			WaveletError("Low pass subband is too small\n");
		}
		
		// Transform
		WaveletTransformStep(wavelet, data[currentIndex], data[1-currentIndex], 
			lowSize, symExt);
		
		highSize = lowSize/2;
		lowSize = (lowSize+1)/2;
		
		// Copy high-pass data to output signal
		copy (data[1-currentIndex] + wavelet->npad + lowSize, 
			output+lowSize, highSize);
		
		// Now pass low-pass data (first 1/2 of signal) back to
		// transform routine
		currentIndex = 1 - currentIndex;
   }
	
	// Copy low-pass data to output signal
   copy (data[currentIndex] + wavelet->npad, output, lowSize);
   
   free(data[1]);
   free(data[0]);
}


/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/
void WaveletInvert1D(WAVELET *wavelet, Real *input, Real *output, int size,
							int nsteps, int symExt)
{
   int i;
   int currentIndex = 0;
   Real *data[2];
	int *lowSize, *highSize;
	
   // If form of extension unspecified, default to symmetric
   // extensions for symmetrical filters and periodic extensions for
   // asymmetrical filters
   if (symExt == -1){
		symExt = wavelet->symmetric;
	}
	
	lowSize = (int *)calloc(nsteps, sizeof(int));
   highSize = (int *)calloc(nsteps, sizeof(int));
	
   lowSize[0] = (size+1)/2;
   highSize[0] = size/2;
	
   for (i = 1; i < nsteps; i++) {
		lowSize[i] = (lowSize[i-1]+1)/2;
		highSize[i] = lowSize[i-1]/2;
   }
	
   data [0] = (Real*)calloc(2*wavelet->npad+size, sizeof(Real));
   data [1] = (Real*)calloc(2*wavelet->npad+size, sizeof(Real));
		
		copy (input, data[currentIndex]+wavelet->npad, lowSize[nsteps-1]);
	
   while (nsteps--)  {
		
		// grab the next high-pass component
		copy (input + lowSize[nsteps], data[currentIndex]+wavelet->npad+lowSize[nsteps], 
			highSize[nsteps]);
		
		// Combine low-pass data (first 1/2^n of signal) with high-pass
		// data (next 1/2^n of signal) to get higher resolution low-pass data
		WaveletInvertStep(wavelet, data[currentIndex], data[1-currentIndex], 
			lowSize[nsteps]+highSize[nsteps], symExt);
		
		// Now pass low-pass data (first 1/2 of signal) back to
		// transform routine
		currentIndex = 1 - currentIndex;
   }
	
   // Copy inverted signal to output signal
   copy (data[currentIndex]+wavelet->npad, output, size);
	
   free(highSize);
   free(lowSize);
	
   free(data[1]);
   free(data[0]);
}


/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/
void WaveletTransform2D(WAVELET *wavelet, Real *input, Real *output, int hsize, int vsize,
								int nsteps, int symExt)
{
   int j;
   int hLowSize = hsize, hHighSize;
   int vLowSize = vsize, vHighSize;
	Real *temp_in, *temp_out;
	
   // If form of extension unspecified, default to symmetric
   // extensions for symmetrical filters and periodic extensions for
   // asymmetrical filters
   if (symExt == -1)
		symExt = wavelet->symmetric;
	
   temp_in = (Real*)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
   temp_out = (Real*)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
	
   copy (input, output, hsize*vsize);
	
   while (nsteps--){
		if ((hLowSize <= 2 || vLowSize <= 2) && symExt == 1) {
			WaveletWarning ("Reduce # of transform steps or increase signal size.\n");
			WaveletWarning ("or switch to periodic extension\n");
			WaveletError ("Low pass subband is too small\n");
      }
		
      // Do a convolution on the low pass portion of each row
      for (j = 0; j < vLowSize; j++)  {
			// Copy row j to data array
			copy (output+(j*hsize), temp_in+wavelet->npad, hLowSize);
			
			// Convolve with low and high pass filters
			WaveletTransformStep (wavelet, temp_in, temp_out, hLowSize, symExt);
			
			// Copy back to image
			copy (temp_out+wavelet->npad, output+(j*hsize), hLowSize);
      }
		
      // Now do a convolution on the low pass portion of  each column
      for (j = 0; j < hLowSize; j++)  {
			// Copy column j to data array
			copy_p1_skip (output+j, hsize, temp_in+wavelet->npad, vLowSize);
			
			// Convolve with low and high pass filters
			WaveletTransformStep (wavelet, temp_in, temp_out, vLowSize, symExt);
			
			// Copy back to image
			copy_p2_skip (temp_out+wavelet->npad, output+j, hsize, vLowSize);
      }
		
      // Now convolve low-pass portion again
      hHighSize = hLowSize/2;
      hLowSize = (hLowSize+1)/2;
      vHighSize = vLowSize/2;
      vLowSize = (vLowSize+1)/2;
   }
	
   free(temp_out);
   free(temp_in);
	
}


/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/
void WaveletInvert2D(WAVELET *wavelet, Real *input, Real *output, int hsize, int vsize,
							int nsteps, int symExt)
{   
	int i, j;
	int *hLowSize, *hHighSize, *vLowSize, *vHighSize;
	Real *temp_in, *temp_out;
	
   // If form of extension unspecified, default to symmetric
   // extensions for symmetrical filters and periodic extensions for
   // asymmetrical filters
   if (symExt == -1)
		symExt = wavelet->symmetric;
	
   hLowSize = (int *) calloc(nsteps, sizeof(int));
   hHighSize = (int *) calloc(nsteps, sizeof(int));
   vLowSize =(int *) calloc(nsteps, sizeof(int));
   vHighSize =(int *) calloc(nsteps, sizeof(int));
	
   hLowSize[0] = (hsize+1)/2;
   hHighSize[0] = hsize/2;
   vLowSize[0] = (vsize+1)/2;
   vHighSize[0] = vsize/2;
	
   for (i = 1; i < nsteps; i++) {
		hLowSize[i] = (hLowSize[i-1]+1)/2;
		hHighSize[i] = hLowSize[i-1]/2;
		vLowSize[i] = (vLowSize[i-1]+1)/2;
		vHighSize[i] = vLowSize[i-1]/2;
   }
	
   temp_in = (Real *)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
   temp_out = (Real *)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
	
   copy (input, output, hsize*vsize);
	
   while (nsteps--)  {
      // Do a reconstruction for each of the columns
      for (j = 0; j < hLowSize[nsteps]+hHighSize[nsteps]; j++)  {
			// Copy column j to data array
			copy_p1_skip (output+j, hsize, temp_in+wavelet->npad, 
				vLowSize[nsteps]+vHighSize[nsteps]);
			
			// Combine low-pass data (first 1/2^n of signal) with high-pass
			// data (next 1/2^n of signal) to get higher resolution low-pass data
			WaveletInvertStep (wavelet, temp_in, temp_out,
				vLowSize[nsteps]+vHighSize[nsteps], symExt);
			
			// Copy back to image
			copy_p2_skip (temp_out+wavelet->npad, output+j, hsize,
				vLowSize[nsteps]+vHighSize[nsteps]);
      }
		
      // Now do a reconstruction pass for each row
      for (j = 0; j < vLowSize[nsteps]+vHighSize[nsteps]; j++)  {
			// Copy row j to data array
			copy (output + (j*hsize), temp_in+wavelet->npad, 
				hLowSize[nsteps]+hHighSize[nsteps]);
			
			// Combine low-pass data (first 1/2^n of signal) with high-pass
			// data (next 1/2^n of signal) to get higher resolution low-pass data