ITK/Examples/SpectralAnalysis/CrossCorrelationInFourierDomain
From KitwarePublic
< ITK | Examples
Jump to navigationJump to search
Revision as of 16:21, 25 January 2011 by Daviddoria (talk | contribs)
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)