Adding upstream version 2.10.36.upstream/2.10.36

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 398 insertions, 0 deletions

diff --git a/plug-ins/flame/rect.c b/plug-ins/flame/rect.c
new file mode 100644
index 0000000..0ff0f5f
--- /dev/null
+++ b/plug-ins/flame/rect.c
@@ -0,0 +1,398 @@
+/*
+   flame - cosmic recursive fractal flames
+   Copyright (C) 1992  Scott Draves <spot@cs.cmu.edu>
+
+   This program 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.
+
+   This program 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 this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include "rect.h"
+
+#include <string.h>
+
+#include "libgimp/gimp.h"
+
+/* for batch
+ *   interpolate
+ *   compute colormap
+ *   for subbatch
+ *     compute samples
+ *     buckets += cmap[samples]
+ *   accum += time_filter[batch] * log(buckets)
+ * image = filter(accum)
+ */
+
+
+typedef short bucket[4];
+
+/* if you use longs instead of shorts, you
+   get higher quality, and spend more memory */
+
+#if 1
+typedef short accum_t;
+#define MAXBUCKET (1<<14)
+#define SUB_BATCH_SIZE 10000
+#else
+typedef long accum_t;
+#define MAXBUCKET (1<<30)
+#define SUB_BATCH_SIZE 10000
+#endif
+
+typedef accum_t abucket[4];
+
+
+
+/* allow this many iterations for settling into attractor */
+#define FUSE 15
+
+/* clamp spatial filter to zero at this std dev (2.5 ~= 0.0125) */
+#define FILTER_CUTOFF 2.5
+
+/* should be MAXBUCKET / (OVERSAMPLE^2) */
+#define PREFILTER_WHITE (MAXBUCKET>>4)
+
+
+#define bump_no_overflow(dest, delta, type) { \
+   type tt_ = dest + delta;            \
+   if (tt_ > dest) dest = tt_;                 \
+}
+
+/* sum of entries of vector to 1 */
+static void
+normalize_vector(double *v,
+                 int     n)
+{
+  double t = 0.0;
+  int i;
+  for (i = 0; i < n; i++)
+    t += v[i];
+  t = 1.0 / t;
+  for (i = 0; i < n; i++)
+    v[i] *= t;
+}
+
+void
+render_rectangle (frame_spec    *spec,
+                  unsigned char *out,
+                  int            out_width,
+                  int            field,
+                  int            nchan,
+                  int progress(double))
+{
+  int      i, j, k, nsamples, nbuckets, batch_size, batch_num, sub_batch;
+  bucket  *buckets;
+  abucket *accumulate;
+  point   *points;
+  double  *filter, *temporal_filter, *temporal_deltas;
+  double   bounds[4], size[2], ppux, ppuy;
+  int      image_width, image_height;    /* size of the image to produce */
+  int      width, height;               /* size of histogram */
+  int      filter_width;
+  int      oversample = spec->cps[0].spatial_oversample;
+  int      nbatches = spec->cps[0].nbatches;
+  bucket   cmap[CMAP_SIZE];
+  int      gutter_width;
+  int      sbc;
+
+  image_width = spec->cps[0].width;
+  if (field)
+    {
+      image_height = spec->cps[0].height / 2;
+      if (field == field_odd)
+        out += nchan * out_width;
+      out_width *= 2;
+    }
+  else
+    image_height = spec->cps[0].height;
+
+  if (1)
+    {
+      filter_width = (2.0 * FILTER_CUTOFF * oversample *
+                      spec->cps[0].spatial_filter_radius);
+      /* make sure it has same parity as oversample */
+      if ((filter_width ^ oversample) & 1)
+        filter_width++;
+
+      filter = g_malloc (sizeof (double) * filter_width * filter_width);
+      /* fill in the coefs */
+      for (i = 0; i < filter_width; i++)
+        for (j = 0; j < filter_width; j++)
+          {
+            double ii = ((2.0 * i + 1.0) / filter_width - 1.0) * FILTER_CUTOFF;
+            double jj = ((2.0 * j + 1.0) / filter_width - 1.0) * FILTER_CUTOFF;
+            if (field)
+              jj *= 2.0;
+            filter[i + j * filter_width] = exp(-2.0 * (ii * ii + jj * jj));
+          }
+      normalize_vector(filter, filter_width * filter_width);
+    }
+  temporal_filter = g_malloc (sizeof (double) * nbatches);
+  temporal_deltas = g_malloc (sizeof (double) * nbatches);
+  if (nbatches > 1)
+    {
+      double t;
+      /* fill in the coefs */
+      for (i = 0; i < nbatches; i++)
+        {
+          t = temporal_deltas[i] = (2.0 * ((double) i / (nbatches - 1)) - 1.0)
+                                   * spec->temporal_filter_radius;
+          temporal_filter[i] = exp(-2.0 * t * t);
+        }
+      normalize_vector(temporal_filter, nbatches);
+    }
+  else
+    {
+      temporal_filter[0] = 1.0;
+      temporal_deltas[0] = 0.0;
+    }
+
+  /* the number of additional rows of buckets we put at the edge so
+     that the filter doesn't go off the edge */
+  gutter_width = (filter_width - oversample) / 2;
+  height = oversample * image_height + 2 * gutter_width;
+  width  = oversample * image_width  + 2 * gutter_width;
+
+  nbuckets = width * height;
+  if (1)
+    {
+      static char *last_block = NULL;
+      static int   last_block_size = 0;
+      int memory_rqd = (sizeof (bucket) * nbuckets +
+                        sizeof (abucket) * nbuckets +
+                        sizeof (point) * SUB_BATCH_SIZE);
+      if (memory_rqd > last_block_size)
+        {
+          if (last_block != NULL)
+            free (last_block);
+          last_block = g_try_malloc (memory_rqd);
+          if (last_block == NULL)
+            {
+              g_printerr ("render_rectangle: cannot malloc %d bytes.\n",
+                          memory_rqd);
+              exit (1);
+            }
+          last_block_size = memory_rqd;
+        }
+      buckets = (bucket *) last_block;
+      accumulate = (abucket *) (last_block + sizeof (bucket) * nbuckets);
+      points = (point *)  (last_block + (sizeof (bucket) + sizeof (abucket)) * nbuckets);
+    }
+
+  memset ((char *) accumulate, 0, sizeof (abucket) * nbuckets);
+  for (batch_num = 0; batch_num < nbatches; batch_num++)
+    {
+      double        batch_time;
+      double        sample_density;
+      control_point cp;
+      memset ((char *) buckets, 0, sizeof (bucket) * nbuckets);
+      batch_time = spec->time + temporal_deltas[batch_num];
+
+      /* interpolate and get a control point */
+      interpolate (spec->cps, spec->ncps, batch_time, &cp);
+
+      /* compute the colormap entries.  the input colormap is 256 long with
+         entries from 0 to 1.0 */
+      for (j = 0; j < CMAP_SIZE; j++)
+        {
+          for (k = 0; k < 3; k++)
+            {
+#if 1
+              cmap[j][k] = (int) (cp.cmap[(j * 256) / CMAP_SIZE][k] *
+                                  cp.white_level);
+#else
+              /* monochrome if you don't have any cmaps */
+              cmap[j][k] = cp.white_level;
+#endif
+            }
+          cmap[j][3] = cp.white_level;
+        }
+      /* compute camera */
+      if (1)
+        {
+          double t0, t1, shift = 0.0, corner0, corner1;
+          double scale;
+
+          scale = pow (2.0, cp.zoom);
+          sample_density = cp.sample_density * scale * scale;
+
+          ppux = cp.pixels_per_unit * scale;
+          ppuy = field ? (ppux / 2.0) : ppux;
+          switch (field)
+            {
+            case field_both:
+              shift =  0.0;
+              break;
+            case field_even:
+              shift = -0.5;
+              break;
+            case field_odd:
+              shift =  0.5;
+              break;
+            }
+          shift = shift / ppux;
+          t0 = (double) gutter_width / (oversample * ppux);
+          t1 = (double) gutter_width / (oversample * ppuy);
+          corner0 = cp.center[0] - image_width / ppux / 2.0;
+          corner1 = cp.center[1] - image_height / ppuy / 2.0;
+          bounds[0] = corner0 - t0;
+          bounds[1] = corner1 - t1 + shift;
+          bounds[2] = corner0 + image_width  / ppux + t0;
+          bounds[3] = corner1 + image_height / ppuy + t1 + shift;
+          size[0] = 1.0 / (bounds[2] - bounds[0]);
+          size[1] = 1.0 / (bounds[3] - bounds[1]);
+        }
+      nsamples = (int) (sample_density * nbuckets /
+                        (oversample * oversample));
+      batch_size = nsamples / cp.nbatches;
+
+      sbc = 0;
+      for (sub_batch = 0;
+           sub_batch < batch_size;
+           sub_batch += SUB_BATCH_SIZE)
+        {
+          if (progress && (sbc++ % 32) == 0)
+              (*progress)(0.5 * sub_batch / (double) batch_size);
+          /* generate a sub_batch_size worth of samples */
+          points[0][0] = random_uniform11 ();
+          points[0][1] = random_uniform11 ();
+          points[0][2] = random_uniform01 ();
+          iterate (&cp, SUB_BATCH_SIZE, FUSE, points);
+
+          /* merge them into buckets, looking up colors */
+          for (j = 0; j < SUB_BATCH_SIZE; j++)
+            {
+              int k, color_index;
+              double *p = points[j];
+              bucket *b;
+
+              /* Note that we must test if p[0] and p[1] is "within"
+               * the valid bounds rather than "not outside", because
+               * p[0] and p[1] might be NaN.
+               */
+              if (p[0] >= bounds[0] &&
+                  p[1] >= bounds[1] &&
+                  p[0] <= bounds[2] &&
+                  p[1] <= bounds[3])
+                {
+                  color_index = (int) (p[2] * CMAP_SIZE);
+
+                  if (color_index < 0)
+                    color_index = 0;
+                  else if (color_index > CMAP_SIZE - 1)
+                    color_index = CMAP_SIZE - 1;
+
+                  b = buckets +
+                      (int) (width * (p[0] - bounds[0]) * size[0]) +
+                      width * (int) (height * (p[1] - bounds[1]) * size[1]);
+
+                  for (k = 0; k < 4; k++)
+                    bump_no_overflow(b[0][k], cmap[color_index][k], short);
+                }
+            }
+        }
+
+      if (1)
+        {
+          double k1 = (cp.contrast * cp.brightness *
+                       PREFILTER_WHITE * 268.0 *
+                       temporal_filter[batch_num]) / 256;
+          double area = image_width * image_height / (ppux * ppuy);
+          double k2 = (oversample * oversample * nbatches) /
+                       (cp.contrast * area * cp.white_level * sample_density);
+
+          /* log intensity in hsv space */
+          for (j = 0; j < height; j++)
+            for (i = 0; i < width; i++)
+              {
+                abucket *a = accumulate + i + j * width;
+                bucket *b = buckets + i + j * width;
+                double c[4], ls;
+                c[0] = (double) b[0][0];
+                c[1] = (double) b[0][1];
+                c[2] = (double) b[0][2];
+                c[3] = (double) b[0][3];
+                if (0.0 == c[3])
+                  continue;
+
+                ls = (k1 * log(1.0 + c[3] * k2))/c[3];
+                c[0] *= ls;
+                c[1] *= ls;
+                c[2] *= ls;
+                c[3] *= ls;
+
+                bump_no_overflow(a[0][0], c[0] + 0.5, accum_t);
+                bump_no_overflow(a[0][1], c[1] + 0.5, accum_t);
+                bump_no_overflow(a[0][2], c[2] + 0.5, accum_t);
+                bump_no_overflow(a[0][3], c[3] + 0.5, accum_t);
+              }
+        }
+    }
+  /*
+   * filter the accumulation buffer down into the image
+   */
+  if (1)
+    {
+      int    x, y;
+      double t[4];
+      double g = 1.0 / spec->cps[0].gamma;
+      y = 0;
+      for (j = 0; j < image_height; j++)
+        {
+          if (progress && (j % 32) == 0)
+            (*progress)(0.5 + 0.5 * j / (double)image_height);
+          x = 0;
+          for (i = 0; i < image_width; i++)
+            {
+              int            ii, jj, a;
+              unsigned char *p;
+              t[0] = t[1] = t[2] = t[3] = 0.0;
+              for (ii = 0; ii < filter_width; ii++)
+                for (jj = 0; jj < filter_width; jj++)
+                  {
+                    double k = filter[ii + jj * filter_width];
+                    abucket *a = accumulate + x + ii + (y + jj) * width;
+
+                    t[0] += k * a[0][0];
+                    t[1] += k * a[0][1];
+                    t[2] += k * a[0][2];
+                    t[3] += k * a[0][3];
+                  }
+              /* FIXME: we should probably use glib facilities to make
+               * this code readable
+               */
+              p = out + nchan * (i + j * out_width);
+              a = 256.0 * pow((double) t[0] / PREFILTER_WHITE, g) + 0.5;
+              if (a < 0) a = 0; else if (a > 255) a = 255;
+              p[0] = a;
+              a = 256.0 * pow((double) t[1] / PREFILTER_WHITE, g) + 0.5;
+              if (a < 0) a = 0; else if (a > 255) a = 255;
+              p[1] = a;
+              a = 256.0 * pow((double) t[2] / PREFILTER_WHITE, g) + 0.5;
+              if (a < 0) a = 0; else if (a > 255) a = 255;
+              p[2] = a;
+              if (nchan > 3)
+                {
+                  a = 256.0 * pow((double) t[3] / PREFILTER_WHITE, g) + 0.5;
+                  if (a < 0) a = 0; else if (a > 255) a = 255;
+                  p[3] = a;
+                }
+              x += oversample;
+            }
+          y += oversample;
+        }
+    }
+
+  free (filter);
+  free (temporal_filter);
+  free (temporal_deltas);
+}