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

[mhg at unizar.es]

Dear Brian,

thanks for the hint. I think that this information should be included in the
wrapper documentation as the fact that fftw only returns half of the
computations should not mean that its wrapper returns half of the 
image, right?

Maybe, for future adaptations of this wrapper there should be a flag to 
let the
algorithm deal with real images considering them as real or complex depending
on the user choice.

Best, Monica.

Quoting Brian Eastwood <beastwoo at cs.unc.edu>:

> Sorry I didn't catch this message earlier--I have recently been 
> working with the fftw Fourier transform implementation.  I believe 
> that when the fftw wrapper returns an image of half the size (for 
> real-to-complex), this is completely intentional.  This is because 
> the "other half" of the image is made up of the complex conjugates of 
> the Fourier coefficients from the returned half.  (After all, the DFT 
> can't create more data than you give it.)  The term for this is 
> Hermitian symmetry.  More details are in the fftw documentation. 
> Similarly, for the complex-to-real transform, fftw expects that a 
> caller will provide half of a Hermitian symmetric data set--this is a 
> requirement for the result to be purely real. (n.b.: the Matlab 
> implementation (which is fftw under the covers) takes care of 
> rescaling your data for you, while the pure fftw libraries (and hence 
> itk implementation) don't.)
>
> I'm writing an itk::FFTWComplexToComplexImageFilter, because I have 
> need for computing the Fourier transform on non-Hermitian symmetric 
> complex data.  One could obtain the full complex result for a DFT of 
> a real image by first converting the input real-valued image to a 
> complex image (with all zero imaginary components) and using this 
> filter.  I was planning to submit the filter to ITK after I have 
> tested it.  (I'm also writing a filter to convert real-valued images 
> to complex images.)
>
> Incidentally, I am using the fftw wrappers as I noticed the same 
> problems with the vnl wrappers you did--two forward transforms do not 
> return the original image.
>
> Cheers,
> Brian Eastwood
>
> Luis Ibanez <luis.ibanez at kitware.com> wrote:
>> 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;
>>> }
>>>
>>>
>>>
>>>
>>
>>
>