// =========================================================================
// @author Leonardo Florez Valencia
// @email florez-l@javeriana.edu.co
// =========================================================================

#include <string>
#include <tclap/CmdLine.h>
#include <itkImage.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkInvertIntensityImageFilter.h>
#include <itkHessianRecursiveGaussianImageFilter.h>
#include <itkHessian3DToVesselnessMeasureImageFilter.h>
#include "Functions.h"

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

// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
  typedef TCLAP::ValueArg< std::string > _TStringArg;
  typedef TCLAP::ValueArg< TScalar > _TRealArg;

  // Parse input line
  _TStringArg in( "i", "input", "Input image", true, "", "file" );
  _TStringArg out( "o", "output", "Output image", true, "", "file" );
  _TRealArg s( "s", "sigma", "Sigma", false, 0.5, "value (0.5)" );
  _TRealArg a1( "a", "alpha1", "Alpha1", false, 0.5, "value (0.5)" );
  _TRealArg a2( "b", "alpha2", "Alpha2", false, 2, "value (2)" );
  try
  {
    TCLAP::CmdLine cmd( "Vesselness computation", ' ', "1.0.0" );
    cmd.add( a2 );
    cmd.add( a1 );
    cmd.add( s );
    cmd.add( out );
    cmd.add( in );
    cmd.parse( argc, argv );
  }
  catch( TCLAP::ArgException& err )
  {
    std::cerr
      << "===============================" << std::endl
      << "Error caught: " << std::endl
      << err.error( ) << " " << err.argId( )
      << "===============================" << std::endl
      << std::endl;
    return( 1 );

  } // yrt

  try
  {
    // Read image
    TImage::Pointer input_image;
    CTBronchi::ReadImage( input_image, in.getValue( ) );

    // Min-max
    typedef itk::MinimumMaximumImageCalculator< TImage > _TMinMax;
    _TMinMax::Pointer minMax = _TMinMax::New( );
    minMax->SetImage( input_image );
    minMax->Compute( );

    // Invert intensities
    typedef itk::InvertIntensityImageFilter< TImage > _TInverter;
    _TInverter::Pointer inverter = _TInverter::New( );
    inverter->SetInput( input_image );
    inverter->SetMaximum( minMax->GetMaximum( ) );
    double t = CTBronchi::MeasureTime( inverter );
    std::cout << "Inversion executed in " << t << " s" << std::endl;

    // Compute hessian image
    typedef itk::HessianRecursiveGaussianImageFilter< TImage > _THessian;
    _THessian::Pointer hessian = _THessian::New( );
    hessian->SetInput( inverter->GetOutput( ) );
    hessian->SetSigma( s.getValue( ) );
    t = CTBronchi::MeasureTime( hessian );
    std::cout << "Hessian executed in " << t << " s" << std::endl;

    // Vesselness
    typedef
      itk::Hessian3DToVesselnessMeasureImageFilter< TScalar >
      _TVesselness;
    _TVesselness::Pointer vesselness = _TVesselness::New( );
    vesselness->SetInput( hessian->GetOutput( ) );
    vesselness->SetAlpha1( a1.getValue( ) );
    vesselness->SetAlpha2( a2.getValue( ) );
    t = CTBronchi::MeasureTime( vesselness );
    std::cout << "Vesselness executed in " << t << " s" << std::endl;

    // Write result
    CTBronchi::WriteImage( vesselness->GetOutput( ), out.getValue( ) );
  }
  catch( std::exception& err )
  {
    std::cerr
      << "===============================" << std::endl
      << "Error caught: " << std::endl
      << err.what( )
      << "===============================" << std::endl
      << std::endl;
    return( 1 );

  } // yrt
  return( 0 );
}

// eof - $RCSfile$