#include <itkImage.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkSignedMaurerDistanceMapImageFilter.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <fpa/Image/SkeletonFilter.h>
#include <cpExtensions/Algorithms/SkeletonToImageFilter.h>

typedef double TScalar;
typedef itk::Image< TScalar, 3 > TImage;
typedef itk::MinimumMaximumImageCalculator< TImage > TMinMax;
typedef itk::SignedMaurerDistanceMapImageFilter< TImage, TImage > TDMap;
typedef fpa::Image::SkeletonFilter< TImage > TSkeletonFilter;
typedef TSkeletonFilter::TSkeleton TSkeleton;
typedef cpExtensions::Algorithms::SkeletonToImageFilter< TSkeleton, TImage > TSkeletonToImage;
typedef itk::ImageFileReader< TImage > TReader;
typedef itk::ImageFileWriter< TImage > TWriter;

int main( int argc, char* argv[] )
{
  if( argc < 4 )
  {
    std::cerr
      << "Usage: " << argv[ 0 ]
      << " image prefix number_of_samples [seed_x seed_y seed_z]"
      << std::endl;
    return( 1 );

  } // fi
  std::string image_name = argv[ 1 ];

  TReader::Pointer reader = TReader::New( );
  reader->SetFileName( image_name );
  reader->Update( );

  TDMap::Pointer dmap = TDMap::New( );
  dmap->SetInput( reader->GetOutput( ) );
  dmap->InsideIsPositiveOn( );
  dmap->Update( );
  const TImage* image = dmap->GetOutput( );

  std::string prefix = argv[ 2 ];
  std::istringstream str( argv[ 3 ] );
  unsigned int number_of_samples;
  str >> number_of_samples;
  TImage::IndexType seed;
  if( argc < 7 )
  {
    TMinMax::Pointer minmax = TMinMax::New( );
    minmax->SetImage( image );
    minmax->Compute( );
    seed = minmax->GetIndexOfMaximum( );
  }
  else
  {
    std::istringstream str_x( argv[ 4 ] );
    std::istringstream str_y( argv[ 5 ] );
    std::istringstream str_z( argv[ 6 ] );
    TImage::PointType pnt;
    str_x >> pnt[ 0 ];
    str_y >> pnt[ 1 ];
    str_z >> pnt[ 2 ];
    image->TransformPhysicalPointToIndex( pnt, seed );

  } // fi

  for( unsigned int i = 0; i < number_of_samples; ++i )
  {
    std::cout << i << std::endl;
    TSkeletonFilter::Pointer skeleton = TSkeletonFilter::New( );
    skeleton->SetInput( image );
    skeleton->AddSeed( seed, 0 );
    skeleton->Update( );

    TSkeletonToImage::Pointer sk2im = TSkeletonToImage::New( );
    sk2im->SetTemplateImage( image );
    sk2im->SetSkeleton( skeleton->GetSkeleton( ) );
    sk2im->SetOutsideValue( 0 );
    sk2im->SetInsideValue( 1 );
    sk2im->Update( );

    TWriter::Pointer writer = TWriter::New( );
    writer->SetInput( sk2im->GetOutput( ) );
    writer->SetFileName( "skeleton.mhd" );
    writer->Update( );

  } // rof
  return( 0 );
}

// eof - $RCSfile$