liftdwt97.cpp

该程序实现各种小波提升分解

// 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;


}