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

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

#include <itkSignedMaurerDistanceMapImageFilter.h>

#include <vtkOutlineSource.h>

#include <fpa/Base/Functors/InvertCostFunction.h>
#include <fpa/Image/DijkstraWithEndPointDetection.h>
#include <fpa/VTK/Image3DObserver.h>
#include <fpa/IO/MinimumSpanningTreeWriter.h>
#include <fpa/IO/UniqueValuesContainerWriter.h>
#include <fpa/IO/MatrixValuesContainerWriter.h>
#include <fpa/VTK/UniqueVerticesToPolyDataFilter.h>

#include "fpa_Utility.h"

// -------------------------------------------------------------------------
const unsigned int Dim = 3;
typedef unsigned char TPixel;
typedef float         TScalar;
typedef itk::Image< TPixel, Dim >  TImage;
typedef itk::Image< TScalar, Dim > TScalarImage;

// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
  if( argc < 7 )
  {
    std::cerr
      << "Usage: " << argv[ 0 ] << std::endl
      << "\tinput_image" << std::endl
      << "\toutput_labels" << std::endl
      << "\toutput_minimum_spanning_tree" << std::endl
      << "\toutput_endpoints" << std::endl
      << "\toutput_bifurcations" << std::endl
      << "\toutput_branches" << std::endl
      << std::endl;
    return( 1 );

  } // fi
  std::string input_image_fn = argv[ 1 ];
  std::string output_labels_fn = argv[ 2 ];
  std::string output_minimum_spanning_tree_fn = argv[ 3 ];
  std::string output_endpoints_fn = argv[ 4 ];
  std::string output_bifurcations_fn = argv[ 5 ];
  std::string output_branches_fn = argv[ 6 ];

  // Read image
  TImage::Pointer input_image;
  std::string err = fpa_Utility::ReadImage( input_image, input_image_fn );
  if( err != "" )
  {
    std::cerr << err << std::endl;
    return( 1 );

  } // fi

  // Show image and wait for, at least, one seed
  itk::ImageToVTKImageFilter< TImage >::Pointer vtk_input_image =
    itk::ImageToVTKImageFilter< TImage >::New( );
  vtk_input_image->SetInput( input_image );
  vtk_input_image->Update( );

  vtkSmartPointer<vtkOutlineSource> outline =
    vtkSmartPointer<vtkOutlineSource>::New();
  outline->SetBounds(vtk_input_image->GetOutput( )->GetBounds());

  vtkSmartPointer<vtkPolyDataMapper> mapper =
    vtkSmartPointer<vtkPolyDataMapper>::New( );
  mapper->SetInputConnection( outline->GetOutputPort( ) );

  vtkSmartPointer<vtkActor> actor =
    vtkSmartPointer<vtkActor>::New( );
  actor->SetMapper( mapper );
 vtkSmartPointer<vtkRenderer> aRenderer =
    vtkSmartPointer<vtkRenderer>::New();
  vtkSmartPointer<vtkRenderWindow> renWin =
    vtkSmartPointer<vtkRenderWindow>::New();
  renWin->AddRenderer(aRenderer);

  vtkSmartPointer<vtkRenderWindowInteractor> iren =
    vtkSmartPointer<vtkRenderWindowInteractor>::New();
  iren->SetRenderWindow(renWin);

  aRenderer->AddActor( actor );
iren->Initialize();
  iren->Start();

 /*
  fpa_Utility::Viewer2DWithSeeds viewer;
  viewer.SetImage( vtk_input_image->GetOutput( ) );
  while( viewer.GetNumberOfSeeds( ) == 0 )
    viewer.Start( );
  */

  // Compute squared distance map
  itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::Pointer
    dmap =
    itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::New( );
  dmap->SetInput( input_image );
  dmap->SetBackgroundValue( TPixel( 0 ) );
  dmap->InsideIsPositiveOn( );
  dmap->SquaredDistanceOn( );
  dmap->UseImageSpacingOn( );

  std::time_t start, end;
  std::time( &start );
  dmap->Update( );
  std::time( &end );
  std::cout
    << "Distance map time = "
    << std::difftime( end, start )
    << "s." << std::endl;

  // Prepare cost conversion function
  typedef fpa::Base::Functors::InvertCostFunction< TScalar > TCostFunction;
  TCostFunction::Pointer cost_f = TCostFunction::New( );

  // Prepare Dijkstra filter
  typedef fpa::Image::
    DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
  TFilter::Pointer filter = TFilter::New( );
  filter->SetInput( dmap->GetOutput( ) );
  filter->SetConversionFunction( cost_f );
  filter->SetNeighborhoodOrder( 2 );
  filter->StopAtOneFrontOff( );

  // Associate seed
  TImage::PointType pnt;
  TImage::IndexType idx;
  /*
    viewer.GetSeed( pnt, 0 );
  */
  pnt[ 0 ] = 0.879066;
  pnt[ 1 ] =  -109.36591;
  pnt[ 2 ] =  1942.480988; 

  if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
    filter->AddSeed( idx, 0 );

  // Prepare graphical debugger
  typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
  TDebugger::Pointer debugger = TDebugger::New( );
  debugger->SetRenderWindow( renWin );
  debugger->SetRenderPercentage( 0.0001 );
  filter->AddObserver( itk::AnyEvent( ), debugger );
  filter->ThrowEventsOn( );

  // Go!
  std::time( &start );
  filter->Update( );
  std::time( &end );
  std::cout
    << "Extraction time = "
    << std::difftime( end, start )
    << "s." << std::endl;

  // Create new actors
  typedef vtkSmartPointer< fpa::VTK::UniqueVerticesToPolyDataFilter< TFilter::TUniqueVertices, TFilter::TInputImage > > TVertices2PD;

  TVertices2PD endpoints2pd = TVertices2PD::New( );
  endpoints2pd->SetInput( filter->GetEndPoints( ) );
  endpoints2pd->SetImage( filter->GetInput( ) );
  endpoints2pd->Update( );

  vtkSmartPointer< vtkPolyDataMapper > endpoints_mapper =
    vtkSmartPointer< vtkPolyDataMapper >::New( );
  endpoints_mapper->SetInputConnection( endpoints2pd->GetOutputPort( ) );

  vtkSmartPointer< vtkActor > endpoints_actor =
    vtkSmartPointer< vtkActor >::New( );
  endpoints_actor->SetMapper( endpoints_mapper );
  endpoints_actor->GetProperty( )->SetColor( 0, 1, 0 );
  endpoints_actor->GetProperty( )->SetPointSize( 5 );
  aRenderer->AddActor( endpoints_actor );

  TVertices2PD bifurcations2pd = TVertices2PD::New( );
  bifurcations2pd->SetInput( filter->GetBifurcations( ) );
  bifurcations2pd->SetImage( filter->GetInput( ) );
  bifurcations2pd->Update( );

  vtkSmartPointer< vtkPolyDataMapper > bifurcations_mapper =
    vtkSmartPointer< vtkPolyDataMapper >::New( );
  bifurcations_mapper->SetInputConnection( bifurcations2pd->GetOutputPort( ) );

  vtkSmartPointer< vtkActor > bifurcations_actor =
    vtkSmartPointer< vtkActor >::New( );
  bifurcations_actor->SetMapper( bifurcations_mapper );
  bifurcations_actor->GetProperty( )->SetColor( 1, 0, 0 );
  bifurcations_actor->GetProperty( )->SetPointSize( 5 );
  aRenderer->AddActor( bifurcations_actor );

  // Some more interaction and finish
iren->Initialize();
  iren->Start();

  // Save results
  err = fpa_Utility::SaveImage( filter->GetLabelImage( ), output_labels_fn );
  if( err != "" )
    std::cerr << err << std::endl;

  fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer
    mst_writer =
    fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( );
  mst_writer->SetInput( filter->GetMinimumSpanningTree( ) );
  mst_writer->SetFileName( output_minimum_spanning_tree_fn );
  mst_writer->Update( );

  fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
    endpoints_writer =
    fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
  endpoints_writer->SetInput( filter->GetEndPoints( ) );
  endpoints_writer->SetFileName( output_endpoints_fn );
  endpoints_writer->Update( );

  fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
    bifurcations_writer =
    fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
  bifurcations_writer->SetInput( filter->GetBifurcations( ) );
  bifurcations_writer->SetFileName( output_bifurcations_fn );
  bifurcations_writer->Update( );

  fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer
    branches_writer =
    fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( );
  branches_writer->SetInput( filter->GetBranches( ) );
  branches_writer->SetFileName( output_branches_fn );
  branches_writer->Update( );

  return( 0 );
}

// eof - $RCSfile$

/* TODO
   vtkSmartPointer< vtkPoints > simple_branches_points =
   vtkSmartPointer< vtkPoints >::New( );
   vtkSmartPointer< vtkCellArray > simple_branches_cells =
   vtkSmartPointer< vtkCellArray >::New( );

   vtkSmartPointer< vtkPoints > detailed_branches_points =
   vtkSmartPointer< vtkPoints >::New( );
   vtkSmartPointer< vtkCellArray > detailed_branches_cells =
   vtkSmartPointer< vtkCellArray >::New( );
   vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
   vtkSmartPointer< vtkFloatArray >::New( );

   TFilter::TBranches::ConstIterator brIt = branches->Begin( );
   for( ; brIt != branches->End( ); ++brIt )
   {
   // Branch's first point
   TImage::PointType first_point;
   input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
   unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
   simple_branches_points->InsertNextPoint(
   first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
   );

   TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt );
   for( ; brRowIt != branches->End( brIt ); ++brRowIt )
   {
   // Branch's second point
   TImage::PointType second_point;
   input_image->
   TransformIndexToPhysicalPoint( brRowIt->first, second_point );
   unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
   simple_branches_points->InsertNextPoint(
   second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
   );
   simple_branches_cells->InsertNextCell( 2 );
   simple_branches_cells->InsertCellPoint( first_id );
   simple_branches_cells->InsertCellPoint( second_id );

   // Detailed path
   double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
   TFilter::TVertices path;
   mst->GetPath( path, brIt->first, brRowIt->first );
   TFilter::TVertices::const_iterator pIt = path.begin( );
   for( ; pIt != path.end( ); ++pIt )
   {
   TImage::PointType path_point;
   input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
   detailed_branches_points->InsertNextPoint(
   path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
   );
   detailed_branches_scalars->InsertNextTuple1( pathId );
   if( pIt != path.begin( ) )
   {
   unsigned long nPoints =
   detailed_branches_points->GetNumberOfPoints( );
   detailed_branches_cells->InsertNextCell( 2 );
   detailed_branches_cells->InsertCellPoint( nPoints - 2 );
   detailed_branches_cells->InsertCellPoint( nPoints - 1 );

   } // fi

   } // rof

   } // rof

   } // rof
   vtkSmartPointer< vtkPolyData > simple_branches_polydata =
   vtkSmartPointer< vtkPolyData >::New( );
   simple_branches_polydata->SetPoints( simple_branches_points );
   simple_branches_polydata->SetLines( simple_branches_cells );
   view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );

   vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
   vtkSmartPointer< vtkPolyData >::New( );
   detailed_branches_polydata->SetPoints( detailed_branches_points );
   detailed_branches_polydata->SetLines( detailed_branches_cells );
   detailed_branches_polydata->
   GetPointData( )->SetScalars( detailed_branches_scalars );
   view.AddPolyData( detailed_branches_polydata, 1 );

   // Let some more interaction
   view.Render( );
   view.Start( );

   return( 0 );
   }
*/

// eof - $RCSfile$