#include <cstdlib>
#include <iostream>
#include <string>

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkSimpleFilterWatcher.h>
#include <itkVector.h>

#include <cpPlugins/Extensions/Algorithms/MultiScaleGaussianImageFilter.h>
#include <cpPlugins/Extensions/Algorithms/GulsunTekMedialness.h>
#include <cpPlugins/Extensions/Algorithms/ImageFunctionFilter.h>

// -------------------------------------------------------------------------
typedef short        TPixel;
typedef float        TScalar;
const   unsigned int Dimension = 3;

typedef itk::Vector< TScalar, Dimension  > TVector;
typedef itk::Image< TPixel,  Dimension >   TImage;
typedef itk::Image< TVector, Dimension >   TGradient;
typedef itk::Image< TScalar,  Dimension >  TScalarImage;

// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
  if( argc < 8 )
  {
    std::cout
      << "Usage: " << argv[ 0 ]
      << " input_image output_medialness"
      << " min_radius max_radius profile_sampling radial_sampling"
      << " s0 s1 ..."
      << std::endl;
    return( 1 );

  } // fi
  std::string image_fn = argv[ 1 ];
  std::string medialness_fn = argv[ 2 ];
  double min_radius = std::atof( argv[ 3 ] );
  double max_radius = std::atof( argv[ 4 ] );
  unsigned int profile_sampling = std::atoi( argv[ 5 ] );
  unsigned int radial_sampling = std::atoi( argv[ 6 ] );

  // Read image
  itk::ImageFileReader< TImage >::Pointer image_reader =
    itk::ImageFileReader< TImage >::New( );
  image_reader->SetFileName( image_fn );
  try
  {
    image_reader->Update( );
  }
  catch( itk::ExceptionObject& err )
  {
    std::cerr << err << std::endl;
    return( 1 );

  } // yrt

  // Prepare filter
  typedef cpPlugins::Extensions::Algorithms::
    MultiScaleGaussianImageFilter< TImage, TGradient > TGradientFilter;

  TGradientFilter::Pointer gradientFilter = TGradientFilter::New( );
  gradientFilter->SetFilterToGradient( );
  gradientFilter->SetInput( image_reader->GetOutput( ) );
  for( int i = 7; i < argc; ++i )
    gradientFilter->AddScale( std::atof( argv[ i ] ) );

  // Execute filter
  itk::SimpleFilterWatcher gradientFilter_watcher( gradientFilter, "Gradient" );
  gradientFilter->Update( );

  // Prepare medialness function
  typedef cpPlugins::Extensions::Algorithms::
    GulsunTekMedialness< TGradient > TMedialness;
  TMedialness::Pointer medialness = TMedialness::New( );
  medialness->SetInputImage( gradientFilter->GetOutput( ) );
  medialness->SetMinRadius( min_radius );
  medialness->SetMaxRadius( max_radius );
  medialness->SetProfileSampling( profile_sampling );
  medialness->SetRadialSampling( radial_sampling );

  // Prepare output filter
  typedef cpPlugins::Extensions::Algorithms::
    ImageFunctionFilter< TGradient, TScalarImage, TMedialness >
    TFunctionApplier;
  TFunctionApplier::Pointer medialnessFilter = TFunctionApplier::New( );
  medialnessFilter->SetInput( gradientFilter->GetOutput( ) );
  medialnessFilter->SetFunction( medialness );

  // Execute filter
  itk::SimpleFilterWatcher medialnessFilter_watcher(
    medialnessFilter, "Medialness"
    );
  medialnessFilter->Update( );

  // Write image
  itk::ImageFileWriter< TScalarImage >::Pointer medialness_writer =
    itk::ImageFileWriter< TScalarImage >::New( );
  medialness_writer->SetInput( medialnessFilter->GetOutput( ) );
  medialness_writer->SetFileName( medialness_fn );
  try
  {
    medialness_writer->Update( );
  }
  catch( itk::ExceptionObject& err )
  {
    std::cerr << err << std::endl;
    return( 1 );

  } // yrt

  return( 0 );
}

// eof - $RCSfile$