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