View Issue Details [ Jump to Notes ] [ Print ]
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 
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.
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.
Resolution Date
Sprint
Sprint Status
Attached Filescpp file icon itk_fft.cpp [^] (1,412 bytes) 2008-02-27 09:36
patch file icon itk_fft.patch [^] (1,052 bytes) 2008-02-27 09:46 [Show Content]

 Relationships
related to 0002161closedBrad Davis FFTDirectInverseTest fails on all platforms 
has duplicate 0003677closedLydia Ng bug @ vnl fft wrapper 

  Notes
(0010654)
Nicolas Savoire (reporter)
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 (developer)
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 (developer)
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 (developer)
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 (developer)
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 (developer)
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 [^]

 Issue History
Date Modified Username Field Change
2008-02-27 09:36 Nicolas Savoire New Issue
2008-02-27 09:36 Nicolas Savoire File Added: itk_fft.cpp
2008-02-27 09:46 Nicolas Savoire Note Added: 0010654
2008-02-27 09:46 Nicolas Savoire File Added: itk_fft.patch
2008-02-27 09:47 Nicolas Savoire Note Edited: 0010654
2008-05-13 04:47 Tom Vercauteren Relationship added related to 0002161
2008-05-20 12:13 kentwilliams Note Added: 0012028
2008-05-22 05:15 Tom Vercauteren Note Added: 0012065
2008-05-22 11:31 Tom Vercauteren Note Added: 0012070
2008-05-22 11:31 Tom Vercauteren Assigned To => Tom Vercauteren
2008-05-22 11:31 Tom Vercauteren Status new => feedback
2008-05-22 11:47 Tom Vercauteren Relationship added has duplicate 0003677
2008-05-22 11:53 kentwilliams Note Added: 0012071
2008-06-02 13:55 Tom Vercauteren Note Added: 0012185
2008-06-03 04:41 Tom Vercauteren Status feedback => resolved
2008-06-03 04:41 Tom Vercauteren Resolution open => fixed
2008-06-04 07:58 Tom Vercauteren Status resolved => closed

Copyright © 2000 - 2018 MantisBT Team