#include <ctime>
#include <iostream>
#include <limits>
#include <set>
#include <string>

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkImageToVTKImageFilter.h>

#include <vtkImageActor.h>
#include <vtkImageMarchingCubes.h>
#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>

#include <fpa/Image/RegionGrowWithMultipleThresholds.h>
#include <fpa/VTK/ImageMPR.h>
#include <fpa/VTK/Image3DObserver.h>

// -------------------------------------------------------------------------
const unsigned int Dim = 3;
typedef short TPixel;
typedef double TScalar;
typedef itk::Image< TPixel, Dim >  TImage;
typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;

typedef itk::ImageFileReader< TImage >   TImageReader;
typedef itk::ImageFileWriter< TImage >   TImageWriter;

typedef
fpa::Image::RegionGrowWithMultipleThresholds< TImage >
TFrontAlgorithm;
typedef
fpa::VTK::Image3DObserver< TFrontAlgorithm, vtkRenderWindow >
TObserver;

// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
  if( argc < 6 )
  {
    std::cerr
      << "Usage: " << argv[ 0 ]
      << " input_image thr_0 thr_1 step output_image"
      << " visual_debug"
      << std::endl;
    return( 1 );

  } // fi
  std::string input_image_fn = argv[ 1 ];
  TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) );
  TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) );
  unsigned int step = std::atoi( argv[ 4 ] );
  std::string output_image_fn = argv[ 5 ];
  bool visual_debug = false;
  if( argc > 6 )
    visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
  
  // Read image
  TImageReader::Pointer input_image_reader = TImageReader::New( );
  input_image_reader->SetFileName( input_image_fn );
  try
  {
    input_image_reader->Update( );
  }
  catch( itk::ExceptionObject& err )
  {
    std::cerr << "Error caught: " << err << std::endl;
    return( 1 );

  } // yrt
  TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  TImage::SpacingType spac = input_image->GetSpacing( );
  double min_spac = spac[ 0 ];
  for( unsigned int d = 1; d < Dim; d++ )
    min_spac = ( spac[ d ] < min_spac )? spac[ d ]: min_spac;

  TVTKImage::Pointer vtk_image = TVTKImage::New( );
  vtk_image->SetInput( input_image );
  vtk_image->Update( );

  fpa::VTK::ImageMPR view;
  view.SetBackground( 0.3, 0.2, 0.8 );
  view.SetSize( 800, 800 );
  view.SetImage( vtk_image->GetOutput( ) );

  // Wait for a seed to be given
  while( view.GetNumberOfSeeds( ) == 0 )
    view.Start( );

  // Convert seed
  double seed[ 3 ];
  view.GetSeed( 0, seed );
  TImage::PointType seed_pnt;
  seed_pnt[ 0 ] = seed[ 0 ];
  seed_pnt[ 1 ] = seed[ 1 ];
  seed_pnt[ 2 ] = seed[ 2 ];
  TImage::IndexType seed_idx;
  input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );

  // Configure algorithm
  TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( );
  algorithm->AddThresholds( thr_0, thr_1, step );
  algorithm->AddSeed( seed_idx, 0 );
  algorithm->SetInput( input_image );
  algorithm->SetNeighborhoodOrder( 1 );
  algorithm->SetDifferenceThreshold( double( 3 ) );

  if( visual_debug )
  {
    // Configure observer
    TObserver::Pointer obs = TObserver::New( );
    obs->SetImage( input_image, view.GetWindow( ) );
    algorithm->AddObserver( itk::AnyEvent( ), obs );
    algorithm->ThrowEventsOn( );
  }
  else
    algorithm->ThrowEventsOff( );

  std::clock_t start = std::clock( );
  algorithm->Update( );
  std::clock_t end = std::clock( );
  double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
  std::cout << "Execution time = " << seconds << std::endl;

  // Show result
  TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
  output_image_vtk->SetInput( algorithm->GetOutput( ) );
  output_image_vtk->Update( );

  vtkSmartPointer< vtkImageMarchingCubes > mc =
    vtkSmartPointer< vtkImageMarchingCubes >::New( );
  mc->SetInputData( output_image_vtk->GetOutput( ) );
  mc->SetValue(
    0,
    double( algorithm->GetInsideValue( ) ) * double( 0.95 )
    );
  mc->Update( );

  // Let some interaction and close program
  view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
  view.Start( );

  // Write resulting image
  TImageWriter::Pointer output_image_writer = TImageWriter::New( );
  output_image_writer->SetInput( algorithm->GetOutput( ) );
  output_image_writer->SetFileName( output_image_fn );
  try
  {
    output_image_writer->Update( );
  }
  catch( itk::ExceptionObject& err )
  {
    std::cerr << "Error caught: " << err << std::endl;
    return( 1 );

  } // yrt
  return( 0 );
}

// eof - $RCSfile$