Examples/Filtering/FFTImageFilter.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
 *
 *         https://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
//
//  In this section we assume that you are familiar with Spectral Analysis, in
//  particular with the concepts of the Fourier Transform and the numerical
//  implementation of the Fast Fourier transform. If you are not familiar with
//  these concepts you may want to consult first any of the many available
//  introductory books to spectral
//  analysis~\cite{Bracewell1999,Bracewell2004}.
//
//  This example illustrates how to use the Fast Fourier Transform filter
//  (FFT) for processing an image in the spectral domain. Given that FFT
//  computation can be CPU intensive, there are multiple hardware specific
//  implementations of FFT. It is convenient in many cases to delegate the
//  actual computation of the transform to local available libraries.
//  Particular examples of those libraries are
//  fftw\footnote{https://www.fftw.org} and the VXL implementation of FFT. For
//  this reason ITK provides a base abstract class that factorizes the
//  interface to multiple specific implementations of FFT. This base class is
//  the \doxygen{ForwardFFTImageFilter}, and two of its derived classes are
//  \doxygen{VnlForwardFFTImageFilter} and
//  \doxygen{FFTWRealToComplexConjugateImageFilter}.
//
//
//  \index{itk::Forward\-FFT\-Image\-Filter}
//  \index{itk::Vnl\-Forward\-FFT\-Image\-Filter}
//  \index{itk::FFTW\-Forward\-FFT\-Image\-Filter}
//
//  Software Guide : EndLatex
 
// Software Guide : BeginLatex
//
// A typical application that uses FFT will need to include the following
// header files.
//
// Software Guide : EndLatex
 
 
// Software Guide : BeginCodeSnippet
#include "itkImage.h"
#include "itkVnlForwardFFTImageFilter.h"
#include "itkComplexToRealImageFilter.h"
#include "itkComplexToImaginaryImageFilter.h"
// Software Guide : EndCodeSnippet
 
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
 
 
int
main(int argc, char * argv[])
{
  if (argc < 5)
  {
    std::cerr << "Usage: " << argv[0]
              << " inputScalarImage  outputRealPartOfComplexImage "
                 "outputRealImaginaryPartOfComplexImage outputComplex"
              << std::endl;
  }
 
 
  // Software Guide : BeginLatex
  //
  // The first decision to make is related to the pixel type and dimension of
  // the images on which we want to compute the Fourier transform.
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using PixelType = float;
  constexpr unsigned int Dimension = 2;
 
  using ImageType = itk::Image<PixelType, Dimension>;
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // We use the same image type in order to instantiate the FFT filter, in
  // this case the \doxygen{VnlForwardFFTImageFilter}. Once the filter type is
  // instantiated, we can use it for creating one object by invoking the
  // \code{New()} method and assigning the result to a SmartPointer.
  //
  // Software Guide : EndLatex
 
 
  // Software Guide : BeginCodeSnippet
  using FFTFilterType = itk::VnlForwardFFTImageFilter<ImageType>;
 
  auto fftFilter = FFTFilterType::New();
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // The input to this filter can be taken from a reader, for example.
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using ReaderType = itk::ImageFileReader<ImageType>;
  auto reader = ReaderType::New();
  reader->SetFileName(argv[1]);
 
  fftFilter->SetInput(reader->GetOutput());
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // The execution of the filter can be triggered by invoking the
  // \code{Update()} method.  Since this invocation can eventually throw an
  // exception, the call must be placed inside a try/catch block.
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  try
  {
    fftFilter->Update();
  }
  catch (const itk::ExceptionObject & excp)
  {
    std::cerr << "Error: " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // In general the output of the FFT filter will be a complex image. We can
  // proceed to save this image in a file for further analysis. This can be
  // done by simply instantiating an \doxygen{ImageFileWriter} using the trait
  // of the output image from the FFT filter. We construct one instance of the
  // writer and pass the output of the FFT filter as the input of the writer.
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using ComplexImageType = FFTFilterType::OutputImageType;
 
  using ComplexWriterType = itk::ImageFileWriter<ComplexImageType>;
 
  auto complexWriter = ComplexWriterType::New();
  complexWriter->SetFileName(argv[4]);
 
  complexWriter->SetInput(fftFilter->GetOutput());
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // Finally we invoke the \code{Update()} method placed inside a try/catch
  // block.
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  try
  {
    complexWriter->Update();
  }
  catch (const itk::ExceptionObject & excp)
  {
    std::cerr << "Error: " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // In addition to saving the complex image into a file, we could also
  // extract its real and imaginary parts for further analysis. This can be
  // done with the \doxygen{ComplexToRealImageFilter} and the
  // \doxygen{ComplexToImaginaryImageFilter}.
  //
  // We instantiate first the ImageFilter that will help us to extract the
  // real part from the complex image.  The \code{ComplexToRealImageFilter}
  // takes as its first template parameter the type of the complex image and
  // as its second template parameter it takes the type of the output image
  // pixel. We create one instance of this filter and connect as its input the
  // output of the FFT filter.
  //
  // \index{itk::ComplexToRealImageFilter}
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using RealFilterType =
    itk::ComplexToRealImageFilter<ComplexImageType, ImageType>;
 
  auto realFilter = RealFilterType::New();
 
  realFilter->SetInput(fftFilter->GetOutput());
  // Software Guide : EndCodeSnippet
 
 
  using WritePixelType = unsigned char;
  using WriteImageType = itk::Image<WritePixelType, Dimension>;
 
 
  // Software Guide : BeginLatex
  //
  // Since the range of intensities in the Fourier domain can be quite
  // concentrated, it is convenient to rescale the image in order to
  // visualize it. For this purpose we instantiate a
  // \doxygen{RescaleIntensityImageFilter} that will rescale the intensities
  // of the \code{real} image into a range suitable for writing in a file. We
  // also set the minimum and maximum values of the output to the range of the
  // pixel type used for writing.
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using RescaleFilterType =
    itk::RescaleIntensityImageFilter<ImageType, WriteImageType>;
 
  auto intensityRescaler = RescaleFilterType::New();
 
  intensityRescaler->SetInput(realFilter->GetOutput());
 
  intensityRescaler->SetOutputMinimum(0);
  intensityRescaler->SetOutputMaximum(255);
  // Software Guide : EndCodeSnippet
 
  using WriterType = itk::ImageFileWriter<WriteImageType>;
 
  auto writer = WriterType::New();
 
  writer->SetFileName(argv[2]);
 
  writer->SetInput(intensityRescaler->GetOutput());
 
  try
  {
    writer->Update();
  }
  catch (const itk::ExceptionObject & excp)
  {
    std::cerr << "Error writing the real image: " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }
 
  // Software Guide : BeginLatex
  //
  // We can now instantiate the ImageFilter that will help us to extract the
  // imaginary part from the complex image.  The filter that we use here is
  // the \doxygen{ComplexToImaginaryImageFilter}. It takes as first template
  // parameter the type of the complex image and as second template parameter
  // it takes the type of the output image pixel. An instance of the filter is
  // created, and its input is connected to the output of the FFT filter.
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using ComplexImageType = FFTFilterType::OutputImageType;
 
  using ImaginaryFilterType =
    itk::ComplexToImaginaryImageFilter<ComplexImageType, ImageType>;
 
  auto imaginaryFilter = ImaginaryFilterType::New();
 
  imaginaryFilter->SetInput(fftFilter->GetOutput());
  // Software Guide : EndCodeSnippet
 
 
  // Software Guide : BeginLatex
  //
  // The Imaginary image can then be rescaled and saved into a file, just as
  // we did with the Real part.
  //
  // Software Guide : EndLatex
 
  intensityRescaler->SetInput(imaginaryFilter->GetOutput());
  writer->SetFileName(argv[3]);
 
  try
  {
    writer->Update();
  }
  catch (const itk::ExceptionObject & excp)
  {
    std::cerr << "Error writing the imaginary image: " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }
 
 
  // Software Guide : BeginLatex
  //
  // For the sake of illustrating the use of a \doxygen{ImageFileReader} on
  // Complex images, here we instantiate a reader that will load the Complex
  // image that we just saved. Note that nothing special is required in this
  // case. The instantiation is done just the same as for any other type of
  // image, which once again illustrates the power of Generic Programming.
  //
  // Software Guide : EndLatex
 
 
  // Software Guide : BeginCodeSnippet
  using ComplexReaderType = itk::ImageFileReader<ComplexImageType>;
 
  auto complexReader = ComplexReaderType::New();
 
  complexReader->SetFileName(argv[4]);
  complexReader->Update();
  // Software Guide : EndCodeSnippet
 
 
  // A way of testing the pixel type of an image in file is to
  // invoke the ImageIO object from the reader and then call
  // \code{GetPixelTypeAsString()}
  complexReader->GetImageIO()->GetPixelTypeAsString(
    complexReader->GetImageIO()->GetPixelType());
 
 
  return EXIT_SUCCESS;
}