Wavelet.cpp

Go to the documentation of this file.

#include "Wavelet.h"

WaveletBase::WaveletBase( double *values_in, unsigned int size_in )
{
  std::cout << "WaveletBase constructor with size_in " << size_in << std::endl;
  std::cout << "values_in[0]= " << values_in[0] << " " << values_in[size_in-1] << std::endl;

  size = getLevels( size_in ); //, nb_level );
  std::cout << " size " << size << " max_level " << max_level << std::endl;

  unsigned int discard_data = size_in - size;
  //need to make the copy of data
  values = new double[ size ];
  //copy discarding data and keeping chronological order
  for ( unsigned int i=0; i< size; i++) {
      // get chronological order ( from latest to new of only the first size data points.
      // inversion of order done in GetCArray
      values[i] = values_in [ i + discard_data ];
  }
  print_values( true ); //option see only first and last
}

unsigned int WaveletBase::getLevels( unsigned int size_in ) //, unsigned int &nb_level )
{
  unsigned int max_size = 0;

  unsigned int v; // compute the next highest power of 2 of 32-bit v, to check with 64-bit just add >> 64
  v = size_in;
  //v--; apply only if not power of 2

  //test first if is a power of 2
  if ( this->isPowerOfTwo( size_in ) == false ) {
      std::cout << "size_in is not a power of 2!" << std::endl;
      v--;

  } else {
      std::cout << " ok power of 2 " << std::endl;
  }

  // http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogObvious
  /*may work but very slow certainly
  while ( isPowerOfTwo ( max_size ) == false ) {
      max_size = max_size - 1
  }*/

  //should be only if not power of 2
  //unsigned int v; // compute the next highest power of 2 of 32-bit v, to check with 64-bit just add >> 64
  //v = size_in;
  //v--;
  v |= v >> 1;
  v |= v >> 2;
  v |= v >> 4;
  v |= v >> 8;
  v |= v >> 16;
  v++;

  std::cout << " size_in " << size_in << " v " << v << " " << v / 2 << std::endl;
  max_size = v / 2;

  // currently signals up to a length of 2^20 supported
  for( unsigned int i = 0; i < 20; ++i) {
      if( max_size == (1 << i)) {
          // here set member max_level
          max_level = i;
          break;
      }
  }
  std::cout << "getLvels return max_size " << max_size << " max_level " << max_level << std::endl;
  return max_size;
}
//
void WaveletBase::print_values( bool only_first_last )
{
  std::cout << "WaveletBase::print_values option " << only_first_last << std::endl;

  if ( only_first_last ) {
      std::cout << "first i=0 " << values[0] << std::endl;
      std::cout << "first i=1 " << values[1] << std::endl;
      std::cout << "last i=" << size-1 << " " << values[ GetSize()-2 ] << std::endl;
      std::cout << "last i=" << size   << " " << values[ GetSize()-1 ] << std::endl;
  } else {
      for ( unsigned int i = 0; i < size; i ++ )
          std::cout << " i= " << i << values[i] << std::endl;
  }
  std::cout << "end print_values \n" << std::endl;
}

int Haar::Transform()
{
  std::cout << "\n\t Haar Transform" << std::endl;
  //std::cout << " size " << size << std::endl;

  // last interface send size/2, retrun a[0]
  haar_value = CalcTransform( values, size / 2, 1 );
 
  std::cout << "End Haar Transform nb_level " << nb_level << std::endl;
  return 0;
}

// need new size ??
double* Haar::InvTransform()
{
  // for return
  double * new_values;
  std::cout << "\n\t Haar Inverse Transform" << std::endl;

  // test entry
  if ( map_coeff_filt.empty() == true ) {
      std::cout << "default coeff use all" << std::endl;
      map_coeff_filt = map_coeff;
  } else
    std::cout << "use modified coefficient " << std::endl;

  /*
  if ( h_value == 0. ) {
      std::cout << "default haar_value " << haar_value << std::endl;
      h_value = haar_value;
  }*/

  //std::cout << " size " << size << std::endl;


  // last interface send size/2, retrun a[0]
  double values[1]; // { haar_value; }
  values[0] = haar_value;
               // tmp_value, size, level to reconstruct
  new_values = CalcInvTransform( values, 2, nb_level - 1 );

  std::cout << "End Haar InvTransform " << nb_level << std::endl;
  return new_values;
}


// need tmp_level to fill the map, to add max_level if want to stop
// first written with tmp_values, finally exist in map average but keep if later want to not store average
double Haar::CalcTransform( double *tmp_values, unsigned int tmp_size, unsigned int tmp_level )
{
  std::cout << "HaarCalc size " << size << " tmp_size " << tmp_size << " tmp_level " << tmp_level << std::endl;
  double retValue;
  
  // temporary array, will be copy after in map
  double a[tmp_size];
  double c[tmp_size];
  
  //correct in c++ j?? compile !? ok !
  for ( unsigned int i = 0, j = 0; i < tmp_size*2; i += 2, j++) {
      a[j] = (tmp_values[i] + tmp_values[i+1])/2.;
      c[j] = (tmp_values[i] - tmp_values[i+1])/2.;
  }

  //insert new coefficient and average
  std::vector<double> tmp_coeff, tmp_average;
  tmp_coeff.reserve(tmp_size);
  tmp_average.reserve(tmp_size);

  for ( unsigned int j = 0; j < tmp_size; j++ ) {
      std::cout << c[j] << " ";
      tmp_coeff.push_back( c[j] );
      tmp_average.push_back( a[j] );
  }
  std::cout << " " << std::endl;
  // copy to the map
  map_coeff[tmp_level] = tmp_coeff;
  map_average[tmp_level] = tmp_average;

  // or max_level reached
  if ( tmp_size == 1 ) { //should be size / 2, not size, change pass tmp_size / 2 in argument
      //std::cout << "last loop sizeof_coeff " << vec_coeff.size() << " a[0] " << a[0] << std::endl;
      retValue = a[0];
      // need to set nb_level data member
      nb_level = tmp_level;

  } else
      retValue = CalcTransform( a, tmp_size / 2, tmp_level+1 );

  return retValue;
}

//
double * Haar::CalcInvTransform ( double *tmp_values, unsigned int tmp_size, unsigned int tmp_level )
{
  double *retValue;
  std::cout << "HaarInvTransform  tmp_size:" << tmp_size << " tmp_level:" << tmp_level << std::endl;

  double new_values[tmp_size];

  for ( unsigned int i = 0, j = 0; i < tmp_size/2; i++, j= j+2 ) {
      //A1 = a1 + c1
      new_values[j]   = tmp_values[i] + map_coeff_filt[tmp_level+1][i];
      //A2 = a1 - c1
      new_values[j+1] = tmp_values[i] - map_coeff_filt[tmp_level+1][i];
  }

  for ( unsigned int i = 0; i < tmp_size; i++ ) {
      std::cout << new_values[i] << " ";
  }
  std::cout << " " << std::endl;

  if ( tmp_level == 0 ) {//or tmp_size
      std::cout << "return recursion " << std::endl;
      retValue = new double[tmp_size];
      // copy data to final array
      for ( unsigned int i=0; i< tmp_size; i++ )
          retValue[i] = new_values[i];
      /*
      for ( unsigned int i = 0; i < tmp_size; i++ ) {
            std::cout << retValue[i] << " ";
        }
      */
      return retValue;
      //return 0.;
  } else
    retValue = CalcInvTransform ( new_values, tmp_size*2, tmp_level - 1);
    //return 0.;


  return retValue;
}

//void Haar::Filter( std::string mode, int numb1, int numb2 )
void Haar::Filter( const char * mode, int numb1, int numb2 )
{
  unsigned int tmp;
  std::cout << "\nEntry Filter nb1:" << numb1 << " nb2:" << numb2 << std::endl;
  std::cout << "mode " << mode << " strcmp " << strcmp( mode,"DEL LEVEL") << std::endl;
  //always make copy of the initial coefficient
  map_coeff_filt = map_coeff;

  // want to delete some level
  //if ( mode == "DEL LEVEL") {
  if ( strcmp( mode,"DEL LEVEL")==0 ) {
      std::cout << "delete level from " << numb1 << " to " << numb2 << std::endl;
      for ( unsigned int i = (unsigned int)numb1; i<= (unsigned int)numb2; i++ ) {
          tmp = map_coeff_filt[i].size();
          std::cout << "tmp " << tmp << std::endl;
          for ( unsigned int j = 0; j<tmp; j++ )
              map_coeff_filt[i][j] = 0.;
      }
  }

  return;
}

void Haar::print_coeff()
{
  std::cout << "Print Coefficient size " << size << std::endl;
  std::cout << "Haar value " << haar_value << std::endl;
  std::cout << "max_level " << max_level << std::endl;
  std::cout << "nb_level " << nb_level << "\n" << std::endl;

  // key map begin at 1
  for ( unsigned int i=1; i <= nb_level; i ++ ) {
      std::cout << "Level " << i << std::endl;
      std::cout << "nb_coeff " << map_coeff[i].size() << std::endl;
      std::cout << "nb_average " << map_average[i].size() << std::endl;
      std::cout << "coeff :";
      for ( unsigned int j=0; j < map_coeff[i].size(); j++ ) {
          std::cout << map_coeff[i][j] << " ";
      }
      std::cout << "\naverage :";
      for ( unsigned int j=0; j < map_average[i].size(); j++ ) {
                std::cout << map_average[i][j] << " ";
      }
      std::cout << " " << std::endl;
  }
}

// may add option, levels, 3d-2d...
void Haar::Draw()
{
  std::cout << "Entry Haar::Draw " << std::endl;
}