--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dmperlin.c	Fri Sep 28 01:54:23 2012 +0300
@@ -0,0 +1,136 @@
+/*
+ * Coherent noise function over 1 or 2 dimensions.
+ */
+#include "dmlib.h"
+#include <math.h>
+
+#define B	(0x100)
+#define BM	(0x0ff)
+#define NP	(12)
+#define N	(0x1000)
+#define NM 	(0x0fff)
+
+
+#define DM_PERLIN_SETUP(i, b0, b1, r0, r1) {	\
+        t = (vec[i] + N);	\
+        b0 = (((int) t) & BM);	\
+        b1 = ((b0 + 1) & BM);	\
+        r0 = (t - ((int) t));	\
+        r1 = (r0 - 1.0f);	\
+        }
+        
+
+#define DM_PERLIN_AT2(rx,ry) ((rx) * q[0] + (ry) * q[1])
+
+
+static int p[B + B + 2];
+static DMFloat g2[B + B + 2][2];
+static DMFloat g1[B + B + 2];
+
+
+static DMFloat dmPerlinDoNoise2(DMFloat vec[2])
+{
+    int bx0, bx1, by0, by1, b00, b10, b01, b11;
+    DMFloat rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
+    int i, j;
+
+    DM_PERLIN_SETUP(0, bx0, bx1, rx0, rx1);
+    DM_PERLIN_SETUP(1, by0, by1, ry0, ry1);
+
+    i = p[bx0];
+    j = p[bx1];
+
+    b00 = p[i + by0];
+    b10 = p[j + by0];
+    b01 = p[i + by1];
+    b11 = p[j + by1];
+
+    sx = DMM_S_CURVE(rx0);
+    sy = DMM_S_CURVE(ry0);
+
+    q = g2[b00];
+    u = DM_PERLIN_AT2(rx0, ry0);
+    q = g2[b10];
+    v = DM_PERLIN_AT2(rx1, ry0);
+    a = DMM_LERP(sx, u, v);
+
+    q = g2[b01];
+    u = DM_PERLIN_AT2(rx0, ry1);
+    q = g2[b11];
+    v = DM_PERLIN_AT2(rx1, ry1);
+    b = DMM_LERP(sx, u, v);
+
+    return DMM_LERP(sy, a, b);
+}
+
+
+static void dmPerlinNormalize2(DMFloat v[2])
+{
+    DMFloat s;
+
+    s = sqrt(v[0] * v[0] + v[1] * v[1]);
+
+    v[0] /= s;
+    v[1] /= s;
+}
+
+
+void dmPerlinInit(void)
+{
+    int i, j, k;
+
+    srand(32);
+    
+    for (i = 0; i < B; i++)
+    {
+        p[i] = i;
+        g1[i] = (DMFloat) ((rand() % (B + B)) - B) / B;
+
+        for (j = 0; j < 2; j++)
+            g2[i][j] = (DMFloat) ((rand() % (B + B)) - B) / B;
+
+        dmPerlinNormalize2(g2[i]);
+    }
+
+    while (--i)
+    {
+        k = p[i];
+        p[i] = p[j = rand() % B];
+        p[j] = k;
+    }
+
+    for (i = 0; i < B + 2; i++)
+    {
+        p[B + i] = p[i];
+        g1[B + i] = g1[i];
+
+        for (j = 0; j < 2; j++)
+            g2[B + i][j] = g2[i][j];
+    }
+}
+
+
+/* Harmonic summing functions - PDB
+ * In what follows "alpha" is the weight when the sum is formed.
+ * Typically it is 2, As this approaches 1 the function is noisier.
+ * "beta" is the harmonic scaling/spacing, typically 2.
+ */
+DMFloat dmPerlinNoise2D(DMFloat x, DMFloat y, DMFloat alpha, DMFloat beta, int n)
+{
+    int i;
+    DMFloat val, sum = 0;
+    DMFloat p[2], scale = 1;
+
+    p[0] = x;
+    p[1] = y;
+    for (i = 0; i < n; i++)
+    {
+        val = dmPerlinDoNoise2(p);
+        sum += val / scale;
+        scale *= alpha;
+        p[0] *= beta;
+        p[1] *= beta;
+    }
+
+    return (sum);
+}