Using Gaussian blurring on heightmaps

In my previous blog post, I talked about a personal project where I am generating heightmaps. I ended that post writing about how blurring has been important to me in generating the maps. Here is more information on the blurring.

First, I found out that rather then generating a lot of points, I could generate a lot less points and apply heavier blurring to get a satisfactory result much more rapidly.

For my blurring, I ended up using a Gaussian blur (sometimes called Gaussian filter). This type of filter is often used with height maps as it gives good natural looking results. Something very similar to erosion.

A Gaussian filter does this by giving more weight to the pixel being blurred rather than giving equal weight to the center pixel and each of it’s neighbors as is the case in the mean filter.

The problem

The problem I was facing was that my heightmaps needed to be blurred. Rather than create a blurring class that allows me to specify various degrees of blurring, work on images or even specify different types of blur filter (ie: gaussian, mean, etc.), I wrote the simplest thing that could solve my problem.

Applying YAGNI here has allowed me to complete this quickly and get on with the rest of the project. Years earlier I would have started a cool blurring library which I wouldn’t had the time to complete.

The solution

The Gaussian filter is a filter where the values of the kernel are calculated using the Gaussian function to produce values falling in a normal distribution. Like other filter (ie: the mean filter), the Gaussian filter works with a kernel which is a matrix.

At it’s simplest, a non-gaussian kernel could look something like this :

0.1, 0.1, 0.1,
0.1, 0.1, 0.1,
0.1, 0.1, 0.1,

We use the kernel to recalculate the value of each pixel of an image (or more simply of an array) by multiplying its actual value with the center position in the kernel and then multiplying the value above that pixel with the value above in the kernel and so on.

After we have we have multiplied all the nine values (for our 3 x 3 kernel), we add them up and store them in a single cell in a second array. We must do this for all points in our image or array.

When doing this, don’t forget to work with a secondary array so as to not modify the original. This is because you will be reading the original again when calculating the other pixels surrounding the first one you calculated.

For a Gaussian filter we would calculate the values of the kernel using a Gaussian distribution formula. I will not go into the details of the mathematical formula. You can find all the nitty gritty on the links at the beginning of the post or the references at the end.

Tips

Rather than create a 2d kernel, it is also possible to use a modified formula to create a 1d kernel. Running this 1d kernel twice, once horizontally and once vertically, produces a result equivalent to running the 2d kernel.

The difference is that running the 1d kernel twice is faster than running the 2d kernel once as it requires less operations. A neat trick.

The degree of blurring (how heavy to blur) is controlled by changing the number of standard deviations in our Gaussian distribution formula.

Rather than calculate a new kernel every time with differing values, another trick is to use the same kernel to blur our image many times in succession to obtain heavier blurring. If this is parallelized, this can also lead to speed optimizations.

My implementation

For my implementation, I have used a precomputed 1d kernel I have found here and run it both horizontally and vertically.

Precomputed 1d kernel: [0.006, 0.061, 0.242, 0.383, 0.242, 0.061, 0.006].

If I need light blurring I will blur once and if I need heavier blurring (my most common scenario) I will blur twice in a row. Also I applied the blurring on an array of int and only created the image once the blurring was done which is much simpler than working with an image file.

Fast and simple to implement.

An issue I had to solve while implementing the blur was what to do when the filter/kernel is working with values on the edge of the array.

In these cases there are many possible solutions but I chose to use the value of the center pixel to fill in for the values of the non existing pixels.

The other choice would have been to wrap around the edge of the array, or to use a predefined value like 0 or 255. Using the center pixel gave pretty good results.

Code

For the following code note that:

1- I am using a 1d array for my map to represent a 2d array which is unnecessary. If you prefer to use a 2d array, you should go with that.

2- The FilterOutOfBoundsSpecification which is in another file has been appended to the end of this code sample.

3- The code could be refactored as both compute_x_filtered_value and compute_y_filtered_value are almost identical and could be one function.

4- Even though the kernel is symmetrical, I store it likewise:

@filter_kernel = [[-3, 0.006], [-2, 0.061], [-1, 0.242],
				 [0, 0.383], [1, 0.242], [2, 0.061], [3, 0.006]]

But I could save space by storing it this way:

@filter_kernel = [[0, 0.383], [1, 0.242], [2, 0.061], [3, 0.006]]

This would necessitate a small change to the code. I decided to keep the longer version of the kernel because I feel it’s more readable and easier to understand. Saving a few characters isn’t worth it if it complicates things.

5- For the sake of brevity I did not include the tests code.

Enough preamble here is the code:

require_relative './filter_out_of_bounds_specification'

class GaussianFilter
	def initialize
		# the filter kernel value pairs. The first element in each pair represent
		# the distance from the central pixel and the second element the
                # weighted filter value
		@filter_kernel = [[-3, 0.006], [-2, 0.061], [-1, 0.242],
						[0, 0.383], [1, 0.242], [2, 0.061], [3, 0.006]]
	end

	# this function assumes the array represents a square 2d matrix
	def filter(array)
		@size = Math.sqrt(array.length)

		intermediate_filtered_array = Array.new(array.length, 0)
		filtered_array = Array.new(array.length, 0)

		# filter horizontally
		(0...@size).each do |y|
			# @size is used since we assume the array is square 2d matrix
			(0...@size).each do |x|
				intermediate_filtered_array[x + y * @size] = compute_x_filtered_value(array, x, y)
			end
		end

		# filter vertically
		(0...@size).each do |x|
			# @size is used since we assume the array is square 2d matrix
			(0...@size).each do |y|
				filtered_array[x + y * @size] = compute_y_filtered_value(intermediate_filtered_array, x, y)
			end
		end

		return filtered_array
	end

private

	def	compute_x_filtered_value(array, x, y)
		computed_value = 0.0

		@filter_kernel.each do |filter_pair|
			if FilterOutOfBoundsSpecification.satisfied_by?(x, @size, filter_pair[0])
				offset = 0
			else
				offset = filter_pair[0]
			end

			computed_value += filter_pair[1] * array[x + offset + y * @size]
		end

		return computed_value.round
	end

	def	compute_y_filtered_value(array, x, y)
		computed_value = 0.0

		@filter_kernel.each do |filter_pair|
			if FilterOutOfBoundsSpecification.satisfied_by?(y, @size, filter_pair[0])
				offset = 0
			else
				offset = filter_pair[0]
			end

			computed_value += filter_pair[1] * array[x + (offset + y) * @size]
		end

		return computed_value.round
	end
end

class FilterOutOfBoundsSpecification
	def self.satisfied_by?(z, size_z, offset)
		return true if z + offset >= size_z

		return true if z + offset < 0

		false
	end
end

Sources

Here are the sources that I have used to write my Gaussian filter and this blog post, in order of relevance to my project and this article.

My post does not talk about the math and science side of the Gaussian distribution and about image processing as both subjects aren’t strengths of mine. So I encourage you to read the following links for more info on these subjects.

Gaussian Smoothing

CS 384G: Image processing

Gaussian filter, or Gaussian blur

Wikipedia Gaussian filter

Wikipedia Gaussian blur

Generating heightmaps using particle deposition

Generating heightmaps using particle deposition

I recently started a small personal project where I wanted to generate heightmaps. More specifically, heightmaps that represent islands.

I remembered seeing something about this in my copy of Game Programming Gems, a book I had purchased years ago. I did not have the cd-rom with the code samples anymore, but this turned out great, as it forced me to work out a particle deposition algorithm by myself with the book’s textual description rather than just having the code handed to me.

The algorithm

Particle deposition is well suited to generate heightmaps that represent volcanic islands but can also be used for other types of maps.

The algorithm works by depositing particles on a flat surface. You can imagine the particles being dropped by a canon from overhead unto the surface. After each particle is shot, the surfaced will be agitated so that the particle will settle in its final location.

A particle is considered stable if it’s exactly one higher than all of its neighbors. Otherwise, when it is agitated it will fall into one of the lower neighboring positions. You can have a fixed or random order to determine which lower neighbor it will fall in.

The stability radius, or how far from a particle the algorithm will look for a lower neighbor, will help control the slope at the drop points. For these values I often use values from 1 to 3. The greater the value, the gentler the slope.

The other variables to control the generation of the heightmaps are the number of drop points, the number of particle per drop points and the number of passes.

The number of drop points is the number of times the particle canon will be moved during the process.

The number of particle per drop is how many particles will be launched at each drop point. I have obtained better results by using a random number that falls between a min and a max value that I can define.

Finally the number of passes. I don’t believe this is in the traditional version of the algorithm so you could probably do away with it. What it does is as the number of passes increase, the drop points get more focused towards the center of the map, the passes include less and less drop points and less and less particles are found in each drop points.

Pseudocode

Here is some pseudocode that gives an overview of the algorithm.

for each pass
    create drops // here I create a random number of particle objects
    set next drop point location // make sure it's not on the edge of the map
    
    for each drop
        for each particle
            drop  // increase value of int array at location
            agitate  // move particle to it's final resting place

     change variables for next pass

blur final results

Blur

After some experimentation, the secret for me has been blurring. Blurring is typically used to smooth out the maps at the end, but I have found out that rather then generating a lot of points, I could generate a lot less points and apply heavier blurring to get a satisfactory result much more rapidly. Blurring, if done correctly will allow to generate the map faster then generating more particles. This has been a real difference maker in my implementation.

Here is an example of how blurring can turn a map that looks sparse and unsatisfactory into something decent and usable.

I will cover how I used blurring in a following blog post which should come out in the following days.