/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2023 by Kay F. Jahnke                    */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/

/// \file roundtrip.cc
///
/// \brief benchmarking and testing code for various vspline capabilities
///
/// load an image, create a b-spline from it, and restore the original data,
/// both by normal evaluation and by convolution with the reconstruction kernel.
/// all of this is done several times each with different boundary conditions,
/// spline degrees and in float and double arithmetic, the processing times
/// and differences between input and restored signal are printed to cout.
///
/// obviously, this is not a useful program, it's to make sure the engine functions as
/// intended and all combinations of float and double as values and coordinates compile
/// and function as intended, also giving an impression of the speed of processing.
///
/// compile:
/// clang++ -std=c++11 -march=native -o roundtrip -O3 -pthread -DUSE_VC roundtrip.cc -lvigraimpex -lVc
/// or: clang++ -std=c++11 -march=native -o roundtrip -O3 -pthread roundtrip.cc -lvigraimpex
///
/// invoke: roundtrip <image file>
///
/// there is no image output.

#include <iostream>

#include <vspline/vspline.h>

#include <vigra/stdimage.hxx>
#include <vigra/imageinfo.hxx>
#include <vigra/impex.hxx>
#include <vigra/accumulator.hxx>
#include <vigra/multi_math.hxx>

#define PRINT_ELAPSED

#ifdef PRINT_ELAPSED
#include <ctime>
#include <chrono>
#endif

using namespace std ;

using namespace vigra ;

/// check for differences between two arrays

template < class view_type >
long double check_diff ( const view_type& a , const view_type& b )
{
  using namespace vigra::multi_math ;
  using namespace vigra::acc;
  
  typedef typename view_type::value_type value_type ;
  typedef typename vigra::ExpandElementResult < value_type > :: type real_type ;
  typedef MultiArray<view_type::actual_dimension,real_type> error_array ;

  error_array ea ( vigra::multi_math::squaredNorm ( b - a ) ) ;
  AccumulatorChain<real_type,Select< Mean, Maximum> > ac ;
  extractFeatures(ea.begin(), ea.end(), ac);
  std::cout << "warped image diff Mean: " << sqrt(get<Mean>(ac)) << std::endl;
  
  long double max_error = sqrt(get<Maximum>(ac)) ;
  std::cout << "warped image diff Maximum: "
            << max_error << std::endl ;
            
//   if ( max_error > .01 )
//   {
//     for ( int y = 0 ; y < a.shape(1) ; y++ )
//     {
//       for ( int x = 0 ; x < a.shape(0) ; x++ )
//       {
//         auto pa = a ( x , y ) ;
//         auto pb = b ( x , y ) ;
//         auto na = pa[0] * pa[0] + pa[1] * pa[1] + pa[2] * pa[2] ;
//         auto nb = pb[0] * pb[0] + pb[1] * pb[1] + pb[2] * pb[2] ;
//         if ( std::abs ( na - nb ) > .01 )
//           std::cout << "excess at ( " << x << " , "
//                     << y << " )" << std::endl ;
//       }
//     }
//   }
  return max_error ;
}

template < class view_type ,
           typename real_type ,
           typename rc_type ,
           int specialize >
long double run_test ( view_type & data ,
                       vspline::bc_code bc ,
                       int DEGREE ,
                       int TIMES = 32 )
{
  typedef typename view_type::value_type pixel_type ;
  typedef typename view_type::difference_type Shape;
  typedef MultiArray < 2 , pixel_type > array_type ;
  typedef int int_type ;
  auto shape = data.shape() ;
  
  long double max_error = 0.0L ;
  long double error ;
  
  // we use simdized types with as many elements as vector_traits
  // considers appropriate for a given real_type, which is the elementary
  // type of the (pixel) data we process:
  
  const int vsize = vspline::vector_traits < real_type > :: size ;
  
  vspline::bcv_type < view_type::actual_dimension > bcv ( bc ) ;
  
  int Nx = data.width() ;
  int Ny = data.height() ;

  vspline::bspline < pixel_type , 2 > bsp ( data.shape() , DEGREE , bcv ) ;
  
  size_t max_extent = shape[0] > shape[1] ? shape[0] : shape[1] ;

  vigra::MultiArray < 1 , rc_type > indexes ( max_extent ) ;
  for ( size_t i = 0 ; i < max_extent ; i++ )
    indexes[i] = i ;

  vigra::MultiArrayView < 1 , rc_type > ix0
    ( vigra::Shape1 ( shape[0] ) , indexes.data() ) ;
  
  vigra::MultiArrayView < 1 , rc_type > ix1
    ( vigra::Shape1 ( shape[1] ) , indexes.data() ) ;
  
  vspline::grid_spec < 2 , rc_type > grid { ix0 , ix1 } ;

  bsp.core = data ;
  
  // first test: time prefilter

#ifdef PRINT_ELAPSED
  std::chrono::system_clock::time_point start = std::chrono::system_clock::now();
  std::chrono::system_clock::time_point end ;
#endif
  
  for ( int times = 0 ; times < TIMES ; times++ )
    bsp.prefilter() ;
  
#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x prefilter:............................ "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif
  
  // to time and test 1D operation, we pretend the data are 1D,
  // prefilter and restore them.

  start = std::chrono::system_clock::now();
  
  // cast data to 1D array
  
  vigra::MultiArrayView < 1 , pixel_type > fake_1d_array
    ( vigra::Shape1 ( prod ( data.shape() ) ) , data.data() ) ;
    
  vigra::TinyVector < vspline::bc_code , 1 > bcv1 ( bcv[0] ) ;
  
  vspline::bspline < pixel_type , 1 > bsp1
    ( fake_1d_array.shape() , DEGREE , bcv1 ) ;
  
  bsp1.core = fake_1d_array ;
  
  for ( int times = 0 ; times < TIMES ; times++ )
    bsp1.prefilter() ;

  #ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x prefilter as fake 1D array:........... "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif

  // use fresh data, data above are useless after TIMES times filtering
  bsp1.core = fake_1d_array ;
  bsp1.prefilter() ;
  
  start = std::chrono::system_clock::now();

  vigra::MultiArray < 1 , pixel_type > fake_1d_target
    ( vigra::Shape1 ( prod ( data.shape() ) ) ) ;
    
  vspline::restore < 1 , pixel_type > ( bsp1 , fake_1d_target ) ;
    
  for ( int times = 1 ; times < TIMES ; times++ )
    vspline::restore < 1 , pixel_type > ( bsp1 , fake_1d_target ) ;
  
#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x restore original data from 1D:........ "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif
  
  cout << "difference original data/restored data:" << endl ;
  error = check_diff<decltype(fake_1d_array)> ( fake_1d_array , fake_1d_target ) ;
  if ( error > max_error )
    max_error = error ; 

  // that's all of the 1D testing we do, back to the 2D data.
  // use fresh data, data above are useless after TIMES times filtering
  bsp.core = data ;
  bsp.prefilter() ;

  // get a view to the core coefficients (those which aren't part of the brace)
  view_type cfview = bsp.core ;

  // set the coordinate type
  typedef vigra::TinyVector < rc_type , 2 > coordinate_type ;
  
  // set the evaluator type
  typedef vspline::evaluator<coordinate_type,pixel_type,vsize,specialize> eval_type ;

  // set the discrete coordinate type
  typedef vigra::TinyVector < int , 2 > dsc_coordinate_type ;
  
  // set the discrete evaluator type
  // currently unused, but is an alternative for unit-spaced
  // index-based transform.
  typedef vspline::evaluator<dsc_coordinate_type,pixel_type,vsize,specialize> dsc_eval_type ;

  // create the evaluator for the b-spline, using plain evaluation (no derivatives)
  
  eval_type raw_ev ( bsp ) ;
  dsc_eval_type dsc_ev ( bsp ) ;
  
  typedef vspline::bspline < pixel_type , 2 > spline_type ;
  
  auto ev = vspline::make_safe_evaluator < spline_type , rc_type > ( bsp ) ;

  // type for coordinate array
  typedef vigra::MultiArray<2, coordinate_type> coordinate_array ;
  
  int Tx = Nx ;
  int Ty = Ny ;

  // now we create a warp array of coordinates at which the spline will be evaluated.
  // Also create a target array to contain the result.

  coordinate_array fwarp ( Shape ( Tx , Ty ) ) ;
  array_type _target ( Shape(Tx,Ty) ) ;
  view_type target ( _target ) ;
  
  rc_type dfx = 0.0 , dfy = 0.0 ; // currently evaluating right at knot point locations
  
  for ( int times = 0 ; times < 1 ; times++ )
  {
    for ( int y = 0 ; y < Ty ; y++ )
    {
      rc_type fy = (rc_type)(y) + dfy ;
      for ( int x = 0 ; x < Tx ; x++ )
      {
        rc_type fx = (rc_type)(x) + dfx ;
        // store the coordinate to fwarp[x,y]
        fwarp [ Shape ( x , y ) ] = coordinate_type ( fx , fy ) ;
      }
    }
  }
 
  // second test. perform a transform using fwarp as warp array. Since fwarp contains
  // the discrete coordinates to the knot points, converted to float, the result
  // should be the same as the input within the given precision

#ifdef PRINT_ELAPSED
  start = std::chrono::system_clock::now();
#endif
  
  for ( int times = 0 ; times < TIMES ; times++ )
    vspline::transform ( ev , fwarp , target ) ;

  
#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x transform with ready-made bspline:.... "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif
       
  error = check_diff<view_type> ( target , data ) ;
  if ( error > max_error )
    max_error = error ; 

  // third test: do the same with the 'classic' remap routine which internally creates
  // a b-spline

#ifdef PRINT_ELAPSED
  start = std::chrono::system_clock::now();
#endif
  
  for ( int times = 0 ; times < TIMES ; times++ )
    vspline::remap ( data , fwarp , target , bcv , DEGREE ) ;

#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x classic remap:........................ "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif

  error = check_diff<view_type> ( target , data ) ;
  if ( error > max_error )
    max_error = error ; 

  // fourth test: perform an transform() directly using the b-spline evaluator
  // as the transform's functor. This is, yet again, the same, because
  // it evaluates at all discrete positions, but now without the warp array:
  // the index-based transform feeds the evaluator with the discrete coordinates.

#ifdef PRINT_ELAPSED
  start = std::chrono::system_clock::now();
#endif
  
  for ( int times = 0 ; times < TIMES ; times++ )
    vspline::transform ( ev , target ) ;

#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x index-based transform................. "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif

  cout << "difference original data/restored data:" << endl ;
  error = check_diff<view_type> ( target , data ) ;
  if ( error > max_error )
    max_error = error ; 

  // fifth test: use 'restore' which internally uses convolution. This is
  // usually slightly faster than the previous way to restore the original
  // data, but otherwise makes no difference.

#ifdef PRINT_ELAPSED
  start = std::chrono::system_clock::now();
#endif
  
  for ( int times = 0 ; times < TIMES ; times++ )
    vspline::restore < 2 , pixel_type > ( bsp , target ) ;
  
#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x restore original data: .............. "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif
       
  cout << "difference original data/restored data:" << endl ;
  error = check_diff<view_type> ( data , target ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << endl ;
  
#ifdef PRINT_ELAPSED
  start = std::chrono::system_clock::now();
#endif
  
  for ( int times = 0 ; times < TIMES ; times++ )
    vspline::transform ( grid , raw_ev , target ) ;
  
#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x grid eval over canonical indexes: ... "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif
       
  cout << "difference original data/restored data:" << endl ;
  error = check_diff<view_type> ( data , target ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << endl ;
  
#ifdef PRINT_ELAPSED
  start = std::chrono::system_clock::now();
#endif
  
  for ( int times = 0 ; times < TIMES ; times++ )
    vspline::transform ( grid , ev , target ) ;
  
#ifdef PRINT_ELAPSED
  end = std::chrono::system_clock::now();
  cout << "avg " << TIMES << " x gen. grid eval over canonical indexes: "
       << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
       << " ms" << endl ;
#endif
       
  cout << "difference original data/restored data:" << endl ;
  error = check_diff<view_type> ( data , target ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << endl ;
  
  return max_error ;
}

template < class real_type , class rc_type >
long double process_image ( char * name )
{
  long double max_error = 0.0L ;
  long double error ;
  
  cout << fixed << showpoint << setprecision(16) ;
  
  // the import and info-displaying code is taken from vigra:

  vigra::ImageImportInfo imageInfo(name);
  // print some information
  std::cout << "Image information:\n";
  std::cout << "  file format: " << imageInfo.getFileType() << std::endl;
  std::cout << "  width:       " << imageInfo.width() << std::endl;
  std::cout << "  height:      " << imageInfo.height() << std::endl;
  std::cout << "  pixel type:  " << imageInfo.getPixelType() << std::endl;
  std::cout << "  color image: ";
  if (imageInfo.isColor())    std::cout << "yes (";
  else                        std::cout << "no  (";
  std::cout << "number of channels: " << imageInfo.numBands() << ")\n";

  typedef vigra::RGBValue<real_type,0,1,2> pixel_type; 
  typedef vigra::MultiArray<2, pixel_type> array_type ;
  typedef vigra::MultiArrayView<2, pixel_type> view_type ;

  // to test that strided data are processed correctly, we load the image
  // to an inner subarray of containArray
  
//   array_type containArray(imageInfo.shape()+vigra::Shape2(3,5));
//   view_type imageArray = containArray.subarray(vigra::Shape2(1,2),vigra::Shape2(-2,-3)) ;
  
  // alternatively, just use the same for both
  
  array_type containArray ( imageInfo.shape() );
  view_type imageArray ( containArray ) ;
  
  vigra::importImage(imageInfo, imageArray);
  
  // test these boundary conditions:

  vspline::bc_code bcs[] =
  {
    vspline::MIRROR ,
    vspline::REFLECT ,
    vspline::NATURAL ,
    vspline::PERIODIC
  } ;

  for ( int b = 0 ; b < 4 ; b++ )
  {
    vspline::bc_code bc = bcs[b] ;
    for ( int spline_degree = 1 ; spline_degree < 8 ; spline_degree++ )
    {
#if defined USE_VC

      cout << "testing bc code " << vspline::bc_name[bc]
           << " spline degree " << spline_degree
           << " using Vc" << endl ;

#elif defined USE_HWY

      cout << "testing bc code " << vspline::bc_name[bc]
           << " spline degree " << spline_degree
           << " using highway" << endl ;

#elif defined USE_STDSIMD

      cout << "testing bc code " << vspline::bc_name[bc]
           << " spline degree " << spline_degree
           << " using std::simd" << endl ;

#else

      cout << "testing bc code " << vspline::bc_name[bc]
           << " spline degree " << spline_degree
           << " using SIMD emulation" << endl ;
           
#endif
           
      if ( spline_degree == 0 )
      {
        std::cout << "using specialized evaluator" << std::endl ;
        error = run_test < view_type , real_type , rc_type , 0 >
          ( imageArray , bc , spline_degree ) ;
        if ( error > max_error )
          max_error = error ; 
        std::cout << "using unspecialized evaluator" << std::endl ;
        error = run_test < view_type , real_type , rc_type , -1 >
          ( imageArray , bc , spline_degree ) ;
        if ( error > max_error )
          max_error = error ; 
      }
      else if ( spline_degree == 1 )
      {
        std::cout << "using specialized evaluator" << std::endl ;
        error = run_test < view_type , real_type , rc_type , 1 >
          ( imageArray , bc , spline_degree ) ;
        if ( error > max_error )
          max_error = error ; 
        std::cout << "using unspecialized evaluator" << std::endl ;
        error = run_test < view_type , real_type , rc_type , -1 >
          ( imageArray , bc , spline_degree ) ;
        if ( error > max_error )
          max_error = error ; 
      }
      else
      {
        error = run_test < view_type , real_type , rc_type , -1 >
          ( imageArray , bc , spline_degree ) ;
        if ( error > max_error )
          max_error = error ; 
      }
    }
  }
  return max_error ;
}

int main ( int argc , char * argv[] )
{
  long double max_error = 0.0L ;
  long double error ;
  
  cout << "testing float data, float coordinates" << endl ;
  error = process_image<float,float> ( argv[1] ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << "max error of float/float test: " << error << std::endl << std::endl ;

  cout << endl << "testing double data, double coordinates" << endl ;
  error = process_image<double,double> ( argv[1] ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << "max error of double/double test: " << error << std::endl << std::endl ;

  cout << endl << "testing long double data, float coordinates" << endl ;
  error = process_image<long double,float> ( argv[1] ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << "max error of ldouble/float test: " << error << std::endl << std::endl ;
  
  cout << endl << "testing long double data, double coordinates" << endl ;
  error = process_image<long double,double> ( argv[1] ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << "max error of ldouble/double test: " << error << std::endl << std::endl ;
  
  cout << "testing float data, double coordinates" << endl ;
  error = process_image<float,double> ( argv[1] ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << "max error of float/double test: " << error << std::endl << std::endl ;
  
  cout << endl << "testing double data, float coordinates" << endl ;
  error = process_image<double,float> ( argv[1] ) ;
  if ( error > max_error )
    max_error = error ; 
  cout << "max error of double/float test: " << error << std::endl << std::endl ;
  
  cout << "reached end. max error of all tests: " << max_error << std::endl ;
}