#include "itkTransform.h"
#include "list"

namespace itk
{
/** Class for serial transformations.
 * This class represents serial transformations, i.e., the concatenation of several concrete
 * coordinate transformations of arbitrary transformation types.
 *
 * The two constraints on the types of the component transformations are that a) for all of them
 * the dimension of input and output space must be equal (so they can be concatenated), and
 * b) they must all operate on data of the same scalar type.
 *
 * This class can be used to concatenate multiple transformations for use in the 
 * ResampleImageFilter. At present, it can NOT be used directly for registration, as it provides
 * no access to the parameter vectors of the component transformations.
 *
 *\author Torsten Rohlfing, Neuroscience Program, SRI International.
 *
 * Work on this class was funded by the National Institutes on Alcohol Abuse and Alcoholism
 * through the Integrative Neuroscience Initiative on Alcoholism (INIA) under 
 * Grant No. AA013521 (PI: A. Pfefferbaum).
 *
 */
template <class TScalarType,unsigned int NDimensions=3>
class ITK_EXPORT  SerialTransform  : public Transform<TScalarType,NDimensions,NDimensions>
{
public:
  /** Standard class typedefs. */
  typedef SerialTransform  Self;
  typedef Transform<TScalarType,NDimensions,NDimensions> Superclass;
  typedef Transform<TScalarType,NDimensions,NDimensions> TransformType;

  typedef SmartPointer< Self >   Pointer;
  typedef SmartPointer< const Self >  ConstPointer;
  
  typedef SmartPointer<Superclass> TransformPointer;
  typedef std::list<TransformPointer> TransformPointerList;
  typedef typename TransformPointerList::const_iterator TransformPointerListConstIterator;
  typedef std::list<TransformBase::Pointer> TransformBasePointerList;

  typedef typename Superclass::InputPointType InputPointType;
  typedef typename Superclass::OutputPointType OutputPointType;

  typedef typename Superclass::InputVectorType InputVectorType;
  typedef typename Superclass::OutputVectorType OutputVectorType;

  typedef typename Superclass::InputVnlVectorType InputVnlVectorType;
  typedef typename Superclass::OutputVnlVectorType OutputVnlVectorType;

  typedef typename Superclass::InputCovariantVectorType InputCovariantVectorType;
  typedef typename Superclass::OutputCovariantVectorType OutputCovariantVectorType;

  typedef typename Superclass::JacobianType JacobianType;
  
  /** New method for creating an object using a factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro( SerialTransform, Transform );

  /** Add a transform to list. */
  void AddTransform( const TransformPointer& transform )
  {
    this->m_TransformList.push_back( transform );
  }

  /** Add list of transforms to list. */
  void AddTransformList( const TransformBasePointerList& transformList )
  {
    for ( typename TransformBasePointerList::const_iterator it = transformList.begin(); it != transformList.end(); ++it )
      this->m_TransformList.push_back( dynamic_cast<TransformType*>( it->GetPointer() ) );
  }

  /** The serial transform is linear if and only if all component transforms are linear. */
  virtual bool IsLinear() const 
  { 
    for ( TransformPointerListConstIterator it = this->m_TransformList.begin(); it != this->m_TransformList.end(); ++it )
      {
      if ( !(*it)->IsLinear() )
	return false;
      }
    return true; 
  }

  /** Method to transform a point. */
  virtual OutputPointType TransformPoint(const InputPointType& p ) const
  {
    OutputPointType Tp( p );
    for ( TransformPointerListConstIterator it = this->m_TransformList.begin(); it != this->m_TransformList.end(); ++it )
      {
      Tp = (*it)->TransformPoint( Tp );
      }
    return Tp;
  } 
  
  /**  Method to transform a vector. */
  virtual OutputVectorType TransformVector(const InputVectorType& p) const
  { 
    OutputVectorType Tp( p );
    for ( TransformPointerListConstIterator it = this->m_TransformList.begin(); it != this->m_TransformList.end(); ++it )
      {
      Tp = (*it)->TransformVector( Tp );
      }
    return Tp; 
  }
  
  /**  Method to transform a vnl_vector. */
  virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &p ) const
  {
    OutputVnlVectorType Tp( p );
    for ( TransformPointerListConstIterator it = this->m_TransformList.begin(); it != this->m_TransformList.end(); ++it )
      {
      Tp = (*it)->TransformVector( Tp );
      }
    return Tp;
  } 
  
  /**  Method to transform a CovariantVector. */
  virtual OutputCovariantVectorType TransformCovariantVector( const InputCovariantVectorType &p ) const
  { 
    OutputCovariantVectorType Tp( p );
    for ( TransformPointerListConstIterator it = this->m_TransformList.begin(); it != this->m_TransformList.end(); ++it )
      {
      Tp = (*it)->TransformCovariantVector( Tp );
      }
    return Tp;
  } 

  /** Compute Jacobian (product of component Jacobians). Not thread safe!! */
  virtual const JacobianType & GetJacobian(const InputPointType &p ) const
  { 
    OutputPointType Tp( p );

    this->m_Jacobian.Fill( 0.0 );
    for(unsigned int dim=0; dim < NDimensions; dim++ ) 
      {
      this->m_Jacobian( dim , dim ) = 1.0;
      }
    
    for ( TransformPointerListConstIterator it = this->m_TransformList.begin(); it != this->m_TransformList.end(); ++it )
      {
      // Q: is this the correct multiplication order?
      this->m_Jacobian *= (*it)->GetJacobian( Tp );
      Tp = (*it)->TransformPoint( Tp );
      }
    return this->m_Jacobian;
  }

protected:
    /** Default constructor. Otherwise we get a run time warning from itkTransform. */
  SerialTransform() : Superclass( NDimensions, 0 ) {}

private:
  /** List of transformations. */
  TransformPointerList m_TransformList;

  /** Temporary storage for transformation Jacobian. */
  mutable JacobianType m_Jacobian;
};

}