From 3c57dd931145d43f2b0aef96c4d178135956bf91 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Fri, 19 Apr 2024 05:13:10 +0200 Subject: Adding upstream version 2.10.36. Signed-off-by: Daniel Baumann --- plug-ins/flame/rect.c | 398 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 398 insertions(+) create mode 100644 plug-ins/flame/rect.c (limited to 'plug-ins/flame/rect.c') 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 + + 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 . +*/ + +#include "rect.h" + +#include + +#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); +} -- cgit v1.2.3