Erosion

erosion.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>

#define RADIUS (1)
#define DIAMETER (RADIUS * 2 + 1)
#define NNEIGHBORS (DIAMETER * DIAMETER)
#define ITERATIONS 64

#define PI 3.14159265359f

#include "image.h"

#define max(x, y) ((x > y) ? x : y)
#define min(x, y) ((x < y) ? x : y)
#define SATURATE(x) (max(min(255, x), 0))

float noise(int x, int y, int o, int seed)
{
	uint32_t hash = 0;
	uint32_t key[4] = {seed, x, y, o};
	int i = 0;

	for (i = 0; i < 4; i++) {
		hash += key[i];
		hash += (hash << 10);
		hash ^= (hash >> 6);
	}

	hash += (hash << 3);
	hash ^= (hash >> 11);
	hash += (hash << 15);

	return (float) hash / UINT32_MAX;
}

inline float lerp(float x1, float x2, float y1, float y2, float x)
{
	return (x - x1) / (x2 - x1) * (y2 - y1) + y1;
}

float perlin(int x, int y, int firstOctave, int lastOctaves, int size, int seed)
{
	float sum = 0.0f;
	int o;
	int x1, x2, y1, y2;
	float m1, m2, m;
	float mult = 0.5f;

	size >>= firstOctave;

	for (o = firstOctave; o <= lastOctaves; ++o) {
		size >>= 1;
		mult *= 0.5f;

		x1 = x - x % size;
		y1 = y - y % size;
		x2 = x1 + size;
		y2 = y1 + size;

		m1 = lerp(x1, x2, noise(x1, y1, o, seed), noise(x2, y1, o, seed), x);
		m2 = lerp(x1, x2, noise(x1, y2, o, seed), noise(x2, y2, o, seed), x);
		m = lerp(y1, y2, m1, m2, y);

		sum += m * mult;
	}

	return sum;
}

int main(int argc, char *argv[])
{
	Image *img = Image_new(256, 256);
	float *height, *water, *sediment;
	int x, y, iteration;
	int offset[NNEIGHBORS];
	float eroded, evaporated;
	int dx, dy, i, j, k = 0;
	int used[NNEIGHBORS];
	int nb_used;
	float di[NNEIGHBORS];
	float dtotal;
	float a, dmi;
	float avg_a;
	float ai[NNEIGHBORS];
	float da;
	float w;
	float dwi[NNEIGHBORS];
	float dmtotal;

	float minnoise = FLT_MAX;
	float maxnoise = FLT_MIN;

	float rain_amount = 0.04f;
	float solubility = 0.2f;
	float evaporation = 0.5f;
	float capacity = 0.01f;

	height = (float *) malloc(sizeof(float) * img->width * img->height);
	water = (float *) malloc(sizeof(float) * img->width * img->height);
	sediment = (float *) malloc(sizeof(float) * img->width * img->height);

	for (y = 0; y < img->height; ++y) {
		for (x = 0; x < img->width; ++x) {
			i = y * img->width + x;
			height[i] = (float) img->data[i];
			water[i] = 0.0f;
			sediment[i] = 0.0f;
			height[i] = perlin(x, y, 2, 7, 256, 9411456) * 255.0f;
			height[i] *= powf(sinf((float) y / img->height * PI) * sinf((float) x / img->width * PI), 2.0f);
			minnoise = min(minnoise, height[i]);
			maxnoise = max(maxnoise, height[i]);
		}
	}

	for (y = -RADIUS; y <= RADIUS; ++y) {
		for (x = -RADIUS; x <= RADIUS; ++x) {
			offset[k++] = y * img->width + x;
		}
	}

	float normalization = 1.0f / (maxnoise - minnoise);

	for (iteration = 0; iteration < ITERATIONS; ++iteration) {
		for (y = RADIUS; y < img->height - RADIUS; ++y) {
			for (x = RADIUS; x < img->width - RADIUS; ++x) {
				i = y * img->width + x;

				if (iteration == 0) {
					height[i] = (height[i] - minnoise) * normalization * 255.0;
				}

				water[i] += rain_amount;

				eroded = solubility * water[i];
				height[i] -= eroded;
				sediment[i] += eroded;

				w = water[i];
				a = w + height[i];
				dtotal = 0.0f;
				avg_a = 0.0f;
				nb_used = 0;
				dmtotal = 0.0f;

				for (k = 0; k < NNEIGHBORS; ++k) {
					j = i + offset[k];

					ai[k] = height[j] + water[j];

					if (ai[k] > a && j != i) {
						used[k] = 0;
						continue;
					}

					nb_used++;
					used[k] = 1;
					avg_a += ai[k];
					di[k] = a - ai[k];
					dtotal += di[k];
				}

				avg_a /= nb_used;
				da = a - avg_a;

				for (k = 0; k < NNEIGHBORS; ++k) {
					j = i + offset[k];
					dmi = 0.0f;

					if (!used[k])
						continue;

					if (dtotal == 0.0f) {
						dwi[k] = 0.0f;
					} else {
						dwi[k] = min(w, da) * di[k] / dtotal;

						water[j] += dwi[k];
						water[i] -= dwi[k];

						dmi = sediment[i] * dwi[k] / w;
						sediment[j] += dmi;
						dmtotal += dmi;
					}
				}

				sediment[i] -= dmtotal;
				water[i] -= evaporation * water[i];

				if (sediment[i] > capacity * water[i]) {
					height[i] += sediment[i] - capacity * water[i];
					sediment[i] = capacity * water[i];
				}
			}
		}
	}

	for (y = 0; y < img->height; ++y) {
		for (x = 0; x < img->width; ++x) {
			int i = y * img->width + x;
			img->data[i] = (uint8_t) SATURATE(height[i]);
		}
	}

	Image_save(img, "test.pgm");
	Image_delete(img);

	return 0;
}