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:
https://i.stack.imgur.com/xBOjN.png
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).
--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
