ITK  5.3.0
Insight Toolkit
Examples/Numerics/FourierDescriptors1.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
// 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) \right)
// \end{equation}
//
// We take now the functions $X(t)$, $Y(t)$ and interpret them as the
// components of a complex number for which 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 <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
using FFTCalculator = vnl_fft_1d<double>;
// 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
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);
using PointIterator = PointsContainer::Iterator;
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 and 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
using FFTCoefficientType = std::complex<double>;
using FFTSpectrumType = std::vector<FFTCoefficientType>;
// 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 auto powerOfTwo = (unsigned int)std::ceil(
std::log((double)(numberOfPoints)) / std::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 = std::sqrt(real * real + imag * imag);
std::cout << k << " " << magnitude << std::endl;
}
// Software Guide : EndCodeSnippet
return EXIT_SUCCESS;
}
itkPoint.h
itk::GTest::TypedefsAndConstructors::Dimension2::PointType
ImageBaseType::PointType PointType
Definition: itkGTestTypedefsAndConstructors.h:51
itkVectorContainer.h
itk::Point< double, 2 >
itk::VectorContainer
Define a front-end to the STL "vector" container that conforms to the IndexedContainerInterface.
Definition: itkVectorContainer.h:48