#include "Dwt.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CDwt::CDwt()
{
   FilterSet(0);
}

CDwt::~CDwt()
{
   delete []h;
   delete []h1;
   delete []g;
   delete []g1;
}

void CDwt::Dwt1D(double *buffer, int buflen, int hBen, int hEnd, int h1Ben, int h1End, double *h, double *g)
{
	double *p;
	int m,n,tempn;
	
	p = new double [buflen];
	for (m=0; m<buflen; m++){
		*(p+m) = *(buffer+m);
		*(buffer+m) = 0;
	}
	
	for (m=0; m<(buflen>>1); m++){
		for (n=hBen+2*m; n<=hEnd+2*m; n++){
			
			tempn=n;
			/* 恒等延拓  */
			/*if ( tempn<0 ) tempn=0;
			if ( tempn>= buflen ) tempn=buflen-1;*/
			
			/*  周期延拓/ *
			if ( tempn<0 ) tempn=tempn+buflen;
			if ( tempn >buflen-1 ) tempn=tempn-buflen;
			
			/*   對稱延拓 */
			/**/if ( tempn<0 ) tempn=-tempn;
			if ( tempn >buflen-1 ) tempn=2*buflen-tempn-2;
			
			/*   0 填充  */
			/*if ( tempn<0 ) tempn=0;
			if ( tempn >buflen-1 ) tempn=0;*/
			
			
			*(buffer+m) = *(buffer+m)+*(p+tempn)*(*(h+n-2*m-hBen));
		}// for n
		
		for (n=2*m+1-h1End; n<=2*m+1-h1Ben; n++){
			tempn=n;
			/* 恒等延拓  */
			/*if ( tempn<0 ) tempn=0;
			if ( tempn>= buflen ) tempn=buflen-1;*/
			
			/* 周期延拓 /*
			if ( tempn<0 ) tempn=tempn+buflen;
			if ( tempn >buflen-1 ) tempn=tempn-buflen;
			
			/*   對稱延拓 */
			/**/if ( tempn<0 ) tempn=-tempn;
			if ( tempn >buflen-1 ) tempn=2*buflen-tempn-2;
			
			/*   0 填充  */
			/*if ( tempn<0 ) tempn=0;
			if ( tempn >buflen-1 ) tempn=0;*/
			
			*(buffer+m+(buflen>>1)) = *(buffer+m+(buflen>>1))+*(p+tempn)*(*(g+n-2*m-1+h1End));
		}// for n
	}// for m
    delete []p;
}

void CDwt::IDwt1D(double buffer[], int buflen, int hBen, int hEnd, int h1Ben, int h1End, double *h, double *g)
{
	double *p;
	int m,n,tempn;
	p = new double [buflen];
	for (m=0; m<buflen; m++){
		*(p+m) = buffer[m];
		buffer[m] = 0;
	}

	for (m=0; m<buflen; m++){
		for(n = int(floor((m-h1End+1)/2.)); n<=(m-h1Ben)/2.; n++){
			tempn=n;
				/* 恒等延拓  */
		/*	if (tempn <0 ) tempn=0;
			if ( tempn > (buflen>>1)-1 ) tempn = (buflen>>1) -1;*/

			/*  周期延拓 /*
			if ( tempn<0 ) tempn=tempn+(buflen>>1);
			if ( tempn >(buflen>>1)-1 ) tempn=tempn-(buflen>>1);

			/*   對稱延拓 */
			/**/if ( tempn<0 ) tempn=-tempn;
			if ( tempn >(buflen>>1)-1 ) tempn=2*(buflen>>1)-tempn-1;

			/*   0 填充  */
			/*if ( tempn<0 ) tempn=0;
			if ( tempn >(buflen>>1)-1 ) tempn=0;*/

			buffer[m] = buffer[m] + *(p+tempn)*(*(h+m-2*n-h1Ben));
		}
		for(n = int(floor((m+hBen)/2.)); n<=(m+hEnd-1)/2.; n++){
			tempn=n;
				/* 恒等延拓  */
			/*if (tempn<0) tempn=0;
			if (tempn>(buflen>>1)-1) tempn=(buflen>>1)-1;*/

			/*  周期延拓 /*
			if ( tempn<0 ) tempn=tempn+(buflen>>1);
			if ( tempn >(buflen>>1)-1 ) tempn=tempn-(buflen>>1);

			/*   對稱延拓 */
			/**/if ( tempn<0 ) tempn=-tempn-1;
			if ( tempn >(buflen>>1)-1 ) tempn=2*(buflen>>1)-tempn-2;

			/*   0 填充  */
			/*if ( tempn<0 ) tempn=0;
			if ( tempn >(buflen>>1)-1 ) tempn=0;*/

			buffer[m]=buffer[m]+*(p+tempn+(buflen>>1))*(*(g+m-2*n-1+hEnd));
		}
	}

	delete []p;
}

void CDwt::Dwt2D(double *buffer, int height, int width, int lv)
{
	int i,j,k;
	int nheight,nwidth;

	for ( k=0; k<lv; k++){

		nheight=height>>k;
		nwidth=width>>k;
		
		double *p;
		p = new double [nwidth];
		for (i=0; i<nheight; i++){
			for(j=0; j<nwidth; j++)
				*(p+j) = *(buffer+i*width+j);
			Dwt1D(p,nwidth,hBen,hEnd,h1Ben,h1End,h,g);
			for(j=0; j<nwidth; j++)
				*(buffer+i*width+j) = *(p+j);
		}
		delete []p;

		double *p1;
		p1 = new double [nheight];
		for(j=0; j<nwidth; j++){
			for (i=0; i<nheight; i++)
				*(p1+i) = *(buffer+i*width+j);
			Dwt1D(p1,nheight,hBen,hEnd,h1Ben,h1End,h,g);
			for(i=0; i<nwidth; i++)
				*(buffer+i*width+j) = *(p1+i);
		}
		delete []p1;
	}

}

void CDwt::IDwt2D(double *buffer, int height, int width, int lv)
{
	int i,j,k;
	int nheight,nwidth;
	
	for (k=0; k<lv; k++){
		nheight=height>>(lv-k-1);
		nwidth=width>>(lv-k-1);
		
		double *p;
		p = new double [nheight];
		for(j=0; j<nwidth; j++){
			for (i=0; i<nheight; i++)
				*(p+i) = *(buffer+i*width+j);
			IDwt1D(p,nheight,hBen,hEnd,h1Ben,h1End,h1,g1);
			for(i=0; i<nwidth; i++)
				*(buffer+i*width+j) = *(p+i);
		}
		delete []p;

		double *p1;
		p1 = new double [nwidth];		
		for (i=0; i<nheight; i++){
			for(j=0; j<nwidth; j++)
				*(p1+j) = *(buffer+i*width+j);
			IDwt1D(p1,nwidth,hBen,hEnd,h1Ben,h1End,h1,g1);
			for(j=0; j<nwidth; j++)
				*(buffer+i*width+j) = *(p1+j);
		}
		delete []p1;
	}
	
}

void CDwt::FilterSet(float t)
{
	hLen = 9;
	h1Len = 7;
	hBen = -4;
	hEnd = 4;
	h1Ben = -3;
	h1End = 3;
	h = new double [hLen];
	h1 = new double [h1Len];
	g = new double [h1Len];
	g1 = new double [hLen];

    *h = (8*t*t*t-6*t*t+3*t)/(1+2*t)*(1/32.)*sqrt(2);
	*(h+1) = (-16*t*t*t+20*t*t-12*t+3)/(1+2*t)*(1/32.)*sqrt(2);
	*(h+2) = (2*t-3)/(1+2*t)*(1/8.)*sqrt(2);
	*(h+3) = (16*t*t*t-20*t*t+28*t+5)/(1+2*t)*(1/32.)*sqrt(2);
	*(h+4) = (-8*t*t*t+6*t*t+5*t+20)/(1+2*t)*(1/16.)*sqrt(2);
	*(h+5) = *(h+3);
	*(h+6) = *(h+2);
	*(h+7) = *(h+1);
	*(h+8) = *(h+0);

	*h1 = -t/16.*sqrt(2);
	*(h1+1) = (1-2*t)/16.*sqrt(2);
	*(h1+2) = (t+4)/16.*sqrt(2);
	*(h1+3) = (3+2*t)/8.*sqrt(2);
	*(h1+4) = *(h1+2);
	*(h1+5) = *(h1+1);
	*(h1+6) = *(h1+0);

	
	*g1 = *(h+8);
	*(g1+1) = -*(h+7);
	*(g1+2) =  *(h+6);
	*(g1+3) = -*(h+5);
	*(g1+4) =  *(h+4);
	*(g1+5) = -*(h+3);
	*(g1+6) =  *(h+2);
	*(g1+7) = -*(h+1);
	*(g1+8) =  *h;
	
	
	*g = -*(h1+6);
	*(g+1) =  *(h1+5);
	*(g+2) = -*(h1+4);
	*(g+3) =  *(h1+3);
	*(g+4) = -*(h1+2);
	*(g+5) =  *(h1+1);
	*(g+6) = -*h1;

}