« back to all changes in this revision
Viewing changes to flam3/rect.c
- browse files at revision 3
- compare with another revision
- download diff
- download tarball
- view history from revision 3
- Committer: Bazaar Package Importer
- Author(s): Alexis Sukrieh
- Date: 2005-06-17 18:28:17 UTC
- mfrom: (1.1.1 upstream)
- Revision ID: james.westby@ubuntu.com-20050617182817-3ju56pnktawusjt8
Tags: 2.6.2-1
* New upstream release: 2.6.2.
(closes: #309852)
* Added myself as maintainer (Sander Smeenk offered it to me).
* Added dpatch for handling changes applied on upstream sources.
* Added all missing manpages.
(closes: #309852)
* Added myself as maintainer (Sander Smeenk offered it to me).
* Added dpatch for handling changes applied on upstream sources.
* Added all missing manpages.
- files added:
- compile
- debian/manpages
- debian/patches
- flam3
- files removed:
- debian/electricsheep.1
- flame
- files modified:
- Makefile.am
added
removed
flam3/rect.c
1
/*
2
flame - cosmic recursive fractal flames
3
Copyright (C) 1992-2003 Scott Draves <source@flam3.com>
4
5
This program is free software; you can redistribute it and/or modify
6
it under the terms of the GNU General Public License as published by
7
the Free Software Foundation; either version 2 of the License, or
8
(at your option) any later version.
9
10
This program is distributed in the hope that it will be useful,
11
but WITHOUT ANY WARRANTY; without even the implied warranty of
12
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
GNU General Public License for more details.
14
15
You should have received a copy of the GNU General Public License
16
along with this program; if not, write to the Free Software
17
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18
*/
19
20
21
22
23
24
/* for batch
25
* interpolate
26
* compute colormap
27
* for subbatch
28
* compute samples
29
* buckets += cmap[samples]
30
* accum += time_filter[batch] * log(buckets)
31
* image = filter(accum)
32
*/
33
34
35
36
/* allow this many iterations for settling into attractor */
37
#define FUSE 15
38
39
/* clamp spatial filter to zero at this std dev (2.5 ~= 0.0125) */
40
#define FILTER_CUTOFF 1.8
41
42
#define PREFILTER_WHITE 255
43
#define WHITE_LEVEL 255
44
#define SUB_BATCH_SIZE 10000
45
46
47
48
static void render_rectangle(spec, out, out_width, field, nchan)
49
flam3_frame *spec;
50
unsigned char *out;
51
int out_width;
52
int field;
53
int nchan;
54
{
55
int i, j, k, nbuckets, batch_num;
56
double nsamples, batch_size, sub_batch;
57
bucket *buckets;
58
abucket *accumulate;
59
double *points;
60
double *filter, *temporal_filter, *temporal_deltas;
61
double bounds[4], size[2], ppux, ppuy;
62
double rot[2][2];
63
int image_width, image_height; /* size of the image to produce */
64
int width, height; /* size of histogram */
65
int filter_width;
66
int oversample = spec->genomes[0].spatial_oversample;
67
int nbatches = spec->genomes[0].nbatches;
68
bucket cmap[CMAP_SIZE];
69
unsigned char *cmap2;
70
int gutter_width;
71
int sbc;
72
double vibrancy = 0.0;
73
double gamma = 0.0;
74
double background[3];
75
int vib_gam_n = 0;
76
time_t progress_timer = 0, progress_began,
77
progress_timer_history[64] = {0};
78
double progress_history[64] = {0};
79
int progress_history_mark = 0;
80
int verbose = spec->verbose;
81
char *fname = getenv("image");
82
83
if (nbatches < 1) {
84
fprintf(stderr, "nbatches must be positive," " not %d.\n", nbatches);
85
exit(1);
86
}
87
88
if (oversample < 1) {
89
fprintf(stderr, "oversample must be positive," " not %d.\n", oversample);
90
exit(1);
91
}
92
93
image_width = spec->genomes[0].width;
94
if (field) {
95
image_height = spec->genomes[0].height / 2;
96
if (field == flam3_field_odd)
97
out += nchan * out_width;
98
out_width *= 2;
99
} else
100
image_height = spec->genomes[0].height;
101
102
if (1) {
103
double fw = (2.0 * FILTER_CUTOFF * oversample *
104
spec->genomes[0].spatial_filter_radius /
105
spec->pixel_aspect_ratio);
106
double adjust;
107
filter_width = ((int) fw) + 1;
108
/* make sure it has same parity as oversample */
109
if ((filter_width ^ oversample) & 1)
110
filter_width++;
111
if (fw > 0.0)
112
adjust = FILTER_CUTOFF * filter_width / fw;
113
else
114
adjust = 1.0;
115
116
#if 0
117
fprintf(stderr, "filter_width = %d adjust=%g\n",
118
filter_width, adjust);
119
#endif
120
121
filter = (double *) malloc(sizeof(double) * filter_width * filter_width);
122
/* fill in the coefs */
123
for (i = 0; i < filter_width; i++)
124
for (j = 0; j < filter_width; j++) {
125
double ii = ((2.0 * i + 1.0) / filter_width - 1.0) * adjust;
126
double jj = ((2.0 * j + 1.0) / filter_width - 1.0) * adjust;
127
if (field)
128
jj *= 2.0;
129
jj /= spec->pixel_aspect_ratio;
130
filter[i + j * filter_width] =
131
exp(-2.0 * (ii * ii + jj * jj));
132
}
133
134
normalize_vector(filter, filter_width * filter_width);
135
#if 0
136
printf("vvvvvvvvvvvvvvvvvvvvvvvvvvvv\n");
137
for (j = 0; j < filter_width; j++) {
138
for (i = 0; i < filter_width; i++) {
139
printf(" %5d", (int)(10000 * filter[i + j * filter_width]));
140
}
141
printf("\n");
142
}
143
printf("^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n");
144
fflush(stdout);
145
#endif
146
}
147
temporal_filter = (double *) malloc(sizeof(double) * nbatches);
148
temporal_deltas = (double *) malloc(sizeof(double) * nbatches);
149
if (nbatches > 1) {
150
double t;
151
/* fill in the coefs */
152
for (i = 0; i < nbatches; i++) {
153
t = temporal_deltas[i] = (2.0 * ((double) i / (nbatches - 1)) - 1.0)
154
* spec->temporal_filter_radius;
155
temporal_filter[i] = exp(-2.0 * t * t);
156
}
157
158
normalize_vector(temporal_filter, nbatches);
159
} else {
160
temporal_filter[0] = 1.0;
161
temporal_deltas[0] = 0.0;
162
}
163
164
#if 0
165
for (j = 0; j < nbatches; j++)
166
fprintf(stderr, "%4f %4f\n", temporal_deltas[j], temporal_filter[j]);
167
fprintf(stderr, "\n");
168
#endif
169
170
/* the number of additional rows of buckets we put at the edge so
171
that the filter doesn't go off the edge */
172
173
gutter_width = (filter_width - oversample) / 2;
174
height = oversample * image_height + 2 * gutter_width;
175
width = oversample * image_width + 2 * gutter_width;
176
177
nbuckets = width * height;
178
if (1) {
179
char *last_block = NULL;
180
int memory_rqd = (sizeof(bucket) * nbuckets +
181
sizeof(abucket) * nbuckets +
182
4 * sizeof(double) * SUB_BATCH_SIZE);
183
last_block = (char *) malloc(memory_rqd);
184
if (NULL == last_block) {
185
fprintf(stderr, "render_rectangle: cannot malloc %d bytes.\n", memory_rqd);
186
fprintf(stderr, "render_rectangle: h=%d w=%d nb=%d.\n", width, height, nbuckets);
187
exit(1);
188
}
189
/* else fprintf(stderr, "render_rectangle: mallocked %dMb.\n", Mb); */
190
191
buckets = (bucket *) last_block;
192
accumulate = (abucket *) (last_block + sizeof(bucket) * nbuckets);
193
points = (double *) (last_block + (sizeof(bucket) + sizeof(abucket)) * nbuckets);
194
}
195
196
if (verbose) {
197
fprintf(stderr, "chaos: ");
198
progress_began = time(NULL);
199
}
200
201
if (fname) {
202
int len = strlen(fname);
203
FILE *fin = fopen(fname, "rb");
204
int w, h;
205
cmap2 = NULL;
206
if (len > 4) {
207
char *ending = fname + len - 4;
208
if (!strcmp(ending, ".png")) {
209
cmap2 = read_png(fin, &w, &h);
210
} else if (!strcmp(ending, ".jpg")) {
211
cmap2 = read_jpeg(fin, &w, &h);
212
}
213
}
214
if (NULL == cmap2) {
215
perror(fname);
216
exit(1);
217
}
218
if (256 != w || 256 != h) {
219
fprintf(stderr, "image must be 256 by 256, not %d by %d.\n", w, h);
220
exit(1);
221
}
222
}
223
224
background[0] = background[1] = background[2] = 0.0;
225
memset((char *) accumulate, 0, sizeof(abucket) * nbuckets);
226
for (batch_num = 0; batch_num < nbatches; batch_num++) {
227
double batch_time;
228
double sample_density;
229
flam3_genome cp;
230
memset((char *) buckets, 0, sizeof(bucket) * nbuckets);
231
batch_time = spec->time + temporal_deltas[batch_num];
232
233
/* interpolate and get a control point */
234
flam3_interpolate(spec->genomes, spec->ngenomes, batch_time, &cp);
235
236
/* compute the colormap entries. the input colormap is 256 long with
237
entries from 0 to 1.0 */
238
if (!fname) {
239
for (j = 0; j < CMAP_SIZE; j++) {
240
for (k = 0; k < 3; k++) {
241
#if 1
242
cmap[j][k] = (int) (cp.palette[(j * 256) / CMAP_SIZE][k] * WHITE_LEVEL);
243
#else
244
/* monochrome if you don't have any cmaps */
245
cmap[j][k] = WHITE_LEVEL;
246
#endif
247
}
248
cmap[j][3] = WHITE_LEVEL;
249
}
250
}
251
252
/* compute camera */
253
if (1) {
254
double t0, t1, shift, corner0, corner1;
255
double scale;
256
257
if (cp.sample_density <= 0.0) {
258
fprintf(stderr,
259
"sample density (quality) must be greater than zero,"
260
" not %g.\n", cp.sample_density);
261
exit(1);
262
}
263
264
scale = pow(2.0, cp.zoom);
265
sample_density = cp.sample_density * scale * scale;
266
267
268
ppux = cp.pixels_per_unit * scale;
269
ppuy = field ? (ppux / 2.0) : ppux;
270
ppux /= spec->pixel_aspect_ratio;
271
switch (field) {
272
case flam3_field_both: shift = 0.0; break;
273
case flam3_field_even: shift = -0.5; break;
274
case flam3_field_odd: shift = 0.5; break;
275
}
276
shift = shift / ppux;
277
t0 = (double) gutter_width / (oversample * ppux);
278
t1 = (double) gutter_width / (oversample * ppuy);
279
corner0 = cp.center[0] - image_width / ppux / 2.0;
280
corner1 = cp.center[1] - image_height / ppuy / 2.0;
281
bounds[0] = corner0 - t0;
282
bounds[1] = corner1 - t1 + shift;
283
bounds[2] = corner0 + image_width / ppux + t0;
284
bounds[3] = corner1 + image_height / ppuy + t1 + shift;
285
size[0] = 1.0 / (bounds[2] - bounds[0]);
286
size[1] = 1.0 / (bounds[3] - bounds[1]);
287
rot[0][0] = cos(cp.rotate * 2 * M_PI / 360.0);
288
rot[0][1] = -sin(cp.rotate * 2 * M_PI / 360.0);
289
rot[1][0] = -rot[0][1];
290
rot[1][1] = rot[0][0];
291
}
292
293
nsamples = (sample_density * (double) nbuckets /
294
(oversample * oversample));
295
#if 0
296
fprintf(stderr, "sample_density=%g nsamples=%g nbuckets=%d\n",
297
sample_density, nsamples, nbuckets);
298
#endif
299
300
batch_size = nsamples / cp.nbatches;
301
302
sbc = 0;
303
for (sub_batch = 0;
304
sub_batch < batch_size;
305
sub_batch += SUB_BATCH_SIZE) {
306
307
if (spec->progress&&!(sbc++&7))
308
if ((*spec->progress)(spec->progress_parameter, 0.5*sub_batch/(double)batch_size))
309
return;
310
311
if (verbose && time(NULL) != progress_timer) {
312
double percent = 100.0 * ((sub_batch / (double) batch_size) + batch_num) / nbatches;
313
double eta;
314
int old_mark = 0;
315
316
fprintf(stderr, "\rchaos: %5.1f%%", percent);
317
progress_timer = time(NULL);
318
if (progress_timer_history[progress_history_mark] &&
319
progress_history[progress_history_mark] < percent)
320
old_mark = progress_history_mark;
321
322
if (percent > 0) {
323
eta = (100 - percent) * (progress_timer - progress_timer_history[old_mark])
324
/ (percent - progress_history[old_mark]);
325
326
if (eta > 100)
327
fprintf(stderr, " ETA: %.1f minutes", eta / 60);
328
else
329
fprintf(stderr, " ETA: %ld seconds ", (long) ceil(eta));
330
fprintf(stderr, " \r");
331
}
332
333
progress_timer_history[progress_history_mark] = progress_timer;
334
progress_history[progress_history_mark] = percent;
335
progress_history_mark = (progress_history_mark + 1) % 64;
336
}
337
338
339
/* generate a sub_batch_size worth of samples */
340
points[0] = flam3_random11();
341
points[1] = flam3_random11();
342
points[2] = flam3_random01();
343
points[3] = flam3_random01();
344
flam3_iterate(&cp, SUB_BATCH_SIZE, FUSE, points);
345
346
347
/* merge them into buckets, looking up colors */
348
for (j = 0; j < SUB_BATCH_SIZE; j++) {
349
double p0, p1, p00, p11;
350
int k, color_index0, color_index1;
351
double *p = &points[4*j];
352
bucket *b;
353
354
if (cp.rotate != 0.0) {
355
p00 = p[0] - cp.rot_center[0];
356
p11 = p[1] - cp.rot_center[1];
357
p0 = p00 * rot[0][0] + p11 * rot[0][1] + cp.rot_center[0];
358
p1 = p00 * rot[1][0] + p11 * rot[1][1] + cp.rot_center[1];
359
} else {
360
p0 = p[0];
361
p1 = p[1];
362
}
363
364
if (p0 < bounds[0] || p1 < bounds[1] ||
365
p0 > bounds[2] || p1 > bounds[3])
366
continue;
367
if (fname) {
368
int ci;
369
color_index0 = (int) (p[2] * CMAP_SIZE);
370
color_index1 = (int) (p[3] * CMAP_SIZE);
371
if (color_index0 < 0) color_index0 = 0;
372
else if (color_index0 > (CMAP_SIZE-1))
373
color_index0 = CMAP_SIZE-1;
374
if (color_index1 < 0) color_index1 = 0;
375
else if (color_index1 > (CMAP_SIZE-1))
376
color_index1 = CMAP_SIZE-1;
377
b = buckets +
378
(int) (width * (p0 - bounds[0]) * size[0]) +
379
width * (int) (height * (p1 - bounds[1]) * size[1]);
380
ci = 4 * (CMAP_SIZE * color_index1 + color_index0);
381
for (k = 0; k < 4; k++)
382
// b[0][k] += cmap2[ci+k];
383
bump_no_overflow(b[0][k], cmap2[ci+k], ACCUM_T);
384
} else {
385
color_index0 = (int) (p[2] * CMAP_SIZE);
386
if (color_index0 < 0) color_index0 = 0;
387
else if (color_index0 > (CMAP_SIZE-1))
388
color_index0 = CMAP_SIZE-1;
389
390
b = buckets +
391
(int) (width * (p0 - bounds[0]) * size[0]) +
392
width * (int) (height * (p1 - bounds[1]) * size[1]);
393
for (k = 0; k < 4; k++)
394
// b[0][k] += cmap[color_index0][k];
395
bump_no_overflow(b[0][k], cmap[color_index0][k], ACCUM_T);
396
}
397
}
398
}
399
400
if (1) {
401
double k1 =(cp.contrast * cp.brightness *
402
PREFILTER_WHITE * 268.0 *
403
temporal_filter[batch_num]) / 256;
404
double area = image_width * image_height / (ppux * ppuy);
405
double k2 = (oversample * oversample * nbatches) /
406
(cp.contrast * area * WHITE_LEVEL * sample_density);
407
408
vibrancy += cp.vibrancy;
409
gamma += cp.gamma;
410
background[0] += cp.background[0];
411
background[1] += cp.background[1];
412
background[2] += cp.background[2];
413
vib_gam_n++;
414
415
/* log intensity */
416
for (j = 0; j < height; j++) {
417
for (i = 0; i < width; i++) {
418
abucket *a = accumulate + i + j * width;
419
bucket *b = buckets + i + j * width;
420
double c[4], ls;
421
422
c[0] = (double) b[0][0];
423
c[1] = (double) b[0][1];
424
c[2] = (double) b[0][2];
425
c[3] = (double) b[0][3];
426
if (0.0 == c[3])
427
continue;
428
429
ls = (k1 * log(1.0 + c[3] * k2))/c[3];
430
c[0] *= ls;
431
c[1] *= ls;
432
c[2] *= ls;
433
c[3] *= ls;
434
435
bump_no_overflow(a[0][0], c[0], ACCUM_T);
436
bump_no_overflow(a[0][1], c[1], ACCUM_T);
437
bump_no_overflow(a[0][2], c[2], ACCUM_T);
438
bump_no_overflow(a[0][3], c[3], ACCUM_T);
439
}
440
}
441
}
442
}
443
444
if (verbose) {
445
fprintf(stderr, "\rchaos: 100.0%% took: %ld seconds \n",
446
time(NULL) - progress_began);
447
fprintf(stderr, "filtering...");
448
}
449
450
/*
451
* filter the accumulation buffer down into the image
452
*/
453
if (1) {
454
int x, y;
455
double t[4];
456
double g = 1.0 / (gamma / vib_gam_n);
457
vibrancy /= vib_gam_n;
458
background[0] /= vib_gam_n/256.0;
459
background[1] /= vib_gam_n/256.0;
460
background[2] /= vib_gam_n/256.0;
461
y = 0;
462
463
for (j = 0; j < image_height; j++) {
464
if (spec->progress && !(j&15))
465
if ((*spec->progress)(spec->progress_parameter, 0.5+0.5*j/(double)image_height))
466
return;
467
x = 0;
468
for (i = 0; i < image_width; i++) {
469
int ii, jj;
470
unsigned char *p;
471
double ls, a;
472
double alpha;
473
t[0] = t[1] = t[2] = t[3] = 0.0;
474
for (ii = 0; ii < filter_width; ii++)
475
for (jj = 0; jj < filter_width; jj++) {
476
double k = filter[ii + jj * filter_width];
477
abucket *a = accumulate + x + ii + (y + jj) * width;
478
479
t[0] += k * a[0][0];
480
t[1] += k * a[0][1];
481
t[2] += k * a[0][2];
482
t[3] += k * a[0][3];
483
}
484
p = out + nchan * (i + j * out_width);
485
486
alpha = t[3];
487
if (alpha > 0.0) {
488
ls = vibrancy * 256.0 * pow(alpha / PREFILTER_WHITE, g)
489
/ (alpha / PREFILTER_WHITE);
490
alpha = pow(alpha / PREFILTER_WHITE, g);
491
if (alpha < 0.0) alpha = 0.0;
492
else if (alpha > 1.0) alpha = 1.0;
493
} else
494
ls = 0.0;
495
496
/* apply to rgb channels the relative scale from gamma of alpha channel */
497
498
a = ls * ((double) t[0] / PREFILTER_WHITE);
499
a += (1.0-vibrancy) * 256.0 * pow((double) t[0] / PREFILTER_WHITE, g);
500
a += ((1.0 - alpha) * background[0]);
501
if (a < 0) a = 0; else if (a > 255) a = 255;
502
p[0] = (unsigned char) a;
503
504
a = ls * ((double) t[1] / PREFILTER_WHITE);
505
a += (1.0-vibrancy) * 256.0 * pow((double) t[1] / PREFILTER_WHITE, g);
506
a += ((1.0 - alpha) * background[1]);
507
if (a < 0) a = 0; else if (a > 255) a = 255;
508
p[1] = (unsigned char) a;
509
510
a = ls * ((double) t[2] / PREFILTER_WHITE);
511
a += (1.0-vibrancy) * 256.0 * pow((double) t[2] / PREFILTER_WHITE, g);
512
a += ((1.0 - alpha) * background[2]);
513
if (a < 0) a = 0; else if (a > 255) a = 255;
514
p[2] = (unsigned char) a;
515
516
if (nchan > 3)
517
p[3] = (unsigned char) (alpha * 255.999);
518
519
x += oversample;
520
}
521
y += oversample;
522
}
523
}
524
525
free(filter);
526
free(buckets);
527
if (fname) free(cmap2);
528
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1