libs/din: add a perlin noise implementation

This is a 3D noise field addressed as a unit cube.

The caller supplies the resolution of the noise field in three
dimensions.

I've just pulled in my v3f.h here, but it probably makes sense to
later on move vector headers into libs/ and share them.  Later.

It's called din as in noise, because it's shorter than perlin and
noise.

1 files changed, 147 insertions, 0 deletions

diff --git a/src/libs/din/din.c b/src/libs/din/din.c
new file mode 100644
index 0000000..2293c30
--- /dev/null
+++ b/src/libs/din/din.c
@@ -0,0 +1,147 @@
+/* implements a classical perlin noise function */
+/* https://en.wikipedia.org/wiki/Perlin_noise */
+
+#include <stdlib.h>
+
+#include "din.h"
+#include "v3f.h"
+
+typedef struct din_t {
+	int	width, height, depth;
+	v3f_t	grid[];
+} din_t;
+
+
+/* return random number between -1 and +1 */
+static inline float randf(void)
+{
+	return 2.f / RAND_MAX * rand() - 1.f;
+}
+
+
+void din_randomize(din_t *din)
+{
+	int	x, y, z;
+
+	for (z = 0; z < din->depth; z++) {
+		for (y = 0; y < din->height; y++) {
+			for (x = 0; x < din->width; x++) {
+				v3f_t	r;
+
+				r.x = randf();
+				r.y = randf();
+				r.z = randf();
+
+				din->grid[z * din->width * din->height + y * din->width + x] = v3f_normalize(&r);
+			}
+		}
+	}
+}
+
+
+din_t * din_new(int width, int height, int depth)
+{
+	din_t	*din;
+
+	din = calloc(1, sizeof(din_t) + (sizeof(v3f_t) * (width * height * depth)));
+	if (!din)
+		return NULL;
+
+	din->width = width;
+	din->height = height;
+	din->depth = depth;
+
+	din_randomize(din);
+
+	return din;
+}
+
+
+void din_free(din_t *din)
+{
+	free(din);
+}
+
+
+static inline float dotgradient(const din_t *din, int x, int y, int z, const v3f_t coordinate)
+{
+	v3f_t	distance = v3f_sub(&coordinate, &(v3f_t){.x = x, .y = y, .z = z});
+
+	return v3f_dot(&din->grid[z * din->width * din->height + y * din->width + x], &distance);
+}
+
+
+static inline float lerp(float a, float b, float t)
+{
+	return (1.0f - t) * a + t * b;
+}
+
+
+static inline float clamp(float x, float lowerlimit, float upperlimit) {
+	if (x < lowerlimit)
+		x = lowerlimit;
+
+	if (x > upperlimit)
+		x = upperlimit;
+
+	return x;
+}
+
+
+static inline float smootherstep(float edge0, float edge1, float x) {
+	x = clamp((x - edge0) / (edge1 - edge0), 0.f, 1.f);
+
+	return x * x * x * (x * (x * 6.f - 15.f) + 10.f);
+}
+
+
+/* coordinate is in a unit cube of -1...+1 */
+float din(din_t *din, v3f_t coordinate)
+{
+	int	x0, y0, z0, x1, y1, z1;
+	float	i1, i2, ii1, ii2;
+	float	tx, ty, tz;
+	float	n0, n1;
+
+	coordinate.x = 1.f + (coordinate.x * .5f + .5f) * (float)(din->width - 2);
+	coordinate.y = 1.f + (coordinate.y * .5f + .5f) * (float)(din->height - 2);
+	coordinate.z = 1.f + (coordinate.z * .5f + .5f) * (float)(din->depth - 2);
+
+	x0 = floorf(coordinate.x);
+	y0 = floorf(coordinate.y);
+	z0 = floorf(coordinate.z);
+
+	x1 = x0 + 1.f;
+	y1 = y0 + 1.f;
+	z1 = z0 + 1.f;
+
+	tx = coordinate.x - (float)x0;
+	ty = coordinate.y - (float)y0;
+	tz = coordinate.z - (float)z0;
+
+	n0 = dotgradient(din, x0, y0, z0, coordinate);
+	n1 = dotgradient(din, x1, y0, z0, coordinate);
+	tx = smootherstep(0.f, 1.f, tx);
+	i1 = lerp(n0, n1, tx);
+
+	n0 = dotgradient(din, x0, y1, z0, coordinate);
+	n1 = dotgradient(din, x1, y1, z0, coordinate);
+	i2 = lerp(n0, n1, tx);
+
+	ty = smootherstep(0.f, 1.f, ty);
+	ii1 = lerp(i1, i2, ty);
+
+	n0 = dotgradient(din, x0, y0, z1, coordinate);
+	n1 = dotgradient(din, x1, y0, z1, coordinate);
+	i1 = lerp(n0, n1, tx);
+
+	n0 = dotgradient(din, x0, y1, z1, coordinate);
+	n1 = dotgradient(din, x1, y1, z1, coordinate);
+	i2 = lerp(n0, n1, tx);
+
+	ii2 = lerp(i1, i2, ty);
+
+	tz = smootherstep(0.f, 1.f, tz);
+
+	return lerp(ii1, ii2, tz);
+}