ITK/Examples/SpectralAnalysis/CrossCorrelationInFourierDomain: Difference between revisions

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...")
 
No edit summary
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

<source lang="cpp">

  1. include "itkImage.h"
  2. include "itkImageFileReader.h"
  3. include "itkImageFileWriter.h"
  1. include "itkVnlFFTRealToComplexConjugateImageFilter.h"
  2. include "itkVnlFFTComplexConjugateToRealImageFilter.h"
  3. include "itkComplexToRealImageFilter.h"
  4. include "itkComplexToImaginaryImageFilter.h"
  5. include "itkRealAndImaginaryToComplexImageFilter.h"
  6. include "itkMultiplyByConstantImageFilter.h"
  7. include "itkMultiplyImageFilter.h"
  8. include "itkRescaleIntensityImageFilter.h"
  9. 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;

} </source>

CMakeLists.txt

<source lang="cmake"> 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) </source>

Retrieved from "https://public.kitware.com/Wiki/index.php?title=ITK/Examples/SpectralAnalysis/CrossCorrelationInFourierDomain&oldid=36753"