//
//  Olaf Ronneberger
//  Sample Solution for Exercise 6 (Lecture "3D Image Analysis" Summer 2014) 
//  Spherical Harmonics
//
//  compile with:
//  g++ -Wall -O2 -g transform_world.cc -lgsl -lgslcblas -lblitz -o transform_world -lhdf5



#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <complex>
#include <gsl/gsl_sf_legendre.h>

// HDF5
#include "BlitzHDF5Helper.hh"
// Blitz++
#include <blitz/array.h>


int main( int argc, char** argv)
{
  /*-----------------------------------------------------------------------
   *  some constants
   *-----------------------------------------------------------------------*/
  std::complex<double> com_i(0,1);
  int ellMax = 30;
  
  /*-----------------------------------------------------------------------
   *  Load world elevation data
   *-----------------------------------------------------------------------*/
  blitz::Array<double,2> elevation;
  readHDF5toBlitz( "world_elevation.h5", "elevation", elevation);
  int nTheta = elevation.extent(0);
  int nPhi = elevation.extent(1);
  
  blitz::Array<std::complex<double>,2> reconstruct( nTheta, nPhi);
  reconstruct = 0;
  blitz::Array<double,3> allReconstructed( ellMax, nTheta, nPhi);
  

  for( int ell = 0; ell < ellMax; ++ell)
  {
    for( int m = -ell; m <= ell; ++m)
    {
      /*-------------------------------------------------------------------
       *  Compute the SH coefficient for current ell,m (projection on
       *  basis function)  
       *-------------------------------------------------------------------*/
      std::complex<double> coef = 0;
      for( int thetaQuant = 0; thetaQuant < nTheta; ++thetaQuant)
      {
        double theta = double( thetaQuant +0.5)/nTheta * M_PI;

        for( int phiQuant = 0; phiQuant < nPhi; ++phiQuant)
        {
         double phi = double(phiQuant+0.5)/nPhi * 2*M_PI;
         std::complex<double> Y;
          if( m >= 0)
          {
            Y = gsl_sf_legendre_sphPlm(ell,m,cos(theta)) 
                * exp(  m * phi * com_i);
          }
          else
          {
            Y = pow(-1,m) * gsl_sf_legendre_sphPlm(ell,-m,cos(theta))
                * exp(  m * phi * com_i);
          }

          coef += elevation(thetaQuant, phiQuant) 
              * conj( Y) 
              * sin(theta) * (2*M_PI/nPhi) * (M_PI/nTheta);
        }
      }
      std::cout << "l = " << ell << ", m = " << m << ": " << coef << std::endl;
      
      /*-------------------------------------------------------------------
       *  Reconstruct (scale basis function with found coefficient)
       *-------------------------------------------------------------------*/
      for( int thetaQuant = 0; thetaQuant < nTheta; ++thetaQuant)
      {
        double theta = double( thetaQuant +0.5)/nTheta * M_PI;

        for( int phiQuant = 0; phiQuant < nPhi; ++phiQuant)
        {
         double phi = double(phiQuant+0.5)/nPhi * 2*M_PI;

         std::complex<double> Y;
          if( m >= 0)
          {
            Y = gsl_sf_legendre_sphPlm(ell,m,cos(theta)) 
                * exp(  m * phi * com_i);
          }
          else
          {
            Y = pow(-1,m) * gsl_sf_legendre_sphPlm(ell,-m,cos(theta))
                * exp(  m * phi * com_i);
          }

          reconstruct(thetaQuant, phiQuant) += coef * Y;
        }
      }
    } 
    allReconstructed( ell, blitz::Range::all(), blitz::Range::all()) = 
        blitz::real(reconstruct);
    
  }
  
  /*-----------------------------------------------------------------------
   *  save the world
   *-----------------------------------------------------------------------*/
  writeBlitzToHDF5(allReconstructed, "elevation", "smooth_world.h5");

  return 0;
}