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

[Brian Eastwood]

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;
>> }
>>
>>
>>
>>     
>
>