ITK/Release 4/Refactor Numerical Libraries/Inventory/Fourier Transforms
The Problem
Functionalities
ITK obtains Fast Fourier Transform (FFT) functionalities from
- VXL - as software included as third party library
- FFTW - as a configuration time option
When using the VXL code that is shipped with ITK, the FFT operations are very slow. In addition, the VXL FFT produces the full complex output image when operating on a real image instead of producing only the non-redundant half of the image as FFTW does. The ITK filter providing access to VXL's forward FFT can crop the output image to produce output consistent with FFTW, but this is computationally wasteful.
FFTW is a lot faster, but it is licensed under GPL and can not be included as part of the ITK distribution.
To choose between these two options you use, at CMake configuration time, the following variables
- USE_FFTWF to use the float version of FFTW
- USE_FFTWD to use the double version of FFTW
- USE_SYSTEM_FFTW to use an FFTW library installed in your system
Architecture
- The current code is *NOT* multi-threaded...
Where is the Code
ITK classes that provide FFT filters can be found in the Module:
ITK/Modules/Filtering/FFT
They include
- itkFFTShiftImageFilter.h
- itkFFTWForwardFFTImageFilter.h
- itkFFTWInverseFFTImageFilter.h
- itkForwardFFTImageFilter.h
- itkInverseFFTImageFilter.h
- itkVnlForwardFFTImageFilter.h
- itkVnlInverseFFTImageFilter.h
- vnl_fft_3d.h
VXL
Code
The FFT code in VXL is in the directory
ITK/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo
and include the following files
- vnl_fft_1d.h
- vnl_fft_1d.txx
- vnl_fft_2d.h
- vnl_fft_2d.txx
- vnl_fft_base.h
- vnl_fft_base.txx
- vnl_fft.cxx
- vnl_fft.h
- vnl_fft_prime_factors.h
- vnl_fft_prime_factors.txx
It seems that VXL ultimately calls FFT functions from FORTRAN
CALL GPFA(A,B,TRIGS,INC,JUMP,N,LOT,ISIGN,NIPQ,INFO)
in the file
vnl_fft.h
Conditions
Image Size
The code accepts images whose number of pixels along every dimension is a multiple of powers of the following prime numbers:
- 2
- 3
- 5
That is, it should be an expression such as
2^P * 3^Q * 5^R
Image Dimensions
- Implementations are available for 1D and 2D in VXL.
- ITK added a 3D instantiation outside of VXL
Benchmark Test
Quick Tests
Quick Test that verify correctness are available here:
ITK/Release/Modules/Filtering/FFT/test
They can run the following
Test #1: itkVnlFFTTest Test #2: itkFFTShiftImageFilterTestOdd0 Test #3: itkFFTShiftImageFilterTestOdd1 Test #4: itkFFTShiftImageFilterTestEven0 Test #5: itkFFTShiftImageFilterTestEven1
Longer Tests
More demanding tests are available here
ITK/Release/Examples/Filtering/test
and can be run with
ctest -R FFT -VV
they are
Test #67: FFTImageFilterTest Test #68: FFTDirectInverseTest Test #69: FFTImageFilterTest2 Test #70: FFTImageFilterTest3
Best Tests
The most convenient code for benchmarking is in
ITK/Examples/Filtering/
They are
- FFTDirectInverse2.cxx : This one uses FFTW
- FFTDirectInverse.cxx : This one uses VXL
Notice that you still have to change the CMake variables at configuration time to switch from one to the other.
The executables of these tests end up in the binary directory:
ITK/bin
with the executable name
FFTDirectInverse
You can run it as
FFTDirectInverse inputImage outputImage
For example
FFTDirectInverse mosaic_im100010_05.png FFT_mosaic_im100010_05.png
The output image is the result of doing Forward Fourier Transform followed by an Inverse Fourier Transform, so in principle it the output image should look just like the input image.
Large Images
Moderate
Very Large
Large images for testing can be found here:
In particular:
- mosaic_im100010_05.png
This is an image of size
8768 x 8640
That FFT expand to
16384 x 16384
Doing Direct + Inverst FFT in this image takes
262.38s user 4.14s system 102% cpu 4:20.20 total
in a computer with a processor:
Intel(R) Xeon(R) CPU E5520 @ 2.27GHz
This machine is a quad-core, but... the code is using only one core. This particular machine also has 24 GB of RAM
Options
Options that have been considered:
- INTEL MKL Library
- AMD Core Math Library (ACML)
MKL
Requires commercial license for distribution
ACML
Does not require commercial license for distribution