FFTW benchmarking, single vs. double precision transforms.

[pipe]

I did some benchmarking with FFTW for another project. I got tired of messing with the graphs, but on my somewhat standard Core2 Quad, switching from double precision FFTW to single precision FFTW would get a speed increase of around 20-40%.

The main benefit of switching to single precision is however the reduced memory requirements. Both because the DFT obviously takes half the size, but also because the transform can be performed directly on the image data if ImageMagick is compiled with HDRI. I mention this because I recently had to take the DFT of very large images, and I couldn't do it with ImageMagick because it had to create a lot of unnecessary buffers to perform the transform. The loss of precision involved in the move from doubles to floats are likely to be negligible for image operations.

Even better gains would be had switching from FFTW_ESTIMATE to FFTW_MEASURE, but that has two serious drawbacks: It will destroy the input data, and it will take a lot of time to perform (roughly 2-3 times it actually takes to execute the plan). Switching to single precision combined with FFTW_MEASURE would improve the speed by 20-150% depending on the size of the transform. There is however a way out of this, since FFTW can save these measured plans and restore them next time. Since it is likely that a user will perform the same size transform many times during the lifetime of the program, taking advantage of this would get most of the benefits without much of the drawbacks. It doesn't take much room, a couple of kB. There is also the standardized concept of a systemwide "wisdom" file for FFTW, that should be used. It's only one function call.

There's a small chance that I could possibly want to implement some or all of these changes in IM7 or IM6 in the future if they aren't there already.

[fmw42]

I know little about the inner workings of FFTW. So any one who can contribute would certainly be welcome, provided the code changes are approved by Magick.

The one major issue that has never been fixed is a limitation in Imagemagick of having square images. This arises because the initial developer (not Magick), who dropped out, had coded the FFT unfolding for square images as a start and test case, and never got back to modifying it for non-square images. Therefore if you or anyone else wants to improve the IM FFT, then this should be one of the first things done.

[pipe]

Interesting, I did not know that. I always wondered why it had to be square, because nothing in the actual transforms mandate this. I'm going to investigate this.

One minor "problem" is that all textbooks I've read about this (very very few since I'm an electrical engineer and hardly even that) almost exclusively work with square images. I suppose it makes everything easier to understand.

[fmw42]

One minor "problem" is that all textbooks I've read about this (very very few since I'm an electrical engineer and hardly even that) almost exclusively work with square images. I suppose it makes everything easier to understand.

The original FFT required power of two images. I am not sure if it needed to be square. It has been too long for me. The FFTW now by using other transforms has eliminated that restriction. So they really only need to be even (and if not, they get automatically padded with an extra row or column of black).

The IM limitation is strictly artificial, due to the original implementation as a test and was not rectified when formally folded into IM.

see
http://www.fmwconcepts.com/misc_tests/F ... index.html

where you can see the very first attempts and the phase show a lack of proper unfolding even for square images.

The current unfolding to put the center of the FFT in the center of the image is shown with this example as per the phase image

square31.png


convert square31.png -fft square31_fft_mp.png

Magnitude (amplified):
convert square31_fft_mp-0.png -auto-level -evaluate log 10000 square31_fft_mp-0_al_log10000.png


Phase:



see the more current documentation at http://www.fmwconcepts.com/imagemagick/ ... urier.html

The unfolding concept is demonstrated at http://www.mathworks.com/help/matlab/ref/fftshift.html.

The FFTW code for an MxN image only produces the right "half" transform (M/2+1)xN and the unfolding is used to fill out the FFT so that the DC (zero frequency) point is in the middle of an MxN image at (M/2),(N/2). The problem with the IM coding is that it was done only for a MXM image with center at (M/2),(M/2).

This may not be apparent in a simple round trip. But when doing deconvolutions with filters, it shows up with the image wrapped and parts lost.