/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFourierSeriesPath.txx,v $
  Language:  C++
  Date:      $Date: 2006/03/19 04:36:58 $
  Version:   $Revision: 1.7 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkFourierSeriesPath_txx
#define __itkFourierSeriesPath_txx

#include "itkFourierSeriesPath.h"
#include <math.h>

// not all versions of math.h seem to define M_PI:
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif



namespace itk
{

template<unsigned int VDimension>
typename FourierSeriesPath<VDimension>::OutputType
FourierSeriesPath<VDimension>
::Evaluate( const InputType & input ) const
{
  InputType   theta;
  OutputType  output;
  int         numHarmonics;
  
  numHarmonics = m_CosCoefficients->Size(); 
  output.Fill(0);
  
  if( numHarmonics > 0 ) { output += m_CosCoefficients->ElementAt(0); }
  
  for(int n=1; n<numHarmonics; n++)
    {
    // input defined over [0,1] maps to theta defined over [0,2pi * n]
    theta = M_PI*2.0*n*input;
    output += ( m_CosCoefficients->ElementAt(n) * vcl_cos(theta) +
                m_SinCoefficients->ElementAt(n) * vcl_sin(theta) ) * 2.0;
    }
  
  return output;
}



template<unsigned int VDimension>
typename FourierSeriesPath<VDimension>::VectorType
FourierSeriesPath<VDimension>
::EvaluateDerivative(const InputType & input) const
{
  InputType   theta;
  VectorType  output;
  int         numHarmonics;
  
  numHarmonics = m_CosCoefficients->Size(); 
  output.Fill(0);
  
  for(int n=1; n<numHarmonics; n++)
    {
    // input defined over [0,1] maps to theta defined over [0,2pi * n]
    theta = M_PI*2.0*n*input;
    output += ( m_SinCoefficients->ElementAt(n) * vcl_cos(theta) -
                m_CosCoefficients->ElementAt(n) * vcl_sin(theta) ) * (2.0 * n);
    }
  
  return output;
}



template<unsigned int VDimension>
void
FourierSeriesPath<VDimension>
::AddHarmonic( const VectorType & CosCoefficients,
               const VectorType & SinCoefficients )
{
  unsigned int numHarmonics = m_CosCoefficients->Size();
  
  m_CosCoefficients->InsertElement(numHarmonics, CosCoefficients);
  m_SinCoefficients->InsertElement(numHarmonics, SinCoefficients);
  this->Modified();
}



/**
 * Constructor
 */
template <unsigned int VDimension>
FourierSeriesPath<VDimension>
::FourierSeriesPath()
{
  this->SetDefaultInputStepSize( 1.0/50.0 );
  m_CosCoefficients = CoefficientsType::New();
  m_SinCoefficients = CoefficientsType::New();
}



/**
 * Standard "PrintSelf" method
 */
template <unsigned int VDimension>
void
FourierSeriesPath<VDimension>
::PrintSelf( std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf( os, indent );
  os << indent << "Cos Coefficients:  " << m_CosCoefficients << std::endl;
  os << indent << "Sin Coefficients:  " << m_SinCoefficients << std::endl;
}

} // end namespaceitk

#endif
