/*-------------------------------------------------------
 * Package:     Wavelet Analysis Tool
 * File name:   wgauss.cc
 *-------------------------------------------------------
*/
//_______________________________________________________
// WaveletG is a class for a continious wavelet 
// transform that uses the family of Gaussian wavelets.
// 


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

/* minimum scale value */
#define ZERO 1e-5
/* minimum value of unused relative square */
#define ZEROUNUSED 1e-10

#include "wgauss.h"
#include "gsbmp.h"
#include <stddef.h>
#include <math.h>

double gauss (double x)
{
   return exp (- x * x / 2);
}

/* Gaussian wavelet g1 (x) */
double g1wavelet (double x)
{
   return - x * exp (- x * x / 2);
}

/* Gaussian wavelet g2 (x) */
double g2wavelet (double x)
{
   double xsqr = x * x;
   return (1 - xsqr) * exp (- xsqr / 2);
}

/* Gaussian wavelet g3 (x) */
double g3wavelet (double x)
{
   double xsqr = x * x;
   return (3 - xsqr) * x * exp (- xsqr / 2);
}

/* Gaussian wavelet g4 (x) */
double g4wavelet (double x)
{
   double xsqr = x * x;
   return ((6 - xsqr) * xsqr - 3) * exp (- xsqr / 2);
}

/* Gaussian wavelet g5 (x) */
double g5wavelet (double x)
{
   double xsqr = x * x;
   return ((10 - xsqr) * xsqr - 15) * x * exp (- xsqr / 2);
}

/* Gaussian wavelet g6 (x) */
double g6wavelet (double x)
{
   double xsqr = x * x;
   return (((15 - xsqr) * xsqr - 45) * xsqr + 15) * exp (- xsqr /2);
}


 WaveletG::WaveletG(int ngw)
{
// Constructor. Create wavelet with ngw vanishing moment.

	pWDC = new WSeries(1);
	pWDC->N =1;

	   switch (ngw)
	   {
		 case 1: wavenum = wg1; _TrFunc = g1waveletFast; break;
		 case 2: wavenum = wg2; _TrFunc = g2waveletFast; break;
		 case 3: wavenum = wg3; _TrFunc = g3waveletFast; break;
		 case 4: wavenum = wg4; _TrFunc = g4waveletFast; break;
		 case 5: wavenum = wg5; _TrFunc = g5waveletFast; break;
		 case 6: wavenum = wg6; _TrFunc = g6waveletFast; break;
		default: wavenum = wg2; _TrFunc = g2waveletFast; break;
	   }

	/* Initial values of information tables */

	histTrInfo.Left=0;         // left margin
	histTrInfo.Right=0;        // right margin
	histTrInfo.Delta=0;        // step == (Right - Left) / Bins
	histTrInfo.Bins=0;         // number of samples
	histTrInfo.Points=0;       // == Bins + 1

	shiftTrInfo.Left=0;         // left margin
	shiftTrInfo.Right=0;        // right margin
	shiftTrInfo.Delta=0;        // step == (Right - Left) / Bins
	shiftTrInfo.Bins=0;         // number of samples
	shiftTrInfo.Points=0;       // == Bins + 1

	scaleTrInfo.Left=0;         // left margin
	scaleTrInfo.Right=0;        // right margin
	scaleTrInfo.Delta=0;        // step == (Right - Left) / Bins
	scaleTrInfo.Bins=0;         // number of samples
	scaleTrInfo.Points=0;       // == Bins + 1
	
	/* Initial value of auxiliary table */

	_TrTable.Min=0;           // minimum value of index
	_TrTable.Max=0;           // maximum number of index
	_TrTable.Size=0;          // range of it
	_TrTable.Rate=0;          // how many shift steps per sample
	_TrTable.Hist=NULL;       // source data samples
	_TrTable.Wave=NULL;       // table of wavelet values
	
	/* Used relative square of the wavelet */
	RelSquare = 0.95;
	/* Half-width of transformation window */
	halfwindow = 1;

	scaleN=0;
	scale0=scale1=0.;
}

 WaveletG::~WaveletG()
{
// Dectructor.
//
	delete pWDC;
}

/* Clears information table */
TR_INFO& WaveletG::clearTrInfo (TR_INFO& ti)
{
// Clears information table.
//
	ti.Left = 0;
	ti.Right = 0;
	ti.Delta = 0;
	ti.Bins = 0;
	ti.Points = 0;

	return ti;
}

/* Creates information table for the histogram of samples */
TR_INFO& WaveletG::buildHistTrInfo()
{
// Creates information table for the histogram of samples.
//
	if (!histTrInfo.Bins) histTrInfo.Bins = histTrInfo.Points - 1;
	else histTrInfo.Points = histTrInfo.Bins + 1;
	histTrInfo.Delta = (histTrInfo.Right - histTrInfo.Left) / histTrInfo.Bins;
	
	return histTrInfo;
}

/* Builds the table for shift selecting */
TR_INFO& WaveletG::buildShiftTrInfo()
{
//  Builds the table for shift selecting.
//
	if (!shiftTrInfo.Bins) shiftTrInfo.Bins = shiftTrInfo.Points - 1;
	else shiftTrInfo.Points = shiftTrInfo.Bins + 1;
	shiftTrInfo.Right = shiftTrInfo.Left + shiftTrInfo.Delta * shiftTrInfo.Bins;
	
	return shiftTrInfo;
}

/* Native algorithm. Computes single wavelet coefficient. */
 wavereal WaveletG::nativeTr (wavereal shift, wavereal scale)
{
// "Native" algorithm. Computes single wavelet coefficient.
//
	/* return value */
	wavereal ret = 0;
	/* current point */
	wavereal x = histTrInfo.Left;
	/* step */
	wavereal width = histTrInfo.Delta;
	/* for all data samples */
	for (unsigned int bin = 0; bin < histTrInfo.Bins; bin++, x += width)
	   ret += (_TrHist[bin] * (_TrFunc ((x - shift) / scale) -
				_TrFunc ((x + width - shift) / scale)));
	return scale * ret;
}

/* Conventional algorithm. Computes single wavelet coefficient. */
 wavereal WaveletG::convTr (wavereal shift, wavereal scale)
{
//  Conventional algorithm. Computes single wavelet coefficient.
//
	wavereal beta = (histTrInfo.Left - shift) / scale;
	/* Return value.                                */
	/* The first and the last samples are involved. */
	wavereal ret = _TrHist [0] * _TrFunc (beta) -
			_TrHist [histTrInfo.Bins - 1] *
			_TrFunc ((histTrInfo.Right - shift) / scale);
	
	wavereal delta = histTrInfo.Delta / scale;
	/* starting point */
	wavereal x = beta + delta;
	/* for all other cells of the histogram */
	for (unsigned int bin = 0; bin < histTrInfo.Bins - 1; bin++, x += delta)
	{
	   ret += (_TrHist [bin + 1] - _TrHist [bin]) * _TrFunc (x);
	}
	return scale * ret;
}

/* Fast algorithm. Computes single wavelet coefficient. */
 wavereal WaveletG::fastTr (int shiftj, wavereal scale)
{
//  Fast algorithm. Computes single wavelet coefficient.
//
	/* return value */
	wavereal ret = 0;
	
	/* starting point in the auxiliary table */
	int index = -shiftj - _TrTable.Min;
	/* compute the convolution */
	for (unsigned int c = 0; c < histTrInfo.Points; c++, index += _TrTable.Rate)
	   ret += _TrTable.Hist[c] * _TrTable.Wave[index];
	
	return scale * ret;
}

/* Ultra method. Computes single wavelet coefficient. */
 wavereal WaveletG::ultrafastTr (int shiftj, wavereal scale)
{
//  Ultra method. Computes single wavelet coefficient.
//
	/* return value */
	wavereal ret = 0;
	
	/* first index in the table of values of the wavelet */
	int k = (int)(bxdelta + shiftj / _TrTable.Rate);
	
	/* selecting boundaries according to window width */
	int binmin = k - halfwindow;
	int binmax = k + halfwindow;
	if (binmin < 0) binmin = 0;
	if (binmax > (int)histTrInfo.Points) binmax = histTrInfo.Points;

	int index = binmin * _TrTable.Rate - shiftj - _TrTable.Min;
	
	/* evaluating the convolution */
	for (int c = binmin; c < binmax; c++, index += _TrTable.Rate)
	   ret += _TrTable.Hist[c] * _TrTable.Wave[index];
	
	return scale * ret;
}

/* Native algorithm for evaluating the line in wavelet spectrum. */
 void WaveletG::nativeTr (wavereal* result, wavereal scale)
{
//  Native algorithm. Computes a  line in wavelet spectrum.
//
	/* initial value of shift */
	wavereal shift = shiftTrInfo.Left;
	/* computing single coefficients at all shifts */
	for (unsigned int point = 0; point < shiftTrInfo.Points;
		point++, shift += shiftTrInfo.Delta)
	   result[point] = nativeTr (shift, scale);
}
	
/* Conventional algorithm for evaluating the line in wavelet spectrum. */
 void WaveletG::convTr (wavereal* result, wavereal scale)
{
//  Conventional algorithm. Computes a line in wavelet spectrum.
//
	/* starting (minimum) shift */
	wavereal shift = shiftTrInfo.Left;
	/* run over all shifts */
	for (unsigned int point = 0; point < shiftTrInfo.Points;
		point++, shift += shiftTrInfo.Delta)
	result[point] = convTr (shift, scale);
}

/* Fast algorithm for evaluating the line in wavelet spectrum. */
 void WaveletG::fastTr (wavereal* result, wavereal scale)
{
//  Fast algorithm. Computes a the line in wavelet spectrum.
//
	/* just go through all points in the array */
	for (unsigned int point = 0; point < shiftTrInfo.Points; point++)
	   result[point] = fastTr (point, scale);
}

/* Ultra algorithm for evaluating the line in wavelet spectrum. */
 void WaveletG::ultrafastTr (wavereal* result, wavereal scale)
{
// Ultra algorithm for evaluating the line in wavelet spectrum.
//
	/* run over all shifts */
	for (unsigned int point = 0; point < shiftTrInfo.Points; point++)
	   result[point] = ultrafastTr (point, scale);
}

/* Builds the transformation tables for given parameters */
 int WaveletG::buildTrTableFast (wavereal scale, wavereal rate)
{
// Builds the transformation tables for given parameters.
//
	_TrTable.Rate = (int)rate;

	/* filling the table of differences */
	_TrTable.Hist = new wavereal[histTrInfo.Points];
	if (!_TrTable.Hist) return 0;
	_TrTable.Hist[0] = _TrHist[0];
	_TrTable.Hist[histTrInfo.Bins] = -_TrHist[histTrInfo.Bins - 1];
	for (unsigned bin = 1; bin < histTrInfo.Bins; bin++)
	   _TrTable.Hist[bin] = _TrHist[bin] - _TrHist[bin - 1];
	
	_TrTable.Min = -shiftTrInfo.Bins;
	_TrTable.Max = (int)(histTrInfo.Bins * rate);
	_TrTable.Size = _TrTable.Max - _TrTable.Min + 1;
	_TrTable.Wave = new wavereal [_TrTable.Size];
	if (!_TrTable.Wave)
	{
	   delete _TrTable.Hist;
	   return 0;
	}
	
	wavereal coef = shiftTrInfo.Delta / scale;
	wavereal add = (histTrInfo.Left - shiftTrInfo.Left) / scale;
	
	/* creating the table of values of the wavelet */
	int id = _TrTable.Min;
//cout<<"id="<<id<<endl;
	for (int index = 0; index < (int)_TrTable.Size; index++)
	{
	   wavereal x = coef * (index + id) + add;
	   _TrTable.Wave[index] = _TrFunc (x);
	}
	
	return 1;
}
 void WaveletG::t2w(WaveData &t, int _scaleN, wavereal _scale0, wavereal _scale1, int method, int norm, wavereal unused)
{
// Time domain to Wavelet domain transform.
// Given parameters t - data for transformation.
// scale - Width of analizing wavelet.
// method - should be 0,1,2,or 3. 
// 0 - Ultra method,
// 1 - Fast method,
// 2 - Conventional method, 
// 3 - Native method.
// norm - should be equal 1, if normalization is needed.
// unused - for Ultra method. Shows what relative square of wavelet you don't want to use.
//
   scaleN=_scaleN;
   scale0=_scale0;
   scale1=_scale1;

   if(method<0 || method>3) method=0;
   if (unused < ZEROUNUSED) unused = ZEROUNUSED;
   if (unused > 1) unused = 1.0 - ZEROUNUSED;
   RelSquare = 1 - unused;
   int size=t.N;
   setTrHist (t.data);
   wavereal rate=1.;
   wavereal shiftmin=0.;
	
   histTrInfo.Left = 0;
   histTrInfo.Right = size;
   histTrInfo.Bins = size;
   buildHistTrInfo();

   shiftTrInfo.Left = shiftmin;
   shiftTrInfo.Delta = histTrInfo.Delta / rate;
   shiftTrInfo.Points = (int) (size * rate);	
   buildShiftTrInfo();
   
   int wpoints=(int)(shiftTrInfo.Points)*scaleN;
   if(pWDC->N!=wpoints) pWDC->Resize(wpoints);

   wavereal *wdata;
   wavereal scale;

   for(int iscale=0;iscale<scaleN;iscale++)
   {
      wdata=pWDC->data+shiftTrInfo.Points*iscale;
      scale=scale0+iscale*(scale1-scale0)/scaleN;

	/* Select which method to use */
	if (method == 3)
	{
	   /* Native */
	   /* the transformation itself */
	   nativeTr(wdata, scale);
	}

	if (method == 2)
	{
	   /* Conventional */

	   /* transform */
	   convTr(wdata, scale);
	}

	if (method == 1)
	{
	   /* Fast */

	   /* build auxilary tables */
	   buildTrTableFast (scale, rate);

	   /* transform */
	   fastTr (wdata, scale);
	   destroyTrTable();
	}

	if (method == 0)
	{
	   /* Ultra */

	   buildTrTableFast (scale, rate);

	   /* select window width according to given unused relative square */
	   bxdelta = (shiftTrInfo.Left - histTrInfo.Left) / histTrInfo.Delta;
	   wavereal halfv = scale * sqrt (-2 * log (1 - RelSquare));
	   halfwindow = (int) ((wavereal) ((wavereal) halfv / histTrInfo.Delta) + 1.0);

	   /* transform */
	   ultrafastTr (wdata, scale);	
	   destroyTrTable();
	}

	/* normalize if necessary */
	if (norm)
	{
	   wavereal wavecoef=g2coefsq;
	   switch (wavenum)
	   {
	      case wg1: wavecoef = g1coefsq; break;
	      case wg2: wavecoef = g2coefsq; break;
	      case wg3: wavecoef = g3coefsq; break;
	      case wg4: wavecoef = g4coefsq; break;
	      case wg5: wavecoef = g5coefsq; break;
	      case wg6: wavecoef = g6coefsq; break;
	      case wempty: break;
	   }
	   wavereal ncoef = 1.0 / (wavecoef * sqrt (fabs (scale)));
	   for (unsigned int c = 0; c < shiftTrInfo.Points; c++)
		 wdata[c]*=ncoef;
	}
   }
}

 void WaveletG::DumpWTAsBitmap(const char* fn,int palette)
{
   if(scaleN==0){ cout<<"DumpWTAsBitmap(const char*,int=0) warning: Nothing to do.n"; return; }

   unsigned char *bits=new unsigned char[pWDC->N];
   
   wavereal mx,mn;
   mx=mn=pWDC->data[0];
   for(long i=1;i<pWDC->N;i++)
   {
      if(pWDC->data[i]<mn)mn=pWDC->data[i];
      if(pWDC->data[i]>mx)mx=pWDC->data[i];
   }
   mx=255./(mx-mn);

   for(long i=0;i<pWDC->N;i++)bits[i]=(unsigned char)(mx*(pWDC->data[i]-mn));

   gs_bm_write(fn,bits,shiftTrInfo.Points,scaleN,palette);

   delete bits;

}

#ifdef _USE_ROOT
  ClassImp(WaveletG)
#endif