/*********************

  I. 9 multiplies for each pel, for each level
		this could be reduced by going explicility two-dimensional
		(since you can merge multiplies) but that's a real hassle
	9 * (10 cycles) * (256 * 256 dimensions) * log2(4 levels) = 10 mill = 0.06 sec , is half the time

 II. unnecessary mem-copies :
		A. transposing
		B. copy in from ints to doubles
		C. copy to temp-work space x_alloc for each line

III. speeds : time to do 256x256x1 , levels=4
		dwt: 0.15 sec
		dct: 0.10 sec
		spt: 0.05 sec
		haar:0.025 sec

***********************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <crblib/inc.h>
#include <crblib/memutil.h>

#define WT_ROUND_IN		0
#define WT_ROUND_OUT	0.5

#define F_EXTPAD 4
#define D_EXTPAD 2

static double *x_alloc = NULL;	//** temp work

void forwardDWT(double *x_in, int N)
{
int i, n, half;
double *x,*r,*d;

	x = x_alloc + F_EXTPAD;
	memcpy(x,x_in,sizeof(double)*N);
	for(i=1;i<=F_EXTPAD;i++) {
		x[-i] = x[i];
		x[(N-1) + i] = x[(N-1) - i];
	}

	half = N>>1;  

	r = x_in; d = x_in + half;
	for (n=half;n--;) {

		*r++ = 0.0378285 * (x[ 4] + x[-4]) - 0.0238495 * (x[ 3] + x[-3]) - 0.110624 * (x[ 2] + x[-2]) + 
					0.377403 * (x[ 1] + x[-1]) + 0.852699 * x[0];

		x++;

		*d++ =  - 0.064539 * (x[ 3] + x[-3]) + 0.040689 * (x[ 2] + x[-2]) + 
				0.418092 * (x[ 1] + x[-1]) - 0.788486 * x[ 0] ;

		x++;
    }
}

void inverseDWT(double *x, int N)
{
int i, n, half;
double *r, *d;

	half = N/2;

	r = x_alloc + D_EXTPAD;  
	d = x_alloc + D_EXTPAD + half+D_EXTPAD+D_EXTPAD;  
	memcpy(r,x,half*sizeof(double));
	memcpy(d,x+half,half*sizeof(double));

	for(i=1;i<=D_EXTPAD;i++) {
		r[-i] = r[i];
		r[(half-1)+i] = r[half - i];
		d[-i] = d[i-1];
		d[(half-1)+i] = d[(half-1) - i];
	}
 
	for (n = half; n--;) {

		*x++ = 0.788486 * r[0] - 0.0406894 * ( r[1] + r[-1] )
		     - 0.023849 * (d[1] + d[-2]) + 0.377403 * (d[0] + d[-1]);

		*x++ = 0.418092 * (r[1] + r[0]) - 0.0645389 * (r[2] + r[-1])
		     - 0.037829 * (d[2] + d[-2]) + 0.110624 * (d[1] + d[-1]) - 0.852699 * d[0];

		d++; r++;
	}
}

void waveletTransform2D(int **raw,double **wave, int width, int height, int levels,bool inverse)
{
int x, y, w, yw, h, l;
int *raw_ptr;
double *buffer,*wave_ptr;

    /* Check the dimensions for compatability. */

    if (width%(1 << levels) || height%(1 << levels)) {
		errputs("width and height must be divisible by 2^levels");
		exit(10);
	}

	if ( (x_alloc = malloc(sizeof(double)*(width+height+F_EXTPAD+F_EXTPAD))) == NULL )  {
		errputs("malloc failed"); exit(10);
	}
  
    /* Allocate a work array (for transposing columns) */
    
   	if ( (buffer = newarray(double,height)) == NULL ) {
		errputs("malloc failed"); exit(10);
	}

    /* Compute the wavelet transform. */
  
	if ( !inverse ) { /* forward transform. */

		/** raw -> wave **/

		wave_ptr = wave[0]; raw_ptr = raw[0];
		for(x=width*height;x--;) *wave_ptr++ = (double)(*raw_ptr++) + WT_ROUND_IN;
		
		for (l = 0; l < levels; l++) {
			w = width >> l;
			h = height >> l;
      
			/* Rows. */
	
			for (y = 0; y < h; y++)
				forwardDWT(wave[y], w);
    
			/* Columns. */
	
			wave_ptr = wave[0];
			for (x = 0; x < w; x++) {
				for (y = 0,yw=0; y < h; y++,yw+=width) buffer[y] = wave_ptr[yw];
				forwardDWT(buffer, h);
				for (y = 0,yw=0; y < h; y++,yw+=width) wave_ptr[yw] = buffer[y];
				wave_ptr ++;
			}
		}

    } else {

		for (l = levels-1; l >= 0; l--) {
			w = width >> l;
			h = height >> l;

			/* Columns. */
	
			wave_ptr = wave[0];
			for (x = 0; x < w; x++) {
				for (y = 0,yw=0; y < h; y++,yw+=width) buffer[y] = wave_ptr[yw];
				inverseDWT(buffer, h);
				for (y = 0,yw=0; y < h; y++,yw+=width) wave_ptr[yw] = buffer[y];
				wave_ptr ++;
			}

			/* Rows. */
	
			for (y = 0; y < h; y++)
				inverseDWT(wave[y], w);
    	}

		wave_ptr = wave[0];
		raw_ptr = raw[0];
		for(x=width*height;x--;)
			*raw_ptr++ = (int) ( (*wave_ptr++)  + WT_ROUND_OUT );
	}

	free(x_alloc); x_alloc = NULL;
	free(buffer);
}