// LiftDwt.cpp: implementation of the CLiftDwt class.

//  此程序為完全提升格式的小波變換，所謂完全是指，最后的拉伸變換也進行了提升分解，即分解成了三個提升結構
//  程序中了使用的小波濾波器為帶參數的9-7小波濾波器
//
//  程序包括以下部分： 
//
//  Dwt2D(double buffer[], int height, int width, int lv) 二維提升小波變換，buffer進行變換的數據，height,width數據的尺寸，lv小波分解的層數
//  IDwt2D(double buffer[], int height, int width, int lv) 二維提升小波反變換
//  
//  Dwt1D(double *buffer, int buflen)  一維提升小波分解 buflen 數據的長度
//  IDwt1D(double *buffer, int buflen) 
//   
//  filterset(double t)  設置小波濾波器，小波和對偶小波濾波器之和必須為sqrt(2)，否則編碼效果不好
//
//  數據延拓方式有兩種： 周期延拓和對稱延拓
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
//#include "ImgPro.h"
#include "LiftDwt97.h"
#include "math.h"

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

//#define MAX(x,y) (x>y?x:y)
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CLiftDwt97::CLiftDwt97()
{
	p_t = 0.730174;
	shift = 4;
	filterset(p_t);  // 設置濾波器的系數
}

CLiftDwt97::~CLiftDwt97()
{
	delete []h;
	delete []h1;

}

void CLiftDwt97::Dwt2D(double buffer[], int height, int width, int lv)
{
	int i,j,k;
 	int nheight,nwidth;
    double *pbuffer1,*pbuffer2;
 
	for ( k=0; k<lv; k++)
	{
       nheight=height>>k;
	   nwidth=width>>k;

	   pbuffer1 = new double [nheight];	 
	   pbuffer2 = new double [nwidth];
       //行變換
	   for (i=0; i<nheight; i++)
		{
			for(j=0; j<nwidth; j++)
				*(pbuffer1+j) = *(buffer+i*width+j);
			
			Dwt1D(pbuffer1,nwidth);
			
			for(j=0; j<nwidth; j++)
				 *(buffer+i*width+j) = *(pbuffer1+j);
		}
		delete []pbuffer1;
       
		//列變換
		for(j=0; j<nwidth; j++)
		{
			for (i=0; i<nheight; i++)
				*(pbuffer2+i)= *(buffer+i*width+j);

			Dwt1D(pbuffer2,nheight);
			
			for(i=0; i<nwidth; i++)
				 *(buffer+i*width+j) = *(pbuffer2+i);
		}
		delete []pbuffer2;
	}
}

void CLiftDwt97::IDwt2D(double buffer[], int height, int width, int lv)
{

    int i,j,k;
 	int nheight,nwidth;
	double *pbuffer1,*pbuffer2;
	  
	for (k=0; k<lv; k++)
	{
		nheight=height>>(lv-k-1);
 		nwidth=width>>(lv-k-1);
	    
		pbuffer1 = new double [nheight];	 
	    pbuffer2 = new double [nwidth];
        
		//列反變換
		for(j=0; j<nwidth; j++)
		{
 			for (i=0; i<nheight; i++)
 				*(pbuffer1+i) =  *(buffer+i*width+j);
 		
			IDwt1D(pbuffer1,nheight);
 			
			for(i=0; i<nwidth; i++)
 				 *(buffer+i*width+j) = *(pbuffer1+i);
 		}
 		delete []pbuffer1;
 
        //行反變換
	
 		for (i=0; i<nheight; i++)
		{
 			for(j=0; j<nwidth; j++)
 				*(pbuffer2+j)=*(buffer+i*width+j);
 			IDwt1D(pbuffer2,nwidth);
 			for(j=0; j<nwidth; j++)
 				 *(buffer+i*width+j) = *(pbuffer2+j);
 		}
		delete []pbuffer2;
 	}

}

void CLiftDwt97::Dwt1D(double *buffer, int buflen)
{
	int i;
	int itemp;

/*    使用周期延拓方式	*/
	double *d, *s, *p;
	p = new double [buflen+2*shift];
	d = new double [(buflen>>1)+shift];
	s = new double [(buflen>>1)+shift];

/*	for (i=0; i<shift; i++){
		*(p+i) = *(buffer+buflen-shift+i);
		*(p+buflen+shift+i) = *(buffer+i);
	}
	for (i=0; i<buflen; i++) *(p+i+shift) = *(buffer+i);
	/* 數據周期延拓結束  計算p1 = d1          p2 = s1 */

	/*  對稱延拓 */
	for (i=0; i<shift; i++)
	{
		*(p+i) = *(buffer+shift-i);
		*(p+buflen+shift+i) = *(buffer+buflen-i-2);
	}
	for (i=0; i<buflen; i++) 
		*(p+i+shift) = *(buffer+i);
    /*  結束  */
	for (i=0; i<(buflen>>1)+shift; i++)
	{
		itemp = i-(shift>>1);
		*(d+i) = *(p+shift+2*itemp+1) + p_alfa*( *(p+shift+2*itemp)+*(p+shift+2*itemp+2) );
	}

	for (i=0; i<(buflen>>1)+shift-1; i++)
	{
		itemp = i+1-(shift>>1);
		*(s+i+1) = *(p+shift+2*itemp) + p_beta*( *(d+i+1)+*(d+i+1-1) );
	}
   /*  d1 s1 計算結束 下計算 d2  */
	for (i=0; i<(buflen>>1)+shift-1; i++)
	{
		//itemp = i-(shift>>1);
		*(p+shift+(buflen>>1)+i) =  *(d+i) + p_gama*( *(s+i)+*(s+i+1)); 
	}
	/*  s2 */
    for (i=0; i<(buflen>>1)+shift-1; i++)
	{
	//	itemp = i-(shift>>1);
		*(p+i+1) = *(s+i+1) + p_delta*(*(p+(buflen>>1)+shift+i+1) + *(p+(buflen>>1)+shift+i+1-1) );
	}

	/* d3 */
	for (i=0; i<(buflen>>1); i++)
		*(d+i) = *(p+i+(buflen>>1)+shift+(shift>>1)) +p_kesa*(1-p_kesa)*(*(p+i+(shift>>1)));
	/* s3 */
	for (i=0; i<(buflen>>1); i++)
		*(s+i) = *(p+i+(shift>>1)) -1/p_kesa*(*(d+i));
    /* d4 */
	for (i=0; i<(buflen>>1); i++)
		*(buffer+i+(buflen>>1)) = *(d+i) +(p_kesa-1)*(*(s+i));

	/* s4 */
	for (i=0; i<(buflen>>1); i++)
		*(buffer+i) = *(s+i) + *(buffer+i+(buflen>>1));

/*	for (i=0; i<(buflen>>1); i++){
   *(buffer+i) = p_kesa*(*(p+i+(shift>>1)));  
 	 *(buffer+i+(buflen>>1)) = 1/p_kesa*(*(p+i+(buflen>>1)+shift+(shift>>1)));	
	}
*/
    delete []p;
	delete []d;
	delete []s;
    
}

void CLiftDwt97::IDwt1D(double *buffer, int buflen)
{
	int i;
    double *p1, *p2, *s, *d;
	p1 = new double [(buflen>>1)+shift];
	p2 = new double [(buflen>>1)+shift];
	s  = new double [(buflen>>1)+shift];
	d  = new double [(buflen>>1)+shift];

   /* s3 */
   for (i=0; i<(buflen>>1); i++)
	   *(s+i) = *(buffer+i) - *(buffer+i+(buflen>>1));
   /* d3 */
   for (i=0; i<(buflen>>1); i++)
	   *(d+i) =*(buffer+i+(buflen>>1))- (p_kesa-1)*(*(s+i));
   /*  s2 */
   for (i=0; i<(buflen>>1); i++)
	   *(buffer+i) = *(s+i) + (1/p_kesa)*(*(d+i));
   /*  d2 */
   for (i=0; i<(buflen>>1); i++)
	   *(buffer+i+(buflen>>1)) = *(d+i) - p_kesa*(1-p_kesa)*(*(buffer+i));

/*    延拓方式為周期延拓     /*
	for (i=0; i<(shift>>1); i++){
		p1[i] = buffer[i+(buflen>>1)-(shift>>1)];
		p1[i+(shift>>1)+(buflen>>1)] = buffer[i];
		p2[i] = buffer[i+buflen-(shift>>1)];
		p2[i+(shift>>1)+(buflen>>1)] = buffer[i+(buflen>>1)];
	}
	for (i=0; i<(buflen>>1); i++){
		p1[i+(shift>>1)] = buffer[i];
		p2[i+(shift>>1)] = buffer[i+(buflen>>1)];
	}

/*  延拓結束 */	

	/*  對稱延拓	*/
   	for (i=0; i<(shift>>1); i++){
		p1[i] = buffer[(shift>>1)-i];
		p1[i+(shift>>1)+(buflen>>1)] = buffer[(buflen>>1)-i-1];
		p2[i] = buffer[(buflen>>1)+(shift>>1)-i-1];
		p2[i+(shift>>1)+(buflen>>1)] = buffer[buflen-i-2];
	}
	for (i=0; i<(buflen>>1); i++){
		p1[i+(shift>>1)] = buffer[i];
		p2[i+(shift>>1)] = buffer[i+(buflen>>1)];
	}
/*   對稱延拓結束  */

   /*  s1  */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(s+i+1) = *(p1+i+1) - p_delta*( *(p2+i+1) + *(p2+i+1-1));
   /*  d1 */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(d+i) = *(p2+i) - p_gama*(*(s+i) + *(s+i+1));
    /* p1 = s0 */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(p1+i+1) = *(s+i+1) - p_beta*( *(d+i+1)+*(d+i+1-1));
   /*  p2 = d0 */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(p2+i) = *(d+i) - p_alfa*(*(p1+i) + *(p1+i+1));

   for (i=0; i<(buflen>>1); i++){
	   *(buffer+2*i) = *(p1+i+(shift>>1)) ;
	   *(buffer+2*i+1) = *(p2+i+(shift>>1));
   }


	delete []p1;
	delete []p2;
	delete []s;
	delete []d;
}

void CLiftDwt97::filterset(double t)
{
	h = new double [9];
	h1 = new double [7];
	*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);

	double r0, r1, s0, t0;
	r0 = (*(h+4))-2*(*h)*(*(h+3))/(*(h+1));
	r1 = (*(h+2))-(*h)-(*h)*(*(h+3))/(*(h+1));
	s0 = (*(h+3))-(*(h+1))-(*(h+1))*r0/r1;
	t0 = r0-2*r1;

	p_alfa = (*h)/(*(h+1));
	p_beta = (*(h+1))/r1;
	p_gama = r1/s0;
	p_delta = s0/t0;
	p_kesa = t0;


}