Compute Image Spectral Density

Synopsis

Compute the magnitude of an image’s spectral components.

Results

Input image

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>;

  using ReaderType = itk::ImageFileReader<RealImageType>;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(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>;
  PadFilterType::Pointer padFilter = PadFilterType::New();
  padFilter->SetInput(reader->GetOutput());
  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;
  ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New();
  forwardFFTFilter->SetInput(padFilter->GetOutput());

  using ComplexToModulusFilterType = itk::ComplexToModulusImageFilter<ComplexImageType, RealImageType>;
  ComplexToModulusFilterType::Pointer 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>;
  WindowingFilterType::Pointer windowingFilter = WindowingFilterType::New();
  windowingFilter->SetInput(complexToModulusFilter->GetOutput());
  windowingFilter->SetWindowMinimum(0);
  windowingFilter->SetWindowMaximum(20000);

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

  using WriterType = itk::ImageFileWriter<OutputImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(outputFileName);
  writer->SetInput(fftShiftFilter->GetOutput());
  try
  {
    writer->Update();
  }
  catch (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.