#ifndef __FPA__IMAGE__FUNCTORS__IMAGEABSOLUTEDIFFERENCECOSTFUNCTION__H__
#define __FPA__IMAGE__FUNCTORS__IMAGEABSOLUTEDIFFERENCECOSTFUNCTION__H__

#include <fpa/Image/Functors/ImageCostFunction.h>

namespace fpa
{
  namespace Image
  {
    namespace Functors
    {
      /**
       */
      template< class I, class R >
      class ImageAbsoluteDifferenceCostFunction
        : public fpa::Image::Functors::ImageCostFunction< I, R >
      {
      public:
        /// Type-related and pointers
        typedef ImageAbsoluteDifferenceCostFunction                   Self;
        typedef fpa::Image::Functors::ImageCostFunction< I, R > Superclass;
        typedef itk::SmartPointer< Self >                          Pointer;
        typedef itk::SmartPointer< const Self >               ConstPointer;

        typedef typename Superclass::TInputImage TInputImage;
        typedef typename Superclass::TResult     TResult;
        typedef typename Superclass::TIndex      TIndex;

      public:
        itkNewMacro( Self );
        itkTypeMacro( ImageAbsoluteDifferenceCostFunction, itkObject );

      public:
        virtual R Evaluate( const TIndex& v, const TIndex& p ) const
          {
            typedef typename I::PixelType                    _TPixel;
            typedef typename itk::NumericTraits< _TPixel >   _TTraits;
            typedef typename _TTraits::MeasurementVectorType _TVector;
            typedef typename _TTraits::ValueType             _TValue;

            if( this->m_InputImage.IsNotNull( ) )
            {
              _TPixel pix = this->m_InputImage->GetPixel( v );
              pix        -= this->m_InputImage->GetPixel( p );
              if( typeid( _TPixel ) != typeid( _TValue ) )
              {
                _TVector* a = reinterpret_cast< _TVector* >( &pix );
                unsigned int n =
                  this->m_InputImage->GetNumberOfComponentsPerPixel( );
                double sum = double( 0 );
                for( unsigned int i = 0; i < n; ++i )
                  sum += ( double( ( *a )[ i ] ) * double( ( *a )[ i ] ) );
                return( R( std::sqrt( sum ) ) );
              }
              else
              {
                double pix_val =
                  double( *( reinterpret_cast< _TValue* >( &pix ) ) );
                return( R( std::fabs( pix_val ) ) );
              }
            }
            else
              return( std::numeric_limits< R >::max( ) );
          }

      protected:
        ImageAbsoluteDifferenceCostFunction( )
          : Superclass( )
          { }
        virtual ~ImageAbsoluteDifferenceCostFunction( )
          { }

      private:
        // Purposely not implemented
        ImageAbsoluteDifferenceCostFunction( const Self& );
        void operator=( const Self& );
      };

    } // ecapseman

  } // ecapseman

} // ecapseman

#endif // __FPA__IMAGE__FUNCTORS__IMAGEABSOLUTEDIFFERENCECOSTFUNCTION__H__

// eof - $RCSfile$