wt.cpp

小波变换程序

// WT.cpp: implementation of the CWT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "wt1d.h"
#include "WT.h"
#include "fstream.h"
#include "math.h"

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

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

float CWT::c1[8]={1.5f,1.12f,1.03f,1.01f,1.0f,1.0f,1.0f,1.0f};

float CWT::g[2]={-2.0f,2.0f};

float CWT::h[4]={0.125f,0.375f,0.375f,0.125f};

float CWT::k1[6]={0.0078125f,0.054685f,0.171875f,-0.171875f,-0.054685f,-0.0078125f};

CWT::CWT()
{
	m_data=m_result=NULL;
	m_pn=m_ln=0;
}

CWT::~CWT()
{
	ReleaseData();
}

bool CWT::ReadGrd(CString &infile)
{
	int i,j;
	char head[5]="";
	CString str;
	ReleaseData();
	ifstream ifp(infile);
	if(ifp.fail()){AfxMessageBox("无法打开文件！");return false;}
	ifp>>head;str=head;
	if(str!=_T("DSAA")){AfxMessageBox("文件不是标准的网格化数据文件！");ifp.close();return false;}
	ifp>>m_pn>>m_ln>>m_minx>>m_maxx>>m_miny>>m_maxy>>m_minz>>m_maxz;
	if(ifp.fail()){m_pn=m_ln=0;ifp.close();return false;}
	m_data=(float**)malloc(m_ln*sizeof(float*));
	if(m_data==NULL){m_pn=m_ln=0;ifp.close();return false;}
	m_result=(float**)malloc(m_ln*sizeof(float*));
	if(m_result==NULL){m_pn=m_ln=0;free(m_data);m_data=NULL;ifp.close();return false;}
	for(i=0;i<m_ln;i++){m_data[i]=m_result[i]=NULL;}
	for(i=0;i<m_ln;i++)
	{
		m_data[i]=new float[m_pn];
		m_result[i]=new float[m_pn];
		if(m_data[i]==NULL||m_result[i]==NULL||ifp.fail())
		{
			ReleaseData();
			ifp.close();
			return false;
		}
		for(j=0;j<m_pn;j++){ifp>>m_data[i][j];m_result[i][j]=m_data[i][j];}
	}
	if(ifp.fail())
	{
		ReleaseData();
		ifp.close();
		return false;
	}
	ifp.close();
	return true;
}

bool CWT::WriteGrd(CString &outfile)
{
	int i,j,k;
	ofstream ofp(outfile);
	if(ofp.fail()){AfxMessageBox("无法打开文件！");return false;}
	ofp<<"DSAA"<<endl;
	ofp<<m_pn<<" "<<m_ln<<endl;
	ofp<<m_minx<<" "<<m_maxx<<endl;
	ofp<<m_miny<<" "<<m_maxy<<endl;
	GetLimit(m_result);
	ofp<<m_minz<<" "<<m_maxz<<endl;
	for(i=0,k=1;i<m_ln;i++)
	{
		for(j=0;j<m_pn;j++,k=(k+1)%10){ofp<<m_result[i][j]<<" ";if(k==0)ofp<<endl;}
	}
	ofp.close();
	return true;
}

bool CWT::Read1D(CString &infile)
{
	int j;
	float f;
	FILE* fp=NULL;
	ReleaseData();
	fp=fopen(infile,"r");
	if(fp==NULL){AfxMessageBox("无法打开文件！");return false;}
	m_ln=1;
	while(!feof(fp)){fscanf(fp,"%f",&f);if(!ferror(fp))m_pn++;}
	m_pn-=1;
	fseek(fp,0L,SEEK_SET);
	m_data=(float**)malloc(m_ln*sizeof(float*));
	if(m_data==NULL){ReleaseData();;fclose(fp);return false;}
	m_data[0]=NULL;
	m_result=(float**)malloc(m_ln*sizeof(float*));
	if(m_result==NULL){ReleaseData();fclose(fp);return false;}
	m_result[0]=NULL;
	m_data[0]=new float[m_pn];
	m_result[0]=new float[m_pn];
	if(m_data[0]==NULL||m_result[0]==NULL)
	{
		ReleaseData();
		fclose(fp);
		return false;
	}
	for(j=0;j<m_pn;j++){fscanf(fp,"%f",&m_data[0][j]);m_result[0][j]=m_data[0][j];}
	fclose(fp);
	return true;
}

bool CWT::Write1D(CString &outfile)
{
	int j;
	ofstream ofp(outfile);
	if(ofp.fail()){AfxMessageBox("无法打开文件！");return false;}
	for(j=0;j<m_pn;j++)ofp<<m_result[0][j]<<endl;
	ofp.close();
	return true;
}

void CWT::ReleaseData()
{
	int i;
	if(m_data)
	{
		for(i=0;i<m_ln;i++)if(m_data[i])delete[] m_data[i];
		free(m_data);
		m_data=NULL;
	}
	if(m_result)
	{
		for(i=0;i<m_ln;i++)if(m_result[i])delete[] m_result[i];
		free(m_result);
		m_result=NULL;
	}
	m_pn=m_ln=0;
	m_minx=m_maxx=m_miny=m_maxy=m_minz=m_maxz=0.0f;
}

void CWT::GetLimit(float **data)
{
	int i,j;
	m_minz=m_maxz=0.0f;
	if(data==NULL)return;
	m_minz=m_maxz=data[0][0];
	for(i=0;i<m_ln;i++)
	{
		for(j=0;j<m_pn;j++)
		{
			if(m_minz>data[i][j])m_minz=data[i][j];
			if(m_maxz<data[i][j])m_maxz=data[i][j];
		}
	}
	return;
}

void CWT::DWT1(float *x, int len, int lev, float *a, float *d)
{
	float p,q,c2;
	int i,j,s,c;
	if(x==NULL||a==NULL||d==NULL||lev<1)return;

	c=(int)pow(2,lev-1);
	c2=1.0f/c1[lev-1];
	s=len+2*c+1;
	j=len-1;
	for(i=len;i<=s;i++,j--)x[i]=x[j];
	s=c+c;
	for(i=c/2;i<=s;i++)
	{
		p=q=0.0f;
		if(lev==1)
		{
			for(j=-1;j<=2;j++)p=p+h[j+1]*x[(int)(fabs(i-j+0.5)-0.5)];
			for(j=0;j<=1;j++)q=q+g[j]*x[(int)(fabs(i-j+0.5)-0.5)];
		}
		else
		{
			for(j=-1;j<=2;j++)p=p+h[j+1]*x[(int)(fabs(i-c*j))];
			for(j=0;j<=1;j++)q=q+g[j]*x[(int)(fabs(i-c*j))];
		}
		a[i-c/2]=p;
		d[i-c/2]=-c2*q;
	}
	s=len+c/2;
	for(i=2*c+1;i<=s;i++)
	{
		p=0.0f;
		q=0.0f;
		for(j=-1;j<=2;j++)p=p+h[j+1]*x[i-c*j];
		for(j=0;j<=1;j++)q=q+g[j]*x[i-c*j];
		a[i-c/2]=p;
		d[i-c/2]=-c2*q;
	}
}

void CWT::DWT1(int lev,int pc)
{
	float *x=NULL,*a=NULL,*d=NULL;
	int j,k,c,pl,ll;

	if(m_data==NULL)return;
	c=(int)pow(2,lev-1);
	c = c*3+1;
	if(c<11)c = 11;
	pl = m_pn+c;
	ll = m_ln+c;
	c=ll;if(ll<pl)c=pl;
	x = new float[c];
	a = new float[c];
	d = new float[c];
	if(x==NULL||a==NULL||d==NULL)
	{
		if(x)delete[] x;
		if(a)delete[] a;
		if(d)delete[] d;
		return;
	}
	for(j=0;j<c;j++)x[j]=a[j]=d[j]=0.0f;
	for(k=1;k<=lev;k++)
	{
		for(j=0;j<m_pn;j++)x[j]=m_data[0][j];
		DWT1(x,m_pn,k,a,d);
		for(j=0;j<m_pn;j++)m_data[0][j]=a[j];
	}
	if(pc==0)for(j=0;j<m_pn;j++)m_result[0][j]=a[j];
	else for(j=0;j<m_pn;j++)m_result[0][j]=d[j];
	return;
}

void CWT::DWT2(int lev,int pc)
{
	float *x=NULL,*a=NULL,*d=NULL;
	int i,j,k,c,pl,ll;

	if(m_data==NULL)return;
	c=(int)pow(2,lev-1);
	c = c*3+1;
	if(c<11)c = 11;
	pl = m_pn+c;
	ll = m_ln+c;
	c=ll;if(ll<pl)c=pl;
	x = new float[c];
	a = new float[c];
	d = new float[c];
	if(x==NULL||a==NULL||d==NULL)
	{
		if(x)delete[] x;
		if(a)delete[] a;
		if(d)delete[] d;
		return;
	}
	for(i=0;i<c;i++)x[i]=a[i]=d[i]=0.0f;
	for(k=1;k<lev;k++)
	{
		for(i=0;i<m_ln;i++)
		{
			for(j=0;j<m_pn;j++)x[j]=m_data[i][j];
			DWT1(x,m_pn,k,a,d);
			for(j=0;j<m_pn;j++)m_data[i][j]=a[j];
		}
		for(i=0;i<m_pn;i++)
		{
			for(j=0;j<m_ln;j++)x[j]=m_data[j][i];
			DWT1(x,m_ln,k,a,d);
			for(j=0;j<m_ln;j++)m_data[j][i]=a[j];
		}
	}

	for(i=0;i<m_ln;i++)
	{
		for(j=0;j<m_pn;j++)x[j]=m_data[i][j];
		DWT1(x,m_pn,k,a,d);
		for(j=0;j<m_pn;j++)m_data[i][j]=a[j];
		for(j=0;j<m_pn;j++)m_result[i][j]=d[j];
	}
	for(i=0;i<m_pn&&(pc==0||pc==1);i++)
	{
		for(j=0;j<m_ln;j++)x[j]=m_data[j][i];
		DWT1(x,m_ln,k,a,d);
		for(j=0;j<m_ln;j++)m_data[j][i]=a[j];
		if(pc==0)for(j=0;j<m_ln;j++)m_result[j][i]=a[j];
		if(pc==1)for(j=0;j<m_ln;j++)m_result[j][i]=d[j];
	}
	for(i=0;i<m_pn&&(pc==2||pc==3);i++)
	{
		for(j=0;j<m_ln;j++)x[j]=m_result[j][i];
		DWT1(x,m_ln,k,a,d);
		if(pc==2)for(j=0;j<m_ln;j++)m_result[j][i]=a[j];
		if(pc==3)for(j=0;j<m_ln;j++)m_result[j][i]=d[j];
	}
	return;
}

void CWT::WT2(CString& infile,CString &outfile,int lev,int pc,int rec)
{
	int i,j;
	if(!ReadGrd(infile))return;
	DWT2(lev,pc);
	if(rec)
	{
		for(i=0;i<m_ln;i++)for(j=0;j<m_pn;j++)m_data[i][j]=m_result[i][j];
		IDWT2(lev,pc);
	}
	if(!WriteGrd(outfile))return;
	return;
}

void CWT::IDWT1(float *a,float *d,int len,int lev,float *x)
{
	int i,j,s,c=(int)pow(2,lev-1);
	float p,q,c2;
	c2=c1[lev-1];
	s=len+c+c+c;
	for(i=len,j=len-1;i<=s;i++,j--)
	{
		d[i]=-d[j];
		a[i]=a[j];
	}
	s=s-len;
	for(i=-c/2;i<=s;i++)
	{
		p=q=0.0f;
		for(j=-1;j<=2;j++)p=p+h[j+1]*a[(int)(fabs(i+c*j))];
		for(j=-3;j<=2;j++)
			if((i-c*j)<0)q=q-k1[j+3]*d[(int)(fabs(i-c*j))];
			else q=q+k1[j+3]*d[i-c*j];
		x[i+c/2]=a[i+c/2]=p-c2*q;
	}
	s=len-c/2;
	for(i=3*c+1;i<=s;i++)
	{
		p=q=0.0f;
		for(j=-1;j<=2;j++)p=p+h[j+1]*a[i+c*j];
		for(j=-3;j<=2;j++)q=q+k1[j+3]*d[i-c*j];
		x[i+c/2]=a[i+c/2]=p-c2*q;
	}
}

void CWT::IDWT1(int lev,int pc)
{
	float *x=NULL,*a=NULL,*d=NULL;
	int j,k,c,pl,ll;

	if(m_data==NULL)return;
	c=(int)pow(2,lev-1);
	c = c*3+1;
	if(c<11)c = 11;
	pl = m_pn+c;
	ll = m_ln+c;
	c=ll;if(ll<pl)c=pl;
	x = new float[c];
	a = new float[c];
	d = new float[c];
	if(x==NULL||a==NULL||d==NULL)
	{
		if(x)delete[] x;
		if(a)delete[] a;
		if(d)delete[] d;
		return;
	}
	for(j=0;j<c;j++)x[j]=a[j]=d[j]=0.0f;
	for(j=0;j<m_pn;j++)
	{
		if(pc==0){a[j]=m_data[0][j];d[j]=0.0f;}
		else {a[j]=0.0f;d[j]=m_data[0][j];}
	}
	IDWT1(a,d,m_pn,lev,x);
	for(j=0;j<m_pn;j++)m_data[0][j]=x[j];

	for(k=lev-1;k>0;k--)
	{
		for(j=0;j<m_pn;j++){a[j]=m_data[0][j];d[j]=0.0f;}
		IDWT1(a,d,m_pn,k,x);
		for(j=0;j<m_pn;j++)m_data[0][j]=x[j];
	}
	for(j=0;j<m_pn;j++)m_result[0][j]=m_data[0][j];
	return;
}

void CWT::IDWT2(int lev,int pc)
{
	float *x=NULL,*a=NULL,*d=NULL;
	int i,j,k,c,pl,ll;

	if(m_data==NULL)return;
	c=(int)pow(2,lev-1);
	c = c*3+1;
	if(c<11)c = 11;
	pl = m_pn+c;
	ll = m_ln+c;
	c=ll;if(ll<pl)c=pl;
	x = new float[c];
	a = new float[c];
	d = new float[c];
	if(x==NULL||a==NULL||d==NULL)
	{
		if(x)delete[] x;
		if(a)delete[] a;
		if(d)delete[] d;
		return;
	}
	for(i=0;i<c;i++)x[i]=a[i]=d[i]=0.0f;
	for(i=0;i<m_pn;i++)
	{
		for(j=0;j<m_ln;j++)
		{
			if(pc==0){a[j]=m_data[j][i];d[j]=0.0f;}
			else {a[j]=0.0f;d[j]=m_data[j][i];}
		}
		IDWT1(a,d,m_ln,lev,x);
		for(j=0;j<m_ln;j++)m_data[j][i]=x[j];
	}
	for(i=0;i<m_ln;i++)
	{
		for(j=0;j<m_pn;j++)
		{
			if(pc==0){a[j]=m_data[i][j];d[j]=0.0f;}
			else {a[j]=0.0f;d[j]=m_data[i][j];}
		}
		IDWT1(a,d,m_pn,lev,x);
		for(j=0;j<m_pn;j++)m_data[i][j]=x[j];
	}

	for(k=lev-1;k>0;k--)
	{
		for(i=0;i<m_pn;i++)
		{
			for(j=0;j<m_ln;j++){a[j]=m_data[j][i];d[j]=0.0f;}
			IDWT1(a,d,m_ln,k,x);
			for(j=0;j<m_ln;j++)m_data[j][i]=x[j];
		}
		for(i=0;i<m_ln;i++)
		{
			for(j=0;j<m_pn;j++){a[j]=m_data[i][j];d[j]=0.0f;}
			IDWT1(a,d,m_pn,k,x);
			for(j=0;j<m_pn;j++)m_data[i][j]=x[j];
		}
	}
	for(i=0;i<m_ln;i++)for(j=0;j<m_pn;j++)m_result[i][j]=m_data[i][j];

	return;
}

void CWT::WT1(CString& infile,CString &outfile,int lev,int pc,int rec)
{
	int j;
	if(!Read1D(infile))return;
	DWT1(lev,pc);
	if(rec)
	{
		for(j=0;j<m_pn;j++)m_data[0][j]=m_result[0][j];
		IDWT1(lev,pc);
	}
	if(!Write1D(outfile))return;
	return;
}