home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / grafik / vifs / ifsshow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-08-22  |  11.3 KB  |  461 lines
  1. /*
  2.  *    ifsshow.c  --  show an ifs.
  3.  *
  4.  *    4 june 1989  Olle Olsson
  5.  */
  6.  
  7. #include <stdlib.h>
  8. #include <conio.h>
  9. #include <math.h>
  10. #include "ifs.h"
  11.  
  12. /* *  /
  13. #define PRONLY 1        /* for test */
  14. /* */
  15.  
  16. #define FSCL    0x10000        /* fixed point scale */
  17.  
  18. typedef union
  19.     {
  20.     long    l;        /* fixed point long */
  21.     struct
  22.         {
  23.         int    lw;    /* low word, fraction part */
  24.         int    hw;    /* high word, integer part */
  25.         } w;
  26.     } fpl;
  27.  
  28. #define FINT( x )    (((fpl *)&x)->w.hw)    /* fixed point integer part */
  29.  
  30. typedef struct            /* fixed point version */
  31.     {
  32.     long a11, a12;        /* one transformation matrix */
  33.     long a21, a22;
  34.     long b1, b2;        /* translation vector */
  35.     int color;        /* color index */
  36.     double prob;        /* probability */
  37.     int iprob;        /* probability fitted to rand() */
  38.     int group;        /* group index */
  39.     } ftrans;
  40.  
  41. typedef struct            /* group member spec. */
  42.     {
  43.     int nels;        /* number of elements */
  44.     long tcount;        /* total number of iterations to do */
  45.     long dcount;        /* delta count (for mixing of groups) */
  46.     double frcnt;        /* fraction part of dcount */
  47.     double sumfrcnt;    /* accumulator for fraction part of dcount */
  48.     long x, y;        /* current point */
  49.     int color;        /* current color */
  50.     } gr_el;
  51.     
  52.  
  53. /* local routines */
  54. static void fixptr( ifsdes *dp, ftrans *tp0 );
  55. static int gcomp( const void *tr1, const void *tr2 );
  56. static int no_stop( void );
  57. static void far no_wpix( int x, int y, int color );
  58. static unsigned far no_rpix( int x, int y );
  59. static long dtl( double d );
  60. static double ltd( long l );
  61. static void probint( ftrans ftp[], int size, gr_el groups[], int ngroups,
  62.             long gt_count, long delta_count );
  63. static int riterate( ftrans *ftp, gr_el *gp, int im, long count, long *ecountp,
  64.             int (*stop)( void ),
  65.             void (far *wpix)( int x, int y, int color ),
  66.             unsigned (far *rpix)( int x, int y )
  67.             );
  68.  
  69. /* C or assembly routine */
  70. /*#define mul4f( a, b )     cmul4f( a, b ); /* */
  71. static long cmul4f( long a, long b );
  72.  
  73.  
  74. void ifsshow( dp, stop, estop, wpix, rpix )
  75. ifsdes *dp;            /* the transformations */
  76. int (*stop)( void );        /* stop function */
  77. int estop;            /* early miss stop flag */
  78. void (far *wpix)( int x, int y, int color );    /* write pixel */
  79. unsigned (far *rpix)( int x, int y );        /* read pixel, returns color */
  80. {
  81. int i;
  82. register int j;
  83. int size = dp -> size;        /* number of transforms */
  84. ftrans *tp0;            /* fixed point base pointer */
  85. register ftrans *tp;        /* tmp */
  86. gr_el *grps;            /* group specifier */
  87. gr_el *gp;            /* tmp */
  88. int ngrps = dp -> maxgroup + 1;    /* number of groups */
  89. int running;            /* flag */
  90. long count;            /* iteration count */
  91. long mcnt;            /* miss count */
  92. long dcnt;            /* tmp */
  93.  
  94. /* normalize probabilities */
  95. normprob( dp -> tp, size );
  96.  
  97. /* allocate a fixed point matrix array */
  98. if (!(tp0 = (ftrans *) malloc( size * sizeof (ftrans) + 100 )))
  99.     error( "Out of memory for transformations (ifsshow)" );
  100.  
  101. /* allocate group specifiers */
  102. if (!(grps = (gr_el *) malloc( ngrps * sizeof (gr_el) + 100 )))
  103.     error( "Out of memory for groups (ifsshow)" );
  104.  
  105. /* recompute the transforms to fixed point */
  106. fixptr( dp, tp0 );
  107.  
  108. /* sort groups */
  109. qsort( tp0, size, sizeof (ftrans), gcomp );
  110.  
  111. /* find the total number of iterations */
  112. count = DENS2ITER( dp -> density );
  113.  
  114. /* set the group probabilities and such things */
  115. probint( tp0, size, grps, ngrps, count, 100L );
  116.  
  117. /* the fifty initial iterations are not shown */
  118. for (gp = grps, i = 0, tp = tp0; (i < ngrps) && (tp < &tp0[size]); ++gp, ++i)
  119.     {
  120.     /* skip if it is an empty group */
  121.     if (gp -> nels <= 0)
  122.         continue;
  123.  
  124.     /* set current point & set color to background */
  125.     gp -> x = gp -> y = gp -> color = 0;
  126.  
  127.     /* iterate */
  128.     (void) riterate( tp, gp, dp -> im, 50L, &mcnt,
  129.                         no_stop, no_wpix, no_rpix );
  130.  
  131.         /* next group */
  132.     tp += gp -> nels;
  133.     }
  134.  
  135. /* max missing count is controlled by the early stop parameter */
  136. mcnt = estop ? (count / 5) : count;
  137.  
  138. /* now show the fractals */
  139. for (running = 1; running && (count > 0);)
  140. for (gp = grps, i = 0, tp = tp0; (i < ngrps) && (tp < &tp0[size]); ++gp, ++i)
  141.     {
  142.     /* skip if it is an empty group */
  143.     if (gp -> nels <= 0)
  144.         continue;
  145.  
  146.     /* how many iterations for this group? */
  147.     /* accumulate fraction count (for groups with very low prob.) */
  148.     j = (gp -> sumfrcnt += gp -> frcnt);
  149.     gp -> sumfrcnt -= j;
  150.     count -= (dcnt = gp -> dcount + j);
  151.  
  152.     /* iterate */
  153.     if (!riterate( tp, gp, dp -> im, dcnt, &mcnt, stop, wpix, rpix ))
  154.         {
  155.         /* the stop() functions has sayed stop, so stop! */
  156.         running = 0;
  157.         break;
  158.         }
  159.     else if (mcnt < 0)
  160.         {
  161.         /* to many misses stop! */
  162.         running = 0;
  163.         break;
  164.         }
  165.  
  166.         /* next group */
  167.     tp += gp -> nels;
  168.     }
  169.  
  170. /* remember the density */
  171. dp -> prrdens = dp -> density - ITER2DENS( count );
  172.  
  173. /* free() the vectors */
  174. free( tp0 );
  175. free( grps );
  176. }
  177.  
  178. /* local routines */
  179.  
  180. static void fixptr( dp, tp0 )
  181. ifsdes *dp;            /* the transformations */
  182. ftrans *tp0;            /* fixed point base pointer */
  183. {
  184. int i;                /* tmp */
  185. area *ap = &dp -> area;        /* the area */
  186. register transform *dtp;    /* tmp */
  187. register ftrans *tp;        /* tmp */
  188. double sx, sy;            /* tmp */
  189.  
  190. /* recompute the transformation matrix to fixed point, scaled to the screen */
  191. /*
  192.  * One transformation is:  Ax + b, x is the position,
  193.  * g(x) takes this to screen coordinates.
  194.  * Recompute the transformation to screen coordinates!
  195.  * I.e. solve H and k from: g(Ax + b) = Hg(x) + k.
  196.  * g(x) is the ap parameter and the screen size combined.
  197.  */
  198. if (fabs( ap -> xlen ) > 1e-20)
  199.     sx = MAX_COL / ap -> xlen;
  200. else
  201.     sx = 1e20;
  202.  
  203. if (fabs( ap -> ylen ) > 1e-20)
  204.     sy = MAX_ROW / ap -> ylen;
  205. else
  206.     sy = 1e20;
  207.  
  208. for (i = dp -> size, tp = tp0, dtp = dp -> tp; i; --i, ++tp, ++dtp)
  209.     {
  210.     double z;
  211.  
  212.     tp -> a11 = dtl( dtp -> a11 );
  213.     tp -> a12 = dtl( dtp -> a12 * sx / sy );
  214.     tp -> a21 = dtl( dtp -> a21 * sy / sx );
  215.     tp -> a22 = dtl( dtp -> a22 );
  216.     z = dtp -> b1 + ap -> xlow * (dtp -> a11 - 1) + dtp -> a12 * ap -> ylow;
  217.     tp -> b1  = dtl( sx * z );
  218.     z = dtp -> b2 + ap -> ylow * (dtp -> a22 - 1) + dtp -> a21 * ap -> xlow;
  219.     tp -> b2  = dtl( sy * z );
  220.  
  221.     /* copy the probability, the color and the group */
  222.     tp -> prob = dtp -> prob;
  223.     tp -> color = dtp -> color;
  224.     tp -> group = dtp -> group;
  225.     }
  226. }
  227.  
  228. static int gcomp( const void *tr1, const void *tr2 )
  229. {
  230. /* group comparing function for sorting */
  231. return (((ftrans *)tr1) -> group - ((ftrans *)tr2) -> group);
  232. }
  233.  
  234. static int no_stop( void )
  235. {
  236. /* fake function for the initial iterations */
  237. return (0);
  238. }
  239.  
  240. static void far no_wpix( int x, int y, int color )
  241. {
  242. /* fake function for the initial iterations */
  243. /* fool the compiler to not warn about unused arguments */
  244. x = y; y = color; color = x;
  245. return;
  246. }
  247.  
  248. static unsigned far no_rpix( int x, int y )
  249. {
  250. /* fake function for the initial iterations */
  251. /* fool the compiler to not warn about unused arguments */
  252. x = y; y = x;
  253. return (0);
  254. }
  255.  
  256.  
  257. static int riterate( ftp, gp, im, count, ecountp, stop, wpix, rpix )
  258. ftrans *ftp;            /* the transformations */
  259. gr_el *gp;            /* current group */
  260. int im;                /* invariant measure mode flag */
  261. long count;            /* number of iterations */
  262. long *ecountp;            /* error max count */
  263. int (*stop)( void );        /* stop function */
  264. void (far *wpix)( int x, int y, int color );    /* write pixel */
  265. unsigned (far *rpix)( int x, int y );        /* read pixel, returns color */
  266. {
  267. register ftrans *tp;        /* current transform */
  268. register int j;            /* tmp */
  269. int r;
  270. int size1;            /* one less than the number of transforms */
  271. long x, y;            /* current point */
  272. long xt;
  273. unsigned int ix, iy;
  274. int prevclr;            /* previous color */
  275.  
  276. /* get number of transforms - 1, in this group */
  277. size1 = gp -> nels - 1;
  278.  
  279. /* pick up where we left off */
  280. x = gp -> x;
  281. y = gp -> y;
  282. prevclr = gp -> color;
  283.  
  284. for (; count > 0; --count)
  285.     {
  286.     if ((*stop)())
  287.         {
  288.         /* remember the stop and then stop */
  289.         count = -1;
  290.         break;
  291.         }
  292.  
  293.     /* select */
  294.     r = rand();
  295.  
  296.     /* use size1 to not fall out of the group */
  297.     for (j = 0, tp = ftp; j < size1; ++j, ++tp)
  298.         if (r < tp -> iprob) break;
  299.  
  300.     /* transform */
  301.     /*xt = tp -> a11 * x + tp -> a12 * y + tp -> b1;*/
  302.     xt = cmul4f( tp -> a11, x ) + cmul4f( tp -> a12, y ) + tp -> b1;
  303.  
  304.     /*y  = tp -> a21 * x + tp -> a22 * y + tp -> b2;*/
  305.     y  = cmul4f( tp -> a21, x ) + cmul4f( tp -> a22, y ) + tp -> b2;
  306.  
  307.     x = xt;
  308.  
  309.     /* show */
  310. #ifdef PRONLY
  311.     printf( "%d: tp0 %d tp %d delta %d ", j, tp0, tp, tp - tp0 );
  312.     printf( "x %ld, y %ld  (%d,%d)\n", x, y, FINT( x ), FINT( y ) );
  313.  
  314. /*    printf( "%10ld%10ld", x, y );*/
  315.     continue;
  316. #endif
  317.     ix = FINT( x );
  318.     iy = FINT( y );
  319.  
  320.     if (ix >= MAX_COL || iy >= MAX_ROW)
  321.         {
  322.         /* decrease the error max */
  323.         --(*ecountp);
  324.         continue;
  325.         }
  326.  
  327.     if (im)
  328.         {
  329.         /* invariant measure mode */
  330.         /* check if max color index is reached */
  331.         if ((j = (*rpix)( ix, MAX_ROW - iy )) >= im)
  332.             continue;
  333.         else ++j;
  334.         }
  335.     else if ((j = tp -> color) < 0)        /* -1 means prev color */
  336.         j = prevclr;
  337.  
  338.     /* show it */
  339.     wpix( ix, MAX_ROW - iy, prevclr = j );
  340.     }
  341.  
  342. /* remember the status */
  343. gp -> x = x;
  344. gp -> y = y;
  345. gp -> color = prevclr;
  346.  
  347. /* return 0 if (*stop)() or 1 if ok */
  348. return (count == 0);
  349. }
  350.  
  351. static long dtl( d )
  352. double d;
  353. {
  354. /* convert to fixed point */
  355. /*printf( " %lf --> %ld ", d, (long)(d*FSCL) );*/
  356. return ((long)(d*FSCL));
  357. }
  358.  
  359. static double ltd( l )
  360. long l;
  361. {
  362. /* convert from fixed point */
  363. /*printf( " %ld --> %lf ", d, ((double) l)/FSCL );*/
  364. return ((double)(l/FSCL));
  365. }
  366.  
  367. /* Test versions */
  368. static long cmul4f( a, b )
  369. long a, b;
  370. {
  371. double d = a;
  372.  
  373. d *= b;
  374. d /= FSCL;
  375. #ifdef PRONLY
  376. /*printf( " cmul4f:%lf ", d );*/
  377. #endif
  378. return ((long) d);
  379. }
  380.  
  381. /* * /
  382. static long csq4f( a )
  383. long a;
  384. {
  385. return (cmul4f( a, a ));
  386. }
  387. /* */
  388.  
  389. static void probint( ftp, size, grps, ngrps, count, dcount )
  390. ftrans ftp[];            /* the transformations */
  391. int size;            /* number of transformations */
  392. gr_el grps[];            /* the group specs to be filled in */
  393. int ngrps;            /* number of groups */
  394. long count;            /* grand total count to distribute */
  395. long dcount;            /* delta count wanted */
  396. {
  397. register int i, ig;
  398. gr_el *gp;            /* current group */
  399. int prevp;
  400. ftrans *tp;
  401. double v;
  402.  
  403. /* the ftp[] vector is sorted on groups */
  404.  
  405. /* set up the group sizes and sum each groups prob. (temp use of sumfrcnt) */
  406. for (gp = grps, ig = 0, tp = ftp; ig < ngrps; ++gp, ++ig)
  407.     {
  408.     gp -> sumfrcnt = 0;
  409.     gp -> dcount = 0;
  410.     gp -> frcnt = 0;
  411.     gp -> nels = 0;
  412.     gp -> tcount = 0;
  413.  
  414.     /* out of transforms means empty group */
  415.     if (tp) 
  416.         {
  417.         /* count and sum up prob. */
  418.         for (; tp -> group == ig; ++gp -> nels, ++tp)
  419.             if (tp >= &ftp[size])
  420.                 {
  421.                 /* out of transforms! */
  422.                 tp = 0;
  423.                 break;
  424.                 }
  425.             else gp -> sumfrcnt += tp -> prob;
  426.  
  427.         /* set the count for this group */
  428.         gp -> tcount = count * gp -> sumfrcnt;
  429.         v = dcount * gp -> sumfrcnt;
  430.         gp -> dcount = v;            /* integer part */
  431.         gp -> frcnt = v - gp -> dcount;        /* fraction */
  432.  
  433.         /* sumfrcnt should not be zero, see below */
  434.         if (gp -> sumfrcnt <= 0) gp -> sumfrcnt = MINPROB;
  435.         }
  436.     }
  437.  
  438. /* */
  439. /* recalculate the probabilities to suit rand() (integers) */
  440. /* make a sequence withint each group */
  441. for (tp = ftp, prevp = 0, i = 0, ig = tp -> group; i < size; ++tp, ++i)
  442.     {
  443.     /* new group? */
  444.     if (tp -> group != ig)
  445.         {
  446.         prevp = 0;
  447.         ig = tp -> group;
  448.         }
  449.  
  450.     /* set and sum (sumfrcnt is not zero!) */
  451.     prevp = (tp -> iprob = prevp +
  452.             MAXRAND * tp -> prob / grps[ig].sumfrcnt);
  453.     }
  454.  
  455. /* clear sumfrcnt (usage as temp storage is over) */
  456. for (gp = grps, ig = 0; ig < ngrps; ++gp, ++ig)
  457.     gp -> sumfrcnt = 0;
  458. }
  459.  
  460.  
  461.