operations: add a cell noise operation

Added a cell noise operation implementing the texture function described in

    Steven Worley. 1996. A cellular texture basis function.
    In Proceedings of the 23rd annual conference on Computer
    graphics and interactive techniques (SIGGRAPH '96).

2 files changed, 390 insertions, 0 deletions

diff --git a/opencl/noise-cell.cl b/opencl/noise-cell.cl
new file mode 100644
index 00000000..1e555853
--- /dev/null
+++ b/opencl/noise-cell.cl
@@ -0,0 +1,194 @@
+#define MAX_RANK 3
+
+/* Random feature counts following the Poisson distribution with
+   lambda equal to 7. */
+
+static const __constant char poisson[256] = {
+  7, 9, 12, 12, 8, 7, 5, 5, 6, 7, 8, 6, 10, 7, 6, 2, 8, 3, 9, 5, 13, 10, 9,
+  8, 8, 9, 3, 8, 9, 6, 8, 7, 4, 9, 6, 3, 10, 7, 7, 7, 6, 7, 4, 14, 7, 6, 11,
+  7, 7, 7, 12, 7, 10, 6, 8, 11, 3, 5, 7, 7, 8, 7, 9, 8, 5, 8, 11, 3, 4, 5, 8,
+  8, 7, 8, 9, 2, 7, 8, 12, 4, 8, 2, 11, 8, 14, 7, 8, 2, 3, 10, 4, 6, 9, 5, 8,
+  7, 10, 10, 10, 14, 5, 7, 6, 4, 5, 6, 11, 8, 7, 3, 11, 5, 5, 2, 9, 7, 7, 7,
+  9, 2, 7, 6, 9, 7, 6, 5, 12, 5, 3, 11, 9, 12, 8, 6, 8, 6, 8, 5, 5, 7, 5, 2,
+  9, 5, 5, 8, 11, 8, 8, 10, 6, 4, 7, 14, 7, 3, 10, 7, 7, 4, 9, 10, 10, 9, 8,
+  8, 7, 6, 5, 10, 10, 5, 10, 7, 7, 10, 7, 4, 9, 9, 6, 8, 5, 10, 7, 3, 9, 9,
+  7, 8, 9, 7, 5, 7, 6, 5, 5, 12, 4, 7, 5, 5, 4, 5, 7, 10, 8, 7, 9, 4, 6, 11,
+  6, 3, 7, 8, 9, 5, 8, 6, 7, 8, 7, 7, 3, 7, 7, 9, 4, 5, 5, 6, 9, 7, 6, 12, 4,
+  9, 10, 8, 8, 6, 4, 9, 9, 8, 11, 6, 8, 13, 8, 9, 12, 6, 9, 8
+};
+
+static uint
+philox (uint s,
+        uint t,
+        uint k)
+{
+  ulong p;
+  int   i;
+
+  for (i = 0 ; i < 3 ; i += 1)
+    {
+      p = s * 0xcd9e8d57ul;
+
+      s = ((uint)(p >> 32)) ^ t ^ k;
+      t = (uint)p;
+
+      k += 0x9e3779b9u;
+    }
+
+  return s;
+}
+
+static float
+lcg (uint *hash)
+{
+  return (*hash = *hash * 1664525u + 1013904223u) / 4294967296.0f;
+}
+
+static void
+search_box (float              *closest,
+            uint               *feature,
+            int                 s,
+            int                 t,
+            float               x,
+            float               y,
+            float               shape,
+            uint                rank,
+            uint                seed)
+{
+  uint hash;
+  int  i, n;
+
+  hash = philox ((uint)s, (uint)t, seed);
+  n = poisson[hash >> 24];
+
+  for (i = 0 ; i < n ; i += 1)
+    {
+      float delta_x, delta_y, d;
+      int   j, k;
+
+      /* Calculate the distance to each feature point. */
+
+      delta_x = s + lcg (&hash) - x;
+      delta_y = t + lcg (&hash) - y;
+
+      if (shape == 2)
+        d = delta_x * delta_x + delta_y * delta_y;
+      else if (shape == 1)
+        d = fabs (delta_x) + fabs (delta_y);
+      else
+        d = pow (fabs (delta_x), shape) +
+          pow (fabs (delta_y), shape);
+
+      /* Insert it into the list of n closest distances if needed. */
+
+      for (j = 0 ; j < rank && d > closest[j] ; j += 1);
+
+      if (j < rank)
+        {
+          for (k = rank - 1 ; k > j ; k -= 1)
+            {
+              closest[k] = closest[k - 1];
+            }
+
+          closest[j] = d;
+
+          if (j == rank - 1)
+            *feature = hash;
+        }
+    }
+}
+
+__kernel void kernel_noise (__global float *out,
+                            const int       x_0,
+                            const int       y_0,
+                            const uint      iterations,
+                            float           scale,
+                            float           shape,
+                            uint            rank,
+                            uint            seed,
+                            int             palettize)
+{
+  const int gidx = get_global_id(0);
+  const int gidy = get_global_id(1);
+
+  float c, d, n, closest[MAX_RANK], *d_0 = &closest[MAX_RANK - 1];
+  uint  feature;
+  int   i, j;
+
+  for (j = 0, n = 0, c = 1, d = scale;
+       j < iterations;
+       c *= 2, d *= 2, j += 1)
+    {
+      float d_l, d_r, d_t, d_b;
+      float x = (float)(gidx + x_0) * d;
+      float y = (float)(gidy + y_0) * d;
+      int   s = (int)floor(x);
+      int   t = (int)floor(y);
+
+      for (i = 0 ; i < rank ; closest[i] = 1.0 / 0.0, i += 1);
+
+      /* Search the box the point is in. */
+
+      search_box (closest, &feature, s, t, x, y, shape, rank, seed);
+
+      d_0 = &closest[rank - 1];
+      d_l = x - s; d_l *= d_l;
+      d_r = 1.0 - x + s; d_r *= d_r;
+      d_b = y - t; d_b *= d_b;
+      d_t = 1.0 - y + t; d_t *= d_t;
+
+      /* Search adjacent boxes if it is possible for them to contain a
+       * nearby feature point. */
+
+      if (d_l < *d_0)
+        {
+          if (d_l + d_b < *d_0)
+            search_box (closest, &feature, s - 1, t - 1, x, y,
+                        shape, rank, seed);
+
+          search_box (closest, &feature, s - 1, t, x, y,
+                      shape, rank, seed);
+
+          if (d_l + d_t < *d_0)
+              search_box (closest, &feature, s - 1, t + 1, x, y,
+                          shape, rank, seed);
+        }
+
+      if (d_b < *d_0)
+          search_box (closest, &feature, s, t - 1, x, y,
+                      shape, rank, seed);
+
+      if (d_t < *d_0)
+          search_box (closest, &feature, s, t + 1, x, y,
+                      shape, rank, seed);
+
+      if (d_r < *d_0)
+        {
+          if (d_r + d_b < *d_0)
+              search_box (closest, &feature, s + 1, t - 1, x, y,
+                          shape, rank, seed);
+
+          search_box (closest, &feature, s + 1, t, x, y,
+                      shape, rank, seed);
+
+          if (d_r + d_t < *d_0)
+              search_box (closest, &feature, s + 1, t + 1, x, y,
+                          shape, rank, seed);
+        }
+
+      /* If palettized output is requested return the normalized hash of
+       * the closest feature point, otherwise return the closest
+       * distance. */
+
+      if(palettize)
+        {
+          n += feature / 4294967295.0f / c;
+        }
+      else
+        {
+          n += pow(closest[rank - 1], 1 / shape) / c;
+        }
+    }
+
+  out[gidy * get_global_size(0) + gidx] = n;
+}
diff --git a/opencl/noise-cell.cl.h b/opencl/noise-cell.cl.h
new file mode 100644
index 00000000..775f5911
--- /dev/null
+++ b/opencl/noise-cell.cl.h
@@ -0,0 +1,196 @@
+static const char* noise_cell_cl_source =
+"#define MAX_RANK 3                                                            \n"
+"                                                                              \n"
+"/* Random feature counts following the Poisson distribution with              \n"
+"   lambda equal to 7. */                                                      \n"
+"                                                                              \n"
+"static const __constant char poisson[256] = {                                 \n"
+"  7, 9, 12, 12, 8, 7, 5, 5, 6, 7, 8, 6, 10, 7, 6, 2, 8, 3, 9, 5, 13, 10, 9,   \n"
+"  8, 8, 9, 3, 8, 9, 6, 8, 7, 4, 9, 6, 3, 10, 7, 7, 7, 6, 7, 4, 14, 7, 6, 11,  \n"
+"  7, 7, 7, 12, 7, 10, 6, 8, 11, 3, 5, 7, 7, 8, 7, 9, 8, 5, 8, 11, 3, 4, 5, 8, \n"
+"  8, 7, 8, 9, 2, 7, 8, 12, 4, 8, 2, 11, 8, 14, 7, 8, 2, 3, 10, 4, 6, 9, 5, 8, \n"
+"  7, 10, 10, 10, 14, 5, 7, 6, 4, 5, 6, 11, 8, 7, 3, 11, 5, 5, 2, 9, 7, 7, 7,  \n"
+"  9, 2, 7, 6, 9, 7, 6, 5, 12, 5, 3, 11, 9, 12, 8, 6, 8, 6, 8, 5, 5, 7, 5, 2,  \n"
+"  9, 5, 5, 8, 11, 8, 8, 10, 6, 4, 7, 14, 7, 3, 10, 7, 7, 4, 9, 10, 10, 9, 8,  \n"
+"  8, 7, 6, 5, 10, 10, 5, 10, 7, 7, 10, 7, 4, 9, 9, 6, 8, 5, 10, 7, 3, 9, 9,   \n"
+"  7, 8, 9, 7, 5, 7, 6, 5, 5, 12, 4, 7, 5, 5, 4, 5, 7, 10, 8, 7, 9, 4, 6, 11,  \n"
+"  6, 3, 7, 8, 9, 5, 8, 6, 7, 8, 7, 7, 3, 7, 7, 9, 4, 5, 5, 6, 9, 7, 6, 12, 4, \n"
+"  9, 10, 8, 8, 6, 4, 9, 9, 8, 11, 6, 8, 13, 8, 9, 12, 6, 9, 8                 \n"
+"};                                                                            \n"
+"                                                                              \n"
+"static uint                                                                   \n"
+"philox (uint s,                                                               \n"
+"        uint t,                                                               \n"
+"        uint k)                                                               \n"
+"{                                                                             \n"
+"  ulong p;                                                                    \n"
+"  int   i;                                                                    \n"
+"                                                                              \n"
+"  for (i = 0 ; i < 3 ; i += 1)                                                \n"
+"    {                                                                         \n"
+"      p = s * 0xcd9e8d57ul;                                                   \n"
+"                                                                              \n"
+"      s = ((uint)(p >> 32)) ^ t ^ k;                                          \n"
+"      t = (uint)p;                                                            \n"
+"                                                                              \n"
+"      k += 0x9e3779b9u;                                                       \n"
+"    }                                                                         \n"
+"                                                                              \n"
+"  return s;                                                                   \n"
+"}                                                                             \n"
+"                                                                              \n"
+"static float                                                                  \n"
+"lcg (uint *hash)                                                              \n"
+"{                                                                             \n"
+"  return (*hash = *hash * 1664525u + 1013904223u) / 4294967296.0f;            \n"
+"}                                                                             \n"
+"                                                                              \n"
+"static void                                                                   \n"
+"search_box (float              *closest,                                      \n"
+"            uint               *feature,                                      \n"
+"            int                 s,                                            \n"
+"            int                 t,                                            \n"
+"            float               x,                                            \n"
+"            float               y,                                            \n"
+"            float               shape,                                        \n"
+"            uint                rank,                                         \n"
+"            uint                seed)                                         \n"
+"{                                                                             \n"
+"  uint hash;                                                                  \n"
+"  int  i, n;                                                                  \n"
+"                                                                              \n"
+"  hash = philox ((uint)s, (uint)t, seed);                                     \n"
+"  n = poisson[hash >> 24];                                                    \n"
+"                                                                              \n"
+"  for (i = 0 ; i < n ; i += 1)                                                \n"
+"    {                                                                         \n"
+"      float delta_x, delta_y, d;                                              \n"
+"      int   j, k;                                                             \n"
+"                                                                              \n"
+"      /* Calculate the distance to each feature point. */                     \n"
+"                                                                              \n"
+"      delta_x = s + lcg (&hash) - x;                                          \n"
+"      delta_y = t + lcg (&hash) - y;                                          \n"
+"                                                                              \n"
+"      if (shape == 2)                                                         \n"
+"        d = delta_x * delta_x + delta_y * delta_y;                            \n"
+"      else if (shape == 1)                                                    \n"
+"        d = fabs (delta_x) + fabs (delta_y);                                  \n"
+"      else                                                                    \n"
+"        d = pow (fabs (delta_x), shape) +                                     \n"
+"          pow (fabs (delta_y), shape);                                        \n"
+"                                                                              \n"
+"      /* Insert it into the list of n closest distances if needed. */         \n"
+"                                                                              \n"
+"      for (j = 0 ; j < rank && d > closest[j] ; j += 1);                      \n"
+"                                                                              \n"
+"      if (j < rank)                                                           \n"
+"        {                                                                     \n"
+"          for (k = rank - 1 ; k > j ; k -= 1)                                 \n"
+"            {                                                                 \n"
+"              closest[k] = closest[k - 1];                                    \n"
+"            }                                                                 \n"
+"                                                                              \n"
+"          closest[j] = d;                                                     \n"
+"                                                                              \n"
+"          if (j == rank - 1)                                                  \n"
+"            *feature = hash;                                                  \n"
+"        }                                                                     \n"
+"    }                                                                         \n"
+"}                                                                             \n"
+"                                                                              \n"
+"__kernel void kernel_noise (__global float *out,                              \n"
+"                            const int       x_0,                              \n"
+"                            const int       y_0,                              \n"
+"                            const uint      iterations,                       \n"
+"                            float           scale,                            \n"
+"                            float           shape,                            \n"
+"                            uint            rank,                             \n"
+"                            uint            seed,                             \n"
+"                            int             palettize)                        \n"
+"{                                                                             \n"
+"  const int gidx = get_global_id(0);                                          \n"
+"  const int gidy = get_global_id(1);                                          \n"
+"                                                                              \n"
+"  float c, d, n, closest[MAX_RANK], *d_0 = &closest[MAX_RANK - 1];            \n"
+"  uint  feature;                                                              \n"
+"  int   i, j;                                                                 \n"
+"                                                                              \n"
+"  for (j = 0, n = 0, c = 1, d = scale;                                        \n"
+"       j < iterations;                                                        \n"
+"       c *= 2, d *= 2, j += 1)                                                \n"
+"    {                                                                         \n"
+"      float d_l, d_r, d_t, d_b;                                               \n"
+"      float x = (float)(gidx + x_0) * d;                                      \n"
+"      float y = (float)(gidy + y_0) * d;                                      \n"
+"      int   s = (int)floor(x);                                                \n"
+"      int   t = (int)floor(y);                                                \n"
+"                                                                              \n"
+"      for (i = 0 ; i < rank ; closest[i] = 1.0 / 0.0, i += 1);                \n"
+"                                                                              \n"
+"      /* Search the box the point is in. */                                   \n"
+"                                                                              \n"
+"      search_box (closest, &feature, s, t, x, y, shape, rank, seed);          \n"
+"                                                                              \n"
+"      d_0 = &closest[rank - 1];                                               \n"
+"      d_l = x - s; d_l *= d_l;                                                \n"
+"      d_r = 1.0 - x + s; d_r *= d_r;                                          \n"
+"      d_b = y - t; d_b *= d_b;                                                \n"
+"      d_t = 1.0 - y + t; d_t *= d_t;                                          \n"
+"                                                                              \n"
+"      /* Search adjacent boxes if it is possible for them to contain a        \n"
+"       * nearby feature point. */                                             \n"
+"                                                                              \n"
+"      if (d_l < *d_0)                                                         \n"
+"        {                                                                     \n"
+"          if (d_l + d_b < *d_0)                                               \n"
+"            search_box (closest, &feature, s - 1, t - 1, x, y,                \n"
+"                        shape, rank, seed);                                   \n"
+"                                                                              \n"
+"          search_box (closest, &feature, s - 1, t, x, y,                      \n"
+"                      shape, rank, seed);                                     \n"
+"                                                                              \n"
+"          if (d_l + d_t < *d_0)                                               \n"
+"              search_box (closest, &feature, s - 1, t + 1, x, y,              \n"
+"                          shape, rank, seed);                                 \n"
+"        }                                                                     \n"
+"                                                                              \n"
+"      if (d_b < *d_0)                                                         \n"
+"          search_box (closest, &feature, s, t - 1, x, y,                      \n"
+"                      shape, rank, seed);                                     \n"
+"                                                                              \n"
+"      if (d_t < *d_0)                                                         \n"
+"          search_box (closest, &feature, s, t + 1, x, y,                      \n"
+"                      shape, rank, seed);                                     \n"
+"                                                                              \n"
+"      if (d_r < *d_0)                                                         \n"
+"        {                                                                     \n"
+"          if (d_r + d_b < *d_0)                                               \n"
+"              search_box (closest, &feature, s + 1, t - 1, x, y,              \n"
+"                          shape, rank, seed);                                 \n"
+"                                                                              \n"
+"          search_box (closest, &feature, s + 1, t, x, y,                      \n"
+"                      shape, rank, seed);                                     \n"
+"                                                                              \n"
+"          if (d_r + d_t < *d_0)                                               \n"
+"              search_box (closest, &feature, s + 1, t + 1, x, y,              \n"
+"                          shape, rank, seed);                                 \n"
+"        }                                                                     \n"
+"                                                                              \n"
+"      /* If palettized output is requested return the normalized hash of      \n"
+"       * the closest feature point, otherwise return the closest              \n"
+"       * distance. */                                                         \n"
+"                                                                              \n"
+"      if(palettize)                                                           \n"
+"        {                                                                     \n"
+"          n += feature / 4294967295.0f / c;                                   \n"
+"        }                                                                     \n"
+"      else                                                                    \n"
+"        {                                                                     \n"
+"          n += pow(closest[rank - 1], 1 / shape) / c;                         \n"
+"        }                                                                     \n"
+"    }                                                                         \n"
+"                                                                              \n"
+"  out[gidy * get_global_size(0) + gidx] = n;                                  \n"
+"}                                                                             \n"
+;