ITK  4.13.0
Insight Segmentation and Registration Toolkit
WikiExamples/SpectralAnalysis/VnlFFTRealToComplexConjugateImageFilter.cxx
#include "itkImage.h"
#if ITK_VERSION_MAJOR < 4
#include "itkVnlFFTRealToComplexConjugateImageFilter.h"
#else
#endif
#include <itksys/SystemTools.hxx>
#include "vnl/vnl_sample.h"
#include <math.h>
#include "QuickView.h"
int main(int argc, char*argv[])
{
// Verify input
if(argc < 2)
{
std::cerr << "Required: filename" << std::endl;
return EXIT_FAILURE;
}
// Define some types
typedef itk::Image<float, 2> FloatImageType;
// Read the image
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(argv[1]);
reader->Update();
// Compute the smallest power of two that is bigger than the image
unsigned int powerOfTwo = 0;
for(unsigned int i = 0; i < 20; i++)
{
if(pow(2.0, static_cast<double>(i)) >= reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0] &&
pow(2.0, static_cast<double>(i)) >= reader->GetOutput()->GetLargestPossibleRegion().GetSize()[1])
{
powerOfTwo = i;
break;
}
}
// Create an image bigger than the input image and that has
// dimensions which are powers of two
start.Fill(0);
size.Fill(pow(2.0,static_cast<double>(powerOfTwo)));
itk::ImageRegion<2> region(start, size);
FloatImageType::Pointer image = FloatImageType::New();
image->SetRegions(region);
image->Allocate();
image->FillBuffer(0);
// The image dimensions must be powers of two
PasteImageFilterType;
itk::Index<2> destinationIndex;
destinationIndex.Fill(0);
PasteImageFilterType::Pointer pasteFilter
= PasteImageFilterType::New ();
pasteFilter->SetSourceImage(reader->GetOutput());
pasteFilter->SetDestinationImage(image);
pasteFilter->SetSourceRegion(reader->GetOutput()->GetLargestPossibleRegion());
pasteFilter->SetDestinationIndex(destinationIndex);
pasteFilter->Update();
image->Graft(pasteFilter->GetOutput());
// Compute the FFT
#if ITK_VERSION_MAJOR < 4
typedef itk::VnlFFTRealToComplexConjugateImageFilter<FloatImageType> FFTType;
#else
#endif
FFTType::Pointer fftFilter = FFTType::New();
fftFilter->SetInput(image);
// Extract the real part
RealFilterType::Pointer realFilter = RealFilterType::New();
realFilter->SetInput(fftFilter->GetOutput());
// Extract the complex part
ImaginaryFilterType::Pointer imaginaryFilter = ImaginaryFilterType::New();
imaginaryFilter->SetInput(fftFilter->GetOutput());
// Compute the magnitude
ModulusFilterType::Pointer modulusFilter = ModulusFilterType::New();
modulusFilter->SetInput(fftFilter->GetOutput());
QuickView viewer;
viewer.AddImage(image.GetPointer(), true,
itksys::SystemTools::GetFilenameName(argv[1]));
viewer.AddImage(realFilter->GetOutput(), true,
"Real");
viewer.AddImage(imaginaryFilter->GetOutput(), true,
"Imaginary");
viewer.AddImage(modulusFilter->GetOutput(), true,
"Magnitude");
viewer.Visualize();
return EXIT_SUCCESS;
}