#include <chrono>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <fpa/Image/Mori.h>

// -------------------------------------------------------------------------
const unsigned int Dim = 3;
typedef short          TPixel;
typedef unsigned short TLabel;

typedef itk::Image< TPixel, Dim > TInputImage;
typedef itk::Image< TLabel, Dim > TLabelImage;
typedef fpa::Image::Mori< TInputImage, TLabelImage > TFilter;

// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
  // Get arguments
  if( argc < 17 )
  {
    std::cerr
      << "Usage: " << argv[ 0 ] << std::endl
      << "   input_image output_image output_signal" << std::endl
      << "   init_threshold(-1024) end_threshold(0) delta(1)" << std::endl
      << "   minimum_threshold(-850)" << std::endl
      << "   inside_value(255) outside_value(0)" << std::endl
      << "   signal_kernel_size(20) signal_threshold(500) signal_influence(0.5)"
      << std::endl
      << "   [index/point] seed_x seed_y seed_z"
      << std::endl;
    return( 1 );

  } // fi
  std::string input_image_filename = argv[ 1 ];
  std::string output_image_filename = argv[ 2 ];
  std::string output_signal_filename = argv[ 3 ];
  TPixel init_threshold = std::atoi( argv[ 4 ] );
  TPixel end_threshold = std::atoi( argv[ 5 ] );
  TPixel delta = std::atoi( argv[ 6 ] );
  TPixel minimum_threshold = std::atoi( argv[ 7 ] );
  TLabel inside_value = std::atoi( argv[ 8 ] );
  TLabel outside_value = std::atoi( argv[ 9 ] );
  unsigned long signal_kernel_size = std::atoi( argv[ 10 ] );
  double signal_threshold = std::atof( argv[ 11 ] );
  double signal_influence = std::atof( argv[ 12 ] );
  std::string seed_type = argv[ 13 ];
  TInputImage::IndexType iseed;
  TInputImage::PointType pseed;
  for( unsigned int i = 0; i < Dim; ++i )
  {
    if( seed_type == "index" )
      iseed[ i ] = std::atoi( argv[ 14 + i ] );
    else
      pseed[ i ] = std::atof( argv[ 14 + i ] );

  } // rof

  // Read image
  itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
    itk::ImageFileReader< TInputImage >::New( );
  input_image_reader->SetFileName( input_image_filename );

  // Prepare filter
  TFilter::Pointer filter = TFilter::New( );
  filter->SetInput( input_image_reader->GetOutput( ) );
  if( seed_type == "index" )
    filter->SetSeed( iseed );
  else
    filter->SetSeed( pseed );
  filter->SetThresholds( init_threshold, end_threshold, delta );
  filter->SetMinimumThreshold( minimum_threshold );
  filter->SetInsideValue( inside_value );
  filter->SetOutsideValue( outside_value );
  filter->SetSignalKernelSize( signal_kernel_size );
  filter->SetSignalThreshold( signal_threshold );
  filter->SetSignalInfluence( signal_influence );

  // Show some information
  std::cout << "----------------------------------------------" << std::endl;
  std::cout << "Image: " << input_image_filename << std::endl;

  // Execute pipeline
  std::chrono::time_point< std::chrono::high_resolution_clock > ts, te;
  std::chrono::duration< double > tr;
  try
  {
    ts = std::chrono::high_resolution_clock::now( );
    input_image_reader->Update( );
    te = std::chrono::high_resolution_clock::now( );
    tr = te - ts;
    std::cout << "Read time: " << tr.count( ) << " s" << std::endl;

    ts = std::chrono::high_resolution_clock::now( );
    filter->Update( );
    te = std::chrono::high_resolution_clock::now( );
    tr = te - ts;
    std::cout
      << "Mori time: " << tr.count( ) << " s" << std::endl
      << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
      << "Number of thresholds: "
      << filter->GetNumberOfEvaluatedThresholds( ) << std::endl;
  }
  catch( std::exception& err )
  {
    std::cerr << "Error caught: " << err.what( ) << std::endl;
    return( 1 );

  } // yrt

  // Save output image
  itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer =
    itk::ImageFileWriter< TLabelImage >::New( );
  output_image_writer->SetInput( filter->GetThresholdedOutput( ) );
  output_image_writer->SetFileName( output_image_filename );
  try
  {
    ts = std::chrono::high_resolution_clock::now( );
    output_image_writer->Update( );
    te = std::chrono::high_resolution_clock::now( );
    tr = te - ts;
    std::cout << "Write time: " << tr.count( ) << " s" << std::endl;
  }
  catch( std::exception& err )
  {
    std::cerr << "Error caught: " << err.what( ) << std::endl;
    return( 1 );

  } // yrt

  // Save output signal
  std::ofstream osignal( output_signal_filename.c_str( ) );
  for(
    unsigned long i = 0; i < filter->GetNumberOfEvaluatedThresholds( ); ++i
    )
  {
    double x, y;
    TFilter::TPeak p;
    filter->GetSignalValues( i, x, y, p );
    osignal << x << " " << y << std::endl;

  } // rof
  osignal.close( );
  std::cout << "----------------------------------------------" << std::endl;

  return( 0 );
}

// eof - $RCSfile$