Using FFTW just to convert to frequency domain and back gives result quite different from the source

Question

In my testcase I'm using FFTW just to convert an image to the frequency domain and back to the spatial domain. But what I end up with is a result quite different from what I started with.

To get some questions out of the way:

• The reason I'm creating a FFTW plan on every invocation of fft() and ifft(), instead of caching them, is for simplicity. That shouldn't be a problem
• The reason I convert the real-valued array to a complex-valued array before passing it to fftw is for simplicity. I find this way is harder to get wrong than using the real-to-complex DFT functions.

My code (with utility functions, macros and classes included) follows. It uses libcinder (an old version of its) for vector classes and graphics, but even if you don't know libcinder, it will be self-explanatory. Excuse the code length, I've worked to get it pretty minimal but some utility things like my Array2D class take up quite some space.

main.cpp:

#include "precompiled.h"
typedef std::complex<float> Complexf;
using namespace ci;

void createConsole()
{
    AllocConsole();
    std::fstream* fs = new std::fstream("CONOUT$");
    std::cout.rdbuf(fs->rdbuf());
}

#define forxy(image) \
    for(Vec2i p(0, 0); p.x < image.w; p.x++) \
        for(p.y = 0; p.y < image.h; p.y++)
template<class T>
class ArrayDeleter
{
public:
    ArrayDeleter(T* arrayPtr)
    {
        refcountPtr = new int(0);
        (*refcountPtr)++;

        this->arrayPtr = arrayPtr;
    }

    ArrayDeleter(ArrayDeleter const& other)
    {
        arrayPtr = other.arrayPtr;
        refcountPtr = other.refcountPtr;
        (*refcountPtr)++;
    }

    ArrayDeleter const& operator=(ArrayDeleter const& other)
    {
        reduceRefcount();

        arrayPtr = other.arrayPtr;
        refcountPtr = other.refcountPtr;
        (*refcountPtr)++;

        return *this;
    }

    ~ArrayDeleter()
    {
        reduceRefcount();
    }

private:
    void reduceRefcount()
    {
        (*refcountPtr)--;
        if(*refcountPtr == 0)
        {
            delete refcountPtr;
            fftwf_free(arrayPtr);
        }
    }

    int* refcountPtr;
    T* arrayPtr;
};

template<class T>
struct Array2D
{
    T* data;
    int area;
    int w, h;
    ci::Vec2i Size() const { return ci::Vec2i(w, h); }
    ArrayDeleter<T> deleter;

    Array2D(Vec2i s) : deleter(Init(s.x, s.y)) { }
    Array2D() : deleter(Init(0, 0)) { }

    T* begin() { return data; }
    T* end() { return data+w*h; }

    T& operator()(Vec2i const& v) { return data[v.y*w+v.x]; }

private:
    T* Init(int w, int h) {
        data = (T*)fftwf_malloc(w * h * sizeof(T));
        area = w * h;
        this->w = w;
        this->h = h;
        return data;
    }
};

Array2D<Complexf> fft(Array2D<float> in)
{
    Array2D<Complexf> in_complex(in.Size());
    forxy(in)
    {
        in_complex(p) = Complexf(in(p));
    }

    Array2D<Complexf> result(in.Size());

    auto plan = fftwf_plan_dft_2d(in.h, in.w,
        (fftwf_complex*)in_complex.data, (fftwf_complex*)result.data, FFTW_FORWARD, FFTW_MEASURE);
    fftwf_execute(plan);
    auto mul = 1.0f / sqrt((float)result.area);
    forxy(result)
    {
        result(p) *= mul;
    }
    return result;
}

Array2D<float> ifft(Array2D<Complexf> in)
{
    Array2D<Complexf> result(in.Size());
    auto plan = fftwf_plan_dft_2d(in.h, in.w,
        (fftwf_complex*)in.data, (fftwf_complex*)result.data, FFTW_BACKWARD, FFTW_MEASURE);
    fftwf_execute(plan);

    Array2D<float> out_real(in.Size());
    forxy(in)
    {
        out_real(p) = result(p).real();
    }
    auto mul = 1.0f / sqrt((float)out_real.area);
    forxy(out_real)
    {
        out_real(p) *= mul;
    }
    return out_real;
}

gl::Texture uploadTex(Array2D<float> a)
{
    gl::Texture tex(a.w, a.h);
    tex.bind();
    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, a.w, a.h, GL_LUMINANCE, GL_FLOAT, a.data);
    return tex;
}

struct SApp : ci::app::AppBasic {
    gl::Texture texSource;
    gl::Texture texbackAndForthed;
    void setup()
    {
        createConsole();

        Array2D<float> source(Vec2i(400, 400));
        forxy(source) {
            float dist = p.distance(source.Size() / 2);
            if(dist < 100)
                source(p) = 1.0f;
            else
                source(p) = 0.0f;
        }
        printSum("source", source);
        texSource = uploadTex(source);
        setWindowSize(source.w, source.h);
        auto backAndForthed = backAndForth(source);
        printSum("backAndForthed", backAndForthed);
        //backAndForthed = backAndForth(loaded);
        texbackAndForthed = uploadTex(backAndForthed);
    }
    void printSum(std::string description, Array2D<float> arr) {
        float sum = std::accumulate(arr.begin(), arr.end(), 0.0f);
        std::cout << "array '" << description << "' has sum " << sum << std::endl;
    }
    void draw()
    {
        gl::clear(Color(0, 0, 0));
        gl::draw(texSource, Rectf(0,0, getWindowWidth() / 2, getWindowHeight() /2));
        gl::draw(texbackAndForthed, Rectf(getWindowWidth() / 2, getWindowWidth() / 2, getWindowWidth(), getWindowHeight()));
    }
    Array2D<float> backAndForth(Array2D<float> in) {
        auto inFd = fft(in);
        auto inResult = ifft(inFd);
        return inResult;
    }
};

CINDER_APP_BASIC(SApp, ci::app::RendererGl)

precompiled.h:

#include <complex>
#include <cinder/app/AppBasic.h>
#include <cinder/gl/Texture.h>
#include <cinder/gl/gl.h>
#include <cinder/Vector.h>
#include <fftw3.h>
#include <numeric>

Console output:

array 'source' has sum 31397
array 'backAndForthed' has sum -0.110077

Graphical output:

As you can see, the lower-right circle is sort of darker and gradienty.

Note: If you uncomment the second backAndForthed = backAndForth(loaded); line, the result is correct (so, the result is wrong only the first time).

Answer 1

The problem is here:

auto plan = fftwf_plan_dft_2d(in.h, in.w,
    (fftwf_complex*)in_complex.data, (fftwf_complex*)result.data,
     FFTW_FORWARD, FFTW_MEASURE);
fftwf_execute(plan);

and there:

auto plan = fftwf_plan_dft_2d(in.h, in.w,
    (fftwf_complex*)in.data, (fftwf_complex*)result.data,
    FFTW_BACKWARD, FFTW_MEASURE);
fftwf_execute(plan);

Using the flag FFTW_MEASURE implies that FFTW tries many DFT algorithms to choose the fastest one. The problem is that the input array in.data is overwritten on the way. Hence, after fftwf_execute(plan);, result.data is not the DFT of the in.data as it was before the creation of the plan. According to the documentation of FFTW on planner flags:

Important: the planner overwrites the input array during planning unless a saved plan (see Wisdom) is available for that problem, so you should initialize your input data after creating the plan. The only exceptions to this are the FFTW_ESTIMATE and FFTW_WISDOM_ONLY flags, as mentioned below.

The documentation provides a solution: using the flag FFTW_ESTIMATE ensures that the input/output arrays are not overwritten during planning.

In this particular case, using FFTW_ESTIMATE instead of FFTW_MEASURE is not expected to trigger a large increase of the computation time. Indeed, as many DFTs are computed during the plan creation, the creation of the fftw_plan using FFTW_MEASURE will be much slower than the creation using FFTW_ESTIMATE. Nevertheless, if many DFTs of the same size are to be performed, the plan must be created once for all using FFTW_MEASURE and stored. If necessary, new-array execution functions can be applied.

My guess is than you are already aware of r2c and c2r transforms, designed to cut storage space and computation time by almost 2.