399 lines
13 KiB
C
399 lines
13 KiB
C
/*
|
|
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 "libpika/pika.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);
|
|
}
|
