home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 3 / 3406 / mb.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-05-24  |  8.2 KB  |  318 lines
  1.  
  2. /* This module does the most essential part of every mandelbrot
  3.  * program: it iterates points (pixels). I am trying to make this
  4.  * part fast using 16-bit or 32-bit fixed point math whenever possible.
  5.  * Unfortunately, this had to be done in somewhat unportable manner.
  6.  * If your machine has problems with these, #undefine USE_INTEGER_MATH
  7.  * from mb.h.
  8.  */
  9.  
  10. /* And right :) You guessed it; because it is me that is writing this
  11.  * program, the contour crawling algorithm is used, too.
  12.  */
  13.  
  14. /* iterate() is the main callable function of this module:
  15.  
  16. int iterate(x0, y0, dx, dy, cx, cy, nx, ny, iter, flags, dest)
  17. double x0, y0, dx, dy, cx, cy;
  18. int nx, ny;
  19. int iter;
  20. int flags;
  21. u_short **dest;
  22.  
  23.  * Where
  24.  o x0+iy0 is the upper lefthand corner of the rectangular area,
  25.  o dx and idy are deltas in two orthogonal directions,
  26.  o cx+icy is the constant used in iterating Julia sets,
  27.  o (nx, ny) is the size of the rectangle,
  28.  o iter is the maximum number of iterations,
  29.  o flags can be set to JULIA and
  30.  
  31.  o dest is an array of u_shorts that returns the result to the caller.
  32.  o iterate() returns the number of usable u_shorts in the array.
  33.  
  34.  * dest will be allocated by iterate(), but it is callers responsibility
  35.  * to free it. dest has the following format:
  36.  
  37.  o 16 bits for color of the first pixel (as returned by first_point(),
  38.  o 16 bits for the number of entries in the chain of directions and
  39.  o 2 bits for each direction code (see mb.h), rounded to next 16 bits.
  40.  * NOTE: to further reduce the bandwidth needed, if color has highest
  41.  * (bit 15 set), the number of directions is zero (0), that is, the
  42.  * area consists of one (1) pixel.
  43.  */
  44.  
  45. #include <sys/types.h>
  46. #include <math.h>
  47. #include <stdio.h>
  48. #include "mb.h"
  49. #include "pix.h"
  50. #define PACKSIZE 512
  51.  
  52. #ifdef USE_INTEGER_MATH
  53. /* First, the 16-bit version */
  54. #define multi16(x,y) (((x)*(y)) >> BIT16)
  55.  
  56. int mb16(xl, yl, cx, cy, iter)
  57. double xl, yl;
  58. double cx, cy;
  59. int iter;
  60. {
  61. register short x= xl * (1 << BIT16), y= yl * (1 << BIT16);
  62. register short x0= cx * (1 << BIT16), y0= cy * (1 << BIT16);
  63. register short x2, y2;
  64. register short tmp;
  65.  
  66. /* It is painful for assembler programmer such as me to write
  67.  * unefficient code such as this. And, worse still, this code
  68.  * is not portable. (It counts on signed bit-shifts and 2's
  69.  * complement math.
  70.  */
  71. while (iter--) {
  72.   x2= multi16(x,x);
  73.   y2= multi16(y,y);
  74.   tmp= x2 + y2;
  75.   if (x < -(2 << BIT16) || x > (2 << BIT16) ||
  76.       y < -(2 << BIT16) || y > (2 << BIT16) ||
  77.       tmp < 0 || tmp > ((4 << BIT16)-1) )
  78.     return(iter);
  79.   y= multi16(x,y);
  80.   tmp= y + y;
  81.   if ((tmp > 0 && y <= 0) || (tmp < 0 && y >= 0))
  82.     return(iter?iter-1:0);
  83.   y= tmp + y0;
  84.   if ((tmp > 0 && y0 > 0 && y <= 0) || (tmp < 0 && y0 < 0 && y >= 0))
  85.     return(iter?iter-1:0);
  86.   tmp= x2 - y2;
  87.   x= tmp + x0;
  88.   if ((tmp > 0 && x0 > 0 && x <= 0) || (tmp < 0 && x0 < 0 && x >= 0))
  89.     return(iter?iter-1:0);
  90.   }
  91. return(0);
  92. }
  93.  
  94.  
  95. /* Then a 32-bit one with even more complicated multiplication */
  96. long multi32(x, y)
  97. register long x,y;
  98. {
  99. register short sign= 1;
  100. register unsigned short xl, yl, xh, yh;
  101.  
  102. if (x < 0) {sign= -1; x= -x;}
  103. if (y < 0) {sign= -sign; y= -y;}
  104. xl= x & 0xffff; xh= x >> 16;
  105. yl= y & 0xffff; yh= y >> 16;
  106. return(sign*(long)( ((xh*yh) << (32-BIT32)) + ((xl*yh) >> (BIT32-16)) +
  107.         ((xh*yl) >> (BIT32-16)) + ((xl*yl) >> (BIT32))) );
  108. }
  109.  
  110. int mb32(xl, yl, cx, cy, iter)
  111. double xl, yl;
  112. double cx, cy;
  113. int iter;
  114. {
  115. long x= xl * (1 << BIT32), y= yl * (1 << BIT32);
  116. long x0= cx * (1 << BIT32), y0= cy * (1 << BIT32);
  117. long x2, y2;
  118. long tmp;
  119.  
  120. while (iter--) {
  121.   x2= multi32(x,x);
  122.   y2= multi32(y,y);
  123.   tmp= x2 + y2;
  124.   if (x < -(2 << BIT32) || x > (2 << BIT32) ||
  125.       y < -(2 << BIT32) || y > (2 << BIT32) ||
  126.       tmp < 0 || tmp > ((4 << BIT32)-1) )
  127.     return(iter);
  128.   y= multi32(x,y);
  129.   tmp= y + y;
  130.   if ((tmp > 0 && y <= 0) || (tmp < 0 && y >= 0))
  131.     return(iter?iter-1:0);
  132.   y= tmp + y0;
  133.   if ((tmp > 0 && y0 > 0 && y <= 0) || (tmp < 0 && y0 < 0 && y >= 0))
  134.     return(iter?iter-1:0);
  135.   tmp= x2 - y2;
  136.   x= tmp + x0;
  137.   if ((tmp > 0 && x0 > 0 && x <= 0 ) || (tmp < 0 && x0 < 0 && x >= 0))
  138.     return(iter?iter-1:0);
  139.   }
  140. return(0);
  141. }
  142. #endif
  143.  
  144.  
  145. /* Then, the simplest version of them all. (Slowest too.) */
  146. int mbf(x0, y0, cx, cy, iter)
  147. double x0, y0;
  148. double cx, cy;
  149. int iter;
  150. {
  151. double x= x0, y= y0;
  152. double x2, y2;
  153.  
  154. x0= cx; y0= cy;
  155. while (iter--) {
  156.   x2= x*x;
  157.   y2= y*y;
  158.   if (x2 + y2 >= 4)
  159.     return(iter);
  160.   y= 2 * x * y + y0;
  161.   x= x2 - y2 + x0;
  162.   }
  163. return(0);
  164. }
  165.  
  166. int iterate(x0, y0, dx, dy, cx, cy, nx, ny, iter, flags, dest)
  167. double x0, y0, dx, dy, cx, cy;
  168. int nx, ny;
  169. int iter;
  170. int flags;
  171. u_short **dest;
  172. {
  173. int px, py;
  174. int **pixel;
  175. int (*f)();
  176. long index= 0;
  177. long colorindex;
  178. long len;
  179. int bitindex= 7;
  180.  
  181. long size= PACKSIZE;
  182.  
  183. *dest= (u_short *)malloc(PACKSIZE * sizeof(u_short));
  184.  
  185. /* Determine which iterating function to use... */
  186. f= mbf; /* default to highest precision */
  187. #ifdef USE_INTEGER_MATH
  188. if (flags & JULIA) {
  189.   if (fabs(dx) > LIMIT_J32 && fabs(dy) > LIMIT_J32)
  190.     f= mb32;
  191.   if (fabs(dx) > LIMIT_J16 && fabs(dy) > LIMIT_J16)
  192.     f= mb16;
  193.   }
  194. else {
  195.   if (fabs(dx) > LIMIT_M32 && fabs(dy) > LIMIT_M32)
  196.     f= mb32;
  197.   if (fabs(dx) > LIMIT_M16 && fabs(dy) > LIMIT_M16)
  198.     f= mb16;
  199.   }
  200. #endif
  201.  
  202. /* alloc space for map of iterated pixels... */
  203. pixel= (int **)malloc(ny * sizeof(int *));
  204. for (py= 0; py < ny; py++) {
  205.   pixel[py]= (int *)malloc(nx*sizeof(int));
  206.   for (px= 0; px < nx; px++)
  207.     pixel[py][px]= -1;
  208.   }
  209.  
  210. /* Move pixel forward (relative to dir) */
  211. #define forward(dir) (ox= x, oy= y, qx= px,qy= py,\
  212.          (dir&1 ? (dir&2?(x-= dx, --px) : (x+= dx, ++px)) : \
  213.              dir&2?(y-= dy,++py) : (y+= dy,--py)))
  214. /* This nullifies the effect of previous forward() */
  215. #define stepback (x= ox, y= oy,px= qx,py= qy)
  216. /* And these are used to rotate dir 90 deg clockwise or anticlockwice */
  217. #define clockwise(dir) (dir==3 ? dir= 0 : ++dir)
  218. #define anticlockwise(dir) (dir==0 ? dir= 3 : --dir)
  219.  
  220.  
  221. /* First check if x and y are within limits; if not, return -1.
  222.  * Then check if pixel already calculated; if so, return its value.
  223.  * Otherwise calculate the pixel, store and return it.
  224.  */
  225. #define getpixel(px,py) (px < 0 || py < 0 || px >= nx || py >= ny ? -1 : \
  226.       (pixel[py][px] >= 0 ? pixel[py][px] : \
  227.       (pixel[py][px]= (flags & JULIA ? \
  228.       f(x, y, cx, cy, iter) : f(x, y, x, y, iter)))))
  229.  
  230. initialize(nx, ny); /* pixel recording module (pix.c) */
  231. for (;;) {
  232.   double x, y, ox, oy;
  233.   int qx, qy;
  234.   int startx, starty;
  235.   int prevpixel, newpixel;
  236.   int dir;
  237.   u_short color;
  238.  
  239. /* This tells us where to start from. */
  240.   if (first_point(&px, &py) == 0)
  241.     break;
  242.  
  243.   if (flags & FAST) {
  244.     if (index+2 >= size) *dest= (u_short *)
  245.       realloc(*dest, (size += PACKSIZE) * sizeof(u_short));
  246.     (*dest)[index++]= px;
  247.     (*dest)[index++]= py;
  248.     }
  249.  
  250.   x= x0+dx*px; y= y0-dy*py; dir= RIGHT;
  251.   startx= px; starty= py;
  252.   color= prevpixel= getpixel(px,py);
  253.   colorindex= index;
  254.   index += 2; /* space for color & len */
  255.   if (index >= size) *dest= (u_short *)
  256.     realloc(*dest, (size += PACKSIZE) * sizeof(u_short));
  257.   len= 0;
  258.  
  259. /* Another macro... the complexity of the crawling algorithm is hidden behind
  260.  * the facade of simple-looking macros ;) (Just kidding)
  261.  */
  262. #define record_dir(dir) { \
  263. if (bitindex <0 ) {bitindex= 7; index++;} \
  264. if (index >= size) *dest= (u_short *) \
  265.   realloc(*dest, (size += PACKSIZE) * sizeof(u_short)); \
  266. if (bitindex == 7) (*dest)[index]= 0; \
  267. (*dest)[index] |= ((u_short)dir) << (2*bitindex); \
  268. bitindex--; \
  269. }
  270.  
  271. /* This is the actual crawling algorithm. See how inherintly simple it is.
  272.  */
  273.   do {
  274.     forward(dir);
  275.     if ((newpixel= getpixel(px,py)) != prevpixel) {
  276.       stepback;
  277.       clockwise(dir);
  278.       }
  279.     else {
  280.       prevpixel= newpixel;
  281.       next_point(dir); /* every single pixel must be recorded */
  282.       record_dir(dir); /* ...twice */
  283.       len++;
  284.       anticlockwise(dir);
  285.       }
  286.     }
  287.   while (startx != px || starty != py || dir != RIGHT);
  288.  
  289.   if (bitindex != 7) {
  290.     bitindex= 7;
  291.     index++;
  292.     }
  293.   if (len==0) {
  294.     (*dest)[colorindex]=  color | 0x8000;
  295.     index--;
  296.     }
  297.   else {
  298.     (*dest)[colorindex]= color & 0x7fff;
  299.     (*dest)[colorindex+1]= len;
  300.     }
  301.  
  302.   } /* for */
  303.  
  304. for (py= 0; py < ny; py++) {
  305.   free(pixel[py]);
  306.   }
  307. free(pixel);
  308. deinit(); /* free some memory (pix.c) */
  309.  
  310. if (flags & FAST) {
  311.   if (index >= size) *dest= (u_short *)
  312.     realloc(*dest, (size += PACKSIZE) * sizeof(u_short));
  313.   (*dest)[index++]= LASTPIX;
  314.   }
  315. return(index);
  316. }
  317.  
  318.