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

#define MULT(x,y)	((((x)*(y))+16383)>>15)

/**********

#define ONE	(1<<15)
#define	C_A	((int)(ONE* 1.586134342  + 0.5))
#define C_B	((int)(ONE* 0.05298011854+ 0.5))
#define C_C ((int)(ONE* 0.8829110762 + 0.5))
#define C_D	((int)(ONE* 0.4435068522 + 0.5))

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

#define	C_A	51974
#define C_B	1736
#define C_C	28931
#define C_D	14533

static void do_tdec_line(int * to,int *from,int len)
{
int x,*ptr,*low,*high,half;

	half = len>>1;

	assert(len & 1 == 0);
	assert(half >= 2);

	low  = to; high = to + half; ptr = from;

		*high = ptr[1] - MULT(C_A,( ptr[0] + ptr[2] ));
		*low  = ptr[0] - MULT(C_B,( high[0] + high[0]));
		low++; high++; ptr += 2;

	for(x=half-2;x--;) {
		*high = ptr[1] - MULT(C_A,( ptr[0] + ptr[2] ));
		*low  = ptr[0] - MULT(C_B,( high[0] + high[-1]));
		low++; high++; ptr += 2;
	}

		*high = ptr[1] - MULT(C_A,( ptr[0] + ptr[0] ));
		*low  = ptr[0] - MULT(C_B,( high[0] + high[-1]));

	low  = to; high = to + half;

		*high += MULT(C_C,( low[0] + low[1] ));
		*low  += MULT(C_D,( high[0] + high[0]));
		low++; high++;

	for(x=half-2;x--;) {
		*high += MULT(C_C,( low[0] + low[1] ));
		*low  += MULT(C_D,( high[0] + high[-1]));
		low++; high++;
	}

		*high += MULT(C_C,( low[0] + low[0] ));
		*low  += MULT(C_D,( high[0] + high[-1]));

}

static void un_tdec_line(int *to,int *from,int len)
{
int x,*ptr,*low,*high,half;

	half = len>>1;

	low  = from; high = from + half;
	for(x=0;x<half;x++) {
		if ( x == 0 )		*low  -= MULT(C_D,( high[0] + high[0]));
		else				*low  -= MULT(C_D,( high[0] + high[-1]));
		low++; high++;
	}

	low  = from; high = from + half;
	for(x=0;x<half;x++) {
		if ( x == half-1 )	*high -= MULT(C_C,( low[0] + low[0] ));
		else				*high -= MULT(C_C,( low[0] + low[1] ));
		low++; high++;
	}

	low  = from; high = from + half; ptr = to;
	for(x=0;x<half;x++) {
		if ( x == 0 )		ptr[0] = *low  + MULT(C_B,( high[0] + high[0]));
		else				ptr[0] = *low  + MULT(C_B,( high[0] + high[-1]));
		low++; high++; ptr += 2;
	}

	low  = from; high = from + half; ptr = to;
	for(x=0;x<half;x++) {
		if ( x == half-1 )	ptr[1] = *high + MULT(C_A,( ptr[0] + ptr[0] ));
		else				ptr[1] = *high + MULT(C_A,( ptr[0] + ptr[2] ));
		low++; high++; ptr += 2;
	}

}

void l97_2D(int **rows, int width, int height, int levels,bool inverse)
{
int x, y, w, h, l;
int *buffer,*tempbuf,*temprow;

    if (width%(1 << (levels+1)) || height%(1 << (levels+1)))
		errexit("width and height must be divisible by 2^(levels+1)");
  
    /* Allocate a work array (for transposing columns) */
    
   	if ( (buffer = newarray(int,height+max(width,height)+height)) == NULL )
		errexit("malloc failed");
	temprow = buffer+height;
	tempbuf = buffer+height+height;

	if ( !inverse ) {

		for (l = 0; l < levels; l++) {
			w = width >> l;
			h = height >> l;
      
			/* Rows */
	
			do_tdec_line(temprow,rows[h-1],w);
			for (y = h-2; y >=0; y--) {
				do_tdec_line(rows[y+1],rows[y],w);
			}
	
			/* Columns */
	
			for (x = 0; x < w; x++) {
					for (y = 1; y < h; y++) buffer[y-1] = rows[y][x];
					buffer[h-1] = temprow[x];
				do_tdec_line(tempbuf,buffer,h);
					for (y = 0; y < h; y++) rows[y][x] = tempbuf[y];
			}
		}

    } else {

		for (l = levels-1; l >= 0; l--) { /** backwards in scale **/
			w = width >> l;
			h = height >> l;

			/* Columns */
	
			for (x = 0; x < w; x++) {
					for (y = 0; y < h; y++) buffer[y] = rows[y][x];
				un_tdec_line(tempbuf,buffer,h);
					for (y = 0; y < h-1; y++) rows[y+1][x] = tempbuf[y];
					temprow[x] = tempbuf[h-1];
			}

			/* Rows */
			for (y = 0; y < h-1; y++) {
				un_tdec_line(rows[y],rows[y+1],w);
			} 
			un_tdec_line(rows[h-1],temprow,w);
		}
	}

	free(buffer);
}


void l97Quad(int *band,int w,int h,int fullw,bool inverse)
{
int x, y;
int *buffer,*tempbuf,*bptr,*temprow;

   	if ( (buffer = newarray(int,h+h+max(w,h))) == NULL ) {
		errputs("malloc failed"); exit(10);
	}
	temprow = buffer+h;
	tempbuf = buffer+h+h;
  
	if ( !inverse ) { /* forward transform. */
	
		bptr = band + (h-1)*fullw;
		do_tdec_line(temprow,bptr,w);
		for (y = (h-1); y--;) {
			bptr -= fullw;
			do_tdec_line(bptr+fullw,bptr,w);
		}
    
		for (x = 0; x < w; x++) {
			bptr = band + x + fullw;
			for (y = 0; y < (h-1); y++) { buffer[y] = *bptr; bptr += fullw; }
			buffer[h-1] = temprow[x];
			do_tdec_line(tempbuf,buffer,h);
			bptr = band + x;
			for (y = 0; y < h; y++) { *bptr = tempbuf[y]; bptr += fullw; }
		}

    } else {

		for (x = 0; x < w; x++) {
			bptr = band + x;
			for (y = 0; y < h; y++) { buffer[y] = *bptr; bptr += fullw; }
			un_tdec_line(tempbuf,buffer,h);
			bptr = band + x + fullw;
			for (y = 0; y < h-1; y++) { *bptr = tempbuf[y]; bptr += fullw; }
			temprow[x] = tempbuf[h-1];
		}

		bptr = band;
		for (y = (h-1); y--; ){
			un_tdec_line(bptr,bptr+fullw,w);
			bptr += fullw;
		}
		un_tdec_line(bptr,temprow,w);

	}

	free(buffer);
}