Compute Image Spectral Density#

Synopsis#

Compute the magnitude of an image's spectral components.

Results#

Output image of spectral density. The DC component is shifted to the center of the image.#

Code#

C++#

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkWrapPadImageFilter.h"
#include "itkForwardFFTImageFilter.h"
#include "itkComplexToModulusImageFilter.h"
#include "itkIntensityWindowingImageFilter.h"
#include "itkFFTShiftImageFilter.h"

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

  const char * inputFileName = argv[1];
  const char * outputFileName = argv[2];

  constexpr unsigned int Dimension = 3;

  using PixelType = float;
  using RealImageType = itk::Image<PixelType, Dimension>;

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

  // Some FFT filter implementations, like VNL's, need the image size to be a
  // multiple of small prime numbers.
  using PadFilterType = itk::WrapPadImageFilter<RealImageType, RealImageType>;
  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 ForwardFFTFilterType = itk::ForwardFFTImageFilter<RealImageType>;
  using ComplexImageType = ForwardFFTFilterType::OutputImageType;
  auto forwardFFTFilter = ForwardFFTFilterType::New();
  forwardFFTFilter->SetInput(padFilter->GetOutput());

  using ComplexToModulusFilterType = itk::ComplexToModulusImageFilter<ComplexImageType, RealImageType>;
  auto complexToModulusFilter = ComplexToModulusFilterType::New();
  complexToModulusFilter->SetInput(forwardFFTFilter->GetOutput());

  // Window and shift the output for visualization.
  using OutputPixelType = unsigned short;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;
  using WindowingFilterType = itk::IntensityWindowingImageFilter<RealImageType, OutputImageType>;
  auto windowingFilter = WindowingFilterType::New();
  windowingFilter->SetInput(complexToModulusFilter->GetOutput());
  windowingFilter->SetWindowMinimum(0);
  windowingFilter->SetWindowMaximum(20000);

  using FFTShiftFilterType = itk::FFTShiftImageFilter<OutputImageType, OutputImageType>;
  auto fftShiftFilter = FFTShiftFilterType::New();
  fftShiftFilter->SetInput(windowingFilter->GetOutput());

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

  return EXIT_SUCCESS;
}

Classes demonstrated#

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 FFTShiftImageFilter : public itk::CyclicShiftImageFilter<TInputImage, TOutputImage>

Shift the zero-frequency components of a Fourier transform to the center of the image.

The Fourier transform produces an image where the zero frequency components are in the corner of the image, making it difficult to understand. This filter shifts the component to the center of the image.

https://www.insight-journal.org/browse/publication/125

Note: For images with an odd-sized dimension, applying this filter twice will not produce the same image as the original one without using SetInverse(true) on one (and only one) of the two filters.

Author: Gaetan Lehmann. Biologie du Developpement et de la Reproduction, INRA de Jouy-en-Josas, France.
See: ForwardFFTImageFilter, InverseFFTImageFilter

See itk::FFTShiftImageFilter for additional documentation.