SolSpider – New version

Share

Check out the update for solspider in the Google Play store.

It now supports Android 4+ and has smoother movement.

SolSpider is open-source. Get the source here.

Google Play Store Link here.

Share
Posted in Linux | Tagged , , , | Comments Off on SolSpider – New version

ZFS Monitor – ZFS status monitoring for Android

Share

ZFS Monitor

admindojo proudly presents: ZFS Monitor.

Easily monitor your ZFS pools with this free Android App.

Market link: ZFS Monitor on Android Market

Direct QR:

Share
Posted in Linux | Tagged , , , | Comments Off on ZFS Monitor – ZFS status monitoring for Android

SolSpider – Android Spider clone

Share


I’ve been experimenting with the android-sdk for a while now and have now released my first app to the Android-Market.

SolSpider is free, open and you can get the source-code here.

The app does not use a SurfaceView with a seperate thread to draw the game at full steam, but rather leaves it to the main thread to issue a redraw. This design choice can make the drag-animation appear a bit sluggish, but it saves a lot of power in the long run.

Links:
SolSpider in Android-Market
SolSpider gitorious project
SolSpider source tree

Share
Posted in Uncategorized | Tagged , , , | Comments Off on SolSpider – Android Spider clone

Discrete Fourier Transform in C++ with fftw

Share

This tutorial will give you a short introduction to using libfftw in C++ under Linux.
The discrete Fourier-Transform (DFT) applied on an Image will translate the Image-Pixels into frequency and magnitude components.
Some knowledge of complex numbers is required to understand the transform itself but we’ll keep the math to the bare minimum.
The final result should yield a set of images that describe the original image and should look like this:

Input image

Magnitude image

Phase image

The magnitude appears to be black but isn’t. To make the information visible we scale the image logarithmically.

Magnitude scaled to visibility

We can recombine the magnitude and phase information with the inverse DFT and get:

Output image

As you can see in the image above, the inverse DFT of the magnitude and phase images produces the original image with barely no visible changes.
If you play around with the phase/magnitude you can alter the output image. Keep in mind that we will be saving 16bit images and most editing software under linux only supports 8bit color-spaces.
Re-saving the magnitude or phase image in 8bit color-space and performing the iDFT will return a very bad result:

Output created from 8bit sources

Step 1: Setup

To open images of different filetypes we will use libmagick and of course libfftw for the transform itself.
To be able to compile your program you’ll need (apart from your common g++ compiler) to download libfftw, libmagick++ and their respective development files.
If you’re running Debian/Ubuntu you can get those from your repos by typing:

sudo apt-get install libmagick++3 imagemagick libmagick++-dev libfftw3-3 libfftw3-dev

or get it on the web:
libfftw
Imagemagick

Create a folder and source structure:

mkdir fourier
touch fourier/Makefile
touch fourier/main.cpp

The makefile to our little project will enable us to use make instead of having to type the whole g++ compile-command.

CXXFLAGS=`Magick++-config --cppflags --cxxflags --ldflags --libs` `pkg-config fftw3 --libs` -g -Wall -O0
CXX=g++

all: fourier

fourier: main.o
	$(CXX) -o $@ $(CXXFLAGS) $<

main.o: main.cpp
	$(CXX) -c -o $@ $(CXXFLAGS) $<

clean:
	rm main.o

install:
	cp fourier /usr/bin/

uninstall:
	rm /usr/bin/fourier

Step 2: Programming

We will keep the code very basic and simple, when implementing the code in your own applications you might want to add some sort of exception handling and comments.
We start off with our includes:

#include <fftw3.h>
#include <Magick++.h>
#include <math.h>
#include <complex>
#include <iostream>
#include <string>
using namespace Magick;
using namespace std;

Most of these should be known to you, we add the libfftw and magick++ headers to get the functionality we need and also the std::complex class to perform some basic operations on complex numbers.
We will want our application to act in the following way:

fourier in outMag outPhase

for the forward and

fourier -ift inMag inPhase out

for the inverse transform.
So in our main.cpp we add

int main(int argc, char** argv)
{
	if(argc < 4) {
		cout << "usage:\t" << argv[0] << " in outMag outPhase" << endl
				<< "\t" << argv[0] << " -ift inMag inPhase out" << endl;
		return 0;
	}

	string arg = argv[1];

	if(arg != "-ift") {
		// FORWARD TRANSFORM
	}
	else {
		// BACKWARD TRANSFORM
	}
	return 0;
}

We will open the images with imagemagick. You can do so in C++ directly in the Magick::Image constructor.
To get access to the underlying pixels we use this method:

Image img("myfilename");

// lock image for modification
img.modifyImage();
Pixels pixelCache(img);
PixelPacket* pixels;

// get desired area in image
pixels = pixelCache.get(0, 0, sizeX, sizeY);

// read/write pixel (12,3)
*(pixel + 12 + (3 * sizeX)) = Color("Blue");

// write changes
pixelCache.sync()

Now that we know how to access image pixels, let’s write the prototypes for the forward and backward Fourier-Transform.
Right above the main method add

void fft(int squareSize, PixelPacket* pixels, PixelPacket* outMag, PixelPacket* outPhase)
{
}

void ift(int squareSize, PixelPacket* inMag, PixelPacket* inPhase, PixelPacket* outPixels)
{
}

Our prototypes describe methods that will take pointers to a PixelPacket array that represent our images.
Also we will only work with square Images (we will adjust them to our needs in main)
In the completed main method we prepare the images and call the respective method:

int main(int argc, char** argv)
{
	if(argc < 4) {
		cout << "usage:\t" << argv[0] << " in outMag outPhase" << endl
				<< "\t" << argv[0] << " -ift inMag inPhase out" << endl;
		return 0;
	}

	string arg = argv[1];

	if(arg != "-ift") {
		// FORWARD FOURIER TRANSFORM
		Image img(argv[1]);

		// get the length of the longer side of the image
		int squareSize = img.columns() < img.rows() ? img.rows() : img.columns();

		// the geometry of our padded image
		Geometry padded(squareSize, squareSize);
		padded.aspect(true);

		// make image square
		img.extent(padded);

		// create templates for magnitude and phase
		Image mag(Geometry(squareSize, squareSize), "Black");
		Image phase(Geometry(squareSize, squareSize), "Black");

		// get image pixels
		img.modifyImage();
		Pixels pixelCache(img);
		PixelPacket* pixels;
		pixels = pixelCache.get(0, 0, squareSize, squareSize);

		// get magnitude pixels
		mag.modifyImage();
		Pixels pixelCacheMag(mag);
		PixelPacket* pixelsMag;
		pixelsMag = pixelCacheMag.get(0, 0, squareSize, squareSize);

		// get phase pixels
		phase.modifyImage();
		Pixels pixelCachePhase(phase);
		PixelPacket* pixelsPhase;
		pixelsPhase = pixelCachePhase.get(0, 0, squareSize, squareSize);

		// perform fft
		fft(squareSize, pixels, pixelsMag, pixelsPhase);

		// write changes
		pixelCache.sync();
		pixelCacheMag.sync();
		pixelCachePhase.sync();

		// save files
		mag.write(argv[2]);
		phase.write(argv[3]);
	}
	else {
		// BACKWARD FOURIER TRANSFORM
	}
	return 0;
}

Both our functions still do nothing with the provided pixels, so let’s take a quick look at the libfftw api so we know how to use it.
libfftw provides several methods to do FFT on arrays of complex or real values and in multiple dimensions which we will not need to use. Our methods will only do transforms from real (although stored in a complex array) to complex and back. To do a transform libfftw provides a so called plan on which the operations are executed on. On creation one must supply the plan with an input and output array which have to be allocated first.

fftw_plan plan;
fftw_complex *in, *out;

// allocate memory for input and output arrays
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sizeX * sizeY);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sizeX * sizeY);

// create the plan
plan = fftw_plan_dft_2d(sizeX, sizeY, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

// fill the array
for(int i = 0, i < sizeY; i++)
	for(int j = 0; j < sizeX; j++)
		in[j + (sizeX * i)][0] = randomNumber; // in[][0] is the real part of in[]

// execute the specified transform
fftw_execute(plan);

// use output
[...]

// free allocated memory
fftw_destroy_plan(plan);
fftw_free(out);
fftw_free(in);

We have to create a plan for every color-channel in the image (we will only cover RGB) and create input and output-arrays respectively.
Let’s add just that to our fft method:

void fft(int squareSize, PixelPacket* pixels, PixelPacket* outMag, PixelPacket* outPhase)
{
	fftw_plan planR, planG, planB;
	fftw_complex *inR, *inG, *inB, *outR, *outG, *outB;

	// allocate input arrays
	inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// allocate output arrays
	outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// create plans
	planR = fftw_plan_dft_2d(squareSize, squareSize, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
	planG = fftw_plan_dft_2d(squareSize, squareSize, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);
	planB = fftw_plan_dft_2d(squareSize, squareSize, inB, outB, FFTW_FORWARD, FFTW_ESTIMATE);

	// TODO: assign color-values to input arrays

	// perform FORWARD fft
	fftw_execute(planR);
	fftw_execute(planG);
	fftw_execute(planB);

	// TODO: calculate output mag/phase

	// free memory
	fftw_destroy_plan(planR);
	fftw_destroy_plan(planG);
	fftw_destroy_plan(planB);
	fftw_free(inR); fftw_free(outR);
	fftw_free(inG); fftw_free(outG);
	fftw_free(inB); fftw_free(outB);
}

Putting some values to perform the fft on would be nice. We will use the color-values as they are and put them into the input arrays for each channel.

void fft(int squareSize, PixelPacket* pixels, PixelPacket* outMag, PixelPacket* outPhase)
{
	fftw_plan planR, planG, planB;
	fftw_complex *inR, *inG, *inB, *outR, *outG, *outB;

	// allocate input arrays
	inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// allocate output arrays
	outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// create plans
	planR = fftw_plan_dft_2d(squareSize, squareSize, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
	planG = fftw_plan_dft_2d(squareSize, squareSize, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);
	planB = fftw_plan_dft_2d(squareSize, squareSize, inB, outB, FFTW_FORWARD, FFTW_ESTIMATE);

	// assign values to real parts (values between 0 and MaxRGB)
	for(int i = 0; i < squareSize * squareSize; i++) {
		PixelPacket current = *(pixels + i);
		double red = current.red;
		double green = current.green;
		double blue = current.blue;

		// save as real numbers
		inR[i][0] = red;
		inG[i][0] = green;
		inB[i][0] = blue;
	}

	// perform FORWARD fft
	fftw_execute(planR);
	fftw_execute(planG);
	fftw_execute(planB);

	// TODO: calculate output mag/phase

	// free memory
	fftw_destroy_plan(planR);
	fftw_destroy_plan(planG);
	fftw_destroy_plan(planB);
	fftw_free(inR); fftw_free(outR);
	fftw_free(inG); fftw_free(outG);
	fftw_free(inB); fftw_free(outB);
}

A tiny bit of math

Now that the FFT has been performed we have an output array of complex numbers for each channel (RGB).
To get the information we need, namely magnitude and phase, we will perform some simple math on these numbers.
The magnitude M is calculated as the square root of the scalar product of the components of the complex result z:
z = a + ib
M = |z| = \sqrt{(a^2 + b^2)}
The phase \varphi is calculated with the \arg operator. Both definitions can be found in the polar coordinate system wiki.
\varphi = \arg(z)
To calculate the phase we will use the implementation of std::arg in the complex class.
As we can only save positive values to files (from 0 to MaxRGB) we will have to shift the phase from [-\pi,\pi] to [0,2\pi] followed by a scaling to [0,\mathrm{MaxRGB}] and scale/shift it back when doing the inverse DFT. Also libfftw for performance reasons scales the output by (sizeX * sizeY) which we will correct before calculating phase and magnitude by dividing it by just that.
The code to do all this stuff:

void fft(int squareSize, PixelPacket* pixels, PixelPacket* outMag, PixelPacket* outPhase)
{
	fftw_plan planR, planG, planB;
	fftw_complex *inR, *inG, *inB, *outR, *outG, *outB;

	// allocate input arrays
	inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// allocate output arrays
	outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// create plans
	planR = fftw_plan_dft_2d(squareSize, squareSize, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
	planG = fftw_plan_dft_2d(squareSize, squareSize, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);
	planB = fftw_plan_dft_2d(squareSize, squareSize, inB, outB, FFTW_FORWARD, FFTW_ESTIMATE);

	// assign values to real parts (values between 0 and MaxRGB)
	for(int i = 0; i < squareSize * squareSize; i++) {
		PixelPacket current = *(pixels + i);
		double red = current.red;
		double green = current.green;
		double blue = current.blue;

		// save as real numbers
		inR[i][0] = red;
		inG[i][0] = green;
		inB[i][0] = blue;
	}

	// perform FORWARD fft
	fftw_execute(planR);
	fftw_execute(planG);
	fftw_execute(planB);

	// transform imaginary number to phase and magnitude and save to output
	for(int i = 0; i < squareSize * squareSize; i++) {
		// normalize values
		double realR = outR[i][0] / (double)(squareSize * squareSize);
		double imagR = outR[i][1] / (double)(squareSize * squareSize);

		double realG = outG[i][0] / (double)(squareSize * squareSize);
		double imagG = outG[i][1] / (double)(squareSize * squareSize);

		double realB = outB[i][0] / (double)(squareSize * squareSize);
		double imagB = outB[i][1] / (double)(squareSize * squareSize);

		// magnitude
		double magR = sqrt((realR * realR) + (imagR * imagR));
		double magG = sqrt((realG * realG) + (imagG * imagG));
		double magB = sqrt((realB * realB) + (imagB * imagB));

		// write to output
		(*(outMag + i)).red = magR;
		(*(outMag + i)).green = magG;
		(*(outMag + i)).blue = magB;

		// std::complex for arg()
		complex<double> cR(realR, imagR);
		complex<double> cG(realG, imagG);
		complex<double> cB(realB, imagB);

		// phase
		double phaseR = arg(cR) + M_PI;
		double phaseG = arg(cG) + M_PI;
		double phaseB = arg(cB) + M_PI;

		// scale and write to output
		(*(outPhase + i)).red = (phaseR / (double)(2 * M_PI)) * MaxRGB;
		(*(outPhase + i)).green = (phaseG / (double)(2 * M_PI)) * MaxRGB;
		(*(outPhase + i)).blue = (phaseB / (double)(2 * M_PI)) * MaxRGB;
	}

	// free memory
	fftw_destroy_plan(planR);
	fftw_destroy_plan(planG);
	fftw_destroy_plan(planB);
	fftw_free(inR); fftw_free(outR);
	fftw_free(inG); fftw_free(outG);
	fftw_free(inB); fftw_free(outB);
}

A first Test

If you havent been writing the code up to now here’s what your main.cpp should be looking like:

#include <fftw3.h>
#include <Magick++.h>
#include <math.h>
#include <complex>
#include <iostream>
#include <string>
using namespace Magick;
using namespace std;

void fft(int squareSize, PixelPacket* pixels, PixelPacket* outMag, PixelPacket* outPhase)
{
	fftw_plan planR, planG, planB;
	fftw_complex *inR, *inG, *inB, *outR, *outG, *outB;

	// allocate input arrays
	inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// allocate output arrays
	outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// create plans
	planR = fftw_plan_dft_2d(squareSize, squareSize, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
	planG = fftw_plan_dft_2d(squareSize, squareSize, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);
	planB = fftw_plan_dft_2d(squareSize, squareSize, inB, outB, FFTW_FORWARD, FFTW_ESTIMATE);

	// assign values to real parts (values between 0 and MaxRGB)
	for(int i = 0; i < squareSize * squareSize; i++) {
		PixelPacket current = *(pixels + i);
		double red = current.red;
		double green = current.green;
		double blue = current.blue;
		// save as real numbers
		inR[i][0] = red;
		inG[i][0] = green;
		inB[i][0] = blue;
	}

	// perform FORWARD fft
	fftw_execute(planR);
	fftw_execute(planG);
	fftw_execute(planB);

	// transform imaginary number to phase and magnitude and save to output
	for(int i = 0; i < squareSize * squareSize; i++) {
		// normalize values
		double realR = outR[i][0] / (double)(squareSize * squareSize);
		double imagR = outR[i][1] / (double)(squareSize * squareSize);

		double realG = outG[i][0] / (double)(squareSize * squareSize);
		double imagG = outG[i][1] / (double)(squareSize * squareSize);

		double realB = outB[i][0] / (double)(squareSize * squareSize);
		double imagB = outB[i][1] / (double)(squareSize * squareSize);

		// magnitude
		double magR = sqrt((realR * realR) + (imagR * imagR));
		double magG = sqrt((realG * realG) + (imagG * imagG));
		double magB = sqrt((realB * realB) + (imagB * imagB));

		// write to output
		(*(outMag + i)).red = magR;
		(*(outMag + i)).green = magG;
		(*(outMag + i)).blue = magB;

		// std::complex for arg()
		complex<double> cR(realR, imagR);
		complex<double> cG(realG, imagG);
		complex<double> cB(realB, imagB);

		// phase
		double phaseR = arg(cR) + M_PI;
		double phaseG = arg(cG) + M_PI;
		double phaseB = arg(cB) + M_PI;

		// scale and write to output
		(*(outPhase + i)).red = (phaseR / (double)(2 * M_PI)) * MaxRGB;
		(*(outPhase + i)).green = (phaseG / (double)(2 * M_PI)) * MaxRGB;
		(*(outPhase + i)).blue = (phaseB / (double)(2 * M_PI)) * MaxRGB;
	}

	// free memory
	fftw_destroy_plan(planR);
	fftw_destroy_plan(planG);
	fftw_destroy_plan(planB);
	fftw_free(inR); fftw_free(outR);
	fftw_free(inG); fftw_free(outG);
	fftw_free(inB); fftw_free(outB);
}

void ift(int squareSize, PixelPacket* inMag, PixelPacket* inPhase, PixelPacket* outPixels)
{
}

int main(int argc, char** argv)
{
	if(argc < 4) {
		cout << "usage:\t" << argv[0] << " in outMag outPhase" << endl
				<< "\t" << argv[0] << " -ift inMag inPhase out" << endl;
		return 0;
	}

	string arg = argv[1];

	if(arg != "-ift") {
		// FORWARD FOURIER TRANSFORM
		Image img(argv[1]);

		// get the length of the longer side of the image
		int squareSize = img.columns() < img.rows() ? img.rows() : img.columns();

		// the geometry of our padded image
		Geometry padded(squareSize, squareSize);
		padded.aspect(true);

		// make image square
		img.extent(padded);

		// create templates for magnitude and phase
		Image mag(Geometry(squareSize, squareSize), "Black");
		Image phase(Geometry(squareSize, squareSize), "Black");

		// get image pixels
		img.modifyImage();
		Pixels pixelCache(img);
		PixelPacket* pixels;
		pixels = pixelCache.get(0, 0, squareSize, squareSize);

		// get magnitude pixels
		mag.modifyImage();
		Pixels pixelCacheMag(mag);
		PixelPacket* pixelsMag;
		pixelsMag = pixelCacheMag.get(0, 0, squareSize, squareSize);

		// get phase pixels
		phase.modifyImage();
		Pixels pixelCachePhase(phase);
		PixelPacket* pixelsPhase;
		pixelsPhase = pixelCachePhase.get(0, 0, squareSize, squareSize);

		// perform fft
		fft(squareSize, pixels, pixelsMag, pixelsPhase);

		// write changes
		pixelCache.sync();
		pixelCacheMag.sync();
		pixelCachePhase.sync();

		// save files
		mag.write(argv[2]);
		phase.write(argv[3]);
	}
	else {
		// BACKWARD FOURIER TRANSFORM
	}
	return 0;
}

Open a shell, enter your fourier directory and type make and sudo make install

user@host:~/fourier$ make
g++ -c -o main.o `Magick++-config --cppflags --cxxflags --ldflags --libs` `pkg-config fftw3 --libs` -g -Wall -O0 main.cpp
g++ -o fourier `Magick++-config --cppflags --cxxflags --ldflags --libs` `pkg-config fftw3 --libs` -g -Wall -O0 main.o
user@host:~/fourier$ sudo make install
cp fourier /usr/bin/

Now you can run fourier with the right input parameters to get a magnitude and phase image.
Save the image at the top of this tutorial to a directory of your choice. In this directory execute the following:

user@host:~/fft$ fourier in.png mag.png phase.png
user@host:~/fft$ convert mag.png -contrast-stretch 0 -evaluate log 10000 mag-lognorm.png

The first command will do the FFT and save the images. The second command will scale the magnitude image logarithmically as seen at the top of this page.
You may notice that the images do not match. Indeed our transform so far produces something like this:

Current magnitude image

The pictude above is the correct representation of the magnitudes of our input images’ frequencies although the zero-frequency is not in the center, but in the corners. The zero frequency is often moved to (sizeX/2, sizeY/2) so it converges to the middle. We get the desired result if we swap the four quadrants of the image diagonally.
The method swapQuadrants will do just that:

void swapQuadrants(int squareSize, PixelPacket* pixels)
{
	int half = floor(squareSize / (double)2);

	// swap quadrants diagonally
	for(int i = 0; i < half; i++) {
		for(int j = 0; j < half; j++) {
			int upper = j + (squareSize * i);
			int lower = upper + (squareSize * half) + half;

			PixelPacket cur0 = *(pixels + upper);
			*(pixels + upper) = *(pixels + lower);
			*(pixels + lower) = cur0;

			PixelPacket cur1 = *(pixels + upper + half);
			*(pixels + upper + half) = *(pixels + lower - half);
			*(pixels + lower - half) = cur1;
		}
	}
}

Add this code just after the includes so we can use it within fft.
Calling the swapQuadrants on magnitude and phase output before saving them in our fft method will then give us the desired image

void fft(int squareSize, PixelPacket* pixels, PixelPacket* outMag, PixelPacket* outPhase)
{
	fftw_plan planR, planG, planB;
	fftw_complex *inR, *inG, *inB, *outR, *outG, *outB;

	// allocate input arrays
	inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// allocate output arrays
	outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
	outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);

	// create plans
	planR = fftw_plan_dft_2d(squareSize, squareSize, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
	planG = fftw_plan_dft_2d(squareSize, squareSize, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);
	planB = fftw_plan_dft_2d(squareSize, squareSize, inB, outB, FFTW_FORWARD, FFTW_ESTIMATE);

	// assign values to real parts (values between 0 and MaxRGB)
	for(int i = 0; i < squareSize * squareSize; i++) {
		PixelPacket current = *(pixels + i);
		double red = current.red;
		double green = current.green;
		double blue = current.blue;

		// save as real numbers
		inR[i][0] = red;
		inG[i][0] = green;
		inB[i][0] = blue;
	}

	// perform FORWARD fft
	fftw_execute(planR);
	fftw_execute(planG);
	fftw_execute(planB);

	// transform imaginary number to phase and magnitude and save to output
	for(int i = 0; i < squareSize * squareSize; i++) {
		// normalize values
		double realR = outR[i][0] / (double)(squareSize * squareSize);
		double imagR = outR[i][1] / (double)(squareSize * squareSize);

		double realG = outG[i][0] / (double)(squareSize * squareSize);
		double imagG = outG[i][1] / (double)(squareSize * squareSize);

		double realB = outB[i][0] / (double)(squareSize * squareSize);
		double imagB = outB[i][1] / (double)(squareSize * squareSize);

		// magnitude
		double magR = sqrt((realR * realR) + (imagR * imagR));
		double magG = sqrt((realG * realG) + (imagG * imagG));
		double magB = sqrt((realB * realB) + (imagB * imagB));

		// write to output
		(*(outMag + i)).red = magR;
		(*(outMag + i)).green = magG;
		(*(outMag + i)).blue = magB;

		// std::complex for arg()
		complex<double> cR(realR, imagR);
		complex<double> cG(realG, imagG);
		complex<double> cB(realB, imagB);

		// phase
		double phaseR = arg(cR) + M_PI;
		double phaseG = arg(cG) + M_PI;
		double phaseB = arg(cB) + M_PI;

		// scale and write to output
		(*(outPhase + i)).red = (phaseR / (double)(2 * M_PI)) * MaxRGB;
		(*(outPhase + i)).green = (phaseG / (double)(2 * M_PI)) * MaxRGB;
		(*(outPhase + i)).blue = (phaseB / (double)(2 * M_PI)) * MaxRGB;
	}

	// move zero frequency to (squareSize/2, squareSize/2)
	swapQuadrants(squareSize, outMag);
	swapQuadrants(squareSize, outPhase);

	// free memory
	fftw_destroy_plan(planR);
	fftw_destroy_plan(planG);
	fftw_destroy_plan(planB);
	fftw_free(inR); fftw_free(outR);
	fftw_free(inG); fftw_free(outG);
	fftw_free(inB); fftw_free(outB);
}

After compiling, installing and executing the fourier and convert commands, you should get an image similar to that at the top of the page.

Implementing the iDFT

As this part is very much the same as everything described above we will keep it short.
All the ift method has to do is:

  • Swap the quadrants back into their original place
  • Scale the phase values back to [0,2\pi]
  • Shift the phase values back to [-\pi,\pi]
  • Transform phase and magnitude back to real and imaginary parts
  • Perform backward fft on these values
  • Store real parts in output
  • The final main.cpp doing all this looks like so:

    #include <fftw3.h>
    #include <Magick++.h>
    #include <math.h>
    #include <complex>
    #include <iostream>
    #include <string>
    using namespace Magick;
    using namespace std;
    
    void swapQuadrants(int squareSize, PixelPacket* pixels)
    {
    	int half = floor(squareSize / (double)2);
    	
    	// swap quadrants diagonally
    	for(int i = 0; i < half; i++) {
    		for(int j = 0; j < half; j++) {
    			int upper = j + (squareSize * i);
    			int lower = upper + (squareSize * half) + half;
    			
    			PixelPacket cur0 = *(pixels + upper);
    			*(pixels + upper) = *(pixels + lower);
    			*(pixels + lower) = cur0;
    			
    			PixelPacket cur1 = *(pixels + upper + half);
    			*(pixels + upper + half) = *(pixels + lower - half);
    			*(pixels + lower - half) = cur1;
    		}
    	}
    }
    
    void fft(int squareSize, PixelPacket* pixels, PixelPacket* outMag, PixelPacket* outPhase)
    {
    	fftw_plan planR, planG, planB;
    	fftw_complex *inR, *inG, *inB, *outR, *outG, *outB;
    	
    	// allocate input arrays
    	inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	
    	// allocate output arrays
    	outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	
    	// create plans
    	planR = fftw_plan_dft_2d(squareSize, squareSize, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
    	planG = fftw_plan_dft_2d(squareSize, squareSize, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);
    	planB = fftw_plan_dft_2d(squareSize, squareSize, inB, outB, FFTW_FORWARD, FFTW_ESTIMATE);
    	
    	// assign values to real parts (values between 0 and MaxRGB)
    	for(int i = 0; i < squareSize * squareSize; i++) {
    		PixelPacket current = *(pixels + i);
    		double red = current.red;
    		double green = current.green;
    		double blue = current.blue;
    		
    		// save as real numbers
    		inR[i][0] = red;
    		inG[i][0] = green;
    		inB[i][0] = blue;
    	}
    	
    	// perform FORWARD fft
    	fftw_execute(planR);
    	fftw_execute(planG);
    	fftw_execute(planB);
    	
    	// transform imaginary number to phase and magnitude and save to output
    	for(int i = 0; i < squareSize * squareSize; i++) {
    		// normalize values
    		double realR = outR[i][0] / (double)(squareSize * squareSize);
    		double imagR = outR[i][1] / (double)(squareSize * squareSize);
    				
    		double realG = outG[i][0] / (double)(squareSize * squareSize);
    		double imagG = outG[i][1] / (double)(squareSize * squareSize);
    		
    		double realB = outB[i][0] / (double)(squareSize * squareSize);
    		double imagB = outB[i][1] / (double)(squareSize * squareSize);
    		
    		// magnitude
    		double magR = sqrt((realR * realR) + (imagR * imagR));
    		double magG = sqrt((realG * realG) + (imagG * imagG));
    		double magB = sqrt((realB * realB) + (imagB * imagB));
    		
    		// write to output
    		(*(outMag + i)).red = magR;
    		(*(outMag + i)).green = magG;
    		(*(outMag + i)).blue = magB;
    		
    		// std::complex for arg()
    		complex<double> cR(realR, imagR);
    		complex<double> cG(realG, imagG);
    		complex<double> cB(realB, imagB);
    		
    		// phase
    		double phaseR = arg(cR) + M_PI;
    		double phaseG = arg(cG) + M_PI;
    		double phaseB = arg(cB) + M_PI;
    		
    		// scale and write to output
    		(*(outPhase + i)).red = (phaseR / (double)(2 * M_PI)) * MaxRGB;
    		(*(outPhase + i)).green = (phaseG / (double)(2 * M_PI)) * MaxRGB;
    		(*(outPhase + i)).blue = (phaseB / (double)(2 * M_PI)) * MaxRGB;
    	}
    	
    	// move zero frequency to (squareSize/2, squareSize/2)
    	swapQuadrants(squareSize, outMag);
    	swapQuadrants(squareSize, outPhase);
    	
    	// free memory
    	fftw_destroy_plan(planR);
    	fftw_destroy_plan(planG);
    	fftw_destroy_plan(planB);
    	fftw_free(inR); fftw_free(outR);
    	fftw_free(inG); fftw_free(outG);
    	fftw_free(inB); fftw_free(outB);
    }
    
    void ift(int squareSize, PixelPacket* inMag, PixelPacket* inPhase, PixelPacket* outPixels)
    {
    	// move zero frequency back to corners
    	swapQuadrants(squareSize, inMag);
    	swapQuadrants(squareSize, inPhase);
    	
    	fftw_plan planR, planG, planB;
    	fftw_complex *inR, *inG, *inB, *outR, *outG, *outB;
    	
    	// allocate input arrays
    	inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	
    	// allocate output arrays
    	outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squareSize * squareSize);
    	
    	// create plans
    	planR = fftw_plan_dft_2d(squareSize, squareSize, inR, outR, FFTW_BACKWARD, FFTW_ESTIMATE);
    	planG = fftw_plan_dft_2d(squareSize, squareSize, inG, outG, FFTW_BACKWARD, FFTW_ESTIMATE);
    	planB = fftw_plan_dft_2d(squareSize, squareSize, inB, outB, FFTW_BACKWARD, FFTW_ESTIMATE);
    	
    	// transform magnitude/phase to real/imaginary
    	for(int i = 0; i < squareSize * squareSize; i++) {
    		double magR = inMag[i].red;
    		double phaseR = ((inPhase[i].red / (double)MaxRGB) * 2 * M_PI) - M_PI;
    		inR[i][0] = (magR * cos(phaseR));
    		inR[i][1] = (magR * sin(phaseR));
    		
    		double magB = inMag[i].blue;
    		double phaseB = ((inPhase[i].blue / (double)MaxRGB) * 2 * M_PI) - M_PI;
    		inB[i][0] = (magB * cos(phaseB));
    		inB[i][1] = (magB * sin(phaseB));
    		
    		double magG = inMag[i].green;
    		double phaseG = ((inPhase[i].green / (double)MaxRGB) * 2 * M_PI) - M_PI;
    		inG[i][0] = (magG * cos(phaseG));
    		inG[i][1] = (magG * sin(phaseG));
    	}
    	
    	// perform ift
    	fftw_execute(planR);
    	fftw_execute(planG);
    	fftw_execute(planB);
    	
    	// save real parts to output
    	for(int i = 0; i < squareSize * squareSize; i++) {
    		double magR = outR[i][0];
    		double magG = outG[i][0];
    		double magB = outB[i][0];
    		
    		// make sure it's capped at MaxRGB
    		(*(outPixels + i)).red = magR > MaxRGB ? MaxRGB : magR;
    		(*(outPixels + i)).green = magG > MaxRGB ? MaxRGB : magG;
    		(*(outPixels + i)).blue = magB > MaxRGB ? MaxRGB : magB;
    	}
    	
    	// free memory
    	fftw_destroy_plan(planR);
    	fftw_destroy_plan(planG);
    	fftw_destroy_plan(planB);
    	fftw_free(inR); fftw_free(outR);
    	fftw_free(inG); fftw_free(outG);
    	fftw_free(inB); fftw_free(outB);
    }
    
    int main(int argc, char** argv)
    {
    	if(argc < 4) {
    		cout << "usage:\t" << argv[0] << " in outMag outPhase" << endl
    				<< "\t" << argv[0] << " -ift inMag inPhase out" << endl;
    		return 0;
    	}
    	
    	string arg = argv[1];
    	
    	if(arg != "-ift") {
    		// FORWARD FOURIER TRANSFORM
    		Image img(argv[1]);
    		
    		// get the length of the longer side of the image
    		int squareSize = img.columns() < img.rows() ? img.rows() : img.columns();
    		
    		// the geometry of our padded image		
    		Geometry padded(squareSize, squareSize);
    		padded.aspect(true);
    		
    		// make image square
    		img.extent(padded);
    		
    		// create templates for magnitude and phase
    		Image mag(Geometry(squareSize, squareSize), "Black");
    		Image phase(Geometry(squareSize, squareSize), "Black");
    	
    		// get image pixels
    		img.modifyImage();
    		Pixels pixelCache(img);
    		PixelPacket* pixels;
    		pixels = pixelCache.get(0, 0, squareSize, squareSize);
    		
    		// get magnitude pixels
    		mag.modifyImage();
    		Pixels pixelCacheMag(mag);
    		PixelPacket* pixelsMag;
    		pixelsMag = pixelCacheMag.get(0, 0, squareSize, squareSize);
    		
    		// get phase pixels
    		phase.modifyImage();
    		Pixels pixelCachePhase(phase);
    		PixelPacket* pixelsPhase;
    		pixelsPhase = pixelCachePhase.get(0, 0, squareSize, squareSize);
    		
    		// perform fft
    		fft(squareSize, pixels, pixelsMag, pixelsPhase);
    		
    		// write changes
    		pixelCache.sync();
    		pixelCacheMag.sync();
    		pixelCachePhase.sync();
    		
    		// save files
    		mag.write(argv[2]);
    		phase.write(argv[3]);
    	}
    	else {
    		// BACKWARD FOURIER TRANSFORM
    		Image mag(argv[2]);
    		Image phase(argv[3]);
    		
    		// get size
    		int squareSize = mag.columns();
    		Image img(Geometry(squareSize, squareSize), "Black");
    		
    		// get image pixels
    		img.modifyImage();
    		Pixels pixelCache(img);
    		PixelPacket* pixels;
    		pixels = pixelCache.get(0, 0, squareSize, squareSize);
    		
    		// get magnitude pixels
    		mag.modifyImage();
    		Pixels pixelCacheMag(mag);
    		PixelPacket* pixelsMag;
    		pixelsMag = pixelCacheMag.get(0, 0, squareSize, squareSize);
    		
    		// get phase pixels
    		phase.modifyImage();
    		Pixels pixelCachePhase(phase);
    		PixelPacket* pixelsPhase;
    		pixelsPhase = pixelCachePhase.get(0, 0, squareSize, squareSize);
    		
    		// perform ift
    		ift(squareSize, pixelsMag, pixelsPhase, pixels);
    		
    		// write changes
    		pixelCache.sync();
    		
    		// save file
    		img.write(argv[4]);
    	}
    
    	return 0;
    }
    

    You can get an archive with all the sources here.
    Remember to remove the -g and -O0 flags form the makefile when actually using it. (These are there for debug purposes)

    Executing the command

    user@host:~/fft$ fourier -ift mag.png phase.png out.png
    

    will produce and save the image obtained from the inverse DFT in out.png.

    I take no responsibility in how you use this code and make no guarantee as to its fitness for production environments.
    Have fun experimenting with libfftw in C++.

    Share
    Posted in Linux | Tagged , , , , , , | Comments Off on Discrete Fourier Transform in C++ with fftw

    Setting up ntpd on FreeBSD

    Share

    Ntpd is a daemon that fetches the current time from servers on the internet. It assures correct time and date for your system, which is important for cronjobs and timetags on files.

    Step 1:

    Ntpd is preinstalled on every FreeBSD install, so we just need to edit the config file. It’s located in /etc/ntp.conf. Edit it as root with your preferred editor or use:

    vi /etc/ntp.conf

    You might find a lot of comments in this file, all of which you can safely ignore. Beneath the comments we’ll append our ntp server locations. On the ntp pool project you will find a complete list of servers. We also need to tell our system where our driftfile is, as it stores the difference between our system clock and the external ones. Once you’ve found your servers, create a file that looks like this:

    server 0.ch.pool.ntp.org
    server 1.ch.pool.ntp.org
    server 2.ch.pool.ntp.org
    server 3.ch.pool.ntp.org
    
    driftfile /var/db/ntp.drift

    This way ntpd uses up to four of the servers in the pool and knows were the driftfile is. Be sure to use servers from your own country or at least ones that are close by ping so the delay stays small.

    Step 2:

    Now we need to set the correct timezone. Go to /usr/share/zoneinfo/ and find your closest locality. Then you create a symbolic link as root to our localtime which is done this way:

    cp /usr/share/zoneinfo/Europe/Zurich /etc/localtime

    This sets my timezone to Central Europe. You can use date to check if everything worked properly.

    Step 3.1:

    If your system clock is more than a half hour apart from the current time, we need to set it by using ntpdate since ntpd doesn’t correct for such vast differences. Grab a server close to you from your ntp.conf file and issue this command as root:

    ntpdate ch.pool.ntp.org

    This grabs the correct time from a timeserver, in this case one from Switzerland, and uses it as your system time.

    Step 3.2:

    We want ntpd to start up every time the server is booted so we need to append as root the following to /etc/rc.conf:

    ntpd_enable="YES"

    We can start up ntpd by issuing this command as root:

    /etc/rc.d/ntpd start

    Now top should be showing you the process ntpd and our work is done.

    Share
    Posted in Unix | Tagged , , , | Comments Off on Setting up ntpd on FreeBSD

    Google.com strange symbols

    Share

    Just recently after a fresh install of Ubuntu 10.10 on my laptop I noticed that the default language in Firefox for my Google.com account had somehow been changed to Afrikaans. Even changing the language and saving my settings would in the end result in the buttons being unreadable to me once again.
    So I did a little digging and it turns out Firefox Sync is the culprit in this case. A bug has already been filed so here is a workaround for those of you who can’t wait.

    Workaround:

    Open Firefox and go to

    about:config

    Read the warning and proceed if you dare.
    Search for

    intl.accept_languages

    and set the value to

    en-us, en

    Restart Firefox and everything should be back to normal.

    Share
    Posted in Applications | Tagged , , | Comments Off on Google.com strange symbols

    Patching “fixed channel mon0: -1” in aircrack-ng suite

    Share

    If you’ve recently updated your kernel to 2.6.35 you may have encountered an error while using airodump-ng where the displayed channel would not match the channel you set in airmon-ng when starting your monitor device.
    Before you start patching anything make sure you have tried the following:

    stop networt-manager; stop avahi-daemon; killall wpa_supplicant

    and then

    airmon-ng start wlan0 x

    where x is the channel you wish to monitor.

    If after these steps you still get

    fixed channel mon0: -1

    in the upper right corner of your airodump-ng output, you may proceed to patching compat-wireless.

    Step 1:

    Get the latest compat-wireless source tarball here, or wget it with

    "wget http://wireless.kernel.org/download/compat-wireless-2.6/
    compat-wireless-2.6.tar.bz2"

    Untar it to the current directory with

    tar jxvf compat-wireless-2.6.tar.bz2

    Now enter the directory with

    cd compat-wireless-<tab>

    Step 2:

    You will now need the patch to make things work again. You can either get it here, or wget it using

    wget http://patches.aircrack-ng.org/channel-negative-one-maxim.patch

    and apply it using patch

    patch -Np1 -i channel-negative-one-maxim.patch

    You should get the following output:

    patching file net/wireless/chan.c
    Hunk #1 succeeded at 50 (offset 1 line).
    Hunk #2 succeeded at 80 (offset 1 line).

    Step 3:

    You now need to compile the patched compat-wireless source with:

    make

    and as soon as it’s finished compiling you’ll need to execute the following as root

    make install

    or if you like sudo it would be

    sudo make install

    followed by

    sudo make unload

    and a reboot.

    Note:
    afaik make install will not overwrite your old drivers, so if you need to revert to your old config run

    cd compat-wireless-<tab>
    sudo make uninstall

    to remove all newly installed drivers.

    Share
    Posted in Linux | Tagged , | 3 Comments

    Installing WordPress 3.0.1 on Ubuntu Server 10.10

    Share

    WordPress is a web application for blogging or creating your own website.
    If you want to try it out without installing anything on your server yet, you can create your own blog for free at wordpress.com although you will not be able to upload any plugins or custom themes.

    Now on to setting up your own WordPress server!

    Step 1:

    The wordpress package is available in the ubuntu repositories so all you have to do to get wordpress on your server is:

    sudo apt-get install wordpress

    This will install all needed dependencies like apache2 and mysql. During the installation of mysql apt will ask you to set a root password for your mysql server. You will need this password later on.

    Step 2:

    WordPress needs a mysql database and username to store its data.
    Open a shell and type

    mysql -u root -p

    and supply the password you entered during the mysql-install.
    You should see something like this:

    Welcome to the MySQL monitor.  Commands end with ; or \g.
    Your MySQL connection id is 1093
    Server version: 5.1.49-1ubuntu8 (Ubuntu)
    
    Copyright (c) 2000, 2010, Oracle and/or its affiliates. All rights
    reserved.
    This software comes with ABSOLUTELY NO WARRANTY. This is free
    software,
    and you are welcome to modify and redistribute it under the GPL v2
    license
    
    Type 'help;' or '\h' for help. Type '\c' to clear the current input
    statement.
    
    mysql>

    To create a new database named wordpress enter:

    CREATE DATABASE wordpress;

    Now grant access to wordpress to user wordpress with password wordpressPassword:

    GRANT ALL PRIVILEGES ON wordpress.* TO "wordpress"@"localhost"
    -> IDENTIFIED BY "wordpressPassword";

    Step 3:

    Copy the wp-config.php file to /etc/wordpress with

    cd /usr/share/wordpress
    sudo cp wp-config-sample.php /etc/wordpress/wp-config.php

    and edit the file with your preferred editor

    sudo nano /etc/wordpress/wp-config.php

    You will need to enter your database information. Find the following section and enter your data accordingly:

    /** The name of the database for WordPress */
    define('DB_NAME', 'wordpress');
    
    /** MySQL database username */
    define('DB_USER', 'wordpress');
    
    /** MySQL database password */
    define('DB_PASSWORD', 'wordpressPassword');
    
    /** MySQL hostname */
    define('DB_HOST', 'localhost');

    Step 4:

    You will need to make a symbolic link to your wordpress root in /var/www/ which is your apache2 web root. You can do this by entering:

    sudo ln -s /usr/share/wordpress/ /var/www/

    Step 5:

    WordPress is now ready to run its install script. Open your browser and go to

    http://hostname/wordpress/wp-admin/install.php

    where hostname is the name or IP of your server running wordpress.
    Follow the instructions on screen.

    Step 6:

    For wordpress to be able to install plugins you need to grant it access over FTP. In this tutorial I’ll be using proftpd, but you may use any FTP-Server.
    Install proftpd by entering:

    sudo apt-get install proftpd

    and make it standalone when prompted.
    Next edit the config-file with:

    sudo nano /etc/proftpd/proftpd.conf

    and append the line

    DefaultRoot     /var/www/wordpress

    and finally restart proftpd with:

    sudo /etc/init.d/proftpd restart
    Share
    Posted in Linux | Tagged , | Comments Off on Installing WordPress 3.0.1 on Ubuntu Server 10.10