Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > sci.image.processing > #4217

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

Show all headers | View raw

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).

Back to sci.image.processing | Previous | Next — Next in thread | Find similar | Unroll thread

Thread

Using FFTW just to convert to frequency domain and back gives result quite different from the source Stefan Monov <logixoul@gmail.com> - 2017-01-28 09:03 -0800
  Re: Using FFTW just to convert to frequency domain and back gives result quite different from the source Martin Leese <please@see.Web.for.e-mail.INVALID> - 2017-01-29 12:06 -0700
    Re: Using FFTW just to convert to frequency domain and back gives result quite different from the source Stefan Monov <logixoul@gmail.com> - 2017-01-30 04:46 -0800
      Re: Using FFTW just to convert to frequency domain and back gives result quite different from the source Stefan Monov <logixoul@gmail.com> - 2017-01-30 09:22 -0800

csiph-web