← Back to branch summary

~ubuntu-branches/ubuntu/hoary/gimp/hoary

« back to all changes in this revision

Viewing changes to plug-ins/flame/libifs.c

  • Committer: Bazaar Package Importer
  • Author(s): Sebastien Bacher
  • Date: 2005-04-04 14:51:23 UTC
  • Revision ID: james.westby@ubuntu.com-20050404145123-9py049eeelfymur8
Tags: upstream-2.2.2
Import upstream version 2.2.2

Show diffs side-by-side

added added

removed removed

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