Add simple 2d perlin noise function

Makefile.am

@@ -53,6 +53,8 @@ UTIL_FILES   =		\
     matrix.cpp		\
     matrix.hpp		\
     nocopy.hpp		\
+    perlin.cpp		\
+    perlin.hpp		\
     point.cpp		\
     point.hpp		\
     pool.cpp		\

perlin.cpp

@@ -0,0 +1,105 @@
+#include "perlin.hpp"
+
+#include "debug.hpp"
+#include "maths.hpp"
+
+#include <cmath>
+#include <limits>
+
+
+namespace lerp {
+    double
+    linear (double a, double b, double weight) {
+        check (weight >= 0 && weight <= 1.0);
+        return (1.0 - weight) * a + weight * b;
+    }
+
+
+    double
+    cosine (double a, double b, double weight) {
+        check (weight >= 0 && weight <= 1.0);
+        double ft = weight * PI;
+        double f  = (1.0 - cos (ft)) * 0.5;
+
+        return a * (1.0 - f) + b * f;
+    }
+}
+
+
+// Just a random generator [-1.0, 1.0]
+static float
+generate (intmax_t x, intmax_t y) {
+    uint32_t n = x + 257 * y;
+    n = (n << 13U) ^ n;
+    return (1.0 - ((n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0f);
+}
+
+
+double
+smooth2d (intmax_t x, intmax_t y) {
+    double corners = (generate(x-1, y-1)+generate(x+1, y-1)+generate(x-1, y+1)+generate(x+1, y+1) ) / 16.0;
+    double sides   = (generate(x-1, y)  +generate(x+1, y)  +generate(x, y-1)  +generate(x, y+1) ) /  8.0;
+    double center  =  generate(x, y) / 4.0;
+    
+    return corners + sides + center;
+}
+
+
+double
+curve (double t) {
+    return t * t * (3.0 - 2.0 * t);
+}
+
+
+double
+sample2d (double x, double y) {
+    intmax_t x_int = intmax_t (x);
+    double   x_fac = x - x_int;
+    intmax_t y_int = intmax_t (y);
+    double   y_fac = y - y_int;
+
+    double p0 = smooth2d (x_int,     y_int);
+    double p1 = smooth2d (x_int + 1, y_int);
+    double p2 = smooth2d (x_int,     y_int + 1);
+    double p3 = smooth2d (x_int + 1, y_int + 1);
+
+    double i0 = lerp::cosine (p0, p1, curve (x_fac));
+    double i1 = lerp::cosine (p2, p3, curve (x_fac));
+
+    double value = lerp::cosine (i0, i1, curve (y_fac));
+    return value;
+}
+
+
+double
+sigmoid (double x) {
+    return -1.0 + 2.0 / (1.0 + exp (-2.0 * x));
+}
+
+
+double
+perlin2d (double x, double y) {
+    size_t levels = 16;
+    double persistence = 0.7;
+    double total = 0.0;
+
+    for (size_t i = 0; i < levels; ++i) {
+        double frequency = powf (2.0, i);
+        double amplitude = powf (persistence, i);
+
+        total += sample2d (x * frequency, y * frequency) * amplitude;
+    }
+
+    return sigmoid (1.0 + total);
+}
+
+
+void
+perlin2d (uint8_t *restrict pixels, size_t width, size_t height) {
+    for (size_t y = 0; y < height; ++y)
+        for (size_t x = 0; x < width; ++x) {
+            double p = perlin2d (x / 10.0, y / 10.0);
+            pixels[x + y * width] = p * std::numeric_limits<uint8_t>::max ();
+        }
+}
+

perlin.hpp

@@ -0,0 +1,28 @@
+/*
+ * This file is part of libgim.
+ *
+ * libgim is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ * 
+ * libgim is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ * details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with libgim.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Copyright 2011 Danny Robson <danny@nerdcruft.net>
+ */
+
+#ifndef __UTIL_PERLIN_HPP
+#define __UTIL_PERLIN_HPP
+
+#include <cstdint>
+#include <cstdlib>
+
+void perlin2d (uint8_t *restrict pixels, size_t width, size_t height);
+
+#endif