heightmap - How can I compute a "steepness" value for height map cells?

Here's a question related to height map generation. From a 2D height map I want to produce a second map, I guess you can call it sort of a derivative map, which shows the "steepness" of all the discrete values of the map. I plan to use it to quickly determine what kinds of objects can be placed on different locations.

In practice you would have many areas where rise/run is a fractional value, since many places do not have 45 degree angle slopes. But running through the pixels of the height map in any direction, the value differences are all integer values.

So what I planned to do in order to represent smooth slopes of varying degrees is, iterate through the pixels row by row, and find the lengths of all the distances where height does not change, and for every change in height, divide that by the length of the distance. Do the same thing column by column and then average the two results for each pixel. I do have to wonder how accurate this is in situations where a lower elevation is surrounded by two walls/plateaus since the rise/run is technically infinity, as it represents a completely vertical climb.

Is this a reasonable way for calculating steepness? I haven't found any resources on how to produce a 2D array of steepness values from height maps, just things on topographic (real-life maps) with the basic rise/run formula.

Answer

You almost have the right idea, but actually it is much simpler than you might think. The simplest way of measuring steepness is by finite differencing:

float GetSteepness(float[,] heightmap, int x, int y)
{
   float height = heightmap[x, y];

   // Compute the differentials by stepping over 1 in both directions.
   // TODO: Ensure these are inside the heightmap before sampling.
   float dx = heightmap[x + 1, y] - height;

   float dy = heightmap[x, y + 1] - height;

   // The "steepness" is the magnitude of the gradient vector
   // For a faster but not as accurate computation, you can just use abs(dx) + abs(dy)
   return sqrt(dx * dx + dy * dy);
}

This is sometimes called the "Sobel Filter." There are more complicated ways of doing it (such as a difference of gaussians filter), but finite differencing should give you a pretty good result! You will get some issues with aliasing when the height is changing very rapidly (such as at walls). You will also get issues when the resolution is too low for finite differencing to return a smooth result.

EDIT: By the way, this is what the sobel filter (bottom) looks like on a simple height map (top)