#include <iostream>
#include <fstream>
#include <limits>
#include <string>

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

#include <vtkCellArray.h>
#include <vtkFloatArray.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkSmartPointer.h>

#include <fpa/VTK/ImageMPR.h>

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

typedef itk::Image< TPixel, Dimension > TImage;
 
// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
  if( argc < 3 )
  {
    std::cerr
      << "Usage: " << argv[ 0 ]
      << " input_image input_seeds"
      << std::endl;
      return( 1 );

  } // fi
  std::string input_image_fn = argv[ 1 ];
  std::string input_seeds_fn = argv[ 2 ];

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

  } // yrt
  TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  itk::ImageToVTKImageFilter< TImage >::Pointer input_image_vtk =
    itk::ImageToVTKImageFilter< TImage >::New( );
  input_image_vtk->SetInput( input_image );
  input_image_vtk->Update( );

  // Read seeds
  std::ifstream input_seeds_fs( input_seeds_fn.c_str( ) );
  if( !input_seeds_fs )
  {
    std::cerr
      << "Error opening file \"" << input_seeds_fn << "\"" << std::endl;
    return( 1 );

  } // fi

  // Read seed points
  std::string str_val;
  unsigned int number_of_points;
  input_seeds_fs >> str_val >> number_of_points;
  input_seeds_fs >> str_val; // X
  input_seeds_fs >> str_val; // Y
  input_seeds_fs >> str_val; // Z
  input_seeds_fs >> str_val; // value
  input_seeds_fs >> str_val; // Label

  vtkSmartPointer< vtkPoints > input_seeds_points =
    vtkSmartPointer< vtkPoints >::New( );
  vtkSmartPointer< vtkCellArray > input_seeds_cells =
    vtkSmartPointer< vtkCellArray >::New( );
  vtkSmartPointer< vtkFloatArray > input_seeds_ids =
    vtkSmartPointer< vtkFloatArray >::New( );

  input_seeds_ids->SetNumberOfComponents( 1 );
  input_seeds_ids->SetNumberOfTuples( number_of_points );

  double min_value = std::numeric_limits< double >::max( );
  double max_value = double( 0 );
  for( unsigned int i = 0; i < number_of_points; ++i )
  {
    unsigned int x, y, z;
    double value;
    input_seeds_fs >> x >> y >> z >> value >> str_val;

    min_value = ( value < min_value )? value: min_value;
    max_value = ( value > max_value )? value: max_value;

    TImage::IndexType idx;
    idx[ 0 ] = x;
    idx[ 1 ] = y;
    idx[ 2 ] = z;
    TImage::PointType pnt;
    input_image->TransformIndexToPhysicalPoint( idx, pnt );

    input_seeds_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
    input_seeds_cells->InsertNextCell( 1 );
    input_seeds_cells->InsertCellPoint( i );
    input_seeds_ids->SetTuple1(
      i, double( i ) / double( number_of_points - 1 )
      );


  } // rof
  input_seeds_fs.close( );

  vtkSmartPointer< vtkPolyData > input_seeds =
    vtkSmartPointer< vtkPolyData >::New( );
  input_seeds->SetPoints( input_seeds_points );
  input_seeds->SetVerts( input_seeds_cells );
  input_seeds->GetPointData( )->SetScalars( input_seeds_ids );

  // Show input image and let some interaction
  fpa::VTK::ImageMPR view;
  view.SetBackground( 0.3, 0.2, 0.8 );
  view.SetSize( 600, 600 );
  view.SetImage( input_image_vtk->GetOutput( ) );
  view.SetWindowLevel(
    max_value - min_value,
    ( min_value + max_value ) / double( 2 )
    );
  view.AddPolyData( input_seeds );
  view.Render( );
  view.Start( );

  return( 0 );
}
 
// eof - $RCSfile$