Inverse-FFT of an image

Question

1.00/5 (1 vote)

See more:

I am using the following code to do the Inverse Fast Fourier Transform (I-FFT) of an Image.

C#

public static Complex[,] InverseFourierTransform2d(Complex[,] image)
    {
        int Width = image.GetLength(0);
        int Height = image.GetLength(1);

        FourierTransform.FastFourierTransform2d(image, Direction.Backward);

        for (int y = 0; y < Height; y++)
        {
            for (int x = 0; x < Width; x++)
            {
                if (((x + y) & 0x1) != 0)
                {
                    image[x, y] *= -1;
                }
            }
        }

        return image;
    }

(The DotNetFiddle link of the complete source code)[^]

Output:

http://i.stack.imgur.com/AmSCT.png[^]

enter image description here

As you can see, FFT is working correctly, but, I-FFT isn't?

What could be the problem?

What I have tried:

C#

using System;
using System.Numerics;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;
using System.Windows.Forms;

namespace Fourier_Transform__Test
{
	public enum Direction
	{
		Forward = 1,
		Backward = -1
	}

	;
	public class Tools
	{
		public static int Pow2(int power)
		{
			return ((power >= 0) && (power <= 30)) ? (1 << power) : 0;
		}

		public static int Log2(int x)
		{
			if (x <= 65536)
			{
				if (x <= 256)
				{
					if (x <= 16)
					{
						if (x <= 4)
						{
							if (x <= 2)
							{
								if (x <= 1)
									return 0;
								return 1;
							}

							return 2;
						}

						if (x <= 8)
							return 3;
						return 4;
					}

					if (x <= 64)
					{
						if (x <= 32)
							return 5;
						return 6;
					}

					if (x <= 128)
						return 7;
					return 8;
				}

				if (x <= 4096)
				{
					if (x <= 1024)
					{
						if (x <= 512)
							return 9;
						return 10;
					}

					if (x <= 2048)
						return 11;
					return 12;
				}

				if (x <= 16384)
				{
					if (x <= 8192)
						return 13;
					return 14;
				}

				if (x <= 32768)
					return 15;
				return 16;
			}

			if (x <= 16777216)
			{
				if (x <= 1048576)
				{
					if (x <= 262144)
					{
						if (x <= 131072)
							return 17;
						return 18;
					}

					if (x <= 524288)
						return 19;
					return 20;
				}

				if (x <= 4194304)
				{
					if (x <= 2097152)
						return 21;
					return 22;
				}

				if (x <= 8388608)
					return 23;
				return 24;
			}

			if (x <= 268435456)
			{
				if (x <= 67108864)
				{
					if (x <= 33554432)
						return 25;
					return 26;
				}

				if (x <= 134217728)
					return 27;
				return 28;
			}

			if (x <= 1073741824)
			{
				if (x <= 536870912)
					return 29;
				return 30;
			}

			return 31;
		}

		public static bool IsPowerOf2(int x)
		{
			return (x > 0) ? ((x & (x - 1)) == 0) : false;
		}
	}

	public static class FourierTransform
	{
		private static void FastFourierTransform1d(Complex[] data, Direction direction)
		{
			int n = data.Length;
			int m = Tools.Log2(n);
			// reorder data first
			ReorderData(data);
			// compute FFT
			int tn = 1, tm;
			for (int k = 1; k <= m; k++)
			{
				Complex[] rotation = FourierTransform.GetComplexRotation(k, direction);
				tm = tn;
				tn <<= 1;
				for (int i = 0; i < tm; i++)
				{
					Complex t = rotation[i];
					for (int even = i; even < n; even += tn)
					{
						int odd = even + tm;
						Complex ce = data[even];
						Complex co = data[odd];
						double tr = co.Real * t.Real - co.Imaginary * t.Imaginary;
						double ti = co.Real * t.Imaginary + co.Imaginary * t.Real;
						data[even] += new Complex(tr, ti);
						data[odd] = new Complex(ce.Real - tr, ce.Imaginary - ti);
					}
				}
			}

			if (direction == Direction.Forward)
			{
				for (int i = 0; i < data.Length; i++)
					data[i] /= (double)n;
			}
		}

		public static void FastFourierTransform2d(Complex[, ] data, Direction direction)
		{
			int k = data.GetLength(0);
			int n = data.GetLength(1);
			// check data size
			if (!Tools.IsPowerOf2(k) || !Tools.IsPowerOf2(n))
				throw new ArgumentException("The matrix rows and columns must be a power of 2.");
			if (k < minLength || k > maxLength || n < minLength || n > maxLength)
				throw new ArgumentException("Incorrect cInputFft length.");
			// process rows
			var row = new Complex[n];
			for (int i = 0; i < k; i++)
			{
				// copy row
				for (int j = 0; j < row.Length; j++)
					row[j] = data[i, j];
				// transform it
				FourierTransform.FastFourierTransform1d(row, direction);
				// copy back
				for (int j = 0; j < row.Length; j++)
					data[i, j] = row[j];
			}

			// process columns
			var col = new Complex[k];
			for (int j = 0; j < n; j++)
			{
				// copy column
				for (int i = 0; i < k; i++)
					col[i] = data[i, j];
				// transform it
				FourierTransform.FastFourierTransform1d(col, direction);
				// copy back
				for (int i = 0; i < k; i++)
					data[i, j] = col[i];
			}
		}

#region Private Region
		private const int minLength = 2;
		private const int maxLength = 16384;
		private const int minBits = 1;
		private const int maxBits = 14;
		private static int[][] reversedBits = new int[maxBits][];
		private static Complex[, ][] complexRotation = new Complex[maxBits, 2][];
		// Get array, indicating which data members should be swapped before FFT
		private static int[] GetReversedBits(int numberOfBits)
		{
			if ((numberOfBits < minBits) || (numberOfBits > maxBits))
				throw new ArgumentOutOfRangeException();
			// check if the array is already calculated
			if (reversedBits[numberOfBits - 1] == null)
			{
				int n = Tools.Pow2(numberOfBits);
				int[] rBits = new int[n];
				// calculate the array
				for (int i = 0; i < n; i++)
				{
					int oldBits = i;
					int newBits = 0;
					for (int j = 0; j < numberOfBits; j++)
					{
						newBits = (newBits << 1) | (oldBits & 1);
						oldBits = (oldBits >> 1);
					}

					rBits[i] = newBits;
				}

				reversedBits[numberOfBits - 1] = rBits;
			}

			return reversedBits[numberOfBits - 1];
		}

		// Get rotation of complex number
		private static Complex[] GetComplexRotation(int numberOfBits, Direction direction)
		{
			int directionIndex = (direction == Direction.Forward) ? 0 : 1;
			// check if the array is already calculated
			if (complexRotation[numberOfBits - 1, directionIndex] == null)
			{
				int n = 1 << (numberOfBits - 1);
				double uR = 1.0;
				double uI = 0.0;
				double angle = System.Math.PI / n * (int)direction;
				double wR = System.Math.Cos(angle);
				double wI = System.Math.Sin(angle);
				double t;
				Complex[] rotation = new Complex[n];
				for (int i = 0; i < n; i++)
				{
					rotation[i] = new Complex(uR, uI);
					t = uR * wI + uI * wR;
					uR = uR * wR - uI * wI;
					uI = t;
				}

				complexRotation[numberOfBits - 1, directionIndex] = rotation;
			}

			return complexRotation[numberOfBits - 1, directionIndex];
		}

		// Reorder data for FFT using
		private static void ReorderData(Complex[] data)
		{
			int len = data.Length;
			// check data length
			if ((len < minLength) || (len > maxLength) || (!Tools.IsPowerOf2(len)))
				throw new ArgumentException("Incorrect cInputFft length.");
			int[] rBits = GetReversedBits(Tools.Log2(len));
			for (int i = 0; i < len; i++)
			{
				int s = rBits[i];
				if (s > i)
				{
					Complex t = data[i];
					data[i] = data[s];
					data[s] = t;
				}
			}
		}
#endregion
	}

	public class FourierImageOperation
	{
		public static Complex[, ] FastFourierTransform2d(Complex[, ] image)
		{
			Complex[, ] comp = (Complex[, ])image.Clone();
			int Width = comp.GetLength(0);
			int Height = comp.GetLength(1);
			for (int y = 0; y < Height; y++)
			{
				for (int x = 0; x < Width; x++)
				{
					if (((x + y) & 0x1) != 0)
					{
						comp[x, y] *= -1;
					}
				}
			}

			FourierTransform.FastFourierTransform2d(comp, Direction.Forward);
			return comp;
		}

		public static Complex[, ] InverseFourierTransform2d(Complex[, ] image)
		{
			int Width = image.GetLength(0);
			int Height = image.GetLength(1);
			FourierTransform.FastFourierTransform2d(image, Direction.Backward);
			for (int y = 0; y < Height; y++)
			{
				for (int x = 0; x < Width; x++)
				{
					if (((x + y) & 0x1) != 0)
					{
						image[x, y] *= -1;
					}
				}
			}

			return image;
		}
	}

	public partial class MainForm : Form
	{
		string path = @"lena.png";
		public MainForm()
		{
			InitializeComponent();
		}

		private void loadImageButton_Click(object sender, EventArgs e)
		{
			Bitmap bitmap = Bitmap.FromFile(path) as Bitmap;
			originalImagePictureBox.Image = bitmap as Image;
		}

		private void grayscaleButton_Click(object sender, EventArgs e)
		{
			Bitmap originalImage = originalImagePictureBox.Image as Bitmap;
			Bitmap grayscale = Grayscale.ToGrayscale(originalImage);
			grayscaleImagePictureBox.Image = grayscale as Image;
		}

		private void fftButton_Click(object sender, EventArgs e)
		{
			Bitmap grayscale = grayscaleImagePictureBox.Image as Bitmap;
			Complex[, ] cGrayscale = ImageDataConverter.ToComplex(grayscale);
			Complex[, ] cFFt = FourierImageOperation.FastFourierTransform2d(cGrayscale);
			fftPictureBox.Image = ImageDataConverter.ToBitmap(cFFt, true);
		}

		private void ifftButton_Click(object sender, EventArgs e)
		{
			Bitmap fftImage = fftPictureBox.Image as Bitmap;
			Complex[, ] image = ImageDataConverter.ToComplex(fftImage);
			Complex[, ] cGrayscale = FourierImageOperation.InverseFourierTransform2d(image);
			ifftPictureBox.Image = ImageDataConverter.ToBitmap(image, false);
		}
	}
}

Posted 20-Jul-16 14:13pm

User 6976447

Updated 20-Jul-16 21:41pm

Patrice T

v2

Add a Solution

Comments

[no name] 21-Jul-16 3:38am

1. Why not use OpenCV?
2. Where are originalImagePictureBox, ImageDataConverter &etc (in dotnetfiddle)?

[no name] 21-Jul-16 11:13am

1. Learning curve is steep. Lack of coherence with the philosophy of C# and .NET.

2. https://dotnetfiddle.net/sakyrj

https://dotnetfiddle.net/E4kCMe

[no name] 21-Jul-16 20:44pm

"Learning curve is steep." - Not really for simple jobs. It is well tested and plenty of examples on the web. http://www.codeproject.com/Articles/257502/Creating-Your-First-EMGU-Image-Processing-Project
"Lack of coherence with the philosophy of C# and .NET." - Meaning? Emgu CV is used in many professional applications.
If you are planning to go into image processing this is a good path.

2 solutions

Add a Solution

Add your solution here

Treat my content as plain text, not as HTML

Preview 0

…

Existing Members

Sign in to your account

...or Join us

Download, Vote, Comment, Publish.

Your Email
Password
Forgot your password?

Your Email
This email is in use. Do you need your password?
Optional Password

I have read and agree to the Terms of Service and Privacy Policy
Please subscribe me to the CodeProject newsletters

When answering a question please:

Read the question carefully.
Understand that English isn't everyone's first language so be lenient of bad spelling and grammar.
If a question is poorly phrased then either ask for clarification, ignore it, or edit the question and fix the problem. Insults are not welcome.
Don't tell someone to read the manual. Chances are they have and don't get it. Provide an answer or move on to the next question.

Let's work to help developers, not make them feel stupid.

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Patrice T · Answer 1 · 2016-07-20T19:05:00

Solution 1

Your Log2 function look weird.
Last time I checked, Log2(9) was not 4.

Posted 20-Jul-16 19:05pm

Patrice T

Comments

[no name] 21-Jul-16 11:14am

I copied it from Accord.NET Framework.

Jochen Arndt · Answer 2 · 2016-07-20T21:41:00

I don't want to analyse the whole FFT code. But it seems that the author of the code has not understand the forward, backward, and inverse terms when seeing this code part in the function FastFourierTransform1d:

C#

if (direction == Direction.Forward)
{
    for (int i = 0; i < data.Length; i++)
        data[i] /= (double)n;
}

Dividing the vector by n is the operation for an inverse FFT after having performed a forward FFT. So this implementatione seems to perform an inverse FFT when passing Direction.Forward.

To perform an inverse FFT, do a forward FFT and divide the vector by n.

A backward FFT is the unscaled version of an inverse FFT. It can be used when the scale does not care because it is faster (dividing the vector is not necessary).

This link provides a good explantion of the differences: 12 Fast Fourier Transforms[^]

So you have two options:

Let the FFT code unchanged, pass Direction.Forward and you should get the inverse FTT.
Remove the above code from the FFT code, pass Direction.Forward, and perform the division by n to get the inverse FFT.

Overall I would not trust such code posted somewhere by anonymous.