ITK/Examples/SpectralAnalysis/CrossCorrelationInFourierDomain

CrossCorrelationInFourierDomain.cxx

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkVnlFFTRealToComplexConjugateImageFilter.h"
#include "itkVnlFFTComplexConjugateToRealImageFilter.h"
#include "itkComplexToRealImageFilter.h"
#include "itkComplexToImaginaryImageFilter.h"
#include "itkRealAndImaginaryToComplexImageFilter.h"
#include "itkMultiplyByConstantImageFilter.h"
#include "itkMultiplyImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkFFTShiftImageFilter.h"

int main(int argc, char*argv[])
{
  const    unsigned int    Dimension = 2;
  typedef  float           PixelType;
  typedef itk::Image< PixelType, Dimension >  ImageType;

  if( argc < 3 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " FixedImage MovingImage"<< std::endl;;
    return EXIT_FAILURE;
    }
  std::string fixedImageFilename = argv[1];
  std::string movingImageFilename = argv[2];

  // Read the input images
  typedef itk::ImageFileReader< ImageType  > ImageReaderType;
  ImageReaderType::Pointer  fixedImageReader  = ImageReaderType::New();
  fixedImageReader->SetFileName( fixedImageFilename );
  fixedImageReader->Update();
  
  ImageReaderType::Pointer movingImageReader = ImageReaderType::New();
  movingImageReader->SetFileName( movingImageFilename );
  movingImageReader->Update();

  // Shift the input images
  typedef itk::FFTShiftImageFilter< ImageType, ImageType > FFTShiftFilterType;
  FFTShiftFilterType::Pointer fixedFFTShiftFilter = FFTShiftFilterType::New();
  fixedFFTShiftFilter->SetInput(fixedImageReader->GetOutput());
  fixedFFTShiftFilter->Update();

  FFTShiftFilterType::Pointer movingFFTShiftFilter = FFTShiftFilterType::New();
  movingFFTShiftFilter->SetInput(movingImageReader->GetOutput());
  movingFFTShiftFilter->Update();
  
  // Compute the FFT of the input
  typedef itk::VnlFFTRealToComplexConjugateImageFilter<
                                PixelType, Dimension >  FFTFilterType;
  FFTFilterType::Pointer fixedFFTFilter = FFTFilterType::New();
  fixedFFTFilter->SetInput( fixedFFTShiftFilter->GetOutput() );
  fixedFFTFilter->Update();
  
  FFTFilterType::Pointer movingFFTFilter = FFTFilterType::New();
  movingFFTFilter->SetInput( movingFFTShiftFilter->GetOutput() );
  
  typedef FFTFilterType::OutputImageType    SpectralImageType;

  // Take the conjugate of the fftFilterMoving
  // Extract the real part
  typedef itk::ComplexToRealImageFilter<SpectralImageType, ImageType> RealFilterType;
  RealFilterType::Pointer realFilter = RealFilterType::New();
  realFilter->SetInput(movingFFTFilter->GetOutput());

  // Extract the imaginary part
  typedef itk::ComplexToImaginaryImageFilter<SpectralImageType, ImageType> ImaginaryFilterType;
  ImaginaryFilterType::Pointer imaginaryFilter = ImaginaryFilterType::New();
  imaginaryFilter->SetInput(movingFFTFilter->GetOutput());

  // Flip the sign of the imaginary and combine with the real part again
  typedef itk::MultiplyByConstantImageFilter<ImageType,PixelType,ImageType> MultiplyConstantFilterType;
  MultiplyConstantFilterType::Pointer flipSignFilter = MultiplyConstantFilterType::New();
  flipSignFilter->SetConstant(-1);
  flipSignFilter->SetInput(imaginaryFilter->GetOutput());
  typedef itk::RealAndImaginaryToComplexImageFilter<PixelType,PixelType,PixelType,2> RealImagToComplexFilterType;
  RealImagToComplexFilterType::Pointer conjugateFilter = RealImagToComplexFilterType::New();
  conjugateFilter->SetInput1(realFilter->GetOutput());
  conjugateFilter->SetInput2(flipSignFilter->GetOutput());

  // The conjugate product of the spectrum
  typedef itk::MultiplyImageFilter< SpectralImageType,
                                  SpectralImageType,
                                  SpectralImageType >  MultiplyFilterType;
  MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New();
  multiplyFilter->SetInput1( fixedFFTFilter->GetOutput() );
  multiplyFilter->SetInput2( conjugateFilter->GetOutput() );

  // IFFT
  typedef itk::VnlFFTComplexConjugateToRealImageFilter<
        PixelType, Dimension >  IFFTFilterType;
  IFFTFilterType::Pointer fftInverseFilter = IFFTFilterType::New();
  fftInverseFilter->SetInput( multiplyFilter->GetOutput() );

  // Write the spectrum
  typedef unsigned char OutputPixelType;
  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
  typedef itk::RescaleIntensityImageFilter< ImageType,  OutputImageType > RescaleFilterType;
  typedef itk::ImageFileWriter< OutputImageType >  WriterType;
  RescaleFilterType::Pointer  rescaler =  RescaleFilterType::New();
  WriterType::Pointer writer =  WriterType::New();
  writer->SetFileName( "CrossCorr.png" );
  rescaler->SetInput( fftInverseFilter->GetOutput() );
  rescaler->SetOutputMinimum(0);
  rescaler->SetOutputMaximum(255);
  writer->SetInput( rescaler->GetOutput() );
  writer->Update();

  return EXIT_SUCCESS;
}

CMakeLists.txt

cmake_minimum_required(VERSION 2.6)

project(CrossCorrelationInFourierDomain)

find_package(ITK REQUIRED)
include(${ITK_USE_FILE})

add_executable(CrossCorrelationInFourierDomain CrossCorrelationInFourierDomain.cxx )
target_link_libraries(CrossCorrelationInFourierDomain ITKIO ITKCommon)