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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    FourierDescriptors1.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  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.

=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

//  Software Guide : BeginLatex
//
//  Fourier Descriptors provide a mechanism for representing a closed curve in
//  space.  The represented curve has infinite continuiity because the
//  parametric coordinate of its points are computed from a Fourier Series.
//
//  In this example we illustrate how a curve that is initially defined by a
//  set of points in space can be represented in terms for Fourier Descriptors.
//  This representation is useful for several purposes. For example, it
//  provides a mechanmism for interpolating values among the points, it
//  provides a way of analyzing the smoothness of the curve.  In this particular
//  example we will focus on this second application of the Fourier Descriptors.
//  
//  The first operation that we should use in this context is the computation
//  of the discrete fourier transform of the point coordinates. The coordinates
//  of the points are considered to be independent functions and each one is
//  decomposed in a Fourier Series. In this section we will use $t$ as the
//  parameter of the curve, and will assume that it goes from $0$ to $1$ and
//  cycles as we go around the closed curve. //
//  \begin{equation}
//  \textbf{V(t)} = \left( X(t), Y(t) \rigth)
//  \end{equation}
//  
//  We take now the functions $X(t)$, $Y(t)$ and interpret them as the
//  components of a complex number for wich we compute its discrete fourier
//  series in the form
//
//  \begin{equation}
//  V(t) = \sum_{k=0}^{N} \exp(-\frac{2 k \pi \textbf{i}}{N}) F_k
//  \end{equation}
//  
//  Where the set of coefficients $F_k$ is the discrete spectrum of the complex
//  function $V(t)$. These coefficients are in general complex numbers and both
//  their magnitude and phase must be considered in any further analysis of the
//  spectrum.
//
//  Software Guide : EndLatex 


//  Software Guide : BeginLatex
//
//  The class \code{vnl_fft_1D} is the VNL class that computes such transform.
//  In order to use it, we should include its header file first. 
//  
//  Software Guide : EndLatex 


// Software Guide : BeginCodeSnippet
#include "vnl/algo/vnl_fft_1d.h"
// Software Guide : EndCodeSnippet


#include "itkPoint.h"
#include "itkVectorContainer.h"


#include <fstream>


int main(int argc, char * argv[] )
{

  if( argc < 2 ) 
    {
    std::cerr << "Missing arguments" << std::endl;
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << "inputFileWithPointCoordinates" << std::endl;
    return EXIT_FAILURE;
    }
    
  //  Software Guide : BeginLatex
  //
  //  We should now instantiate the filter that will compute the Fourier
  //  transform of the set of coordinates.
  //
  //  Software Guide : EndLatex 


  // Software Guide : BeginCodeSnippet
  typedef vnl_fft_1d< double > FFTCalculator;
  // Software Guide : EndCodeSnippet

   
  //  Software Guide : BeginLatex
  //
  //  The points representing the curve are stored in a
  //  \doxygen{VectorContainer} of \doxygen{Point}.
  //
  //  Software Guide : EndLatex 


  // Software Guide : BeginCodeSnippet
  typedef itk::Point< double, 2 >  PointType;

  typedef itk::VectorContainer< unsigned int, PointType >  PointsContainer;

  PointsContainer::Pointer points = PointsContainer::New();
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  In this example we read the set of points from a text file.
  //
  //  Software Guide : EndLatex 


  // Software Guide : BeginCodeSnippet
  std::ifstream inputFile;
  inputFile.open( argv[1] );

  if( inputFile.fail() )
    {
    std::cerr << "Problems opening file " << argv[1] << std::endl;
    }

  unsigned int numberOfPoints;
  inputFile >> numberOfPoints;

  points->Reserve( numberOfPoints );

  typedef PointsContainer::Iterator PointIterator;
  PointIterator pointItr = points->Begin();

  PointType point;
  for( unsigned int pt=0; pt<numberOfPoints; pt++)
    {
    inputFile >> point[0] >> point[1];
    pointItr.Value() = point;
    ++pointItr;
    }
  // Software Guide : EndCodeSnippet

   
  //  Software Guide : BeginLatex
  //
  //  This class will compute the Fast Fourier transform of the input an it will
  //  return it in the same array. We must therefore copy the original data into
  //  an auxiliary array that will in its turn contain the results of the
  //  transform.
  //
  //  Software Guide : EndLatex 


  // Software Guide : BeginCodeSnippet
  typedef vcl_complex<double>   FFTCoefficientType;
  typedef vcl_vector< FFTCoefficientType > FFTSpectrumType;
  // Software Guide : EndCodeSnippet



   
  //  Software Guide : BeginLatex
  //
  // The choice of the spectrum size is very important. Here we select to use
  // the next power of two that is larger than the number of points.
  //
  //  Software Guide : EndLatex 

  // Software Guide : BeginCodeSnippet
  const unsigned int powerOfTwo   = 
    (unsigned int)vcl_ceil( vcl_log( (double)(numberOfPoints)) /
                            vcl_log( (double)(2.0)) );

  const unsigned int spectrumSize = 1 << powerOfTwo;



  //  Software Guide : BeginLatex
  //
  //  The Fourier Transform type can now be used for constructing one of such
  //  filters. Note that this is a VNL class and does not follows ITK notation
  //  for construction and assignment to SmartPointers.
  //
  //  Software Guide : EndLatex 


  // Software Guide : BeginCodeSnippet
  FFTCalculator  fftCalculator( spectrumSize );
  // Software Guide : EndCodeSnippet



  FFTSpectrumType signal( spectrumSize );  

  pointItr = points->Begin();
  for(unsigned int p=0; p<numberOfPoints; p++)
    {
    signal[p] = FFTCoefficientType( pointItr.Value()[0], pointItr.Value()[1] );
    ++pointItr;
    }
  // Software Guide : EndCodeSnippet



   
  //  Software Guide : BeginLatex
  //
  // Fill in the rest of the input with zeros. This padding may have
  // undesirable effects on the spectrum if the signal is not attenuated to
  // zero close to their boundaries. Instead of zero-padding we could have used
  // repetition of the last value or mirroring of the signal.
  //
  //  Software Guide : EndLatex 

  // Software Guide : BeginCodeSnippet
  for(unsigned int pad=numberOfPoints; pad<spectrumSize; pad++)
    {
    signal[pad] = 0.0;
    }
  // Software Guide : EndCodeSnippet




  //  Software Guide : BeginLatex
  //
  //  Now we print out the signal as it is passed to the transform calculator
  //
  //  Software Guide : EndLatex 

  // Software Guide : BeginCodeSnippet
  std::cout << "Input to the FFT transform" << std::endl;
  for(unsigned int s=0; s<spectrumSize; s++)
    {
    std::cout << s << " : ";
    std::cout << signal[s] << std::endl;
    }
  // Software Guide : EndCodeSnippet



  //  Software Guide : BeginLatex
  //
  //  The actual transform is computed by invoking the \code{fwd_transform}
  //  method in the FFT calculator class.
  //
  //  Software Guide : EndLatex 

  // Software Guide : BeginCodeSnippet
  fftCalculator.fwd_transform( signal );
  // Software Guide : EndCodeSnippet
 



  //  Software Guide : BeginLatex
  //
  //  Now we print out the results of the transform.
  //
  //  Software Guide : EndLatex 

  // Software Guide : BeginCodeSnippet
  std::cout << std::endl;
  std::cout << "Result from the FFT transform" << std::endl;
  for(unsigned int k=0; k<spectrumSize; k++)
    {
    const double real = signal[k].real();
    const double imag = signal[k].imag();
    const double magnitude = vcl_sqrt( real * real + imag * imag );
    std::cout << k << "  " << magnitude << std::endl;
    }
  // Software Guide : EndCodeSnippet


  return EXIT_SUCCESS;
}