Compute Forward FFT#

Synopsis#

Compute forward FFT of an image.

Results#

Code#

C++#

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkWrapPadImageFilter.h"
#include "itkForwardFFTImageFilter.h"
#include "itkComplexToRealImageFilter.h"
#include "itkComplexToImaginaryImageFilter.h"
#include "itkComplexToModulusImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc != 5)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " <InputFileName> <Real filename> <Imaginary filename> <Modulus filename>";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputFileName = argv[1];
  const char * realFileName = argv[2];
  const char * imaginaryFileName = argv[3];
  const char * modulusFileName = argv[4];

  constexpr unsigned int Dimension = 3;

  using FloatPixelType = float;
  using FloatImageType = itk::Image<FloatPixelType, Dimension>;

  const auto input = itk::ReadImage<FloatImageType>(inputFileName);

  using UnsignedCharPixelType = unsigned char;
  using UnsignedCharImageType = itk::Image<UnsignedCharPixelType, Dimension>;

  // Some FFT filter implementations, like VNL's, need the image size to be a
  // multiple of small prime numbers.
  using PadFilterType = itk::WrapPadImageFilter<FloatImageType, FloatImageType>;
  auto padFilter = PadFilterType::New();
  padFilter->SetInput(input);
  PadFilterType::SizeType padding;
  // Input size is [48, 62, 42].  Pad to [48, 64, 48].
  padding[0] = 0;
  padding[1] = 2;
  padding[2] = 6;
  padFilter->SetPadUpperBound(padding);

  using FFTType = itk::ForwardFFTImageFilter<FloatImageType>;
  auto fftFilter = FFTType::New();
  fftFilter->SetInput(padFilter->GetOutput());

  using FFTOutputImageType = FFTType::OutputImageType;

  // Extract the real part
  using RealFilterType = itk::ComplexToRealImageFilter<FFTOutputImageType, FloatImageType>;
  auto realFilter = RealFilterType::New();
  realFilter->SetInput(fftFilter->GetOutput());

  using RescaleFilterType = itk::RescaleIntensityImageFilter<FloatImageType, UnsignedCharImageType>;
  auto realRescaleFilter = RescaleFilterType::New();
  realRescaleFilter->SetInput(realFilter->GetOutput());
  realRescaleFilter->SetOutputMinimum(itk::NumericTraits<UnsignedCharPixelType>::min());
  realRescaleFilter->SetOutputMaximum(itk::NumericTraits<UnsignedCharPixelType>::max());

  try
  {
    itk::WriteImage(realRescaleFilter->GetOutput(), realFileName);
  }
  catch (const itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  // Extract the imaginary part
  using ImaginaryFilterType = itk::ComplexToImaginaryImageFilter<FFTOutputImageType, FloatImageType>;
  auto imaginaryFilter = ImaginaryFilterType::New();
  imaginaryFilter->SetInput(fftFilter->GetOutput());

  auto imaginaryRescaleFilter = RescaleFilterType::New();
  imaginaryRescaleFilter->SetInput(imaginaryFilter->GetOutput());
  imaginaryRescaleFilter->SetOutputMinimum(itk::NumericTraits<UnsignedCharPixelType>::min());
  imaginaryRescaleFilter->SetOutputMaximum(itk::NumericTraits<UnsignedCharPixelType>::max());

  try
  {
    itk::WriteImage(imaginaryRescaleFilter->GetOutput(), imaginaryFileName);
  }
  catch (const itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  // Compute the magnitude
  using ModulusFilterType = itk::ComplexToModulusImageFilter<FFTOutputImageType, FloatImageType>;
  auto modulusFilter = ModulusFilterType::New();
  modulusFilter->SetInput(fftFilter->GetOutput());

  auto magnitudeRescaleFilter = RescaleFilterType::New();
  magnitudeRescaleFilter->SetInput(modulusFilter->GetOutput());
  magnitudeRescaleFilter->SetOutputMinimum(itk::NumericTraits<UnsignedCharPixelType>::min());
  magnitudeRescaleFilter->SetOutputMaximum(itk::NumericTraits<UnsignedCharPixelType>::max());

  try
  {
    itk::WriteImage(magnitudeRescaleFilter->GetOutput(), modulusFileName);
  }
  catch (const itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated#

template<typename TInputImage, typename TOutputImage = Image<std::complex<typename TInputImage::PixelType>, TInputImage::ImageDimension>> class ForwardFFTImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>

Base class for forward Fast Fourier Transform.

This is a base class for the “forward” or “direct” discrete Fourier Transform. This is an abstract base class: the actual implementation is provided by the best child class available on the system when the object is created via the object factory system.

This class transforms a real input image into its full complex Fourier transform. The Fourier transform of a real input image has Hermitian symmetry: $f(\mathbf{x}) = f^*(-\mathbf{x})$ . That is, when the result of the transform is split in half along the x-dimension, the values in the second half of the transform are the complex conjugates of values in the first half reflected about the center of the image in each dimension.

This filter works only for real single-component input image types.

The output generated from a ForwardFFTImageFilter is in the dual space or frequency domain. Refer to FrequencyFFTLayoutImageRegionConstIteratorWithIndex for a description of the layout of frequencies generated after a forward FFT. Also see ITKImageFrequency for a set of filters requiring input images in the frequency domain.

See

InverseFFTImageFilter, FFTComplexToComplexImageFilter

ITK Sphinx Examples:

All ITK Sphinx Examples

Compute Forward FFT

Subclassed by itk::FFTWForwardFFTImageFilter< TInputImage, TOutputImage >, itk::VnlForwardFFTImageFilter< TInputImage, TOutputImage >

See itk::ForwardFFTImageFilter for additional documentation.

template<typename TInputImage, typename TOutputImage> class ComplexToRealImageFilter : public itk::UnaryGeneratorImageFilter<TInputImage, TOutputImage>: Computes pixel-wise the real(x) part of a complex image.

See itk::ComplexToRealImageFilter for additional documentation.

template<typename TInputImage, typename TOutputImage> class ComplexToImaginaryImageFilter : public itk::UnaryGeneratorImageFilter<TInputImage, TOutputImage>: Computes pixel-wise the imaginary part of a complex image.

See itk::ComplexToImaginaryImageFilter for additional documentation.

template<typename TInputImage, typename TOutputImage> class ComplexToModulusImageFilter : public itk::UnaryGeneratorImageFilter<TInputImage, TOutputImage>: Computes pixel-wise the Modulus of a complex image.

See itk::ComplexToModulusImageFilter for additional documentation.

template<typename TInputImage, typename TOutputImage> class WrapPadImageFilter : public itk::PadImageFilter<TInputImage, TOutputImage>

Increase the image size by padding with replicants of the input image value.

WrapPadImageFilter changes the image bounds of an image. Added pixels are filled in with a wrapped replica of the input image. For instance, if the output image needs a pixel that is two pixels to the left of the LargestPossibleRegion of the input image, the value assigned will be from the pixel two pixels inside the right boundary of the LargestPossibleRegion. The image bounds of the output must be specified.

../../../../_images/WrapPadImageFilter.png

This filter is implemented as a multithreaded filter. It provides a ThreadedGenerateData() method for its implementation.

See

MirrorPadImageFilter, ConstantPadImageFilter

ITK Sphinx Examples:

All ITK Sphinx Examples

Pad Image By Wrapping

See itk::WrapPadImageFilter for additional documentation.

template<typename TInputImage, typename TOutputImage = TInputImage> class RescaleIntensityImageFilter : public itk::UnaryFunctorImageFilter<TInputImage, TOutputImage, Functor::IntensityLinearTransform<TInputImage::PixelType, TOutputImage::PixelType>>

Applies a linear transformation to the intensity levels of the input Image.

RescaleIntensityImageFilter applies pixel-wise a linear transformation to the intensity values of input image pixels. The linear transformation is defined by the user in terms of the minimum and maximum values that the output image should have.

The following equation gives the mapping of the intensity values

All computations are performed in the precision of the input pixel’s RealType. Before assigning the computed value to the output pixel.

$outputPixel = ( inputPixel - inputMin) \cdot \frac{(outputMax - outputMin )}{(inputMax - inputMin)} + outputMin$

NOTE: In this filter the minimum and maximum values of the input image are computed internally using the MinimumMaximumImageCalculator. Users are not supposed to set those values in this filter. If you need a filter where you can set the minimum and maximum values of the input, please use the IntensityWindowingImageFilter. If you want a filter that can use a user-defined linear transformation for the intensity, then please use the ShiftScaleImageFilter.

See

IntensityWindowingImageFilter

ITK Sphinx Examples:

All ITK Sphinx Examples

Rescale An Image

See itk::RescaleIntensityImageFilter for additional documentation.