ITK/Examples/SpectralAnalysis/CrossCorrelationInFourierDomain: Difference between revisions

From KitwarePublic
< ITK‎ | Examples
Jump to navigationJump to search
No edit summary
No edit summary
Line 14: Line 14:
#include "itkRescaleIntensityImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkFFTShiftImageFilter.h"
#include "itkFFTShiftImageFilter.h"
#include "itkMinimumMaximumImageCalculator.h"


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


   if( argc < 3 )
   if( argc < 3 )
Line 32: Line 34:


   // Read the input images
   // Read the input images
   typedef itk::ImageFileReader< ImageType > ImageReaderType;
   typedef itk::ImageFileReader< FloatImageType > ImageReaderType;
   ImageReaderType::Pointer  fixedImageReader  = ImageReaderType::New();
   ImageReaderType::Pointer  fixedImageReader  = ImageReaderType::New();
   fixedImageReader->SetFileName( fixedImageFilename );
   fixedImageReader->SetFileName( fixedImageFilename );
Line 42: Line 44:


   // Shift the input images
   // Shift the input images
   typedef itk::FFTShiftImageFilter< ImageType, ImageType > FFTShiftFilterType;
   typedef itk::FFTShiftImageFilter< FloatImageType, FloatImageType > FFTShiftFilterType;
   FFTShiftFilterType::Pointer fixedFFTShiftFilter = FFTShiftFilterType::New();
   FFTShiftFilterType::Pointer fixedFFTShiftFilter = FFTShiftFilterType::New();
   fixedFFTShiftFilter->SetInput(fixedImageReader->GetOutput());
   fixedFFTShiftFilter->SetInput(fixedImageReader->GetOutput());
Line 65: Line 67:
   // Take the conjugate of the fftFilterMoving
   // Take the conjugate of the fftFilterMoving
   // Extract the real part
   // Extract the real part
   typedef itk::ComplexToRealImageFilter<SpectralImageType, ImageType> RealFilterType;
   typedef itk::ComplexToRealImageFilter<SpectralImageType, FloatImageType> RealFilterType;
   RealFilterType::Pointer realFilter = RealFilterType::New();
   RealFilterType::Pointer realFilter = RealFilterType::New();
   realFilter->SetInput(movingFFTFilter->GetOutput());
   realFilter->SetInput(movingFFTFilter->GetOutput());


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


   // Flip the sign of the imaginary and combine with the real part again
   // Flip the sign of the imaginary and combine with the real part again
   typedef itk::MultiplyByConstantImageFilter<ImageType,PixelType,ImageType> MultiplyConstantFilterType;
   typedef itk::MultiplyByConstantImageFilter<FloatImageType,PixelType,FloatImageType> MultiplyConstantFilterType;
   MultiplyConstantFilterType::Pointer flipSignFilter = MultiplyConstantFilterType::New();
   MultiplyConstantFilterType::Pointer flipSignFilter = MultiplyConstantFilterType::New();
   flipSignFilter->SetConstant(-1);
   flipSignFilter->SetConstant(-1);
Line 99: Line 101:


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


  typedef itk::MinimumMaximumImageCalculator <UnsignedCharImageType>
          ImageCalculatorFilterType;
  ImageCalculatorFilterType::Pointer imageCalculatorFilter
          = ImageCalculatorFilterType::New ();
  imageCalculatorFilter->SetImage(rescaler->GetOutput());
  imageCalculatorFilter->Compute();
 
  UnsignedCharImageType::IndexType maximumLocation = imageCalculatorFilter->GetIndexOfMaximum();
  std::cout << maximumLocation << std::endl; // should be (17,15)
 
  /*
  if ypeak < size(I,1)/2 ypeak = -(ypeak-1);
  else ypeak = size(I,1) - (ypeak-1);
  end
  if xpeak < size(I,2)/2 xpeak = -(xpeak-1);
  else xpeak = size(I,2) - (xpeak-1);
  end
  */
 
   return EXIT_SUCCESS;
   return EXIT_SUCCESS;
}
}

Revision as of 16:32, 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"
  10. include "itkMinimumMaximumImageCalculator.h"

int main(int argc, char*argv[]) {

 const    unsigned int    Dimension = 2;
 typedef  float           PixelType;
 typedef itk::Image< PixelType, Dimension >  FloatImageType;
 typedef itk::Image< unsigned char, Dimension >  UnsignedCharImageType;
 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< FloatImageType  > 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< FloatImageType, FloatImageType > 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, FloatImageType> RealFilterType;
 RealFilterType::Pointer realFilter = RealFilterType::New();
 realFilter->SetInput(movingFFTFilter->GetOutput());
 // Extract the imaginary part
 typedef itk::ComplexToImaginaryImageFilter<SpectralImageType, FloatImageType> 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<FloatImageType,PixelType,FloatImageType> 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 itk::RescaleIntensityImageFilter< FloatImageType,  UnsignedCharImageType > RescaleFilterType;
 RescaleFilterType::Pointer  rescaler =  RescaleFilterType::New();
 rescaler->SetInput( fftInverseFilter->GetOutput() );
 rescaler->SetOutputMinimum(0);
 rescaler->SetOutputMaximum(255);
 rescaler->Update();
   
 typedef itk::ImageFileWriter< UnsignedCharImageType >  WriterType;
 WriterType::Pointer writer =  WriterType::New();
 writer->SetFileName( "CrossCorr.png" );
 writer->SetInput( rescaler->GetOutput() );
 writer->Update();
 typedef itk::MinimumMaximumImageCalculator <UnsignedCharImageType>
         ImageCalculatorFilterType;
 ImageCalculatorFilterType::Pointer imageCalculatorFilter
         = ImageCalculatorFilterType::New ();
 imageCalculatorFilter->SetImage(rescaler->GetOutput());
 imageCalculatorFilter->Compute();
 
 UnsignedCharImageType::IndexType maximumLocation = imageCalculatorFilter->GetIndexOfMaximum();
 std::cout << maximumLocation << std::endl; // should be (17,15)
 
 /*
 if ypeak < size(I,1)/2 ypeak = -(ypeak-1);
 else ypeak = size(I,1) - (ypeak-1);
 end
 if xpeak < size(I,2)/2 xpeak = -(xpeak-1);
 else xpeak = size(I,2) - (xpeak-1);
 end
 */
 
 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>