// -------------------------------------------------------------------------
// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
// -------------------------------------------------------------------------

#ifndef __CPEXTENSIONS__ALGORITHMS__CPRFILTER__HXX__
#define __CPEXTENSIONS__ALGORITHMS__CPRFILTER__HXX__

#include <itkPoint.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkNearestNeighborInterpolateImageFunction.h>
#include <cpExtensions/Algorithms/BezierCurveFunction.h>

// -------------------------------------------------------------------------
template< class I, class S >
const typename cpExtensions::Algorithms::CPRFilter< I, S >::
TAxis* cpExtensions::Algorithms::CPRFilter< I, S >::
GetAxis( ) const
{
}

// -------------------------------------------------------------------------
template< class I, class S >
void cpExtensions::Algorithms::CPRFilter< I, S >::
SetAxis( const TAxis* axis )
{
}

// -------------------------------------------------------------------------
template< class I, class S >
cpExtensions::Algorithms::CPRFilter< I, S >::
CPRFilter( )
  : Superclass( )
{
  this->m_Interpolator =
    itk::NearestNeighborInterpolateImageFunction< I, S >::New( );
}

// -------------------------------------------------------------------------
template< class I, class S >
cpExtensions::Algorithms::CPRFilter< I, S >::
~CPRFilter( )
{
}

// -------------------------------------------------------------------------
template< class I, class S >
void cpExtensions::Algorithms::CPRFilter< I, S >::
GenerateData( )
{
  typedef itk::Point< S, VDimension >                         _P;
  typedef typename _P::VectorType                             _V;
  typedef cpExtensions::Algorithms::BezierCurveFunction< _V > _Bezier;
  typedef typename _Bezier::TFrame                            _M;
  typedef itk::MinimumMaximumImageCalculator< I >             _MinMax;

  auto image = this->GetInput( );
  auto axis = this->GetAxis( );

  // 0. Get image properties
  typename _MinMax::Pointer minmax = _MinMax::New( );
  minmax->SetImage( image );
  minmax->Compute( );
  auto background = minmax->GetMinimum( );

  // 1. Construct Bezier curve
  typename _Bezier::Pointer bezier = _Bezier::New( );
  auto vertices = axis->GetVertexList( );
  for( unsigned int i = 0; i < vertices->Size( ); ++i )
  {
    auto index = vertices->GetElement( i );
    _P point;
    image->TransformContinuousIndexToPhysicalPoint( index, point );
    bezier->AddPoint( point.GetVectorFromOrigin( ) );

  } // rof

  // 2. Slice image and prepare join filter
  this->m_Join = TJoin::New( );
  unsigned int nSlices = this->m_NumberOfSlices;
  if( nSlices == 0 )
    nSlices = bezier->GetNumberOfPoints( );
  _V vpre;
  _M mpre;
  S len = S( 0 );
  for( unsigned int i = 0; i <= nSlices; ++i )
  {
    // Get geometrical data and correct Frenet frame
    S u = S( i ) / S( nSlices );
    _V vnex = bezier->Evaluate( u );
    _M mnex = bezier->EvaluateFrenetFrame( u );
    if( i > 0 )
    {
      len += ( vpre - vnex ).GetNorm( );

      // TODO: Update Frenet frame orientation

    } // fi

    // Slice image
    typename TSlicer::Pointer slicer = TSlicer::New( );
    slicer->SetInput( image );
    slicer->SetRotation( mnex );
    slicer->SetTranslation( vnex );
    slicer->SetSize( S( this->m_SliceRadius ) );
    slicer->SpacingFromMinimumOn( );
    slicer->SetDefaultValue( background );
    slicer->SetInterpolator( this->m_Interpolator );

    // Add output to join filter and keep track of the current slicer
    this->m_Join->SetInput( i, slicer->GetOutput( ) );
    this->m_Slices.push_back( slicer );

    // Update values
    vpre = vnex;
    mpre = mnex;

  } // rof

  // 3. Finish join filter configuration
  this->m_Join->SetSpacing( double( len / S( nSlices + 1 ) ) );
  this->m_Join->SetOrigin( double( 0 ) );

  // 4. Graft outputs
  this->m_Join->GraftOutput( this->GetOutput( ) );
  this->m_Join->UpdateLargestPossibleRegion( );
  this->GraftOutput( this->m_Join->GetOutput( ) );
}

#endif // __CPEXTENSIONS__ALGORITHMS__CPRFILTER__HXX__

// eof - $RCSfile$