initial commit master
authorNathan Wagner <nw@hydaspes.if.org>
Tue, 25 Sep 2012 04:39:04 +0000 (04:39 +0000)
committerNathan Wagner <nw@hydaspes.if.org>
Tue, 25 Sep 2012 04:39:04 +0000 (04:39 +0000)
Makefile [new file with mode: 0644]
simplex.c [new file with mode: 0644]
simplex.h [new file with mode: 0644]
sp [new file with mode: 0755]

diff --git a/Makefile b/Makefile
new file mode 100644 (file)
index 0000000..51c0203
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,9 @@
+simplex:       simplex.c
+       $(CC) $(CFLAGS) -DMAIN -o $@ $+ -lm
+
+patches: simplex
+       ./simplex 1 | ./sp > sp1.png
+       ./simplex 2 | ./sp > sp2.png
+       ./simplex 4 | ./sp > sp4.png
+       ./simplex 8 | ./sp > sp8.png
+       ./simplex 16 | ./sp > sp16.png
diff --git a/simplex.c b/simplex.c
new file mode 100644 (file)
index 0000000..f492650
--- /dev/null
+++ b/simplex.c
@@ -0,0 +1,811 @@
+#include <math.h>
+
+static double 
+dot3d(int g[], double x, double y, double z)
+{
+       return g[0] * x + g[1] * y + g[2] * z;
+}
+
+static double 
+dot2d(int g[], double x, double y)
+{
+       return g[0] * x + g[1] * y;
+}
+
+static double 
+dot4d(int g[], double x, double y, double z, double w)
+{
+       return g[0] * x + g[1] * y + g[2] * z + g[3] * w;
+}
+
+static double 
+fade(double t)
+{
+       return t * t * t * (t * (t * 6 - 15) + 10);
+}
+
+static double 
+mix(double a, double b, double t)
+{
+       return (1 - t) * a + t * b;
+}
+
+static int      grad3[][3] = {
+       {1, 1, 0}, {-1, 1, 0}, {1, -1, 0}, {-1, -1, 0},
+       {1, 0, 1}, {-1, 0, 1}, {1, 0, -1}, {-1, 0, -1},
+       {0, 1, 1}, {0, -1, 1}, {0, 1, -1}, {0, -1, -1}
+};
+
+static int      grad4[][4] = {{0, 1, 1, 1}, {0, 1, 1, -1}, {0, 1, -1, 1}, {0, 1, -1, -1},
+{0, -1, 1, 1}, {0, -1, 1, -1}, {0, -1, -1, 1}, {0, -1, -1, -1},
+{1, 0, 1, 1}, {1, 0, 1, -1}, {1, 0, -1, 1}, {1, 0, -1, -1},
+{-1, 0, 1, 1}, {-1, 0, 1, -1}, {-1, 0, -1, 1}, {-1, 0, -1, -1},
+{1, 1, 0, 1}, {1, 1, 0, -1}, {1, -1, 0, 1}, {1, -1, 0, -1},
+{-1, 1, 0, 1}, {-1, 1, 0, -1}, {-1, -1, 0, 1}, {-1, -1, 0, -1},
+{1, 1, 1, 0}, {1, 1, -1, 0}, {1, -1, 1, 0}, {1, -1, -1, 0},
+{-1, 1, 1, 0}, {-1, 1, -1, 0}, {-1, -1, 1, 0}, {-1, -1, -1, 0}
+};
+
+/* Classic Perlin noise in 3D, for comparison */
+
+static int      perm[] = {
+       151, 160, 137, 91, 90, 15,
+       131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23,
+       190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
+       88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166,
+       77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244,
+       102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196,
+       135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123,
+       5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42,
+       223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
+       129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228,
+       251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107,
+       49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254,
+       138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180,
+
+       /*
+        * To remove the need for index wrapping, double the permutation
+        * table length
+        */
+
+       151, 160, 137, 91, 90, 15,
+       131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23,
+       190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
+       88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166,
+       77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244,
+       102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196,
+       135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123,
+       5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42,
+       223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
+       129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228,
+       251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107,
+       49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254,
+       138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180
+};
+
+/* Classic Perlin noise, 3D version */
+static double 
+noise3dp(double x, double y, double z)
+{
+       /* Find unit grid cell containing point */
+       int             X = floor(x);
+       int             Y = floor(y);
+       int             Z = floor(z);
+
+       /* Get relative xyz coordinates of point within that cell */
+       x = x - X;
+       y = y - Y;
+       z = z - Z;
+
+       /*
+        * Wrap the integer cells at 255 (smaller integer period can be
+        * introduced here)
+        */
+       X = X & 255;
+       Y = Y & 255;
+       Z = Z & 255;
+
+       /* Calculate a set of eight hashed gradient indices */
+       int             gi000 = perm[X + perm[Y + perm[Z]]] % 12;
+       int             gi001 = perm[X + perm[Y + perm[Z + 1]]] % 12;
+       int             gi010 = perm[X + perm[Y + 1 + perm[Z]]] % 12;
+       int             gi011 = perm[X + perm[Y + 1 + perm[Z + 1]]] % 12;
+       int             gi100 = perm[X + 1 + perm[Y + perm[Z]]] % 12;
+       int             gi101 = perm[X + 1 + perm[Y + perm[Z + 1]]] % 12;
+       int             gi110 = perm[X + 1 + perm[Y + 1 + perm[Z]]] % 12;
+       int             gi111 = perm[X + 1 + perm[Y + 1 + perm[Z + 1]]] % 12;
+
+       /*
+         * The gradients of each corner are now:
+         * g000 = grad3[gi000];
+         * g001 = grad3[gi001];
+         * g010 = grad3[gi010];
+         * g011 = grad3[gi011];
+         * g100 = grad3[gi100];
+         * g101 = grad3[gi101];
+         * g110 = grad3[gi110];
+         * g111 = grad3[gi111];
+         */
+
+       /* Calculate noise contributions from each of the eight corners */
+       double          n000 = dot3d(grad3[gi000], x, y, z);
+       double          n100 = dot3d(grad3[gi100], x - 1, y, z);
+       double          n010 = dot3d(grad3[gi010], x, y - 1, z);
+       double          n110 = dot3d(grad3[gi110], x - 1, y - 1, z);
+       double          n001 = dot3d(grad3[gi001], x, y, z - 1);
+       double          n101 = dot3d(grad3[gi101], x - 1, y, z - 1);
+       double          n011 = dot3d(grad3[gi011], x, y - 1, z - 1);
+       double          n111 = dot3d(grad3[gi111], x - 1, y - 1, z - 1);
+
+       /* Compute the fade curve value for each of x, y, z */
+       double          u = fade(x);
+       double          v = fade(y);
+       double          w = fade(z);
+
+       /* Interpolate along x the contributions from each of the corners */
+       double          nx00 = mix(n000, n100, u);
+       double          nx01 = mix(n001, n101, u);
+       double          nx10 = mix(n010, n110, u);
+       double          nx11 = mix(n011, n111, u);
+
+       /* Interpolate the four results along y */
+       double          nxy0 = mix(nx00, nx10, v);
+       double          nxy1 = mix(nx01, nx11, v);
+
+       /* Interpolate the two last results along z */
+       double          nxyz = mix(nxy0, nxy1, w);
+
+       return nxyz;
+}
+
+/*
+ * A lookup table to traverse the simplex around a given point in 4D.
+ * Details can be found where this table is used, in the 4D noise method.
+ */
+int             simplex[][4] = {
+       {0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 2, 3, 0},
+       {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0},
+       {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+       {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 3, 0, 1}, {2, 3, 1, 0},
+       {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0},
+       {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
+       {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 0, 1, 2}, {3, 0, 2, 1}, {0, 0, 0, 0}, {3, 1, 2, 0},
+{2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}};
+
+double grad1(int hash) {
+       double g;
+       int h = hash & 15;
+       g = 1.0 + (h & 7);   /* Gradient value is one of 1.0, 2.0, ..., 8.0 */
+       if (h&8) g = -g;   /* Make half of the gradients negative */
+       return g;
+}
+
+double noise1d(double x, double *dx)
+{
+       int i0 = floor(x);
+       int i1 = i0 + 1;
+       double x0 = x - i0;
+       double x1 = x0 - 1.0;
+
+       double gx0, gx1;
+       double n0, n1;
+       double t20, t40, t21, t41;
+
+       double x20 = x0*x0;
+       double t0 = 1.0 - x20;
+
+       /*  if(t0 < 0.0f) t0 = 0.0f; /* Never happens for 1D: x0<=1 always */
+       t20 = t0 * t0;
+       t40 = t20 * t20;
+       gx0 = grad1(perm[i0 & 0xff]);
+       n0 = t40 * gx0 * x0;
+
+       double x21 = x1*x1;
+       double t1 = 1.0 - x21;
+       /*  if(t1 < 0.0f) t1 = 0.0f; /* Never happens for 1D: |x1|<=1 always */
+       t21 = t1 * t1;
+       t41 = t21 * t21;
+       gx1 = grad1(perm[i1 & 0xff]);
+       n1 = t41 * gx1 * x1;
+
+       /*
+        * Compute derivative according to:
+        * *dnoise_dx = -8.0f * t20 * t0 * x0 * (gx0 * x0) + t40 * gx0;
+        * *dnoise_dx += -8.0f * t21 * t1 * x1 * (gx1 * x1) + t41 * gx1;
+        */
+
+       if (dx) {
+               *dx = t20 * t0 * gx0 * x20;
+               *dx += t21 * t1 * gx1 * x21;
+               *dx *= -8.0;
+               *dx += t40 * gx0 + t41 * gx1;
+               *dx *= 0.25; /* Scale derivative to match the noise scaling */
+       }
+
+       /*
+        * The maximum value of this noise is 8*(3/4)^4 = 2.53125
+        * A factor of 0.395 would scale to fit exactly within [-1,1], but
+        * to better match classic Perlin noise, we scale it down some more.
+        * scale = 1.0 / 2.53125
+        */
+       return (n0 + n1)/2.53125;
+}
+
+/* 2D simplex noise */
+static double 
+noise2d(double xin, double yin)
+{
+       double          n0, n1, n2;     /* Noise contributions from the three
+                                        * corners */
+
+       /* Skew the input space to determine which simplex cell we're in */
+       double          F2 = 0.5 * (sqrt(3.0) - 1.0);
+       double          s = (xin + yin) * F2;   /* Hairy factor for 2D */
+
+       int             i = floor(xin + s);
+       int             j = floor(yin + s);
+
+       double          G2 = (3.0 - sqrt(3.0)) / 6.0;
+       double          t = (i + j) * G2;
+       double          X0 = i - t;     /* Unskew the cell origin back to
+                                        * (x,y) space */
+       double          Y0 = j - t;
+       double          x0 = xin - X0;  /* The x,y distances from the cell
+                                        * origin */
+       double          y0 = yin - Y0;
+
+       /* For the 2D case, the simplex shape is an equilateral triangle. */
+
+       /* Determine which simplex we are in. */
+       int             i1, j1; /* Offsets for second (middle) corner of
+                                * simplex in (i,j) coords */
+       if (x0 > y0) {
+               /* lower triangle, XY order: (0,0)->(1,0)->(1,1) */
+               i1 = 1;
+               j1 = 0;
+       } else {
+               /* upper triangle, YX order: (0,0)->(0,1)->(1,1) */
+               i1 = 0;
+               j1 = 1;
+       }
+       /*
+         * A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
+         * a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
+         * c = (3-sqrt(3))/6
+         */
+       double          x1 = x0 - i1 + G2;      /* Offsets for middle corner
+                                                * in (x,y) unskewed coords */
+       double          y1 = y0 - j1 + G2;
+       double          x2 = x0 - 1.0 + 2.0 * G2;       /* Offsets for last
+                                                        * corner in (x,y)
+                                                        * unskewed coords */
+       double          y2 = y0 - 1.0 + 2.0 * G2;
+
+       /* Work out the hashed gradient indices of the three simplex corners */
+       int             ii = i & 255;
+       int             jj = j & 255;
+       int             gi0 = perm[ii + perm[jj]] % 12;
+       int             gi1 = perm[ii + i1 + perm[jj + j1]] % 12;
+       int             gi2 = perm[ii + 1 + perm[jj + 1]] % 12;
+
+       /* Calculate the contribution from the three corners */
+       double          t0 = 0.5 - x0 * x0 - y0 * y0;
+       if (t0 < 0) {
+               n0 = 0.0;
+       } else {
+               t0 *= t0;
+               n0 = t0 * t0 * dot2d(grad3[gi0], x0, y0);       /* (x,y) of grad3 used
+                                                                * for 2D gradient */
+       }
+
+       double          t1 = 0.5 - x1 * x1 - y1 * y1;
+       if (t1 < 0) {
+               n1 = 0.0;
+       } else {
+               t1 *= t1;
+               n1 = t1 * t1 * dot2d(grad3[gi1], x1, y1);
+       }
+
+       double          t2 = 0.5 - x2 * x2 - y2 * y2;
+       if (t2 < 0) {
+               n2 = 0.0;
+       } else {
+               t2 *= t2;
+               n2 = t2 * t2 * dot2d(grad3[gi2], x2, y2);
+       }
+
+       /*
+        * Add contributions from each corner to get the final noise value.
+        * The result is scaled to return values in the interval [-1,1].
+        */
+       return 70.0 * (n0 + n1 + n2);
+}
+
+/* 3D simplex noise */
+static double 
+noise3d(double xin, double yin, double zin)
+{
+       double          n0, n1, n2, n3; /* Noise contributions from the four
+                                        * corners */
+       /* Skew the input space to determine which simplex cell we're in */
+       double          F3 = 1.0 / 3.0;
+       double          s = (xin + yin + zin) * F3;     /* Very nice and simple
+                                                        * skew factor for 3D */
+       int             i = floor(xin + s);
+       int             j = floor(yin + s);
+       int             k = floor(zin + s);
+       double          G3 = 1.0 / 6.0; /* Very nice and simple unskew
+                                        * factor, too */
+       double          t = (i + j + k) * G3;
+       double          X0 = i - t;     /* Unskew the cell origin back to
+                                        * (x,y,z) space */
+       double          Y0 = j - t;
+       double          Z0 = k - t;
+       double          x0 = xin - X0;  /* The x,y,z distances from the cell
+                                        * origin */
+       double          y0 = yin - Y0;
+       double          z0 = zin - Z0;
+
+       /*
+        * For the 3D case, the simplex shape is a slightly irregular
+        * tetrahedron.
+        */
+
+       /* Determine which simplex we are in. */
+
+       int             i1, j1, k1;     /* Offsets for second corner of
+                                        * simplex in (i,j,k) coords */
+
+       int             i2, j2, k2;     /* Offsets for third corner of
+                                        * simplex in (i,j,k) coords */
+       if (x0 >= y0) {
+               if (y0 >= z0) {
+                       i1 = 1;
+                       j1 = 0;
+                       k1 = 0;
+                       i2 = 1;
+                       j2 = 1;
+                       k2 = 0;
+               }
+                /* X Y Z order */ 
+               else if (x0 >= z0) {
+                       i1 = 1;
+                       j1 = 0;
+                       k1 = 0;
+                       i2 = 1;
+                       j2 = 0;
+                       k2 = 1;
+               }
+                /* X Z Y order */ 
+               else {
+                       i1 = 0;
+                       j1 = 0;
+                       k1 = 1;
+                       i2 = 1;
+                       j2 = 0;
+                       k2 = 1;
+               }               /* Z X Y order */
+       } else {                /* x0<y0 */
+               if (y0 < z0) {
+                       i1 = 0;
+                       j1 = 0;
+                       k1 = 1;
+                       i2 = 0;
+                       j2 = 1;
+                       k2 = 1;
+               }
+                /* Z Y X order */ 
+               else if (x0 < z0) {
+                       i1 = 0;
+                       j1 = 1;
+                       k1 = 0;
+                       i2 = 0;
+                       j2 = 1;
+                       k2 = 1;
+               }
+                /* Y Z X order */ 
+               else {
+                       i1 = 0;
+                       j1 = 1;
+                       k1 = 0;
+                       i2 = 1;
+                       j2 = 1;
+                       k2 = 0;
+               }               /* Y X Z order */
+       }
+
+       /*
+        * A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in
+        * (x,y,z),
+        */
+
+       /*
+        * a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in
+        * (x,y,z), and
+        */
+
+       /*
+        * a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in
+        * (x,y,z), where
+        */
+
+       /* c = 1/6. */
+       double          x1 = x0 - i1 + G3;      /* Offsets for second corner
+                                                * in (x,y,z) coords */
+       double          y1 = y0 - j1 + G3;
+       double          z1 = z0 - k1 + G3;
+       double          x2 = x0 - i2 + 2.0 * G3;        /* Offsets for third
+                                                        * corner in (x,y,z)
+                                                        * coords */
+       double          y2 = y0 - j2 + 2.0 * G3;
+       double          z2 = z0 - k2 + 2.0 * G3;
+       double          x3 = x0 - 1.0 + 3.0 * G3;       /* Offsets for last
+                                                        * corner in (x,y,z)
+                                                        * coords */
+       double          y3 = y0 - 1.0 + 3.0 * G3;
+       double          z3 = z0 - 1.0 + 3.0 * G3;
+
+       /* Work out the hashed gradient indices of the four simplex corners */
+       int             ii = i & 255;
+       int             jj = j & 255;
+       int             kk = k & 255;
+       int             gi0 = perm[ii + perm[jj + perm[kk]]] % 12;
+       int             gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]] % 12;
+       int             gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]] % 12;
+       int             gi3 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]] % 12;
+
+       /* Calculate the contribution from the four corners */
+       double          t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
+       if (t0 < 0)
+               n0 = 0.0;
+       else {
+               t0 *= t0;
+               n0 = t0 * t0 * dot3d(grad3[gi0], x0, y0, z0);
+       }
+
+       double          t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
+       if (t1 < 0) {
+               n1 = 0.0;
+       }
+       else {
+               t1 *= t1;
+               n1 = t1 * t1 * dot3d(grad3[gi1], x1, y1, z1);
+       }
+
+       double          t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
+       if (t2 < 0) {
+               n2 = 0.0;
+       }
+       else {
+               t2 *= t2;
+               n2 = t2 * t2 * dot3d(grad3[gi2], x2, y2, z2);
+       }
+
+       double          t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
+       if (t3 < 0) {
+               n3 = 0.0;
+       }
+       else {
+               t3 *= t3;
+               n3 = t3 * t3 * dot3d(grad3[gi3], x3, y3, z3);
+       }
+       /* Add contributions from each corner to get the final noise value. */
+       /* The result is scaled to stay just inside [-1,1] */
+       return 32.0 * (n0 + n1 + n2 + n3);
+}
+
+/* 4D simplex noise */
+double 
+noise4d(double x, double y, double z, double w)
+{
+
+       /* The skewing and unskewing factors are hairy again for the 4D case */
+       double          F4 = (sqrt(5.0) - 1.0) / 4.0;
+       double          G4 = (5.0 - sqrt(5.0)) / 20.0;
+       double          n0, n1, n2, n3, n4;     /* Noise contributions from
+                                                * the five corners */
+       /*
+        * Skew the (x,y,z,w) space to determine which cell of 24 simplices
+        * we're in
+        */
+       double          s = (x + y + z + w) * F4;       /* Factor for 4D skewing */
+       int             i = floor(x + s);
+       int             j = floor(y + s);
+       int             k = floor(z + s);
+       int             l = floor(w + s);
+       double          t = (i + j + k + l) * G4;       /* Factor for 4D
+                                                        * unskewing */
+       double          X0 = i - t;     /* Unskew the cell origin back to
+                                        * (x,y,z,w) space */
+       double          Y0 = j - t;
+       double          Z0 = k - t;
+       double          W0 = l - t;
+       double          x0 = x - X0;    /* The x,y,z,w distances from the
+                                        * cell origin */
+       double          y0 = y - Y0;
+       double          z0 = z - Z0;
+       double          w0 = w - W0;
+       /*
+        * For the 4D case, the simplex is a 4D shape I won't even try to
+        * describe.
+        */
+       /*
+        * To find out which of the 24 possible simplices we're in, we need
+        * to
+        */
+       /* determine the magnitude ordering of x0, y0, z0 and w0. */
+       /*
+        * The method below is a good way of finding the ordering of x,y,z,w
+        * and
+        */
+       /* then find the correct traversal order for the simplex we’re in. */
+       /*
+        * First, six pair-wise comparisons are performed between each
+        * possible pair
+        */
+       /*
+        * of the four coordinates, and the results are used to add up binary
+        * bits
+        */
+       /* for an integer index. */
+       int             c1 = (x0 > y0) ? 32 : 0;
+       int             c2 = (x0 > z0) ? 16 : 0;
+       int             c3 = (y0 > z0) ? 8 : 0;
+       int             c4 = (x0 > w0) ? 4 : 0;
+       int             c5 = (y0 > w0) ? 2 : 0;
+       int             c6 = (z0 > w0) ? 1 : 0;
+       int             c = c1 + c2 + c3 + c4 + c5 + c6;
+       int             i1, j1, k1, l1; /* The integer offsets for the second
+                                        * simplex corner */
+       int             i2, j2, k2, l2; /* The integer offsets for the third
+                                        * simplex corner */
+       int             i3, j3, k3, l3; /* The integer offsets for the fourth
+                                        * simplex corner */
+       /*
+        * simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some
+        * order.
+        */
+       /*
+        * Many values of c will never occur, since e.g. x>y>z>w makes x<z,
+        * y<w and x<w
+        */
+       /*
+        * impossible. Only the 24 indices which have non-zero entries make
+        * any sense.
+        */
+       /*
+        * We use a thresholding to set the coordinates in turn from the
+        * largest magnitude.
+        */
+       /*
+        * The number 3 in the "simplex" array is at the position of the
+        * largest coordinate.
+        */
+       i1 = simplex[c][0] >= 3 ? 1 : 0;
+       j1 = simplex[c][1] >= 3 ? 1 : 0;
+       k1 = simplex[c][2] >= 3 ? 1 : 0;
+       l1 = simplex[c][3] >= 3 ? 1 : 0;
+       /*
+        * The number 2 in the "simplex" array is at the second largest
+        * coordinate.
+        */
+       i2 = simplex[c][0] >= 2 ? 1 : 0;
+       j2 = simplex[c][1] >= 2 ? 1 : 0;
+       k2 = simplex[c][2] >= 2 ? 1 : 0;
+       l2 = simplex[c][3] >= 2 ? 1 : 0;
+       /*
+        * The number 1 in the "simplex" array is at the second smallest
+        * coordinate.
+        */
+       i3 = simplex[c][0] >= 1 ? 1 : 0;
+       j3 = simplex[c][1] >= 1 ? 1 : 0;
+       k3 = simplex[c][2] >= 1 ? 1 : 0;
+       l3 = simplex[c][3] >= 1 ? 1 : 0;
+       /*
+        * The fifth corner has all coordinate offsets = 1, so no need to
+        * look that up.
+        */
+       double          x1 = x0 - i1 + G4;      /* Offsets for second corner
+                                                * in (x,y,z,w) coords */
+       double          y1 = y0 - j1 + G4;
+       double          z1 = z0 - k1 + G4;
+       double          w1 = w0 - l1 + G4;
+       double          x2 = x0 - i2 + 2.0 * G4;        /* Offsets for third
+                                                        * corner in (x,y,z,w)
+                                                        * coords */
+       double          y2 = y0 - j2 + 2.0 * G4;
+       double          z2 = z0 - k2 + 2.0 * G4;
+       double          w2 = w0 - l2 + 2.0 * G4;
+       double          x3 = x0 - i3 + 3.0 * G4;        /* Offsets for fourth
+                                                        * corner in (x,y,z,w)
+                                                        * coords */
+       double          y3 = y0 - j3 + 3.0 * G4;
+       double          z3 = z0 - k3 + 3.0 * G4;
+       double          w3 = w0 - l3 + 3.0 * G4;
+       double          x4 = x0 - 1.0 + 4.0 * G4;       /* Offsets for last
+                                                        * corner in (x,y,z,w)
+                                                        * coords */
+       double          y4 = y0 - 1.0 + 4.0 * G4;
+       double          z4 = z0 - 1.0 + 4.0 * G4;
+       double          w4 = w0 - 1.0 + 4.0 * G4;
+       /* Work out the hashed gradient indices of the five simplex corners */
+       int             ii = i & 255;
+       int             jj = j & 255;
+       int             kk = k & 255;
+       int             ll = l & 255;
+       int             gi0 = perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32;
+       int             gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32;
+       int             gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32;
+       int             gi3 = perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32;
+       int             gi4 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32;
+       /* Calculate the contribution from the five corners */
+       double          t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
+       if (t0 < 0) {
+               n0 = 0.0;
+       }
+       else {
+               t0 *= t0;
+               n0 = t0 * t0 * dot4d(grad4[gi0], x0, y0, z0, w0);
+       }
+       double          t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
+       if (t1 < 0) {
+               n1 = 0.0;
+       }
+       else {
+               t1 *= t1;
+               n1 = t1 * t1 * dot4d(grad4[gi1], x1, y1, z1, w1);
+       }
+       double          t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
+       if (t2 < 0) {
+               n2 = 0.0;
+       }
+       else {
+               t2 *= t2;
+               n2 = t2 * t2 * dot4d(grad4[gi2], x2, y2, z2, w2);
+       }
+       double          t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
+       if (t3 < 0) {
+               n3 = 0.0;
+       }
+       else {
+               t3 *= t3;
+               n3 = t3 * t3 * dot4d(grad4[gi3], x3, y3, z3, w3);
+       }
+       double          t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
+       if (t4 < 0) {
+               n4 = 0.0;
+       }
+       else {
+               t4 *= t4;
+               n4 = t4 * t4 * dot4d(grad4[gi4], x4, y4, z4, w4);
+       }
+
+       /* Sum up and scale the result to cover the range [-1,1] */
+       return 27.0 * (n0 + n1 + n2 + n3 + n4);
+}
+
+double noise_simplex(int dim, double scale, double *pt) {
+       switch (dim) {
+               case 1:
+                       return noise1d(pt[0]*scale,0);
+                       break;
+               case 2:
+                       return noise2d(pt[0]*scale,pt[1]*scale);
+                       break;
+               case 3:
+                       return noise3d(pt[0]*scale,pt[1]*scale,pt[2]*scale);
+                       break;
+               case 4:
+                       return noise4d(pt[0]*scale,pt[1]*scale,pt[2]*scale,pt[3]*scale);
+                       break;
+               default:
+                       break;
+       }
+       return 0.0;
+}
+
+double noise_fractal(int dim, double *pt, int octaves, double lacunarity, double persistence) {
+       double noise,scale,lac;
+       noise = 0.0;
+       scale = 1.0;
+       lac = 1.0;
+
+       for (scale=1.0; octaves--;scale*=persistence) {
+               noise += noise_simplex(dim, lac, pt) * scale;
+               scale *= persistence;
+               lac *= lacunarity;
+       }
+       return noise;
+}
+
+#ifdef MAIN
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+
+int 
+main(int argc, char *argv[])
+{
+       double  r;
+       double  base;
+       double  s[2];
+
+       int x,y;
+
+       int octaves;
+       int pos = 0, neg = 0;
+       double max, min;
+#ifdef INFINITY
+       max = -INFINITY;
+       min = +INFINITY;
+#else
+       max = DBL_MIN;
+       min = DBL_MAX;
+#endif
+
+       octaves = argc > 1 ? atoi(argv[1]) : 1;
+       base = 40.0;
+
+       for (x=0;x<256;x++) {
+               s[0] = x/base;
+               for (y=0;y<256;y++) {
+                       s[1] = y/base;
+                       r = noise_fractal(2, s, octaves, 2.0, 0.7);
+                       r /= 1.5;
+
+                       /* grain */
+                       r = noise_simplex(2,0.20,s)*10;
+                       r = (r - (int)r);
+                       /* fine */
+                       r += noise_simplex(2, 75.0, s) * 0.20;
+#if 1
+                       /* streak */
+                       r += noise2d(s[0]*0.45, s[1]*200.0) * 0.3;
+#endif
+
+                       r /= 2.5;
+                       printf("%f\n", r);
+                       if (r > 1.0) {
+                               pos++;
+                       } else if (r < -1.0) {
+                               neg++;
+                       }
+                       if (r > max) max = r;
+                       if (r < min) min = r;
+
+               }
+       }
+       fprintf(stderr, "max value: %f\n", max);
+       fprintf(stderr, "min value: %f\n", min);
+       if (pos) fprintf(stderr, "%d values > +1.0\n", pos);
+       if (neg) fprintf(stderr, "%d values < -1.0\n", neg);
+       
+       return 0;
+
+       switch (--argc) {
+               case 1:
+               r = noise1d(atof(argv[1]),0);
+               break;
+
+       case 2:
+               r = noise2d(atof(argv[1]), atof(argv[2]));
+               break;
+       case 3:
+               r = noise3d(atof(argv[1]), atof(argv[2]), atof(argv[3]));
+               break;
+       case 4:
+               r = noise4d(atof(argv[1]), atof(argv[2]), atof(argv[3]), atof(argv[4]));
+               break;
+       default:
+               fprintf(stderr, "Wrong number of arguments: %d, must be 2, 3, or 4\n", argc);
+               exit(EXIT_FAILURE);
+               break;
+       };
+
+       printf("%f\n", r);
+
+       return 0;
+}
+
+#endif
diff --git a/simplex.h b/simplex.h
new file mode 100644 (file)
index 0000000..f5ba1e1
--- /dev/null
+++ b/simplex.h
@@ -0,0 +1,40 @@
+#ifndef LIBNOISE_H_
+#define LIBNOISE_H_
+
+enum {
+       NOISE_PERLIN, NOISE_SIMPLEX, NOISE_WHITE, NOISE_BLUE
+};
+
+struct noise_generator {
+       int type;
+       int seed;
+       int perm[512];
+};
+
+double noise_simplex(int dim, double *vec);
+double noise_perlin(int dim, double *vec);
+double noise_gen(struct noise_generator *gen, int dim, double *vec);
+int noise_gen_vec(struct noise_generator *gen, int dim, double *vec, int npoints, double *result);
+double noise_fbm(struct noise_generator *gen, int dim, double *vec, double persistence);
+
+struct noise_generator *noise_gen_init(struct noise_generator *g, int type, int dim, int seed);
+
+/*
+ * takes a start and end point and generates a noisy line between them
+ *
+ * line is constrained to be in a quadrilateral defined by the two end points
+ * and two points equidistant from ... hmm, not sure how to specify dimensions
+ * less than two don't make any sense.
+ *
+ * linevec is filled in with computed points
+ */
+
+int noise_line(struct noise_generator *gen, int dim, double *start, double *end, double width, int segments, double *linevec);
+
+/*
+ * noise_polygon?
+ * noise_patch?
+ * noise_mesh?
+ */
+
+#endif
diff --git a/sp b/sp
new file mode 100755 (executable)
index 0000000..e0c30ae
--- /dev/null
+++ b/sp
@@ -0,0 +1,24 @@
+#!/usr/bin/perl
+
+use warnings;
+use strict;
+
+use GD;
+
+GD::Image->trueColor(1);
+
+my $png = GD::Image->new(256,256);
+
+for my $x (0..255) {
+       for my $y (0..255) {
+               my $s = <>;
+               chomp($s);
+               my $dx = $x/255;
+               my $dy = $y/255;
+               $s = ($s + 1.0) / 2.0;
+               $png->setPixel($x, $y, $png->colorResolve($s *180, $s*130, $s*80));
+       }
+}
+
+binmode STDOUT;
+print $png->png;