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

#define F_EXTPAD 4
#define D_EXTPAD 2

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

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

	x = x_alloc + F_EXTPAD;
	memcpy(x,x_in,sizeof(int)*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 inversef97(int *x, int N)
{
int i, n, half;
int *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(int));
	memcpy(d,x+half,half*sizeof(int));

	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 f97_2D(int **rows, int width, int height, int levels,bool inverse)
{
int x, y, w, yw, h, l;
int *buffer,*rows_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(int)*(width+height+F_EXTPAD+F_EXTPAD))) == NULL )  {
		errputs("malloc failed"); exit(10);
	}
  
    /* Allocate a work array (for transposing columns) */
    
   	if ( (buffer = newarray(int,height)) == NULL ) {
		errputs("malloc failed"); exit(10);
	}

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

		for (l = 0; l < levels; l++) {
			w = width >> l;
			h = height >> l;
      
			/* Rows. */
	
			for (y = 0; y < h; y++)
				forwardf97(rows[y], w);
    
			/* Columns. */
	
			rows_ptr = rows[0];
			for (x = 0; x < w; x++) {
				for (y = 0,yw=0; y < h; y++,yw+=width) buffer[y] = rows_ptr[yw];
				forwardf97(buffer, h);
				for (y = 0,yw=0; y < h; y++,yw+=width) rows_ptr[yw] = buffer[y];
				rows_ptr ++;
			}
		}

    } else {

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

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

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

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


void f97Quad(int *wave,int width,int height,int fullw,bool inverse)
{
int x, y, yw;
int *buffer,*wave_ptr;

	if ( (x_alloc = malloc(sizeof(int)*(width+height+F_EXTPAD+F_EXTPAD))) == NULL )  {
		errputs("malloc failed"); exit(10);
	}
 
   	if ( (buffer = newarray(int,height)) == NULL ) {
		errputs("malloc failed"); exit(10);
	}

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

		wave_ptr = wave;
		for (y = 0; y < height; y++) {
			forwardf97(wave_ptr, width);
			wave_ptr += fullw;
		}
    
		wave_ptr = wave;
		for (x = 0; x < width; x++) {
			for (y = 0,yw=0; y < height; y++,yw+=fullw) buffer[y] = wave_ptr[yw];
			forwardf97(buffer, height);
			for (y = 0,yw=0; y < height; y++,yw+=fullw) wave_ptr[yw] = buffer[y];
			wave_ptr ++;
		}

    } else {

		wave_ptr = wave;
		for (x = 0; x < width; x++) {
			for (y = 0,yw=0; y < height; y++,yw+=fullw) buffer[y] = wave_ptr[yw];
			inversef97(buffer, height);
			for (y = 0,yw=0; y < height; y++,yw+=fullw) wave_ptr[yw] = buffer[y];
			wave_ptr ++;
		}

		wave_ptr = wave;
		for (y = 0; y < height; y++) {
			inversef97(wave_ptr, width);
			wave_ptr += fullw;
		}
    }

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