Using Fourier synthesis to generate a fractional Brownian motion surface

In this post we will discuss generating fractal terrain. In the previous post we implemented our own fast Fourier transform in order to simulate an ocean surface. In that post we implemented an unnormalized inverse transform. It seemed logical to employ our FFT object as a means of generating terrain, but in order to do so we will need to add a method to compute the forward transform. This will be a relatively brief post. Most of the foundation has been laid with our ocean simulation, so we will simply retool our FFT object for the purpose of generating terrain.

The process we will implement begins with generating a two-dimensional Gaussian white noise. We then take the Fourier transform of the white noise and apply a filter in the frequency domain. This filter has the effect of curbing the amplitudes of the higher frequencies. We then apply the inverse transform and interpret the result as a height field. Below are two screen captures, one for and the other for . Periodicity of the Fourier transform allows us to tile our maps. The images below have been tiled.

Terrain generated with a beta of 2.4 and tiled.

Terrain generated with a beta of 2.4 and tiled.

Terrain generated with a beta of 2.0 and tiled.

Terrain generated with a beta of 2.0 and tiled.

In the previous post we found our twiddle factors for the inverse transform to be,

\begin{align*}
T_N^x &= \exp\left(\frac{i2\pi x}{N}\right) \\
\end{align*}

If we follow that analysis for the case of the forward transform, we find the twiddle factors,

\begin{align*}
{T_N^x}' &= \exp\left(\frac{-i2\pi x}{N}\right) \\
\end{align*}

So the only modifications we need to make to our FFT object is to precompute a second set of twiddle factors and to apply those in the case of the forward transform. See src/fft.h and src/fft.cc in the project download.

In order to generate terrain using the method of Fourier synthesis our first step will be to generate a Gaussian white noise (see src/terrain.h and src/terrain.cc for the complete implementation).

// generate noise in the spatial domain
for (int i = 0; i < N * N; i++) {
	c[i] = gaussianRandomVariable();
}

A rendering of the white noise is below.

A rendering of the Gaussian white noise.

A rendering of the Gaussian white noise.

We then apply the forward transform to transform the white noise into the frequency domain.

// transform to the frequency domain
for (int j = 0; j < N; j++) fft.fft_(c, c, 1, j * N);
for (int i = 0; i < N; i++) fft.fft_(c, c, N, i);

Below we have mapped the real component of the result of the forward transform to the height of the terrain.

Applying the forward transform and interpreting the real values as the height map.

Applying the forward transform and interpreting the real values as the height map.

In the third step we apply the filter to curb the higher frequencies.

	// apply 1/(f^p) filter
	double f;
	for (int j = 0; j < N; j++) {
		for (int i = 0; i < N; i++) {
			index = i + j * N;
			f = sqrt((i-N/2)/float(N) * (i-N/2)/float(N) + (j-N/2)/float(N) * (j-N/2)/float(N));
			f = f < 1.0/N ? 1.0/N : f;
			c[index] = c[index]*(1.0f/ pow(f, power));
		}
	}

The result of applying the low-pass filter is below. Note that we have four spikes because we have tiled our map in the and directions.

Applying the 1/f^beta filter in the frequency domain. Note that we have four spikes because we have tiled the map.

Applying the 1/f^beta filter in the frequency domain. Note that we have four spikes because we have tiled the map.

Lastly, we transform back into the spatial domain by applying the inverse transform.

// transform back to spatial domain
for (int j = 0; j < N; j++) {
	fft.fft(c, c, 1, j * N);
	for (int i = 0; i < N; i++) {
		c[i + j * N] = c[i + j * N] * (1.0f/N);
	}
}
for (int i = 0; i < N; i++) {
	fft.fft(c, c, N, i);
	for (int j = 0; j < N; j++) {
		c[i + j * N] = c[i + j * N] * (1.0f/N);
	}
}

We interpret the components of the result of the inverse transform as our height map yielding the rendering below.

In the last step we apply the inverse transform to return to the spatial domain.

In the last step we apply the inverse transform to return to the spatial domain.

In the cTerrain object we have a simple color method to generate a color for a vertex given the height at that vertex.

vector3 cTerrain::color(double h, double min, double max) {
	// c0 is the color at minimum height, c0 + c1 at maximum
	vector3 c0( 66.0/255,  40.0/255,  18.0/255);
	vector3 c1(189.0/255, 185.0/255, 138.0/255);
	vector3 c = c0 + c1 * ((h - min) / (max - min));
	return c;
}

In the project we have also evaluated the normal for each vertex and applied our lighting calculations in the vertex and fragment shaders.

Hopefully, in the next post we will be able to implement a level of detail algorithm in order to generate much larger terrains and display them at interactive frame rates. The source code for this project is available as always. The code is a bit sloppy this time around, but have a look at it. If you have any comments, let me know.

Download this project: terrain.tar.bz2

One thought on “Using Fourier synthesis to generate a fractional Brownian motion surface

  1. I have a question about your method. Some algorithms and definitions of Brownian surfaces have a parameter ɑ that determines the variance of two neighboring points. If one point's height is f(x,y), and we let h and k be increments, then the resulting random values f(x+h, y+k) - f(x,y) have a Gaussian distribution with a variance of (h^2 + k^2)^ɑ. The Hausdorff dimension of the resulting surface is then d = 3-ɑ.

    How does your method account for this ɑ parameter?

Leave a Reply

Your email address will not be published. Required fields are marked *