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

#ifndef __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__
#define __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__

#include <vector>
#include <itkConceptChecking.h>
#include <itkObject.h>
#include <itkObjectFactory.h>

#include <itkCovariantVector.h>
#include <itkMatrix.h>
#include <itkPoint.h>
#include <itkVector.h>

namespace cpExtensions
{
  namespace Algorithms
  {
    /**
     * Recursive formulation taken from:
     * http://lmb.informatik.uni-freiburg.de/lectures/mustererkennung/Englische_Folien/07_c_ME_en.pdf
     */
    template< class _TScalar, unsigned int _VDimension >
    class IterativeGaussianModelEstimator
      : public itk::Object
    {
    public:
      typedef IterativeGaussianModelEstimator Self;
      typedef itk::Object                     Superclass;
      typedef itk::SmartPointer< Self >       Pointer;
      typedef itk::SmartPointer< const Self > ConstPointer;

      typedef _TScalar TScalar;
      itkStaticConstMacro( Dimension, unsigned int, _VDimension );

      // Begin concept checking
#ifdef ITK_USE_CONCEPT_CHECKING
      itkConceptMacro(
        ScalarTypeHasFloatResolution,
        ( itk::Concept::IsFloatingPoint< _TScalar > )
        );
#endif
      // End concept checking

      typedef itk::Matrix< TScalar, Dimension, Dimension > TMatrix;
      typedef itk::Vector< TScalar, Dimension > TVector;

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

      itkGetConstMacro( Mean, TVector );
      itkGetConstMacro( Covariance, TMatrix );
      itkGetConstMacro( UnbiasedCovariance, TMatrix );
      itkGetConstMacro( NumberOfSamples, unsigned long );

      itkSetMacro( Mean, TVector );
      itkSetMacro( Covariance, TMatrix );
      itkSetMacro( NumberOfSamples, unsigned long );

    public:
      void Clear( );

      template< class _TOtherScalar >
      void AddSample( const _TOtherScalar& x1, ... );

      template< class _TOtherScalar >
      void AddSample( const _TOtherScalar* array );

      template< class _TOtherScalar >
      void AddSample( const vnl_vector< _TOtherScalar >& v );

      template< class _TOtherScalar >
      void AddSample(
        const itk::CovariantVector< _TOtherScalar, _VDimension >& v
        );

      template< class _TOtherScalar >
      void AddSample( const itk::Point< _TOtherScalar, _VDimension >& p );

      template< class _TOtherScalar >
      void AddSample( const itk::Vector< _TOtherScalar, _VDimension >& v );

      template< class _TOtherScalar >
      _TScalar SquaredMahalanobis( const _TOtherScalar& x1, ... ) const;

      template< class _TOtherScalar >
      _TScalar SquaredMahalanobis( const _TOtherScalar* array ) const;

      template< class _TOtherScalar >
      _TScalar SquaredMahalanobis(
        const vnl_vector< _TOtherScalar >& v
        ) const;

      template< class _TOtherScalar >
      _TScalar SquaredMahalanobis(
        const itk::CovariantVector< _TOtherScalar, _VDimension >& v
        ) const;

      template< class _TOtherScalar >
      _TScalar SquaredMahalanobis(
        const itk::Point< _TOtherScalar, _VDimension >& p
        ) const;

      template< class _TOtherScalar >
      _TScalar SquaredMahalanobis(
        const itk::Vector< _TOtherScalar, _VDimension >& v
        ) const;

      template< class _TOtherScalar >
      _TScalar Density( const _TOtherScalar& x1, ... ) const;

      template< class _TOtherScalar >
      _TScalar Density( const _TOtherScalar* array ) const;

      template< class _TOtherScalar >
      _TScalar Density(
        const vnl_vector< _TOtherScalar >& v
        ) const;

      template< class _TOtherScalar >
      _TScalar Density(
        const itk::CovariantVector< _TOtherScalar, _VDimension >& v
        ) const;

      template< class _TOtherScalar >
      _TScalar Density(
        const itk::Point< _TOtherScalar, _VDimension >& p
        ) const;

      template< class _TOtherScalar >
      _TScalar Density(
        const itk::Vector< _TOtherScalar, _VDimension >& v
        ) const;

    protected:
      IterativeGaussianModelEstimator( );
      virtual ~IterativeGaussianModelEstimator( );

    protected:
      void _AddSample( const TVector& v ) const;
      _TScalar _SquaredMahalanobis( const TVector& v ) const;
      _TScalar _Density( const TVector& v ) const;

    private:
      // Purposely not implemented
      IterativeGaussianModelEstimator( const Self& other );
      void operator=( const Self& other );

    protected:
      // Recursive avg/cov values
      mutable unsigned long m_NumberOfSamples;
      mutable TVector m_Mean;
      mutable TMatrix m_Covariance;
      mutable TMatrix m_UnbiasedCovariance;
      mutable bool m_InverseUpdated;
      mutable TMatrix m_InverseUnbiasedCovariance;
    };

  } // ecapseman

} // ecapseman

#include <cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx>

#endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__

// eof - $RCSfile$