[Insight-users] Re: Possible bugs in vnl and fftw wrappers

[Luis Ibanez]

Hi Monica,

Thanks for exploring this problem further.

Could you please create a bug entry in the phpBugTracker ?

      http://public.kitware.com/Bug/index.php

Note that the bug report has a space for attaching files.

It will be great is you could add there the code that
you are posting in this email.

That will make easier for us to replicate the bug, and find
a fix for it.


Please let us know if you find any problem creating
the bug entry.


    Thanks


       Luis


==========================
mhg at unizar.es wrote:
> Dear Luis,
> 
> I have been performing some experiments with the wrappers for vnl and 
> fftw and I
> think that there are still some errors in the implementation.
> 
> 1. VNL: I have been comparing the fft results using Matlab, vtk and itk
> implementations and I found that Matlab and vtk perform similar. However,
> the itk implementation provides the CONJUGATE number of the results 
> provided,
> this is
> 
> FFT_Matlab = a+bi
> FFT_VTK = a+bi
> FTTK_ITK = a-bi
> 
> Besides, if you perform fft o fft^-1 using the pipeline
> 1. vnl_fft
> 2. vnl_ifft
> 3. flipping
> you dont get the identity. This problem may be due to the bug mentioned 
> above.
> 
> 2. FFTW: I mentioned in an earlier e-mail that this wrapper returned
> computations in half of the size of the original image. Is this purposely
> implemented?
> 
> Could you tell me if these are real mistakes or I am doing a wrong use 
> of the
> wrappers?
> 
> Thank you, Monica.
> 
> P.S. This is the toy code I am using...
> 
> #include <itkVnlFFTRealToComplexConjugateImageFilter.h>
> #include <itkVnlFFTComplexConjugateToRealImageFilter.h>
> #include <itkComplexToRealImageFilter.h>
> #include <itkComplexToImaginaryImageFilter.h>
> #include <itkFlipImageFilter.h>
> 
> using namespace std;
> 
> typedef FloatPixelType PixelType;
> typedef FloatImageType ImageType;
> typedef FloatIteratorType IteratorType;
> typedef FloatNumericTraits NumericTraits;
> typedef FloatReaderType ReaderType;
> typedef FloatWriterType WriterType;
> 
> typedef itk::VnlFFTRealToComplexConjugateImageFilter< PixelType, DIM >
> FFTFilterType;
> typedef itk::VnlFFTComplexConjugateToRealImageFilter< PixelType, DIM >
> IFFTFilterType;
> typedef itk::FlipImageFilter< ImageType > FlipFilterType;
> 
> typedef FFTFilterType::OutputImageType ComplexImageType;
> typedef ComplexImageType::PixelType ComplexPixelType;
> typedef itk::ImageRegionIterator< ComplexImageType > ComplexIteratorType;
> typedef itk::ComplexToRealImageFilter< ComplexImageType, ImageType >
> RealFilterType;
> typedef itk::ComplexToImaginaryImageFilter< ComplexImageType, ImageType >
> ImaginaryFilterType;
> 
> void usage()
> {
>  cerr << "Usage: fftDemo [input] [output]" << endl;
>  exit(1);
> }
> 
> int main( int argc, char **argv )
> {
>  if( argc < 3 )
>    usage();
> 
>  char *input_name = NULL, *output_name = NULL;
>  input_name  = argv[1];
>  argc--;
>  argv++;
>  output_name = argv[1];
>  argc--;
>  argv++;
> 
>  // Reading the images
>  ReaderType::Pointer sourceR = ReaderType::New();
>  sourceR->SetFileName( input_name );
>  sourceR->Update();
> 
>  // Allocating the output
>  ImageType::Pointer output = ImageType::New();
>  output->SetSpacing( sourceR->GetOutput()->GetSpacing() );
>  output->SetOrigin( sourceR->GetOutput()->GetOrigin() );
>  output->SetRegions( sourceR->GetOutput()->GetRequestedRegion() );
>  output->Allocate();
> 
>  // In order to get the Inverse Fourier Transform you should
>  // apply the direct transform first and then flip the output.
>  //
>  // So, your pipeline will look like:
>  //
>  //   1) A = (image/volume from disk)
>  //   2) B = VnlFFTRealToComplexConjugateImageFilter(A)
>  //   3) [do some operation on B in Fourier space]
>  //   4) C = VnlFFTComplexConjugateToRealImageFilter(B)
>  //   5) D = FlipImageFilter(C) [with all dimensions
>  //                              set to true]
> 
>  // Forward FFT filter
> 
>  FFTFilterType::Pointer fft = FFTFilterType::New();
>  fft->SetInput( sourceR->GetOutput() );
>  try
>  {
>    fft->Update();
>  }
>  catch( itk::ExceptionObject & excp )
>  {
>    std::cerr << "Error: " << std::endl;
>    std::cerr << excp << std::endl;
>  }
> 
>  //Definition of the real and imaginary part
>  RealFilterType::Pointer realFilter = RealFilterType::New();
>  realFilter->SetInput( fft->GetOutput() );
>  realFilter->Update();
> 
>  ImaginaryFilterType::Pointer imFilter = ImaginaryFilterType::New();
>  imFilter->SetInput( fft->GetOutput() );
>  imFilter->Update();
> 
>  output = imFilter->GetOutput();
> 
>  IteratorType out( output, output->GetRequestedRegion() );
> 
>  for( out.GoToBegin(); !out.IsAtEnd(); ++out )
>    cout << out.Get() << endl;
> 
>  return 0;
> }
> 
> 
>