Difference between revisions of "ITK/Examples/SpectralAnalysis/CrossCorrelationInFourierDomain"

From KitwarePublic
< ITK‎ | Examples
Jump to navigationJump to search
(Created page with "==VnlFFTRealToComplexConjugateImageFilter.cxx== <source lang="cpp"> #include "itkImage.h" #include "itkImageFileReader.h" #include "itkImageFileWriter.h" #include "itkVnlFFTRealT...")
 
Line 1: Line 1:
==VnlFFTRealToComplexConjugateImageFilter.cxx==
==CrossCorrelationInFourierDomain.cxx==
<source lang="cpp">
<source lang="cpp">
#include "itkImage.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageFileWriter.h"
#include "itkVnlFFTRealToComplexConjugateImageFilter.h"
#include "itkVnlFFTRealToComplexConjugateImageFilter.h"
#include "itkVnlFFTComplexConjugateToRealImageFilter.h"
#include "itkVnlFFTComplexConjugateToRealImageFilter.h"
#include "itkComplexToRealImageFilter.h"
#include "itkComplexToImaginaryImageFilter.h"
#include "itkRealAndImaginaryToComplexImageFilter.h"
#include "itkMultiplyByConstantImageFilter.h"
#include "itkMultiplyImageFilter.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;


int main( int argc, char * argv [] )
   if( argc < 3 )
{
   if( argc < 4 )
     {
     {
     std::cerr << "Usage: " << argv[0] << " inputScalarImage  inputMaskImage";
    std::cerr << "Missing Parameters " << std::endl;
     std::cerr << " outputFilteredImage" << 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];


   typedef float InputPixelType;
  // Read the input images
   const unsigned int Dimension = 2;
   typedef itk::ImageFileReader< ImageType > ImageReaderType;
 
   ImageReaderType::Pointer  fixedImageReader  = ImageReaderType::New();
   typedef itk::Image< InputPixelType, Dimension > InputImageType;
  fixedImageReader->SetFileName( fixedImageFilename );
 
   fixedImageReader->Update();
   typedef itk::ImageFileReader< InputImageType >   InputReaderType;
 
   ImageReaderType::Pointer movingImageReader = ImageReaderType::New();
  movingImageReader->SetFileName( movingImageFilename );
  movingImageReader->Update();


   InputReaderType::Pointer inputReader1 = InputReaderType::New();
   // Shift the input images
   InputReaderType::Pointer inputReader2 = InputReaderType::New();
  typedef itk::FFTShiftImageFilter< ImageType, ImageType > FFTShiftFilterType;
 
   FFTShiftFilterType::Pointer fixedFFTShiftFilter = FFTShiftFilterType::New();
   inputReader1->SetFileName( argv[1] );
   fixedFFTShiftFilter->SetInput(fixedImageReader->GetOutput());
   inputReader2->SetFileName( argv[2] );
   fixedFFTShiftFilter->Update();


  FFTShiftFilterType::Pointer movingFFTShiftFilter = FFTShiftFilterType::New();
  movingFFTShiftFilter->SetInput(movingImageReader->GetOutput());
  movingFFTShiftFilter->Update();
 
  // Compute the FFT of the input
   typedef itk::VnlFFTRealToComplexConjugateImageFilter<
   typedef itk::VnlFFTRealToComplexConjugateImageFilter<
                                  InputPixelType, Dimension >  FFTFilterType;
                                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;


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


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


   typedef FFTFilterType::OutputImageType    SpectralImageType;
  // 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,
   typedef itk::MultiplyImageFilter< SpectralImageType,
                                    SpectralImageType,
                                  SpectralImageType,
                                    SpectralImageType >  MultiplyFilterType;
                                  SpectralImageType >  MultiplyFilterType;
 
   MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New();
   MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New();
  multiplyFilter->SetInput1( fixedFFTFilter->GetOutput() );
  multiplyFilter->SetInput2( conjugateFilter->GetOutput() );


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


   typedef itk::ImageFileWriter< InputImageType > WriterType;
  // Write the spectrum
   WriterType::Pointer writer = WriterType::New();
  typedef unsigned char OutputPixelType;
   writer->SetFileName( argv[3] );
  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
   writer->SetInput( fftInverseFilter->GetOutput() );
  typedef itk::RescaleIntensityImageFilter< ImageType,  OutputImageType > RescaleFilterType;
 
   typedef itk::ImageFileWriter< OutputImageType > WriterType;
   try
  RescaleFilterType::Pointer  rescaler =  RescaleFilterType::New();
    {
   WriterType::Pointer writer = WriterType::New();
    writer->Update();
   writer->SetFileName( "CrossCorr.png" );
    }
   rescaler->SetInput( fftInverseFilter->GetOutput() );
   catch( itk::ExceptionObject & excp )
   rescaler->SetOutputMinimum(0);
    {
   rescaler->SetOutputMaximum(255);
    std::cerr << "Error writing the real image: " << std::endl;
  writer->SetInput( rescaler->GetOutput() );
    std::cerr << excp << std::endl;
  writer->Update();
    return EXIT_FAILURE;
    }


   return EXIT_SUCCESS;
   return EXIT_SUCCESS;

Revision as of 16:21, 25 January 2011

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)