#define DEBUGing

#include "itkSWaveletTest.h"

#include <algorithm>
#include <fstream>



int main( int argc, char * argv [] )
{

  if (argc != 2) {
    cerr<<"format: itkSWaveletTest 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[1]);
  if (n>=7 || n<0) 
    {
      std::cerr<<"Finest resolution can't be less than 0, and best to be no larger than 6..."<<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 2:
  std::cerr<<"testing task 2:"<<std::endl; 
  std::cerr<<"    step 1.1 we have a scalar function defined on the vertexes of the finest mesh. Here we just use the z-coordinates as the function."<<std::endl;
  std::cerr<<"    step 1.2 Do a wavelet transform of the function and we get wavelet coefficients."<<std::endl;
  std::cerr<<"    step 1.3 Do a wavelet inverse transform from the coefficients to reconstruct the original function."<<std::endl;
  std::cerr<<"    step 1.4 Compare the original functin with the reconstructed function to show the fedelity of the transform."<<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->inverseSWaveletTransform();

  std::vector< double > reF =  mySphereMeshSource->GetReconstructedScalarFunction();

  std::vector< double > dif(f.size());
  std::vector< double >::iterator itDif=dif.begin(), itDifEnd = dif.end(), itReF=reF.begin();
  std::vector< double >::iterator itfEnd = f.end();

  //      std::ofstream file("error.dat");        
  double d = 0;
  for ( itf = f.begin(); itf != itfEnd; ++itf, ++itReF ) {
    d = d>fabs(*itf-*itReF)?d:fabs(*itf-*itReF);
    //        file<<*itf<<"\t\t\t"<<*itReF<<"\t\t\t"<<fabs(*itf-*itReF)<<std::endl;
  }

  std::cerr<<std::endl;
  std::cerr<<"The largest abs error is:..... "<<d<<std::endl;
  std::cerr<<std::endl;

  return EXIT_SUCCESS;
}