MantisBT - ITK
View Issue Details
IDProjectCategoryView StatusDate SubmittedLast Update
0006453ITKpublic2008-02-27 09:362008-06-04 07:59
ReporterNicolas Savoire 
Assigned ToTom Vercauteren 
PrioritynormalSeveritymajorReproducibilityalways
StatusclosedResolutionfixed 
PlatformOSOS Version
Product Version 
Target VersionFixed in Version 
Resolution Date
Sprint
Sprint Status
Summary0006453: VnlFFTRealToComplexConjugateImageFilter result is wrong
DescriptionComputation of N-d fft with VnlFFTRealToComplexConjugateImageFilter returns invalid result if N>1.
There are two issues. First VnlFFTRealToComplexConjugateImageFilter computes the conjugate of the transform (i.e. the "backward" transform).
Second vnl_fft_base expects the array dimensions with the fastest varying dimension last: if an image is nx x ny x nz then the dimensions given to vnl_fft_base (through its variable member vnl_fft_prime_factors<T> factors_[D]) should be [nz,ny,nx]. At the present time dimensions are given with the fastest varying dimension first, which is incorrect.
Steps To Reproduce
Additional InformationI included a small test case (itk_fft.cpp) which computes the 2d fft of the image:
0 1 4
9 16 25
36 49 64
81 100 121
144 169 196
The correct result (computed with matlab) is:
(1015.00,0.00) (-102.50,64.95) (-102.50,-64.95)
(-187.13,526.47) (4.62,-43.96) (40.38,-17.98)
(-307.87,124.28) (18.28,-20.30) (26.72,5.68)
(-307.87,-124.28) (26.72,-5.68) (18.28,20.30)
(-187.13,-526.47) (40.38,17.98) (4.62,43.96)
Instead the test case returns:
(1015,0) (-90.7918,-154.843) (-104.208,-36.5535)
(-104.208,36.5535) (-90.7918,154.843) (-462.5,-411.362)
(7.70046,73.265) (30.4653,33.8351) (44.5347,9.46615)
(67.2995,-29.9637) (-462.5,411.362) (67.2995,29.9637)
(44.5347,-9.46615) (30.4653,-33.8351) (7.70046,-73.265)
TagsNo tags attached.
Relationships
related to 0002161closed Brad Davis FFTDirectInverseTest fails on all platforms 
has duplicate 0003677closed Lydia Ng bug @ vnl fft wrapper 
Attached Filescpp itk_fft.cpp (1,412) 2008-02-27 09:36
https://public.kitware.com/Bug/file/1330/itk_fft.cpp
patch itk_fft.patch (1,052) 2008-02-27 09:46
https://public.kitware.com/Bug/file/1331/itk_fft.patch
Issue History
Date ModifiedUsernameFieldChange
2008-02-27 09:36Nicolas SavoireNew Issue
2008-02-27 09:36Nicolas SavoireFile Added: itk_fft.cpp
2008-02-27 09:46Nicolas SavoireNote Added: 0010654
2008-02-27 09:46Nicolas SavoireFile Added: itk_fft.patch
2008-02-27 09:47Nicolas SavoireNote Edited: 0010654
2008-05-13 04:47Tom VercauterenRelationship addedrelated to 0002161
2008-05-20 12:13kentwilliamsNote Added: 0012028
2008-05-22 05:15Tom VercauterenNote Added: 0012065
2008-05-22 11:31Tom VercauterenNote Added: 0012070
2008-05-22 11:31Tom VercauterenAssigned To => Tom Vercauteren
2008-05-22 11:31Tom VercauterenStatusnew => feedback
2008-05-22 11:47Tom VercauterenRelationship addedhas duplicate 0003677
2008-05-22 11:53kentwilliamsNote Added: 0012071
2008-06-02 13:55Tom VercauterenNote Added: 0012185
2008-06-03 04:41Tom VercauterenStatusfeedback => resolved
2008-06-03 04:41Tom VercauterenResolutionopen => fixed
2008-06-04 07:58Tom VercauterenStatusresolved => closed

Notes
(0010654)
Nicolas Savoire   
2008-02-27 09:46   
(edited on: 2008-02-27 09:47)
The joined patch (itk_fft.patch) for itkVnlFFTRealToComplexConjugateImageFilter.txx gives the correct result for the previous test case.
A similar patch should be applied to itkVnlFFTComplexConjugateToRealImageFilter.txx .
(0012028)
kentwilliams   
2008-05-20 12:13   
I originally wrote those filters, but can claim not much expertise with FFTs. Have you tried putting your patches into the library, and then running ITK testing?

In Insight/Testing/Code/Algorithms/itkFFTTest.cxx the same code (in a template function) is used to create a random image, run the FFT on it, and then the IFFT, and compare the original image with the one piped through FFT->IFFT.

As far as I know that test has been succeeding on all platforms for at least a couple of years. You need to verify that your proposed fix can at least pass this test case.
(0012065)
Tom Vercauteren   
2008-05-22 05:15   
I just submitted an experimental build on CDash. The proposed patch (with the reverse patch for itkVnlFFTComplexConjugateToRealImageFilter.txx) passes the test case.

What should be done to enforce itkFFTTest.cxx is to add dual FFT implementations tests. For example:


if((test_fft<float,1,
      itk::VnlFFTRealToComplexConjugateImageFilter<float,1> ,
      itk::FFTWComplexConjugateToRealImageFilter<float,1> >(SizeOfDimensions2)) != 0)
    {
    std::cerr << "--------------------- Failed!" << std::endl;
    rval++;;
    }


and


if((test_fft<float,1,
      itk::FFTWRealToComplexConjugateImageFilter<float,1> ,
      itk::VnlFFTComplexConjugateToRealImageFilter<float,1> >(SizeOfDimensions2)) != 0)
    {
    std::cerr << "--------------------- Failed!" << std::endl;
    rval++;;
    }
(0012070)
Tom Vercauteren   
2008-05-22 11:31   
A fix has been committed:
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Algorithms/itkVnlFFTRealToComplexConjugateImageFilter.txx?root=Insight&r1=1.10&r2=1.11&sortby=date [^]
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Algorithms/itkVnlFFTComplexConjugateToRealImageFilter.txx?root=Insight&r1=1.10&r2=1.11&sortby=date [^]

A more stringent test is still needed.
(0012071)
kentwilliams   
2008-05-22 11:53   
The cross-library test is a good idea. Unfortunately only the VNL FFT is configured in by default -- there are licensing issues with FFTW and SCSL.

There's no reason not to put together the cross-library tests, as long as they're only compiled in when both libraries are present. Of course, for testing purposes a test system could be configured with all 3 supported libraries...
(0012185)
Tom Vercauteren   
2008-06-02 13:55   
I just committed a cross library test. Unfortunately since Vnl returns the complete FFT and FFTW only half of it, a I only implemented a test to compare the real to complex conjugate case.
http://www.itk.org/cgi-bin/viewcvs.cgi/Testing/Code/Algorithms/itkFFTTest.cxx?root=Insight&r1=1.17&r2=1.18&sortby=date [^]