#define DEBUGing

#include "itkSWaveletTest.h"

#include <algorithm>
#include <fstream>



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

  if (argc != 6) 
    {
    cerr << "format: itkSWaveletTest levelOfSubdivision" << endl;
    cerr << " reconstructedScalarFunctionFile " << endl;
    cerr << " reconstructedScalarFunctionFileForKWMeshVisu1 " << endl;
    cerr << " reconstructedScalarFunctionFileForKWMeshVisu2 " << endl;
    cerr << " MeshFile.meta " << endl;
    return EXIT_FAILURE;
    }


  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;
    return EXIT_FAILURE;
    }

  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 1:
  std::cerr<<"testing task 1:"<<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();

  double d = 0;
  for ( itf = f.begin(); itf != itfEnd; ++itf, ++itReF )
    {
    d = d>fabs(*itf-*itReF)?d:fabs(*itf-*itReF);
    }

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






  mySphereMeshSource->WriteReconstructedScalarFunction( argv[2] );

  mySphereMeshSource->WriteReconstructedScalarFunctionForKWVisu( argv[3] );

  double a[]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  std::vector<double> va(a, a+12);
  mySphereMeshSource->SetAllCoarsestScalingCoefficients( va );
  mySphereMeshSource->SetWaveletCoefficientAtScale( 1,0, 1 );

  mySphereMeshSource->inverseSWaveletTransform();

  mySphereMeshSource->WriteReconstructedScalarFunctionForKWVisu( argv[4] );
  mySphereMeshSource->WriteMeshToMetaFile( argv[5], n );

  MeshType::Pointer myMesh = mySphereMeshSource->GetOutput();
  //     std::vector< std::vector< double > > coeff = mySphereMeshSource->GetCoefficients();
  //     //    mySphereMeshSource->SetAllCoarsestScalingCoefficients(3);
  //     mySphereMeshSource->SetWaveletCoefficientsAtOneScale(1, 1.2);
  //     std::vector< std::vector< double > > newCoeff = mySphereMeshSource->GetCoefficients();
  //    MeshType::Pointer myMesh = mySphereMeshSource->GetMeshAtLevel(1);
  //    int a = myMesh->GetReferenceCount();


  //     PointType  pt;
  //     std::cerr<<"# of vertics = "<<myMesh->GetNumberOfPoints()<<std::endl;
  //     std::cerr<<"# of faces = "<<myMesh->GetNumberOfCells()<<std::endl;
  //     cerr<<"# of edges = V + F -2 = "<<myMesh->GetNumberOfCells() + myMesh->GetNumberOfPoints() -2 <<endl;

  //     std::cerr<<std::endl;
  //     std::cerr<<std::endl;
  //     std::cerr<<"resulotion is: "<<mySphereMeshSource->GetResolution()<<std::endl;
  //     std::cerr<<"size of scalar function: "<<f.size()<<std::endl;
  //     std::cerr<<"size of reconstructed scalar function: "<<reF.size()<<std::endl;



  //    std::cout << "Testing itk::IcosahedronMeshSource "<< std::endl;

  //     for(unsigned int i=0; i<myMesh->GetNumberOfPoints(); i++) 
  //       {
  //         myMesh->GetPoint(i, &pt);
  //         std::cout << "Point[" << i << "]: " << pt << std::endl;
  //       }

  //     std::cout << "Test End "<< std::endl;

  displayITKmesh(myMesh);

  return EXIT_SUCCESS;
}