[Insight-users] Question about vnl_fft_1d

Tue Jul 6 21:02:38 EDT 2004

Hi Jiang,

Thanks for pointing this out.

This seems to be a bug in the example.

The fftCalculator should be constructed with
"spectrumSize" instead of "numberOfPoints".

spectrumSize is the smalles power of two that
is larger than the number of points.

Please find attached the fixed version,
or simply update your CVS checkou where
the fix has been made.

Please let us know if you find any further
problems,

    Thanks

      Luis

--------------
Jiang wrote:

> Dear itk-users,
> 
> I tried to use Fourier decomposition as the new example:
> 
> Insight/Examples/Numerics/
> 
>                  FourierDescriptors1.cxx
> 
> In new CVS checkout.
> 
> When I input my data, for example, the number of points is 53, I will 
> get the following
> 
> error message:
> 
> Assertion failed: !"you probably gave a signal size not of the form 2^p 
> 3^q 5^r",
> 
> file 
> F:\ITK\InsightToolkit-1.4.0\Utilities\vxl\vnl/algo/vnl_fft_prime_factors.txx, 
> line 25
> 
>  
> 
> This error happens at:
> 
>          typedef vnl_fft_1d< double > FFTCalculator;
> 
>          FFTCalculator  fftCalculator( numberOfPoints ); // If 
> numberOfPoints is 53
> 
>  
> 
>  
> 
>  
> 
> How should I do to figure out this problem?
> 
>  
> 
>  
> 
> Thank you very much!
> 
>  
> 
>  
> 
> Chunyan Jiang
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users

-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: FourierDescriptors1.cxx,v $
  Language:  C++
  Date:      $Date: 2004/07/07 00:59:55 $
  Version:   $Revision: 1.6 $

  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.

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

//  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 -1;
    }

  //  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   = ceil( log( (double)(numberOfPoints)) /
                                          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 = sqrt( real * real + imag * imag );
    std::cout << k << "  " << magnitude << std::endl;
    }
  // Software Guide : EndCodeSnippet

  return 0;
}