/*-------------------------------------------------------
 * Package:     Wavelet Analysis Tool
 * File name:   wdwt.cc
 *-------------------------------------------------------
*/
//_______________________________________________________
// WaveDWT is the abstract base class for discrete wavelet
// transform  methods. It includes methods, like getLayer,
// putLayer, t2w, w2t, etc.

#ifndef _STREAM_H
#include <stream.h>
#define _STREAM_H
#endif

#include "wavedwt.h"
#include "gsbmp.h"

#ifdef _USE_ROOT
  ClassImp(WaveDWT)
#endif

WaveDWT::WaveDWT(int m, int ist, int brd)
{
// Constructor which creates object of class WSeries
// with N=m elements of data. Variable 'border'
// controls processing of data array boundaries.
	pWDC = new WSeries(m);
	pWDC->N = m; 

	IsTree=ist;
	Level=0;
	MaxLevel=0;
	border=brd;
}

 WaveDWT::~WaveDWT()
{
// Destructor.
  delete pWDC;
}

 void WaveDWT::t2w(const WaveData& td,int ldeep)
{
// Virtual procedure for forward wavelet transform that 
// take initial data from WaveData "td" and places result
// in WSeries class.
// Only ldeep steps of decomposition will be done.
// By default ldeep=-1, which means do maximum number of allowed steps.

        register int n=td.N;
        register wavereal *in,*wd;
        in = td.data;
	Level=0;

	if(td.N!=pWDC->N) 
	{
		pWDC->Resize(td.N);
	}
	pWDC->Rate=td.Rate;
        wd=pWDC->data;
      
	for(int q=0;q<n;q++,in++,wd++) *wd=*in;

	if (ldeep!=0) t2w(ldeep);
}	

 void WaveDWT::t2w(int ldeep)
{
// Virtual procedure for forward wavelet transform. 
// Real code will appear in
// derivative classes since it depends of the type of wavelet.
// Only ldeep steps of decomposition will be done.
// By default ldeep=1, which means do one step.

	getMaxLevel();

	int levs=Level;
	int levf=Level+ldeep;
	if((ldeep==-1)||(levf>MaxLevel)) levf=MaxLevel;

	int layf;
	for(int level=levs;level<levf;level++)
	{
		layf=IsTree?1<<level:1;
		for(int layer=0;layer<(layf);layer++)
		{
			decompose(level,layer);
		}

		Level=level+1;

	}

	Level=levf;
}

 void WaveDWT::w2t(int ldeep)
{
// Virtual procedure for inverse wavelet transform that
// take initial data from current WSeries class without places
// result into object of  WaveData class.
// Only ldeep steps of backward composition will be done.
// By default ldeep=1, which means do one step.

	int levs=Level;
	int levf=Level-ldeep;
	if((ldeep==-1)||(levf<0)) levf=0;

	int layf;
	for(int level=levs-1;level>=levf;level--)
	{
		layf=IsTree?1<<level:1;
		for(int layer=0;layer<(layf);layer++)
		{
			reconstruct(level,layer);
		}

		Level=level;

	}
	Level=levf;
}
 void WaveDWT::w2t(WaveData &td, int ldeep)
{
// Virtual procedure for inverse wavelet transform that
// take initial data from current WSeries class and places
// result to object td of  WaveData class.
// Only ldeep steps of backward composition will be done.
// By default ldeep=-1, which means do maximum allowed steps.

  w2t(ldeep);
  register int n=pWDC->N;
  register wavereal *out, *d;
  d = pWDC->data;
  
  if(n!=td.N) td.Resize(n);
  out = td.data;

  td.Rate=pWDC->Rate;

  for(int q=0;q<n;q++,out++,d++) *out = *d;
//  memcpy(out,d,n*sizeof(out));

}

 void WaveDWT::Denoise(wavereal cut)
{
// Simple denoising procedure. Sets to zero all elements of
// wavelet spectrum which absolute value less then "cut".
 
  int n = pWDC->N;
  if (cut < 0.) cut = -cut;
  cout<<"Denoising array of "<< n << " elements, cut = "<< cut << "n";
  for (int i = 0; i < n; i++) {
    if (fabs(pWDC->data[i]) < cut) pWDC->data[i] = 0.;
  }
}

 void WaveDWT::getLayer(WaveData &td, int k)
{
// Get layer with given number k and put it into td array.
// 

  register int start, step;

  if(k<0) { k=0; cout<<"getLayer() warning: Selected layer number 0.n";}

  if(IsTree) {
    int max_layers=1<<Level;
    if(k>max_layers-1) k=max_layers-1;
	
    int m = (pWDC->N)>>Level;
    if (td.N!=m) td.Resize(m);

    int k1=getIdA(Level,k);
    for(int i=0;i<m;i++) td.data[i]=pWDC->data[(i<<Level)+k1];
  } 
  else {

    bool is_top=false;

    if(k>MaxLevel){
      k=MaxLevel;
      cout<<"getLayer() warning: Selected layer with maximum number "
      <<MaxLevel<<".n";
    }
 
    if(k==0) { is_top=true; k=Level; }
	
    int m = (pWDC->N)>>k;
    if (td.N!=m) td.Resize(m);

    start = (is_top) ? 0: 1<<(k-1);
    step = 1<<k;

    register wavereal *out, *d;
    out=td.data;
    d=pWDC->data+start;

    for(int i=0;i<m;i++,out++,d+=step) *out=*d;
  }
}

 void WaveDWT::putLayer(WaveData &td, int k)
{
// Replace content of the layer number k with data from array td.
//
  register int start, step;

  if(k<0) { k=0;cout<<"Selected layer number 0.n";}

  if(IsTree){
    int max_layers=1<<Level;
    if(k>max_layers-1){
      k=max_layers-1;
      cout<<"putLayer() warning: Selected layer with maximum number "
          <<max_layers-1<<".n";
    }
	
    int m=(pWDC->N)>>Level;
    if(td.N!=m){
      cout<<
      "putLayer() warning: number of data is not correct! Does nothing.n";
      return;
    }
	
    int k1=getIdA(Level,k);
    for(int i=0;i<m;i++) pWDC->data[(i<<Level)+k1]=td.data[i];
  }
  else {
    char is_top=0;

    if(k>MaxLevel) {
      k=MaxLevel;
      cout<<"putLayer() warning: Selected layer with maximum number "
          <<MaxLevel<<".n";
    }

    if(k==0) { is_top=1; k=Level; }
		
    int m=(pWDC->N)>>k;
    if(td.N!=m) {
      cout<<
      "putLayer() warning: number of data is not correct! Does nothing.n";
      return;
    }
	
    start = (is_top) ? 0: 1<<(k-1);
    step = 1<<k;

    register wavereal *in, *d;
    in=td.data;
    d=pWDC->data+start;

    for(int i=0;i<m;i++,in++,d+=step) *d=*in;
  }
}

 void WaveDWT::getFreqLayer(WaveData &td, int k)
{
// The same as getLayer, but select layers in frequency order.

	int k1=getIdL(Level,k);
	getLayer(td,k1);	
}

 void WaveDWT::putFreqLayer(WaveData &td, int k)
{
// The same as putLayer, but put layers in frequency order.

	putLayer(td,getIdL(Level,k));
}

 int WaveDWT::getMaxLevel()
{
// Calculate maximal level of wavelet decomposition.
// It depends on filter length.

	MaxLevel=0;
	int jc1=pWDC->N;//>>1;

	for(;jc1>=nh;jc1=jc1>>1) MaxLevel++;
	return MaxLevel;
}

 int WaveDWT::getIdA(int level,int layer)
{	
// Gets wavelet layer number for given level and layer. 

	int IdA=0;
	for(int i=0;i<level;i++)
	{
		if((layer>>i)&1) IdA+=1<<(level-1-i);	//some operations with bits
	}
	return IdA;
}


 int WaveDWT::getIdF(int level,int IdL)
{
// Get frequency ID number for given level and layer.
//

	int IdF=IdL;
	for(int i=1;i<level;i++)
	{
		int j=((1<<i)&(IdF));
		if(j)IdF=((1<<i)-1)^(IdF);
	}

	return IdF;
}

 int WaveDWT::getIdL(int level,int IdF)
{
// Get the layer number in binary tree for given level and frequency ID number.
//

	int IdL=IdF;
	for(int i=level-1;i>=1;i--)
	{
		int j=((1<<i)&(IdL));
		if(j)IdL=((1<<i)-1)^(IdL);
	}

	return IdL;
}

 int WaveDWT::DumpWTAsBitmap(const char* fn,unsigned int palette)
{
//
   if(IsTree==0) return 1;
   if(pWDC==NULL) return 1;

   unsigned int ni=1<<Level;
   unsigned int nj=pWDC->N/ni;
   unsigned char *bits=new unsigned char[pWDC->N];

   double mx,mn;
   mx=mn=pWDC->data[0];
   for(int s=0;s<pWDC->N;s++){
      if(pWDC->data[s]>mx) mx=pWDC->data[s];
      if(pWDC->data[s]<mn) mn=pWDC->data[s];
   }

  WaveData l;
  for(unsigned int i=0;i<ni;i++)
  {
	getFreqLayer(l,i);	
	for(unsigned int j=0;j<nj;j++) bits[j+nj*i]=(unsigned char)(255*(l.data[j]-mn)/(mx-mn));
  }

  cout<<"mn="<<mn<<" mx="<<mx<<"n";

//   unsigned char *bits=new unsigned char[256*256];
//   for(int i=0;i<256;i++)for(int j=0;j<256;j++)bits[i*256+j]=(i+j)>255?(i+j-255):(i+j);

   gs_bm_write(fn,bits,nj,ni,palette);

   delete [] bits;

   return 0;
}

#ifdef _USE_DMT
 void WaveDWT::t2w(const TSeries& ts,int ldeep)
{
   unsigned int Ninput = ts.getNSample();
   Interval ti = ts.getTStep();
   wavereal rate = ti.GetSecs();
   if(rate) rate = 1./rate;
   else cout<<"warning: TSeries::GetSecs() returned 0.\n";
   if((int)Ninput != pWDC->N) pWDC->Resize(Ninput);
   pWDC->Rate = rate;
   ts.getData(Ninput,(double *)pWDC->data);

   t2w(ldeep);

}

 void WaveDWT::w2t(TSeries& ts,int ldeep)
{
   double rate = 1.;
   w2t(ldeep);
   if(pWDC->Rate > 0.) rate = 1./pWDC->Rate;
   unsigned int Ninput = ts.getNSample();
   if((int)Ninput!=pWDC->N) ts.ReSize(pWDC->N);
   if(!ts.setData(Time(0.0),Interval(rate),pWDC->data,(unsigned int)pWDC->N))
     cout<<"warning: TSeries::setData() returned 0.\n";   
   return;
}
#endif // use DMT