#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>

#include <cpPlugins/Extensions/Algorithms/InertiaTensorFunction.h>

#include <itkImageToVTKImageFilter.h>

#include <vtkAxesActor.h>
#include <vtkMatrix4x4.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkSmartPointer.h>
#include <cpPlugins/Extensions/Visualization/MPRWithDifferentWindows.h>

#include <cpPlugins/Extensions/Algorithms/InertiaMedialness.h>


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

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

  // Inertia tensor
  /* TODO
  typedef cpPlugins::Extensions::Algorithms::
    InertiaTensorFunction< TScalar, Dimension >
    TInertiaFunction;
  TInertiaFunction::Pointer inertia = TInertiaFunction::New( );

  itk::ImageRegionConstIteratorWithIndex< TImage > it(
    image_reader->GetOutput( ),
    image_reader->GetOutput( )->GetRequestedRegion( )
    );
  for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
  {
    TInertiaFunction::TPoint pnt;
    image_reader->GetOutput( )->
      TransformIndexToPhysicalPoint( it.GetIndex( ), pnt );
    inertia->AddMass( pnt, TScalar( it.Get( ) ) );

  } // rof

  TInertiaFunction::TMatrix I = inertia->GetInertia( );
  TInertiaFunction::TVector m = inertia->GetCenterOfGravity( );

  TInertiaFunction::TVector pv;
  TInertiaFunction::TVector r;
  TInertiaFunction::TMatrix pm;
  inertia->GetEigenAnalysis( pm, pv, r );

  vtkSmartPointer< vtkMatrix4x4 > matrix =
    vtkSmartPointer< vtkMatrix4x4 >::New( );
  matrix->Identity( );
  matrix->SetElement( 0, 0, pm[ 0 ][ 0 ] * r[ 0 ] );
  matrix->SetElement( 0, 1, pm[ 0 ][ 1 ] * r[ 1 ] );
  matrix->SetElement( 1, 0, pm[ 1 ][ 0 ] * r[ 0 ] );
  matrix->SetElement( 1, 1, pm[ 1 ][ 1 ] * r[ 1 ] );
  matrix->SetElement( 0, 3, m[ 0 ] );
  matrix->SetElement( 1, 3, m[ 1 ] );

  vtkSmartPointer< vtkAxesActor > axes =
    vtkSmartPointer< vtkAxesActor >::New( );
  axes->SetUserMatrix( matrix );

  std::cout << "---------------------------------" << std::endl;
  std::cout << I << std::endl;
  std::cout << m << std::endl;
  std::cout << "---------------------------------" << std::endl;
  std::cout << pm << std::endl;
  std::cout << pv << std::endl;
  std::cout << "---------------------------------" << std::endl;

  typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  TVTKImage::Pointer vtk_input_image = TVTKImage::New( );
  vtk_input_image->SetInput( image_reader->GetOutput( ) );
  vtk_input_image->Update( );

  // Visualization window and interactor
  vtkSmartPointer< vtkRenderWindowInteractor > interactor =
    vtkSmartPointer< vtkRenderWindowInteractor >::New( );
  vtkSmartPointer< vtkRenderWindow > window =
    vtkSmartPointer< vtkRenderWindow >::New( );
  window->SetSize( 500, 500 );
  window->SetInteractor( interactor );

  // Viewer
  cpPlugins::Extensions::Visualization::MPRWithDifferentWindows
    view( NULL, NULL, window, NULL );
  view.SetImage( vtk_input_image->GetOutput( ) );
  view.GetRenderer( 2 )->AddActor( axes );
  view.ResetCamera( 2 );
  view.Render( 2 );

  // Start
  interactor->Initialize( );
  interactor->Start( );

  */

  /*
  // 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::
    InertiaMedialness< TImage, TScalar > TMedialness;
  TMedialness::Pointer medialness = TMedialness::New( );
  medialness->SetInputImage( image_reader->GetOutput( ) );
  medialness->SetMaxRadius( max_radius );
  /*
    medialness->SetMinRadius( min_radius );
    medialness->SetProfileSampling( profile_sampling );
    medialness->SetRadialSampling( radial_sampling );

    TGradient::IndexType idx;
    idx[ 0 ] = 288;
    idx[ 1 ] = 690;

    std::cout << "Medialness = " << medialness->EvaluateAtIndex( idx ) << std::endl;
    // return( 0 );
    */

  // Prepare output filter
  typedef cpPlugins::Extensions::Algorithms::
    ImageFunctionFilter< TImage, TScalarImage, TMedialness >
    TFunctionApplier;
  TFunctionApplier::Pointer medialnessFilter = TFunctionApplier::New( );
  medialnessFilter->SetInput( image_reader->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$