wtlib.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* WTLIB: $Revision: 1.3 $ ; $Date: 1998/08/03 21:01:21 $	*/

/*********************** self documentation **********************/
/***************************************************************************
WTLIB - Functions for wavelet transforms

wt_cascade - generate the mother wavelet and scaling function
fhierfromcd - calculates f(x,z) hierarchically from CD (wavelet coefficients)
wto1d - wavelet transform operator, 1D
wto1dset - setup wavelet operators
wt1 - 1D wavelet transform
wtn - n-D wavelet transform

********************************************************************
Function prototypes:
void wt_cascade(WtFilter *wfilt,WtSizes *wtsizes);
void fhierfromcd(float *f,float *cd,
     WtFilter *wfilt,WtSizes *wtsizes,int justDC);
void wto1d(float *cd,int npoints,enum ToCorD tocord,
     WtFilter *wfilt);
void wto1dset(WtFilter *wfilt,WtSizes *wtsizes);
void wt1(float *cd,enum ToCorD tocord,int npoints,
     WtFilter *wfilt,WtSizes *wtsizes,int idim);
void wtn(float *cd,enum ToCorD tocord,
     	WtFilter *wfilt,WtSizes *wtsizes, int dconly);

********************************************************************
wt_cascade:
Input:
WtFilter *wfilt		pointer to wavelet operator (filter)
WtSizes  *wtsizes	pointer to sizes of the wavelet operator

Returns:
wfilt->psi		pointer to mother wavelet
wfilt->phi		pointer to scaling function

fhierfromcd:
Input:
float *cd		pointer to wavelet coefficients
WtFilter *wfilt		pointer to wavelet operator (filters)
WtSizes *wtsizes	pointer to wavelet sizes
int justDC		flag =1 do DC (zero frequency)

Returns:
float *f		pointer to function f(x,z)

wto1d:
Input:
float *cd		pointer to wavelet coefficients, or f(x)
int npoints		size of the input signal
enum ToCorD tocord	=ToC inverse  or =ToD  forward transform
WtFilter *wfilt		pointer to wavelet operator (filters)

Returns:
float *cd		pointer to f(x) or to wavelet coefficients

wto1dset:
Input:
WtFilter *wfilt		pointer to wavelet filter
WtSizes *wtsizes	pointer to wavelet filter sizes

Return:
wfilt->order		order of Daubechies wavelet
wtsizes->sizes		sizes of the wavelet operator

wt1:
Input:
float *cd		pointer to wavelet coeff. or input f(x)
enum ToCorD tocord	=ToC (invers) =ToD (forward) wavelet transform
int npoints		size of signal f(x) or cd
WtFilter *wfilt		pointer to wavelet operator (filters)
WtSizes *wtsizes	pointer to wavelet sizes
int idim		index for dimension

Returns:
float *cd		pointer to input f(x) or to wavelet coefficients

wtn:
Input:
float *cd		pointer to n-D wavelet coeff, or f(x_1,..,x_n)
enum ToCorD tocord	=ToC (inverse) =ToD (forward) wavelet transform
WtFilter *wfilt		pointer to wavelet operator (filters)
WtSizes *wtsizes	pointer to wavelet sizes 
int dconly		keep and transform back form mra dc component only

Returns:
float *cd		pointer to f(x_1,...,x_n) or n-D wavelet coeffs.
********************************************************************
Notes:
wt_cascade:
Generate data sets phi[0..lengthphi-1] or psi[0..lengthphi-1].
phi or psi[i] = \phi or \psi (i>>MaxLevel).

fhierfromcd:
Once we have the wavelet coefficients cd, we could applied a	     
inverse wavelet transform to reconstruct the function f, this     
is the ordinary way to view the reconstructed f. A good alternative
is to make the "MRA display", which is to show f hierarchically.  
In "MRA display", the DC is put on the upper-left corner, and higher
frequency rectanges are on the lower-right corner.		     

wto1d:
1D wavelet transform operator, used to apply multi-dimensional  
wavelet transform, in which data of each dimension is extracted 
and 1D wavelet transform (forward: tocord==ToD; backward: tocord 
==ToC) is applied to this 1-D data. Before using wto1d, we must 
call wto1dset to set up the wavelet filters.			   
							   
wto1dset:
This subroutine initialize the wavelet filters for wavelet
transform operator wto1d. According to the order you choose
in your main (wtsizes->order=4 up to 20), wto1dset set up
the filter to be the Daubechies wavelet of this order/length. 
This subroutine will be need before calling wt1 or wtn. 

wt1:
1D wavelet transform. cd[1..npoints] as the discrete values of a 1D
function f(x) will be replaced by its wavelet transform if tocord
==ToD, or the wavelet coefficients cd[1..npoints] will be replaced
by the inverse wavelet transform, thus to obtain the reconstructed
f. The size of f npoints MUST be an integer power of 2.

wtn:
Given a n-D discrete values cd[...] of a n-D function f(x_1,...,x_n)
n-D wavelet transform will be applied if tocord==ToD; or given
a n-D wavelet coefficients of an n-D function, inverse wavelet
transform will be applied to reconstruct the n-D function f. Must
call wto1dset before calling wtn.
********************************************************************
References:
Daubechies, I., 1988, Orthonormal Bases of Compactly Supported
Wavelets, Communications on Pure or Applied Mathematics, Vol. XLI,
909-996.

Z. Meng & J. Scales, Multi-resolution Analysis with wavelet 	  
transform: An application to 2-D tomography, CWP-223, 1996    

W. Press et al., Numberical Recipes in C (second edition),     
Cambridge University Press, 1988				   
********************************************************************
Author:
CWP: Zhaobo Meng, Colorado School of Mines, Sept 1996
Modified: Carlos E. Theodoro, Colorado School of Mines, Jun 97
***************************************************************************/

/**************** end self doc ********************************/

#include "par.h"
#include "wt.h"

void wt_cascade(WtFilter *wfilt,WtSizes *wtsizes)
/***********************************************************************
wt_cascade - generate the mother wavelet and scaling function
***********************************************************************
Function Prototype:
void wt_cascade(WtFilter *wfilt,WtSizes *wtsizes);
***********************************************************************
Input:
WtFilter *wfilt		pointer to wavelet operator (filter)
WtSizes  *wtsizes	pointer to sizes of the wavelet operator

Returns:
wfilt->psi		pointer to mother wavelet
wfilt->phi		pointer to scaling function
***********************************************************************
Notes:
Generate data sets phi[0..lengthphi-1] or psi[0..lengthphi-1].
phi or psi[i] = \phi or \psi (i>>MaxLevel).
***********************************************************************
Reference:
Daubechies, I., 1988, Orthonormal Bases of Compactly Supported
Wavelets, Communications on Pure or Applied Mathematics, Vol. XLI,
909-996.
***********************************************************************
Authors:
CWP: Zhaobo Meng, Colorado School of Mines
Modified:
CWP: Carlos E. Theodoro, Colorado School of Mines (06/25/97)
     Included options for:
     - different level of resolution for each dimension;
     - transform back the lower level of resolution, only.
***********************************************************************/
{
        int   sig,k,j,minusn,n,ktemp;
        int   nphi;        /*length of the scaling function phi*/
	int   interval;    /*grid length in the coarsest version*/
	float *phiplus;    /*temp. phi*/
	float coe;         /*square root of 2 thing*/
	float sum=0;       /*summing varible*/
	int   lengthphi;   /*length of phi*/
	int   MaxLevel;    /*finest level phi is generated*/
  
	MaxLevel=wtsizes->MaxLevel;
	lengthphi=wtsizes->order<<MaxLevel;

	wfilt->psi=ealloc1float(lengthphi);
	wfilt->phi=ealloc1float(lengthphi);
	phiplus=ealloc1float(lengthphi);
	interval=1<<MaxLevel;
	nphi=wtsizes->order*interval;
	minusn=interval-nphi/2; /*minusn!*/
	for (k=0;k<wtsizes->order;k++) wfilt->phi[k]=wfilt->cc[k+1];
	for (k=wtsizes->order;k<nphi;k++) wfilt->phi[k]=0;
	coe=1;
	for (k=0;k<MaxLevel;k++) coe=coe*sqrt(2.0);
	/*coe=pow(2.0,MaxLevel/2.0);*/
	for (j=1;j<MaxLevel;j++) { /*the -jth level*/
	        for (n=0;n<nphi;n++) { /*calculate the nth point*/
		        phiplus[n]=0;
			for (k=0;k<nphi;k++){ /*summing up*/
			        if ((n-k-k>=0) && (n-k-k<wtsizes->order))
				        phiplus[n]+=wfilt->cc[n-k-k+1]*wfilt->phi[k];
			}/*next k*/
		}/*next n*/
		for (k=0;k<nphi;k++) wfilt->phi[k]=phiplus[k];
	}/*next j*/
	for (n=0;n<nphi;n++) {wfilt->phi[n]*=coe;sum+=wfilt->phi[n];}
	sum/=interval;
	coe=sqrt(2.0);
	for (k=0;k<nphi;k++){
	        wfilt->psi[k]=0;
		sig=-1;
		for (n=0;n<wtsizes->order;n++){
		        ktemp=2*(k+minusn)+(n-1)*interval;
			if ((ktemp>=0) && (ktemp<nphi))
			        wfilt->psi[k]+=sig*wfilt->cc[n+1]*coe*wfilt->phi[ktemp];
			sig=-sig;
		}/*next n*/
	}/*next k*/
	free1float(phiplus);
}/*end of wt_cascade*/

void fhierfromcd(float *f,float *cd,
     WtFilter *wfilt,WtSizes *wtsizes,int justDC)
/******************************************************************************
fhierfromcd - calculates f(x,z) hierarchically from CD (wavelet coefficients)
*******************************************************************************
Function Prototype:
void fhierfromcd(float *f,float *cd,
     WtFilter *wfilt,WtSizes *wtsizes,int justDC);
*******************************************************************************
Input:
float *cd		pointer to wavelet coefficients
WtFilter *wfilt		pointer to wavelet operator (filters)
WtSizes *wtsizes	pointer to wavelet sizes
int justDC		flag =1 do DC (zero frequency)

Returns:
float *f		pointer to function f(x,z)
*******************************************************************************
Notes:
Once we have the wavelet coefficients cd, we could applied a	     
inverse wavelet transform to reconstruct the function f, this     
is the ordinary way to view the reconstructed f. A good alternative
is to make the "MRA display", which is to show f hierarchically.  
In "MRA display", the DC is put on the upper-left corner, and higher
frequency rectanges are on the lower-right corner.		     
*******************************************************************************
Reference:                                                      
Z. Meng & J. Scales, Multi-resolution Analysis with wavelet 	  
transform: An application to 2-D tomography, CWP-223, 1996    
*******************************************************************************
 Author:   Zhaobo Meng, 11/25/95,  Colorado School of Mines         
 Modified: Carlos E. Theodoro, 06/25/97, Colorado School of Mines
           - modified to allow:
	     - different levels of mra for each dimension; and
	     - transform back only the lower level component, if required.
*******************************************************************************/
{
        int   level;	   /* level of decomposition */
	int   iz,jx,i,j; /* work varibles */
	int   StLevel;   /* starting level for MRA */
	int   StL1;      /* starting level for MRA */
	int   StL2;      /* starting level for MRA */
	float *workf;    /* work varible for signal f */
	float *workcd;   /* work varible for wavelet transform of signal f */
	int   n1;        /* signal size in the fastest direction */
	int   n2;        /* signal size in the slowest direction */
	int   sizez1;    /* signal size in the fastest direction */
	int   sizex2;    /* signal size in the slowest direction */
	int   leap;      /* leap for coarse image */
	float sum,leap2;
	long  int ntot;  /* 1-D length of the multi-D signal */

	n1=wtsizes->NPointsn[1];
	n2=wtsizes->NPointsn[2];
	ntot=n2*n1;
	workf=ealloc1float(ntot+1);
	workcd=ealloc1float(ntot+1);

	StL1=wtsizes->Mraleveln[1];
	StL2=wtsizes->Mraleveln[2];
	if (StL1<=StL2) {
	        StLevel=StL1;
	}
	else {
	        StLevel=StL2;
	}

	for (i=0;i<=ntot;i++) f[i]=0;

	for (level=0;level<StLevel;level++){
    
	  /* fprintf(stderr,"level=%d\n",level); */
    
		sizez1=n1>>level;
		sizex2=n2>>level;
		leap=2<<level;
		leap2=leap*leap;
              
		for (i=0;i<ntot+1;i++) workcd[i]=0;
		for (iz=sizez1/2;iz<sizez1;iz++)
		        for (jx=sizex2/2;jx<sizex2;jx++)
			        workcd[iz+jx*n1]=cd[iz+jx*n1]; 
		wtn(workcd,ToC,wfilt,wtsizes,justDC);
		for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
		for (iz=sizez1/2;iz<sizez1;iz++){
		        for (jx=sizex2/2;jx<sizex2;jx++){
			        sum=0;
				for (i=0;i<leap;i++)
				        for (j=0;j<leap;j++)
					        sum+=workf[iz*leap+i-n1+(jx*leap-n2+j)*n1];
				f[iz+jx*n1]=sum/leap2;
			}
		}
	    
		for (i=0;i<ntot+1;i++) workcd[i]=0;
		for (iz=sizez1/2;iz<sizez1;iz++)
		        for (jx=0;jx<sizex2/2;jx++)
			        workcd[iz+jx*n1]=cd[iz+jx*n1];
		wtn(workcd,ToC,wfilt,wtsizes,justDC);
		for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
		for (iz=sizez1/2;iz<sizez1;iz++){
		        for (jx=0;jx<sizex2/2;jx++){
			        sum=0;
				for (i=0;i<leap;i++)
				        for (j=0;j<leap;j++)
					        sum+=workf[iz*leap+i-n1+(jx*leap+j)*n1];
				f[iz+jx*n1]=sum/leap2;
			}
		}

		for (i=0;i<ntot+1;i++) workcd[i]=0;
		for (iz=0;iz<sizez1/2;iz++)
		        for (jx=sizex2/2;jx<sizex2;jx++)
			        workcd[iz+jx*n1]=cd[iz+jx*n1];
		wtn(workcd,ToC,wfilt,wtsizes,justDC);
		for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
		for (iz=0;iz<sizez1/2;iz++){
		        for (jx=sizex2/2;jx<sizex2;jx++){
			        sum=0;
				for (i=0;i<leap;i++)
				        for (j=0;j<leap;j++)
					        sum+=workf[iz*leap+i+(jx*leap+j-n2)*n1];
				f[iz+jx*n1]=sum/leap2;
			}
		}
	}

	if (StL1==StL2) {
	        level=StLevel-1;

		sizez1=n1>>level;
		sizex2=n2>>level;
		leap=2<<level;
		leap2=leap*leap;
              
		/* fprintf(stderr,"level()=%d\n",level); */

		for (i=0;i<ntot+1;i++) workcd[i]=0;
		for (iz=0;iz<sizez1/2;iz++)
		        for (jx=0;jx<sizex2/2;jx++) 
			        workcd[iz+jx*n1]=cd[iz+jx*n1];
		wtn(workcd,ToC,wfilt,wtsizes,justDC);
		for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
		for (iz=0;iz<sizez1/2;iz++){
		        for (jx=0;jx<sizex2/2;jx++){
			        sum=0;
				for (i=0;i<leap;i++)
				        for (j=0;j<leap;j++)
					        sum+=workf[iz*leap+i+(jx*leap+j)*n1];
				f[iz+jx*n1]=sum/leap2;
			}
		}    
	}

	if (StL1>StL2) {
	        for (level=StLevel;level<StL1;level++){
    
		  /* fprintf(stderr,"level=%d\n",level); */

			sizez1=n1>>level;
			sizex2=n2>>StL2;
			leap=2<<level;
			leap2=leap*leap;
              
			for (i=0;i<ntot+1;i++) workcd[i]=0;
			for (iz=sizez1/2;iz<sizez1;iz++)
			        for (jx=0;jx<sizex2;jx++)
				        workcd[iz+jx*n1]=cd[iz+jx*n1]; 
			wtn(workcd,ToC,wfilt,wtsizes,justDC);
			for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
			for (iz=sizez1/2;iz<sizez1;iz++){
			        for (jx=0;jx<sizex2;jx++){
				        sum=0;
					for (i=0;i<leap;i++)
					        sum+=workf[iz*leap+i-n1+jx*n1];
					f[iz+jx*n1]=sum/leap;
				}
			}