1 /*
2 flame - cosmic recursive fractal flames
3 Copyright (C) 1992 Scott Draves <spot@cs.cmu.edu>
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 3 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, see <https://www.gnu.org/licenses/>.
17 */
18
19 #include "config.h"
20
21 #include <stdlib.h>
22 #include <string.h> /* strcmp */
23
24 #include "libgimp/gimp.h"
25
26 #include "libifs.h"
27
28 #define CHOOSE_XFORM_GRAIN 100
29
30 static int flam3_random_bit (void);
31 static double flam3_random01 (void);
32
33 /*
34 * run the function system described by CP forward N generations.
35 * store the n resulting 3 vectors in POINTS. the initial point is passed
36 * in POINTS[0]. ignore the first FUSE iterations.
37 */
38
39 void
iterate(control_point * cp,int n,int fuse,point * points)40 iterate (control_point *cp,
41 int n,
42 int fuse,
43 point *points)
44 {
45 int i, j, count_large = 0, count_nan = 0;
46 int xform_distrib[CHOOSE_XFORM_GRAIN];
47 double p[3], t, r, dr;
48 p[0] = points[0][0];
49 p[1] = points[0][1];
50 p[2] = points[0][2];
51
52 /*
53 * first, set up xform, which is an array that converts a uniform random
54 * variable into one with the distribution dictated by the density
55 * fields
56 */
57 dr = 0.0;
58 for (i = 0; i < NXFORMS; i++)
59 dr += cp->xform[i].density;
60 dr = dr / CHOOSE_XFORM_GRAIN;
61
62 j = 0;
63 t = cp->xform[0].density;
64 r = 0.0;
65 for (i = 0; i < CHOOSE_XFORM_GRAIN; i++)
66 {
67 while (r >= t)
68 {
69 j++;
70 t += cp->xform[j].density;
71 }
72 xform_distrib[i] = j;
73 r += dr;
74 }
75
76 for (i = -fuse; i < n; i++)
77 {
78 /* FIXME: the following is supported only by gcc and c99 */
79 int fn = xform_distrib[g_random_int_range (0, CHOOSE_XFORM_GRAIN)];
80 double tx, ty, v;
81
82 if (p[0] > 100.0 || p[0] < -100.0 ||
83 p[1] > 100.0 || p[1] < -100.0)
84 count_large++;
85 if (p[0] != p[0])
86 count_nan++;
87
88 #define coef cp->xform[fn].c
89 #define vari cp->xform[fn].var
90
91 /* first compute the color coord */
92 p[2] = (p[2] + cp->xform[fn].color) / 2.0;
93
94 /* then apply the affine part of the function */
95 tx = coef[0][0] * p[0] + coef[1][0] * p[1] + coef[2][0];
96 ty = coef[0][1] * p[0] + coef[1][1] * p[1] + coef[2][1];
97
98 p[0] = p[1] = 0.0;
99 /* then add in proportional amounts of each of the variations */
100 v = vari[0];
101 if (v > 0.0)
102 {
103 /* linear */
104 double nx, ny;
105 nx = tx;
106 ny = ty;
107 p[0] += v * nx;
108 p[1] += v * ny;
109 }
110
111 v = vari[1];
112 if (v > 0.0)
113 {
114 /* sinusoidal */
115 double nx, ny;
116 nx = sin (tx);
117 ny = sin (ty);
118 p[0] += v * nx;
119 p[1] += v * ny;
120 }
121
122 v = vari[2];
123 if (v > 0.0)
124 {
125 /* spherical */
126 double nx, ny;
127 double r2 = tx * tx + ty * ty + 1e-6;
128 nx = tx / r2;
129 ny = ty / r2;
130 p[0] += v * nx;
131 p[1] += v * ny;
132 }
133
134 v = vari[3];
135 if (v > 0.0)
136 {
137 /* swirl */
138 double r2 = tx * tx + ty * ty; /* /k here is fun */
139 double c1 = sin (r2);
140 double c2 = cos (r2);
141 double nx = c1 * tx - c2 * ty;
142 double ny = c2 * tx + c1 * ty;
143 p[0] += v * nx;
144 p[1] += v * ny;
145 }
146
147 v = vari[4];
148 if (v > 0.0)
149 {
150 /* horseshoe */
151 double a, c1, c2, nx, ny;
152 if (tx < -EPS || tx > EPS ||
153 ty < -EPS || ty > EPS)
154 a = atan2(tx, ty); /* times k here is fun */
155 else
156 a = 0.0;
157 c1 = sin (a);
158 c2 = cos (a);
159 nx = c1 * tx - c2 * ty;
160 ny = c2 * tx + c1 * ty;
161 p[0] += v * nx;
162 p[1] += v * ny;
163 }
164
165 v = vari[5];
166 if (v > 0.0)
167 {
168 /* polar */
169 double nx, ny;
170 if (tx < -EPS || tx > EPS ||
171 ty < -EPS || ty > EPS)
172 nx = atan2 (tx, ty) / G_PI;
173 else
174 nx = 0.0;
175
176 ny = sqrt (tx * tx + ty * ty) - 1.0;
177 p[0] += v * nx;
178 p[1] += v * ny;
179 }
180
181 v = vari[6];
182 if (v > 0.0)
183 {
184 /* bent */
185 double nx, ny;
186 nx = tx;
187 ny = ty;
188 if (nx < 0.0) nx = nx * 2.0;
189 if (ny < 0.0) ny = ny / 2.0;
190 p[0] += v * nx;
191 p[1] += v * ny;
192 }
193
194 v = vari[7];
195 if (v > 0.0)
196 {
197 /* folded handkerchief */
198 double theta, r2, nx, ny;
199 if (tx < -EPS || tx > EPS ||
200 ty < -EPS || ty > EPS)
201 theta = atan2( tx, ty );
202 else
203 theta = 0.0;
204 r2 = sqrt (tx * tx + ty * ty);
205 nx = sin (theta + r2) * r2;
206 ny = cos (theta - r2) * r2;
207 p[0] += v * nx;
208 p[1] += v * ny;
209 }
210
211 v = vari[8];
212 if (v > 0.0)
213 {
214 /* heart */
215 double theta, r2, nx, ny;
216 if (tx < -EPS || tx > EPS ||
217 ty < -EPS || ty > EPS)
218 theta = atan2( tx, ty );
219 else
220 theta = 0.0;
221 r2 = sqrt (tx * tx + ty * ty);
222 theta *= r2;
223 nx = sin (theta) * r2;
224 ny = cos (theta) * -r2;
225 p[0] += v * nx;
226 p[1] += v * ny;
227 }
228
229 v = vari[9];
230 if (v > 0.0)
231 {
232 /* disc */
233 double theta, r2, nx, ny;
234 if ( tx < -EPS || tx > EPS ||
235 ty < - EPS || ty > EPS)
236 theta = atan2 (tx, ty);
237 else
238 theta = 0.0;
239 nx = tx * G_PI;
240 ny = ty * G_PI;
241 r2 = sqrt (nx * nx * ny * ny);
242 p[0] += v * sin(r2) * theta / G_PI;
243 p[1] += v * cos(r2) * theta / G_PI;
244 }
245
246 v = vari[10];
247 if (v > 0.0)
248 {
249 /* spiral */
250 double theta, r2;
251 if (tx < -EPS || tx > EPS ||
252 ty < -EPS || ty > EPS)
253 theta = atan2( tx, ty );
254 else
255 theta = 0.0;
256 r2 = sqrt (tx * tx + ty * ty) + 1e-6;
257 p[0] += v * (cos (theta) + sin (r2)) / r2;
258 p[1] += v * (cos (theta) + cos (r2)) / r2;
259 }
260
261 v = vari[11];
262 if (v > 0.0)
263 {
264 /* hyperbolic */
265 double theta, r2;
266 if (tx < -EPS || tx > EPS ||
267 ty < -EPS || ty > EPS)
268 theta = atan2 (tx, ty);
269 else
270 theta = 0.0;
271 r2 = sqrt (tx * tx + ty * ty) + 1e-6;
272 p[0] += v * sin (theta) / r2;
273 p[1] += v * cos (theta) * r2;
274 }
275
276 v = vari[12];
277 if (v > 0.0 )
278 {
279 double theta, r2;
280 /* diamond */
281 if ( tx < -EPS || tx > EPS ||
282 ty < -EPS || ty > EPS)
283 theta = atan2 (tx, ty);
284 else
285 theta = 0.0;
286 r2 = sqrt( tx * tx + ty * ty );
287 p[0] += v * sin (theta) * cos (r2);
288 p[1] += v * cos (theta) * sin (r2);
289 }
290
291 v = vari[13];
292 if (v > 0.0)
293 {
294 /* ex */
295 double theta, r2, n0, n1, m0, m1;
296 if ( tx < -EPS || tx > EPS ||
297 ty < -EPS || ty > EPS)
298 theta = atan2 (tx, ty);
299 else
300 theta = 0.0;
301 r2 = sqrt( tx * tx + ty * ty );
302 n0 = sin(theta + r2);
303 n1 = cos(theta - r2);
304 m0 = n0 * n0 * n0 * r2;
305 m1 = n1 * n1 * n1 * r2;
306 p[0] += v * (m0 + m1);
307 p[1] += v * (m0 - m1);
308 }
309
310 v = vari[14];
311 if ( v > 0.0)
312 {
313 double theta, r2, nx, ny;
314 /* julia */
315 if (tx < -EPS || tx > EPS ||
316 ty < -EPS || ty > EPS)
317 theta = atan2 (tx, ty);
318 else
319 theta = 0.0;
320 if (flam3_random_bit ())
321 theta += G_PI;
322 r2 = pow (tx * tx + ty * ty, 0.25);
323 nx = r2 * cos (theta);
324 ny = r2 * sin (theta);
325 p[0] += v * nx;
326 p[1] += v * ny;
327 }
328
329 v = vari[15];
330 if (v > 0.0)
331 {
332 /* waves */
333 double dx, dy, nx, ny;
334 dx = coef[2][0];
335 dy = coef[2][1];
336 nx = tx + coef[1][0] * sin (ty / ((dx * dx) + EPS));
337 ny = ty + coef[1][1] * sin (tx / ((dy * dy) + EPS));
338 p[0] += v * nx;
339 p[1] += v * ny;
340 }
341
342 v = vari[16];
343 if (v > 0.0)
344 {
345 /* fisheye */
346 double theta, r2, nx, ny;
347 if (tx < -EPS || tx > EPS ||
348 ty < -EPS || ty > EPS)
349 theta = atan2 (tx, ty);
350 else
351 theta = 0.0;
352 r2 = sqrt (tx * tx + ty * ty);
353 r2 = 2 * r2 / (r2 + 1);
354 nx = r2 * cos (theta);
355 ny = r2 * sin (theta);
356 p[0] += v * nx;
357 p[1] += v * ny;
358 }
359
360 v = vari[17];
361 if (v > 0.0)
362 {
363 /* popcorn */
364 double dx, dy, nx, ny;
365 dx = tan (3 * ty);
366 dy = tan (3 * tx);
367 nx = tx + coef[2][0] * sin (dx);
368 ny = ty + coef[2][1] * sin (dy);
369 p[0] += v * nx;
370 p[1] += v * ny;
371 }
372
373 v = vari[18];
374 if (v > 0.0)
375 {
376 /* exponential */
377 double dx, dy, nx, ny;
378 dx = exp (tx - 1.0);
379 dy = G_PI * ty;
380 nx = cos (dy) * dx;
381 ny = sin (dy) * dx;
382 p[0] += v * nx;
383 p[1] += v * ny;
384 }
385
386 v = vari[19];
387 if (v > 0.0)
388 {
389 /* power */
390 double theta, r2, tsin, tcos, nx, ny;
391 if (tx < -EPS || tx > EPS ||
392 ty < -EPS || ty > EPS)
393 theta = atan2 (tx, ty);
394 else
395 theta = 0.0;
396 tsin = sin (theta);
397 tcos = cos (theta);
398 r2 = sqrt (tx * tx + ty * ty);
399 r2 = pow (r2, tsin);
400 nx = r2 * tcos;
401 ny = r2 * tsin;
402 p[0] += v * nx;
403 p[1] += v * ny;
404 }
405
406 v = vari[20];
407 if (v > 0.0)
408 {
409 /* cosine */
410 double nx, ny;
411 nx = cos (tx * G_PI) * cosh (ty);
412 ny = -sin (tx * G_PI) * sinh (ty);
413 p[0] += v * nx;
414 p[1] += v * ny;
415 }
416
417 v = vari[21];
418 if (v > 0.0)
419 {
420 /* rings */
421 double theta, r2, dx, nx, ny;
422 if (tx < -EPS || tx > EPS ||
423 ty < -EPS || ty > EPS)
424 theta = atan2 (tx, ty);
425 else
426 theta = 0;
427 dx = coef[2][0];
428 dx = dx * dx + EPS;
429 r2 = sqrt (tx * tx + ty * ty);
430 r2 = fmod (r2 + dx, 2 * dx) - dx + r2 * (1 - dx);
431 nx = cos (theta) * r2;
432 ny = sin (theta) * r2;
433 p[0] += v * nx;
434 p[1] += v * ny;
435 }
436
437 v = vari[22];
438 if (v > 0.0)
439 {
440 /* fan */
441 double theta, r2, dx, dy, dx2, nx, ny;
442 if (tx < -EPS || tx > EPS ||
443 ty < -EPS || ty > EPS)
444 theta = atan2 (tx, ty);
445 else
446 theta = 0.0;
447 dx = coef[2][0];
448 dy = coef[2][1];
449 dx = G_PI * (dx * dx + EPS);
450 dx2 = dx / 2;
451 r2 = sqrt (tx * tx + ty * ty );
452 theta += (fmod (theta + dy, dx) > dx2) ? -dx2: dx2;
453 nx = cos (theta) * r2;
454 ny = sin (theta) * r2;
455 p[0] += v * nx;
456 p[1] += v * ny;
457 }
458
459 v = vari[23];
460 if (v > 0.0)
461 {
462 /* eyefish */
463 double r2;
464 r2 = 2.0 * v / (sqrt(tx * tx + ty * ty) + 1.0);
465 p[0] += r2 * tx;
466 p[1] += r2 * ty;
467 }
468
469 v = vari[24];
470 if (v > 0.0)
471 {
472 /* bubble */
473 double r2;
474 r2 = v / ((tx * tx + ty * ty) / 4 + 1);
475 p[0] += r2 * tx;
476 p[1] += r2 * ty;
477 }
478
479 v = vari[25];
480 if (v > 0.0)
481 {
482 /* cylinder */
483 double nx;
484 nx = sin (tx);
485 p[0] += v * nx;
486 p[1] += v * ty;
487 }
488
489 v = vari[26];
490 if (v > 0.0)
491 {
492 /* noise */
493 double rx, sinr, cosr, nois;
494 rx = flam3_random01 () * 2 * G_PI;
495 sinr = sin (rx);
496 cosr = cos (rx);
497 nois = flam3_random01 ();
498 p[0] += v * nois * tx * cosr;
499 p[1] += v * nois * ty * sinr;
500 }
501
502 v = vari[27];
503 if (v > 0.0)
504 {
505 /* blur */
506 double rx, sinr, cosr, nois;
507 rx = flam3_random01 () * 2 * G_PI;
508 sinr = sin (rx);
509 cosr = cos (rx);
510 nois = flam3_random01 ();
511 p[0] += v * nois * cosr;
512 p[1] += v * nois * sinr;
513 }
514
515 v = vari[28];
516 if (v > 0.0)
517 {
518 /* gaussian */
519 double ang, sina, cosa, r2;
520 ang = flam3_random01 () * 2 * G_PI;
521 sina = sin (ang);
522 cosa = cos (ang);
523 r2 = v * (flam3_random01 () + flam3_random01 () + flam3_random01 () +
524 flam3_random01 () - 2.0);
525 p[0] += r2 * cosa;
526 p[1] += r2 * sina;
527 }
528
529 /* if fuse over, store it */
530 if (i >= 0)
531 {
532 points[i][0] = p[0];
533 points[i][1] = p[1];
534 points[i][2] = p[2];
535 }
536 }
537 }
538
539 /* args must be non-overlapping */
540 void
mult_matrix(double s1[2][2],double s2[2][2],double d[2][2])541 mult_matrix (double s1[2][2],
542 double s2[2][2],
543 double d[2][2])
544 {
545 d[0][0] = s1[0][0] * s2[0][0] + s1[1][0] * s2[0][1];
546 d[1][0] = s1[0][0] * s2[1][0] + s1[1][0] * s2[1][1];
547 d[0][1] = s1[0][1] * s2[0][0] + s1[1][1] * s2[0][1];
548 d[1][1] = s1[0][1] * s2[1][0] + s1[1][1] * s2[1][1];
549 }
550
551 static
det_matrix(double s[2][2])552 double det_matrix (double s[2][2])
553 {
554 return s[0][0] * s[1][1] - s[0][1] * s[1][0];
555 }
556
557 static void
interpolate_angle(double t,double s,double * v1,double * v2,double * v3,int tie,int cross)558 interpolate_angle (double t,
559 double s,
560 double *v1,
561 double *v2,
562 double *v3,
563 int tie,
564 int cross)
565 {
566 double x = *v1;
567 double y = *v2;
568 double d;
569 static double lastx, lasty;
570
571 /* take the shorter way around the circle... */
572 if (x > y)
573 {
574 d = x - y;
575 if (d > G_PI + EPS ||
576 (d > G_PI - EPS && tie))
577 y += 2 * G_PI;
578 }
579 else
580 {
581 d = y - x;
582 if (d > G_PI + EPS ||
583 (d > G_PI - EPS && tie))
584 x += 2 * G_PI;
585 }
586 /* unless we are supposed to avoid crossing */
587 if (cross)
588 {
589 if (lastx > x)
590 {
591 if (lasty < y)
592 y -= 2 * G_PI;
593 }
594 else
595 {
596 if (lasty > y)
597 y += 2 * G_PI;
598 }
599 }
600 else
601 {
602 lastx = x;
603 lasty = y;
604 }
605
606 *v3 = s * x + t * y;
607 }
608
609 static void
interpolate_complex(double t,double s,double * r1,double * r2,double * r3,int flip,int tie,int cross)610 interpolate_complex (double t,
611 double s,
612 double *r1,
613 double *r2,
614 double *r3,
615 int flip,
616 int tie,
617 int cross)
618 {
619 double c1[2], c2[2], c3[2];
620 double a1, a2, a3, d1, d2, d3;
621
622 c1[0] = r1[0];
623 c1[1] = r1[1];
624 c2[0] = r2[0];
625 c2[1] = r2[1];
626 if (flip)
627 {
628 double t = c1[0];
629 c1[0] = c1[1];
630 c1[1] = t;
631 t = c2[0];
632 c2[0] = c2[1];
633 c2[1] = t;
634 }
635
636 /* convert to log space */
637 a1 = atan2 (c1[1], c1[0]);
638 a2 = atan2 (c2[1], c2[0]);
639 d1 = 0.5 * log (c1[0] * c1[0] + c1[1] * c1[1]);
640 d2 = 0.5 * log (c2[0] * c2[0] + c2[1] * c2[1]);
641
642 /* interpolate linearly */
643 interpolate_angle (t, s, &a1, &a2, &a3, tie, cross);
644 d3 = s * d1 + t * d2;
645
646 /* convert back */
647 d3 = exp (d3);
648 c3[0] = cos (a3) * d3;
649 c3[1] = sin (a3) * d3;
650
651 if (flip)
652 {
653 r3[1] = c3[0];
654 r3[0] = c3[1];
655 }
656 else
657 {
658 r3[0] = c3[0];
659 r3[1] = c3[1];
660 }
661 }
662
663 static void
interpolate_matrix(double t,double m1[3][2],double m2[3][2],double m3[3][2])664 interpolate_matrix (double t,
665 double m1[3][2],
666 double m2[3][2],
667 double m3[3][2])
668 {
669 double s = 1.0 - t;
670
671 interpolate_complex (t, s, &m1[0][0], &m2[0][0], &m3[0][0], 0, 0, 0);
672 interpolate_complex (t, s, &m1[1][0], &m2[1][0], &m3[1][0], 1, 1, 0);
673
674 /* handle the translation part of the xform linearly */
675 m3[2][0] = s * m1[2][0] + t * m2[2][0];
676 m3[2][1] = s * m1[2][1] + t * m2[2][1];
677 }
678
679 #define INTERP(x) result->x = c0 * cps[i1].x + c1 * cps[i2].x
680
681 /*
682 * create a control point that interpolates between the control points
683 * passed in CPS. for now just do linear. in the future, add control
684 * point types and other things to the cps. CPS must be sorted by time.
685 */
686 void
interpolate(control_point cps[],int ncps,double time,control_point * result)687 interpolate (control_point cps[],
688 int ncps,
689 double time,
690 control_point *result)
691 {
692 int i, j, i1, i2;
693 double c0, c1, t;
694
695 if (ncps == 1)
696 {
697 *result = cps[0];
698 return;
699 }
700 if (cps[0].time >= time)
701 {
702 i1 = 0;
703 i2 = 1;
704 }
705 else if (cps[ncps - 1].time <= time)
706 {
707 i1 = ncps - 2;
708 i2 = ncps - 1;
709 }
710 else
711 {
712 i1 = 0;
713 while (cps[i1].time < time)
714 i1++;
715 i1--;
716 i2 = i1 + 1;
717 if (time - cps[i1].time > -1e-7 &&
718 time - cps[i1].time < 1e-7)
719 {
720 *result = cps[i1];
721 return;
722 }
723 }
724
725 c0 = (cps[i2].time - time) / (cps[i2].time - cps[i1].time);
726 c1 = 1.0 - c0;
727
728 result->time = time;
729
730 if (cps[i1].cmap_inter)
731 {
732 for (i = 0; i < 256; i++)
733 {
734 double spread = 0.15;
735 double d0, d1, e0, e1, c = 2 * G_PI * i / 256.0;
736 c = cos(c * cps[i1].cmap_inter) + 4.0 * c1 - 2.0;
737 if (c > spread) c = spread;
738 if (c < -spread) c = -spread;
739 d1 = (c + spread) * 0.5 / spread;
740 d0 = 1.0 - d1;
741 e0 = (d0 < 0.5) ? (d0 * 2) : (d1 * 2);
742 e1 = 1.0 - e0;
743 for (j = 0; j < 3; j++)
744 {
745 result->cmap[i][j] = (d0 * cps[i1].cmap[i][j] +
746 d1 * cps[i2].cmap[i][j]);
747 #define bright_peak 2.0
748 result->cmap[i][j] = (e1 * result->cmap[i][j] +
749 e0 * 1.0);
750 }
751 }
752 }
753 else
754 {
755 for (i = 0; i < 256; i++)
756 {
757 double t[3], s[3];
758 rgb2hsv (cps[i1].cmap[i], s);
759 rgb2hsv (cps[i2].cmap[i], t);
760 for (j = 0; j < 3; j++)
761 t[j] = c0 * s[j] + c1 * t[j];
762 hsv2rgb (t, result->cmap[i]);
763 }
764 }
765
766 result->cmap_index = -1;
767 INTERP(brightness);
768 INTERP(contrast);
769 INTERP(gamma);
770 INTERP(width);
771 INTERP(height);
772 INTERP(spatial_oversample);
773 INTERP(center[0]);
774 INTERP(center[1]);
775 INTERP(pixels_per_unit);
776 INTERP(spatial_filter_radius);
777 INTERP(sample_density);
778 INTERP(zoom);
779 INTERP(nbatches);
780 INTERP(white_level);
781 for (i = 0; i < 2; i++)
782 for (j = 0; j < 2; j++)
783 {
784 INTERP(pulse[i][j]);
785 INTERP(wiggle[i][j]);
786 }
787
788 for (i = 0; i < NXFORMS; i++)
789 {
790 double r;
791 double rh_time;
792
793 INTERP(xform[i].density);
794 if (result->xform[i].density > 0)
795 result->xform[i].density = 1.0;
796 INTERP(xform[i].color);
797 for (j = 0; j < NVARS; j++)
798 INTERP(xform[i].var[j]);
799 t = 0.0;
800 for (j = 0; j < NVARS; j++)
801 t += result->xform[i].var[j];
802 t = 1.0 / t;
803 for (j = 0; j < NVARS; j++)
804 result->xform[i].var[j] *= t;
805
806 interpolate_matrix(c1, cps[i1].xform[i].c, cps[i2].xform[i].c,
807 result->xform[i].c);
808
809 rh_time = time * 2 * G_PI / (60.0 * 30.0);
810
811 /* apply pulse factor. */
812 r = 1.0;
813 for (j = 0; j < 2; j++)
814 r += result->pulse[j][0] * sin(result->pulse[j][1] * rh_time);
815 for (j = 0; j < 3; j++)
816 {
817 result->xform[i].c[j][0] *= r;
818 result->xform[i].c[j][1] *= r;
819 }
820
821 /* apply wiggle factor */
822 for (j = 0; j < 2; j++)
823 {
824 double tt = result->wiggle[j][1] * rh_time;
825 double m = result->wiggle[j][0];
826 result->xform[i].c[0][0] += m * cos(tt);
827 result->xform[i].c[1][0] += m * -sin(tt);
828 result->xform[i].c[0][1] += m * sin(tt);
829 result->xform[i].c[1][1] += m * cos(tt);
830 }
831 } /* for i */
832 }
833
834
835
836 /*
837 * split a string passed in ss into tokens on whitespace.
838 * # comments to end of line. ; terminates the record
839 */
840 void
tokenize(char ** ss,char * argv[],int * argc)841 tokenize (char **ss,
842 char *argv[],
843 int *argc)
844 {
845 char *s = *ss;
846 int i = 0, state = 0;
847
848 while (*s != ';')
849 {
850 char c = *s;
851 switch (state)
852 {
853 case 0:
854 if (c == '#')
855 state = 2;
856 else if (!g_ascii_isspace (c))
857 {
858 argv[i] = s;
859 i++;
860 state = 1;
861 }
862 case 1:
863 if (g_ascii_isspace (c))
864 {
865 *s = 0;
866 state = 0;
867 }
868 case 2:
869 if (c == '\n')
870 state = 0;
871 }
872 s++;
873 }
874 *s = 0;
875 *ss = s + 1;
876 *argc = i;
877 }
878
879 static int
compare_xforms(const void * va,const void * vb)880 compare_xforms (const void *va,
881 const void *vb)
882 {
883 double aa[2][2];
884 double bb[2][2];
885 double ad, bd;
886 const xform *a = va;
887 const xform *b = vb;
888
889 aa[0][0] = a->c[0][0];
890 aa[0][1] = a->c[0][1];
891 aa[1][0] = a->c[1][0];
892 aa[1][1] = a->c[1][1];
893 bb[0][0] = b->c[0][0];
894 bb[0][1] = b->c[0][1];
895 bb[1][0] = b->c[1][0];
896 bb[1][1] = b->c[1][1];
897 ad = det_matrix (aa);
898 bd = det_matrix (bb);
899 if (ad < bd)
900 return -1;
901 if (ad > bd)
902 return 1;
903 return 0;
904 }
905
906 #define MAXARGS 1000
907 #define streql(x,y) (!strcmp(x,y))
908
909 /*
910 * given a pointer to a string SS, fill fields of a control point CP.
911 * return a pointer to the first unused char in SS. totally barfucious,
912 * must integrate with tcl soon...
913 */
914
915 void
parse_control_point(char ** ss,control_point * cp)916 parse_control_point (char **ss,
917 control_point *cp)
918 {
919 char *argv[MAXARGS];
920 int argc, i, j;
921 int set_cm = 0, set_image_size = 0, set_nbatches = 0, set_white_level = 0, set_cmap_inter = 0;
922 int set_spatial_oversample = 0;
923 double *slot = NULL, xf, cm, t, nbatches, white_level, spatial_oversample, cmap_inter;
924 double image_size[2];
925
926 for (i = 0; i < NXFORMS; i++)
927 {
928 cp->xform[i].density = 0.0;
929 cp->xform[i].color = (i == 0);
930 cp->xform[i].var[0] = 1.0;
931 for (j = 1; j < NVARS; j++)
932 cp->xform[i].var[j] = 0.0;
933 cp->xform[i].c[0][0] = 1.0;
934 cp->xform[i].c[0][1] = 0.0;
935 cp->xform[i].c[1][0] = 0.0;
936 cp->xform[i].c[1][1] = 1.0;
937 cp->xform[i].c[2][0] = 0.0;
938 cp->xform[i].c[2][1] = 0.0;
939 }
940 for (j = 0; j < 2; j++)
941 {
942 cp->pulse[j][0] = 0.0;
943 cp->pulse[j][1] = 60.0;
944 cp->wiggle[j][0] = 0.0;
945 cp->wiggle[j][1] = 60.0;
946 }
947
948 tokenize (ss, argv, &argc);
949 for (i = 0; i < argc; i++)
950 {
951 if (streql("xform", argv[i]))
952 slot = &xf;
953 else if (streql("time", argv[i]))
954 slot = &cp->time;
955 else if (streql("brightness", argv[i]))
956 slot = &cp->brightness;
957 else if (streql("contrast", argv[i]))
958 slot = &cp->contrast;
959 else if (streql("gamma", argv[i]))
960 slot = &cp->gamma;
961 else if (streql("zoom", argv[i]))
962 slot = &cp->zoom;
963 else if (streql("image_size", argv[i]))
964 {
965 slot = image_size;
966 set_image_size = 1;
967 }
968 else if (streql("center", argv[i]))
969 slot = cp->center;
970 else if (streql("pulse", argv[i]))
971 slot = (double *) cp->pulse;
972 else if (streql("wiggle", argv[i]))
973 slot = (double *) cp->wiggle;
974 else if (streql("pixels_per_unit", argv[i]))
975 slot = &cp->pixels_per_unit;
976 else if (streql("spatial_filter_radius", argv[i]))
977 slot = &cp->spatial_filter_radius;
978 else if (streql("sample_density", argv[i]))
979 slot = &cp->sample_density;
980 else if (streql("nbatches", argv[i]))
981 {
982 slot = &nbatches;
983 set_nbatches = 1;
984 }
985 else if (streql("white_level", argv[i]))
986 {
987 slot = &white_level;
988 set_white_level = 1;
989 }
990 else if (streql("spatial_oversample", argv[i]))
991 {
992 slot = &spatial_oversample;
993 set_spatial_oversample = 1;
994 }
995 else if (streql("cmap", argv[i]))
996 {
997 slot = &cm;
998 set_cm = 1;
999 }
1000 else if (streql("density", argv[i]))
1001 slot = &cp->xform[(int)xf].density;
1002 else if (streql("color", argv[i]))
1003 slot = &cp->xform[(int)xf].color;
1004 else if (streql("coefs", argv[i]))
1005 {
1006 slot = cp->xform[(int)xf].c[0];
1007 cp->xform[(int)xf].density = 1.0;
1008 }
1009 else if (streql("var", argv[i]))
1010 slot = cp->xform[(int)xf].var;
1011 else if (streql("cmap_inter", argv[i]))
1012 {
1013 slot = &cmap_inter;
1014 set_cmap_inter = 1;
1015 }
1016 else
1017 *slot++ = g_strtod(argv[i], NULL);
1018 }
1019 if (set_cm)
1020 {
1021 cp->cmap_index = (int) cm;
1022 get_cmap(cp->cmap_index, cp->cmap, 256);
1023 }
1024 if (set_image_size)
1025 {
1026 cp->width = (int) image_size[0];
1027 cp->height = (int) image_size[1];
1028 }
1029 if (set_cmap_inter)
1030 cp->cmap_inter = (int) cmap_inter;
1031 if (set_nbatches)
1032 cp->nbatches = (int) nbatches;
1033 if (set_spatial_oversample)
1034 cp->spatial_oversample = (int) spatial_oversample;
1035 if (set_white_level)
1036 cp->white_level = (int) white_level;
1037 for (i = 0; i < NXFORMS; i++)
1038 {
1039 t = 0.0;
1040 for (j = 0; j < NVARS; j++)
1041 t += cp->xform[i].var[j];
1042 t = 1.0 / t;
1043 for (j = 0; j < NVARS; j++)
1044 cp->xform[i].var[j] *= t;
1045 }
1046 qsort ((char *) cp->xform, NXFORMS, sizeof(xform), compare_xforms);
1047 }
1048
1049 void
print_control_point(FILE * f,control_point * cp,int quote)1050 print_control_point (FILE *f,
1051 control_point *cp,
1052 int quote)
1053 {
1054 int i, j;
1055 char *q = quote ? "# " : "";
1056 fprintf (f, "%stime %g\n", q, cp->time);
1057 if (cp->cmap_index != -1)
1058 fprintf (f, "%scmap %d\n", q, cp->cmap_index);
1059 fprintf (f, "%simage_size %d %d center %g %g pixels_per_unit %g\n",
1060 q, cp->width, cp->height, cp->center[0], cp->center[1],
1061 cp->pixels_per_unit);
1062 fprintf (f, "%sspatial_oversample %d spatial_filter_radius %g",
1063 q, cp->spatial_oversample, cp->spatial_filter_radius);
1064 fprintf (f, " sample_density %g\n", cp->sample_density);
1065 fprintf (f, "%snbatches %d white_level %d\n",
1066 q, cp->nbatches, cp->white_level);
1067 fprintf (f, "%sbrightness %g gamma %g cmap_inter %d\n",
1068 q, cp->brightness, cp->gamma, cp->cmap_inter);
1069
1070 for (i = 0; i < NXFORMS; i++)
1071 if (cp->xform[i].density > 0.0)
1072 {
1073 fprintf (f, "%sxform %d density %g color %g\n",
1074 q, i, cp->xform[i].density, cp->xform[i].color);
1075 fprintf (f, "%svar", q);
1076 for (j = 0; j < NVARS; j++)
1077 fprintf (f, " %g", cp->xform[i].var[j]);
1078 fprintf (f, "\n%scoefs", q);
1079 for (j = 0; j < 3; j++)
1080 fprintf (f, " %g %g", cp->xform[i].c[j][0], cp->xform[i].c[j][1]);
1081 fprintf (f, "\n");
1082 }
1083 fprintf (f, "%s;\n", q);
1084 }
1085
1086 /* returns a uniform variable from 0 to 1 */
1087 double
random_uniform01(void)1088 random_uniform01 (void)
1089 {
1090 return g_random_double ();
1091 }
1092
random_uniform11(void)1093 double random_uniform11 (void)
1094 {
1095 return g_random_double_range (-1, 1);
1096 }
1097
1098 /* returns a mean 0 variance 1 random variable
1099 see numerical recipes p 217 */
random_gaussian(void)1100 double random_gaussian(void)
1101 {
1102 static int iset = 0;
1103 static double gset;
1104 double fac, r, v1, v2;
1105
1106 if (iset == 0)
1107 {
1108 do
1109 {
1110 v1 = random_uniform11 ();
1111 v2 = random_uniform11 ();
1112 r = v1 * v1 + v2 * v2;
1113 }
1114 while (r >= 1.0 || r == 0.0);
1115 fac = sqrt (-2.0 * log (r) / r);
1116 gset = v1 * fac;
1117 iset = 1;
1118 return v2 * fac;
1119 }
1120 iset = 0;
1121 return gset;
1122 }
1123
1124 void
copy_variation(control_point * cp0,control_point * cp1)1125 copy_variation (control_point *cp0,
1126 control_point *cp1)
1127 {
1128 int i, j;
1129 for (i = 0; i < NXFORMS; i++)
1130 {
1131 for (j = 0; j < NVARS; j++)
1132 cp0->xform[i].var[j] = cp1->xform[i].var[j];
1133 }
1134 }
1135
1136 #define random_distrib(v) ((v)[g_random_int_range (0, vlen(v))])
1137
1138 void
random_control_point(control_point * cp,int ivar)1139 random_control_point (control_point *cp,
1140 int ivar)
1141 {
1142 int i, nxforms, var;
1143 static int xform_distrib[] =
1144 {
1145 2, 2, 2,
1146 3, 3, 3,
1147 4, 4,
1148 5
1149 };
1150 static int var_distrib[] =
1151 {
1152 -1, -1, -1,
1153 0, 0, 0, 0,
1154 1, 1, 1,
1155 2, 2, 2,
1156 3, 3,
1157 4, 4,
1158 5
1159 };
1160
1161 static int mixed_var_distrib[] =
1162 {
1163 0, 0, 0,
1164 1, 1, 1,
1165 2, 2, 2,
1166 3, 3,
1167 4, 4,
1168 5, 5
1169 };
1170
1171 get_cmap (cmap_random, cp->cmap, 256);
1172 cp->time = 0.0;
1173 nxforms = random_distrib (xform_distrib);
1174 var = (0 > ivar) ?
1175 random_distrib(var_distrib) :
1176 ivar;
1177 for (i = 0; i < nxforms; i++)
1178 {
1179 int j, k;
1180 cp->xform[i].density = 1.0 / nxforms;
1181 cp->xform[i].color = i == 0;
1182 for (j = 0; j < 3; j++)
1183 for (k = 0; k < 2; k++)
1184 cp->xform[i].c[j][k] = random_uniform11();
1185 for (j = 0; j < NVARS; j++)
1186 cp->xform[i].var[j] = 0.0;
1187 if (var >= 0)
1188 cp->xform[i].var[var] = 1.0;
1189 else
1190 cp->xform[i].var[random_distrib(mixed_var_distrib)] = 1.0;
1191
1192 }
1193 for (; i < NXFORMS; i++)
1194 cp->xform[i].density = 0.0;
1195 }
1196
1197 /*
1198 * find a 2d bounding box that does not enclose eps of the fractal density
1199 * in each compass direction. works by binary search.
1200 * this is stupid, it should just use the find nth smallest algorithm.
1201 */
1202 void
estimate_bounding_box(control_point * cp,double eps,double * bmin,double * bmax)1203 estimate_bounding_box (control_point *cp,
1204 double eps,
1205 double *bmin,
1206 double *bmax)
1207 {
1208 int i, j, batch = (eps == 0.0) ? 10000 : 10.0/eps;
1209 int low_target = batch * eps;
1210 int high_target = batch - low_target;
1211 point min, max, delta;
1212 point *points = malloc (sizeof (point) * batch);
1213 iterate (cp, batch, 20, points);
1214
1215 min[0] = min[1] = 1e10;
1216 max[0] = max[1] = -1e10;
1217
1218 for (i = 0; i < batch; i++)
1219 {
1220 if (points[i][0] < min[0]) min[0] = points[i][0];
1221 if (points[i][1] < min[1]) min[1] = points[i][1];
1222 if (points[i][0] > max[0]) max[0] = points[i][0];
1223 if (points[i][1] > max[1]) max[1] = points[i][1];
1224 }
1225
1226 if (low_target == 0)
1227 {
1228 bmin[0] = min[0];
1229 bmin[1] = min[1];
1230 bmax[0] = max[0];
1231 bmax[1] = max[1];
1232 return;
1233 }
1234
1235 delta[0] = (max[0] - min[0]) * 0.25;
1236 delta[1] = (max[1] - min[1]) * 0.25;
1237
1238 bmax[0] = bmin[0] = min[0] + 2.0 * delta[0];
1239 bmax[1] = bmin[1] = min[1] + 2.0 * delta[1];
1240
1241 for (i = 0; i < 14; i++)
1242 {
1243 int n, s, e, w;
1244 n = s = e = w = 0;
1245 for (j = 0; j < batch; j++)
1246 {
1247 if (points[j][0] < bmin[0]) n++;
1248 if (points[j][0] > bmax[0]) s++;
1249 if (points[j][1] < bmin[1]) w++;
1250 if (points[j][1] > bmax[1]) e++;
1251 }
1252 bmin[0] += (n < low_target) ? delta[0] : -delta[0];
1253 bmax[0] += (s < high_target) ? delta[0] : -delta[0];
1254 bmin[1] += (w < low_target) ? delta[1] : -delta[1];
1255 bmax[1] += (e < high_target) ? delta[1] : -delta[1];
1256 delta[0] = delta[0] / 2.0;
1257 delta[1] = delta[1] / 2.0;
1258 }
1259 }
1260
1261 /* this has serious flaws in it */
1262
1263 double
standard_metric(control_point * cp1,control_point * cp2)1264 standard_metric (control_point *cp1,
1265 control_point *cp2)
1266 {
1267 int i, j, k;
1268 double t;
1269 double dist = 0.0;
1270
1271 for (i = 0; i < NXFORMS; i++)
1272 {
1273 double var_dist = 0.0;
1274 double coef_dist = 0.0;
1275 for (j = 0; j < NVARS; j++)
1276 {
1277 t = cp1->xform[i].var[j] - cp2->xform[i].var[j];
1278 var_dist += t * t;
1279 }
1280 for (j = 0; j < 3; j++)
1281 for (k = 0; k < 2; k++)
1282 {
1283 t = cp1->xform[i].c[j][k] - cp2->xform[i].c[j][k];
1284 coef_dist += t *t;
1285 }
1286 /* weight them equally for now. */
1287 dist += var_dist + coef_dist;
1288 }
1289 return dist;
1290 }
1291
1292 static int
flam3_random_bit(void)1293 flam3_random_bit (void)
1294 {
1295 static int n = 0;
1296 static int l;
1297 if (n == 0)
1298 {
1299 l = g_random_int ();
1300 n = 20;
1301 }
1302 else
1303 {
1304 l = l >> 1;
1305 n--;
1306 }
1307 return l & 1;
1308 }
1309
1310 static double
flam3_random01(void)1311 flam3_random01 (void)
1312 {
1313 return (g_random_int () & 0xfffffff) / (double) 0xfffffff;
1314 }
1315