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