#define DEBUGing

#include "itkSWaveletTest.h"

#include <algorithm>
#include <fstream>



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

  if (argc != 1) {
    cerr<<"format: itkSWaveletTest"<<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);
//     }

  int n = 4;

  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 5:
  std::cerr<<"testing task 5:"<<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< double > oriF;
  getScalarFn("x10.txt", oriF);
  mySphereMeshSource->SetScalarFunctionWithSphere( "para10.meta", oriF);

  mySphereMeshSource->SWaveletTransform();

  mySphereMeshSource->inverseSWaveletTransform();

  std::vector< double > f =  mySphereMeshSource->GetScalarFunction();
  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;
  std::vector< double >::iterator itf = f.begin();
  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;

  mySphereMeshSource->WriteFinestSubdivisionMeshToMetaFile("hahaha.meta");

  return EXIT_SUCCESS;
}


void getScalarFn(char* scalarFileName, std::vector< double >& scalarFn)
{
  std::ifstream f(scalarFileName);
  double a;

  for (; !f.eof(); ) {
    f>>a;
    scalarFn.push_back(a);
  }
  scalarFn.erase(scalarFn.end()-1, scalarFn.end()); // The last line is a blank line.
} //void getScalarFn()