Generating Realistic Terrain
Simulating wind and erosion
Robert Krten
Robert is a contract programmer in Kanata, Ontario. You can contact him at [email protected].
Realistic landscapes are the bread and butter of the computer graphics used in movies, video games, multimedia applications, and similar simulations. Unfortunately, generating landscapes and terrain that look real isn't always a straightforward process. However, fault generation, the technique I present here, is easy to grasp and implement--and fast. Fault generation is a rough simulation of the way some mountains and other geological features are formed in nature.
For instance, imagine that the landscape that you wish to generate is represented by a two-dimensional array in the computer's memory. The value at any given x, y position within that array represents the height at that point in the landscape. To make the landscape look interesting and realistic, each height value should bear some relationship to its neighbor's height.
To visualize how this can be accomplished, start with a flat terrain (all height values set to 0). Next, draw an imaginary straight line through any given part of the array, such that it passes entirely through the array. Then, change all of the values on one side of that line to be lower by a certain amount, and all of the values on the other side of that line to be higher by a certain amount. This results in a landscape with a big fault running through some part of it. While this may be interesting as a first step, it certainly doesn't offer much variety or realism.
To enhance realism, you need to repeat the process several times, (using a different imaginary line each time), decreasing the amount by which the landscape changes (the height) with each iteration. This causes the landscape to have a large shift in one place (corresponding to the first fault), with smaller "detailing" shifts throughout. You can get creative with the function that you use for the decrease in height, but I've used a 1/x-height reduction with good results. This way, the first fault line causes the landscape to change in height by the same amount as there are iterations, the next fault line changes in height by 1 less, and so on, until the last fault line changes in height by one unit.
This results in a rugged terrain, with sharp corners and sudden changes in certain places. This terrain is usable for some applications without further modification, but could be made more realistic by smoothing out the rough edges.
Smoothing It Out
There are a number of ways to smooth out the terrain. By controlling the fault-decay factor, you can have one large feature, and relatively small subsequent features. Applying many such small features to the landscape makes it statistically likely that the large feature will become "broken down" over time, with fewer sharp edges. This approach may require many iterations, however, consuming a fair amount of CPU time.
An alternative is to apply a digital filter over the generated landscape, smoothing out the rough spots. This allows the entire landscape to be generated and smoothed in one final pass. Even though the digital-filter algorithm is somewhat expensive in CPU time, it is still a good solution because it happens only once.
A more realistic filtering approach is simulation of erosion. Imagine that after the first major land shift, an innocuous filter is passed over the terrain. It is just enough to smooth out a few of the rough spots in the terrain in the neighborhood of the first fault. When the second land shift occurs, the land being shifted is already smoother than it would have been. Again, a filter is run over the newly shifted land. This produces the equivalent of geological wind erosion.
The "low-pass" filter (a single-constant FIR filter) in Listing Two operates by propagating a certain small portion of the previous sample into the current sample. This is repeated for all samples. The two filter types (the two-pass and four-pass, both in Listing Two) differ only in that the two-pass sweeps along the x-axis once and then the y-axis once, whereas the four-pass sweeps along the x-axis in one direction and then the other direction (two passes) and performs the same operation for the y-axis (two more passes, for a total of four). With the two-pass filter, a landscape shift is introduced into the array, whereas with the four-pass this shift is not apparent.
The two-pass filter can be thought of as simulating constant wind erosion, where the wind is always coming from the same direction. The four-pass filter simulates rain erosion, where the rain is (obviously) coming from the top, and eroding particles away in all directions.
Allocating the Landscape Memory
Memory was an interesting side issue during the development of the fault generator. I have put together a 2-D calloc library call (see Listing Two) that allows the landscape size to be determined at run time, rather than compile time. This should be especially helpful for systems where memory is tight and you want to get the biggest array that will fit. The technique used also ensures that the 64-Kbyte segmentation barrier will not be reached (unless your array is bigger than 16x16 Kbytes, in which case you will more than likely first run out of physical memory and processor time). An advantage of allowing the landscape size to be determined at run time is that you can batch up a large number of different-sized landscapes to be processed, then go home for the evening, without having to recompile each time.
There are a number of ways of declaring (and allocating storage for) 2-D arrays in C. For instance, a statement such as Example 1(a) results in a different memory layout than Example 1(b). However, both can be addressed as in Example 1(c). The first statement declares an array-of-array and allocates all of the required storage in a contiguous chunk of memory. This can easily exceed 64 Kbytes. On systems without the 64-Kbyte barrier, it can still be a problem, as it may exceed the largest chunk of contiguous free memory available. The second statement declares a pointer-to-a-pointer and typically allocates four bytes. The terrain generator presented here uses the second statement, in conjunction with a library, for allocating and freeing the structure.
The approach used to allocate memory is a two-stage allocation. Assume that a 200*400*(sizeof int) array is being allocated. In the first stage, the "backbone" is allocated with space for 200 pointers to integers. (This typically consumes 200*sizeof (int*), or 800 bytes.) In the second stage, one 400-integer array is allocated for each backbone pointer and stored in that pointer.
The end result is an array of 200 pointers, each of which points to a different 400-element array of integers. In terms of overhead, this introduces the 800 extra bytes for the backbone array.
C allows the land [x][y] addressing style to work because the compiler is aware of the details of the base type of "land" (that is, it is a pointer-to-a-pointer). The compiler looks at land [x] and generates code to reference the xth pointer (of the 200). It then generates code that indexes into the yth location of that pointer, thus referencing the given array element. By allocating just a little more memory than I actually need, I can store some information at the beginning of the allocated-memory array, then return a pointer to just after that header.
The extra information in the header (the x- and y-array size allocated and the size of the individual element) is stored by the ECalloc2d routine when the array is created, and it is especially useful in EFree2d, ERead2d, and EWrite2d. Without this header, I'd have to pass the information around all of the time or maintain it as a global structure.
Another feature of the pointer-to-a-pointer approach is that the entire array can be assigned to a variable using the C assignment operator, rather than a series of nested for loops. For example, after calling the digital filter in the main loop, you need to swap the input and output arrays. This happens in the procedure fault using just three assignments.
Storing the Landscape to Disk
The easiest way to store the landscape to disk is to write out all of the elements, using two nested for loops. With this approach, the whole first row will be written out (in column order), then the second row, and so on.
Agreeing on a common-output format makes it easy to have other utilities process (or even generate) the terrain data. For example, you can write an alternate filtering or display program. This underscores another advantage of dynamically allocated array sizes. By writing out the array size (xxy) as the first few bytes of the file, you can read in differently sized files for processing in other utilities, without having to recompile all of the utilities for the new size.
The Source Code
The fault-generation source code consists of the makefile (Listing One), common.c (Listing Two), common.h (Listing Three), and fault.c (Listing Four). I've developed this code under QNX 4.2 with the Watcom 9.5 C compiler.
The makefile contains C-compiler flags (that you may change or ignore), linker flags, and dependencies. As defined, fault.c is the main module, with common.c containing the filtering and 2-D manipulation routines. The ANSI C prototypes for common.c are contained in common.h.
The fault.c module contains the main routine, which calls the command-line option parser (optproc), allocates the landscape arrays (ECalloc), and calls the fault generator (fault) with the number of iterations to be performed.
The fault routine calls the single fault-line generator generate_line (with the height of the land shift) and the optional wind-erosion filter. At the end of fault, the optional final filtering is performed, and the file is written to disk.
To use the fault generator, type the name of the executable followed by the name of the output file you wish to generate. For example, fault terrain.xy will invoke the fault generator with the defaults and generate an output file called terrain.xy.
Possible Enhancements
There are a number of enhancements you can add to the characteristics of the terrain generated by this program.
For one thing, you can change how the fault heights are determined by establishing a number of discrete sizes and choosing one of those each time, perhaps with a weighted probability. For example, if you were generating a 200-iteration terrain, you would take 10 100-unit heights, 150 20-unit heights, and the rest (40) in 2-unit heights.
The land movement on each side of the fault need not be constant. You can scale the fault height by the distance from the fault line; for example, the further away a point is from the fault line, the less it is affected. This is also closer to the way that things happen geologically--an earthquake in San Francisco usually has no effect in Pocatello, Idaho, 1000 kilometers away.
As previously mentioned, digital filters, such as FIR or IIR filters, can lead to some quite spectacular effects. For another effect, run the landscape through an FFT one row at a time, chop out a portion of the frequency domain, and then reconstitute it. Regardless of what you do to improve the algorithm, I'd like to hear about it. I can be reached on the Internet at [email protected].
Example 1: Declaring and allocating storage for 2-D arrays.
(a) int land [200][400]; (b) int **land; (c) land [x][y] = stuff;
Listing One
# Makefile for fractal fault routine
# For QNX 4.2, but should be reasonably portable
# 1993 11 25 R. Krten released for publication / public use
CFLAGS = -4 -Oxr -w9 -mf
fault: common.o fault.o
cc $(CFLAGS) -o fault fault.o common.o
fault.o : Makefile fault.c common.h
cc $(CFLAGS) -c fault.c
common.o : Makefile common.c common.h
cc $(CFLAGS) -c common.c
Listing Two
/* common.c
* QNX 4
* (C) Copyright 1993 by Robert Krten, all rights reserved.
* This module contains common utilities for the X * Y fault programs.
* 1993 10 26 R. Krten created
* 1993 11 25 R. Krten released for publication / public use
*/
#include <stdio.h>
#include <stdlib.h>
typedef struct {
int xSize; /* number of backbone entries */
int ySize; /* number of entries in each backbone entry */
int eSize; /* size of each entry */
} MHead;
extern char *progname;
extern int **land;
extern int dimensions [2];
/* Two-dimensional support routines:
* ECalloc2d (x size, y size, size of each member)
* EFree2d (pointer to allocated array to free)
* ERead2d (FILE pointer, pointer to array to read)
* EWrite2d (FILE pointer, pointer to array to write)
* The above routines operate on two dimensional arrays, based upon the
* "pointer to pointer" type. The allocation routine first allocates
* a backbone (consisting of "x size" number of pointers + a header, and
* then fills in the backbone with "y sized" members. The free routine
* frees each backbone member in turn, and then the whole backbone itself,
* including the header. The read and write routines read and write the
* arrays from and to disk.
*/
void **
ECalloc2d (x, y, esize)
int x;
int y;
int esize;
{
void **ptr; /* pointer to allocated memory */
MHead *mptr; /* pointer to memory header */
int i;
if ((ptr = calloc (1, sizeof (MHead) + x * sizeof (void *))) == NULL) {
fprintf (stderr, "%s: out of memory on first allocation\n", progname);
exit (1);
}
mptr = (MHead *) ptr;
mptr -> xSize = x;
mptr -> ySize = y;
mptr -> eSize = esize;
ptr = (void *) (mptr + 1);
for (i = 0; i < x; i++) {
if ((ptr [i] = calloc (y, esize)) == NULL) {
fprintf (stderr, "%s: out of memory (at [%d])!\n", progname, i);
exit (1);
}
}
return (ptr);
}
void
EFree2d (ptr)
void **ptr;
{
MHead *mptr; /* pointer to memory header */
int x, y;
int i;
mptr = (MHead *) ptr - 1;
x = mptr -> xSize;
y = mptr -> ySize;
for (i = 0; i < x; i++) {
free (ptr [i]);
}
free (ptr);
}
ERead2d (fp, l)
FILE *fp;
void **l;
{
MHead *mptr; /* pointer to memory header */
int x, y, esize;
int i;
mptr = (MHead *) l - 1;
x = mptr -> xSize;
y = mptr -> ySize;
esize = mptr -> eSize;
for (i = 0; i < x; i++) {
fread (l [i], y, esize, fp);
}
}
EWrite2d (fp, l)
FILE *fp;
void **l;
{
MHead *mptr; /* pointer to memory header */
int x, y, esize;
int i;
int dim [2];
mptr = (MHead *) l - 1;
dim [0] = x = mptr -> xSize;
dim [1] = y = mptr -> ySize;
esize = mptr -> eSize;
fwrite (dim, 2, sizeof (int), fp);
for (i = 0; i < x; i++) {
fwrite (l [i], y, esize, fp);
}
}
/* The filter algorithm -- This implements a single-constant FIR filter. */
filter (input, output, kx, ky, flag)
int **input;
int **output;
double kx, ky;
int flag;
{
MHead *mptr; /* pointer to memory header */
double acc;
double acckx, accky;
register x, y;
/* we assume that dimensions of "input" == dimensions of "output" */
mptr = (MHead *) input - 1;
/* first pass X direction */
printf ("1"); fflush (stdout);
accky = 1. / (1. - ky);
for (y = 0; y < mptr -> ySize; y++) {
acc = input [0][y] * accky;
for (x = 0; x < mptr -> xSize; x++) {
output [x][y] = acc / accky;
acc = acc * ky + input [x][y];
}
}
/* second pass X direction */
printf ("2"); fflush (stdout);
if (flag == '4') {
for (y = 0; y < mptr -> ySize; y++) {
acc = input [mptr -> xSize - 1][y] * accky;
for (x = mptr -> xSize - 1; x >= 0; x--) {
output [x][y] += acc / accky;
acc = acc * ky + input [x][y];
}
}
}
/* first pass Y direction */
printf ("3"); fflush (stdout);
acckx = 1. / (1. - kx);
for (x = 0; x < mptr -> xSize; x++) {
acc = input [x][0] * acckx;
for (y = 0; y < mptr -> ySize; y++) {
output [x][y] += acc / acckx;
acc = acc * kx + input [x][y];
}
}
/* second pass Y direction */
printf ("4"); fflush (stdout);
if (flag == '4') {
for (x = 0; x < mptr -> xSize; x++) {
acc = input [x][mptr -> ySize - 1] * acckx;
for (y = mptr -> ySize - 1; y >= 0; y--) {
output [x][y] += acc / acckx;
acc = acc * kx + input [x][y];
}
}
}
/* averaging for 2 or 4 passes */
printf ("A"); fflush (stdout);
for (x = 0; x < mptr -> xSize; x++) {
for (y = 0; y < mptr -> ySize; y++) {
output [x][y] /= flag - '0';
}
}
printf ("\r"); fflush (stdout);
}
Listing Three
/* common.h * QNX 4 * (C) Copyright 1993 by Robert Krten, all rights reserved. * This module contains the common utility functions for the fault * handlers. * 1993 10 26 R. Krten created * 1993 11 25 R. Krten released for publication / public use */ /* prototypes */ void **ECalloc2d (int, int, int); void ERead2d (FILE *, void **); void EWrite2d (FILE *, void **); void EFree2d (void **); void filter (int **, int **, double, double, int);
Listing Four
/* fault.c
* QNX 4
* (C) Copyright 1988 by Robert Krten, all rights reserved.
* 1988 05 16 R. Krten created
* 1991 08 29 R. Krten ported to QNX 2
* 1993 02 20 R. Krten ported to QNX 4/Windows
* 1993 10 16 R. Krten expanded array (novelty of 32bit)
* 1993 10 31 R. Krten allow filtering with generator (wind erosion)
* 1993 11 25 R. Krten released for publication / public use
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <time.h>
#include <math.h>
#include "common.h"
static void optproc (int, char **);
char *progname = "fault"; /* for diagnostics */
char outfile [256]; /* for output of fault program */
int **land; /* terrain */
int **outland; /* for filtering */
int dimensions [2]; /* holds size of terrain */
int nit; /* Number of ITerations */
double kx, ky; /* FIR filter constants */
int postFilter; /* flag indicating post-filtering */
int windFilter; /* flag indicating wind-erosion filtering */
int filtering; /* flag indicating any filtering */
int filterType; /* flag indicating type of filter (2 or 4 pass) */
main (argc, argv)
int argc;
char **argv;
{
optproc (argc, argv);
srand (time (NULL));
land = ECalloc2d (dimensions [0], dimensions [1], sizeof (int));
if (filtering || postFilter) {
outland = ECalloc2d (dimensions [0], dimensions [1], sizeof (int));
}
fault (nit);
}
fault (n)
int n;
{
int i;
FILE *fp;
int **tmp;
if ((fp = fopen (outfile, "w")) == NULL) {
fprintf (stderr, "%s: couldn't open %s for w\n", progname, outfile);
exit (1);
}
printf ("Fault /%d [%d x %d]\r", nit, dimensions [0], dimensions [1]);
fflush (stdout);
for (i = 0; i < n; i++) {
printf ("Fault %5d\r", i); fflush (stdout);
generate_line (i + 1);
if (windFilter || (postFilter && i == n - 1)) {
filter (land, outland, kx, ky, filterType);
/* swap input and output arrays */
tmp = land;
land = outland;
outland = tmp;
}
}
printf (" \n");
EWrite2d (fp, land);
EFree2d (land);
fclose (fp);
}
generate_line (height)
int height;
{
int x1, y1, x2, y2;
register x3, y3;
int xintercept;
int sign;
double slope;
int landset;
do {
y1 = random (dimensions [1]);
y2 = random (dimensions [1]);
} while (abs (y2 - y1) <= 2);
do {
x1 = random (dimensions [0]);
x2 = random (dimensions [0]);
} while (abs (x2 - x1) <= 2);
slope = (double) (y2 - y1) / (double) (x2 - x1);
sign = (random (10) > 5) ? -1 : 1;
landset = (int) ((double) nit / (double) height);
for (y3 = 0; y3 < dimensions [1]; y3++) {
xintercept = (int) ((double) (y3 - y1) / slope) + x1;
if (xintercept < 0) {
xintercept = 0;
}
if (xintercept > dimensions [0]) {
xintercept = dimensions [0];
}
for (x3 = 0; x3 < xintercept; x3++) {
land [x3][y3] += sign * landset;
}
for (x3 = xintercept; x3 < dimensions [0]; x3++) {
land [x3][y3] -= sign * landset;
}
}
}
random (n)
int n;
{
return (rand () % n);
}
void
usageError ()
{
fprintf (stderr, "%s: error in command line.\n", progname);
fprintf (stderr,
"\n"
"use: fault [options] outFile\n"
"\n"
"where [options] are optional arguments from:\n"
" -n NIT specify number of iterations (default 100)\n"
" -x xSize specify X size (default 256)\n"
" -y ySize specify Y size (default 256)\n"
" -f const X and Y filtering constant (default 0.5)\n"
" -F type specify filter type (2 or 4 pass, default 4)\n"
" -p filter once at end of run\n"
" -w filter each time through\n"
"\n");
exit (1);
}
void
optproc (argc, argv)
int argc;
char **argv;
{
int opt;
if (!argc) {
usageError ();
}
dimensions [0] = dimensions [1] = 256;
nit = 100;
kx = ky = 0.5;
filterType = '4';
postFilter = 0;
windFilter = 0;
while ((opt = getopt (argc, argv, "f:x:y:n:F:pw")) != -1) {
switch (opt) {
case 'f':
kx = ky = atof (optarg);
break;
case 'x':
dimensions [0] = atoi (optarg);
break;
case 'y':
dimensions [1] = atoi (optarg);
break;
case 'n':
nit = atoi (optarg);
break;
case 'F':
if (*optarg == '4') {
filterType = '4';
} else if (*optarg == '2') {
filterType = '2';
} else {
usageError ();
}
break;
case 'p':
postFilter = 1;
break;
case 'w':
windFilter = 1;
break;
default:
usageError ();
break;
}
}
for (; optind < argc; optind++) {
strcpy (outfile, argv [optind]);
}
filtering = windFilter | postFilter;
}
Copyright © 1994, Dr. Dobb's Journal
