#define DEBUGing

#include "itkSWaveletTest.h"

#include <algorithm>
#include <fstream>


int main( int argc, char * argv [] )
{
  if (argc != 3) {
    cerr<<"format: itkSWaveletTest levelOfWaveletFn levelOfSubdivision"<<endl;
    exit(-1);
  }

  SphereMeshSourceType::Pointer  mySphereMeshSource = SphereMeshSourceType::New();

  PointType center; 
  center.Fill( 0.0 );

  VectorType scale;
  scale.Fill( 1.0 );
  
  int n = atoi(argv[2]);
  int w = atoi(argv[1]);
  if ( n>=7 || n<0 || w<=0 || w>n )
    {
      std::cerr<<"Finest resolution should be from 0 to 6. Level of wavelet function should be greater or equal to 1 and less or equal than finest resolution."<<std::endl;
      exit(-1);
    }

  mySphereMeshSource->SetCenter( center );
  mySphereMeshSource->SetResolution( n );
  mySphereMeshSource->SetScale( scale );

  mySphereMeshSource->Modified();

  try
    {
      mySphereMeshSource->Update();
    }
  catch( itk::ExceptionObject & excp )
    {
      std::cerr << "Error during Update() " << std::endl;
      std::cerr << excp << std::endl;
    }


  // testing task 4:
  std::cerr<<"testing task 4:"<<std::endl; 
  std::cerr<<"    Visualize the Wavelet functions at different levels on the n-th resolution subdivision."<<std::endl;
  std::cerr<<"    the n-th subdivision mesh is written in finestSubdibision.meta."<<std::endl;
  std::cerr<<"    the scalar wavelet function is written in WaveletFnKWVisu.txt."<<std::endl;

  std::vector< PointType > v = mySphereMeshSource->GetVerts();
  std::vector< PointType >::const_iterator itv = v.begin();
  std::vector< PointType >::const_iterator itvEnd = v.end();

  std::vector< double > f( v.size() );
  std::vector< double >::iterator itf = f.begin();

  for ( ; itv != itvEnd; ++itv, ++itf ) {
    *itf = (*itv)[2];
  }

  mySphereMeshSource->SetScalarFunction( f );

  mySphereMeshSource->SWaveletTransform();

  mySphereMeshSource->SetAllCoarsestScalingCoefficients(0);
  mySphereMeshSource->SetAllWaveletCoefficientsOfAllScales(0);
  mySphereMeshSource->SetWaveletCoefficientAtScale(w, 0, 1);

  mySphereMeshSource->inverseSWaveletTransform();

  mySphereMeshSource->WriteReconstructedScalarFunctionForKWVisu("WaveletFnKWVisu.txt");


  return EXIT_SUCCESS;
}