NetNews Usenet Archive 1993 #1/NN_1993_1.iso/spool/sci/fractals/561

Text File | 1993-01-08 | 7.1 KB | 304 lines

Xref: sparky sci.fractals:561 comp.sys.mac.programmer:20963
Newsgroups: sci.fractals,comp.sys.mac.programmer
Path: sparky!uunet!mcsun!sun4nl!dutrun!donau!dutecai!reinoud
From: [e-mail] (R. Lamberts)
Subject: Julia's Dream algorithm + example code
Message-ID: <[e-mail]>
Originator: [e-mail]
Sender: [e-mail] (UseNet News System)
Nntp-Posting-Host: dutecai.et.tudelft.nl
Reply-To: [e-mail] (R. Lamberts)
Organization: Delft University of Technology, Dept. of Electrical Engineering
Date: Fri, 08 Jan 1993 05:27:18 GMT
Lines: 289


Hi there!

This is for all of you who remember Julia's Dream (JD), a Mac
application which calculates Julia sets really fast. Many
people asked me how it works, but I'm not that good at
explaining these things, so a lot of you must have been
dissappointed. I humbly apologise for never having written a
decent piece about the algorithm! However, here's a simplified
implementation of JD written in C, and before that another
attempt at explaining it.


The idea of Julia's Dream is very simple, even though I think
it's not easy to explain. Maybe it's best to have a look at the
source below if you don't like vague stories ;-). This code is
based on the same idea as the Mac version, though the
implementation is quite different.

The idea of JD is that with the escape-time algorithm, each
iteration is essentially independent of its context. It does
not make any difference whether a certain point in an orbit is
an intermediate position or it is a starting position
corresponding to a pixel. This is true for Julia sets, but not
for the Mandelbrot set, where the (c or mu) parameter varies
with the starting point (pixel). Now, because for a Julia set
each intermediate point in an orbit may also be considered as a
starting point, if a pixel is associated with each intermediate
orbit point, the calculation of a pixel value effectively
requires just one iteration.

There are various algorithms possible to implement this idea.
One way is to calculate an orbit for a certain starting point,
and assign pixel locations and values to each point in the
orbit afterwards. This is used in the code below. The original
Mac version of JD did calculations with each orbit point
rounded to a pixel location each iteration, allowing for
significant speed-up if the CPU is slow relative to memory
access.

If you want to see the Mac code, let me know. It's a lot of
hairy 68020 assembly code though...

To improve accuracy around the origin, some pixels are
calculated using the good old escape-time algorithm there.
There are probably more elegant ways to do this, but I didn't
exercise the options.

The code below compiles on a SGI Indigo with entry (8-bit)
graphics (sorry, I didn't even initialise the colormap). It
uses the GL, not X. This is *not* finished code, especially the
user interface is a quick hack. I've included the GL stuff with
the hope that a concrete example may be easier to understand
for some. Sorry about the lack of comments, this was really an
experiment to see how well the JD algorithm would work out in C
with straightforward FP calculations.


Well, I hope this is of any help (or interest, if you didn't
know JD). Many thanks to all who contacted me about this! If
you have ideas for (algorithmic) improvement, spiffy programs,
or further questions, I'd very much like to know.

Have fun,

- Reinoud


/*
 *  Julia's Dream II  (XP hack version)
 *
 *  The nightmare continues...
 *
 *  Copyright (c) 1992 Reinoud Lamberts, all rights reserved
 *
 *  You are free to use this stuff as you like, as long as
 *  you give credit where credit is due.
 *
 *
 *  Reinoud Lamberts
 *  Oldambt 69
 *  3524 BD Utrecht
 *  The Netherlands
 *
 *  Phone (int'l code 31) 30 896775
 *
 *  [e-mail]
 *  [e-mail] (UUCP)
 *  [e-mail]
 *
 */


#include <stdlib.h>
#include <gl/gl.h>
#include <gl/device.h>


#define XRES 256            /* better be power of 2 */
#define YRES 256
#define XLOG 8                /* shift amount i.s.o multiply */


typedef unsigned char byte;        /* should be 8 bits (and int 32) */
typedef double fl;            /* optimise for performance */


byte *ospmap;                /* things happen in here */


main()
{
    int xi, yi, ospindex, it, i, ospwsiz, bhsiz, jinx;
    fl deltax, deltay, idx, idy;
    fl xc, yc, sqre, sqim, sqyc, re, im, cre, cim;
    int orb[10000];
    int bound, empty, hitval;
    int pro, con;

/*
 *  yeah, let's arrange for some accomodation or something
 */

    ospmap=malloc(YRES*XRES);
    if (ospmap==NULL)
      exit(1);

/*
 *  let's do the window thing now
 */

    prefsize(XRES, YRES);
    winopen("Julia's Dream II  XP hack version  (c) 1992 Reinoud Lamberts");
    pixmode(PM_SIZE, 8);        /* 8 bits offscreen*/

/*
 *  now let's start talking business
 */


/*  constants */

    deltax=4.0/(XRES-1);            /* section of c plane is fixed */
    deltay=4.0/(YRES-1);
    idx=1/deltax;
    idy=1/deltay;

    ospwsiz=((YRES*XRES)>>2);
    bhsiz=XRES>>4;

    do
    {

      /* init */

      cre=-2+getvaluator(MOUSEX)*deltax;
      cim=-2+getvaluator(MOUSEY)*deltax;

      for (i=((YRES>>1)*XRES)>>2; i<ospwsiz; i++)
    ((int *) ospmap)[i]=0;

      /* enter the black hole */

      yc=deltay*0.5;

      for (yi=(YRES>>1); yi<(YRES>>1)+bhsiz; yi++)
      {
        xc=-(bhsiz+0.5)*deltax;
    sqyc=yc*yc;
    ospindex=(yi<<XLOG)+(XRES>>1)-1-bhsiz;

    for(xi=(XRES>>1)-1-bhsiz; xi<((XRES>>1)+bhsiz); xi++)
    {
      sqim=sqyc;
      sqre=xc*xc;
      re=xc;
      im=yc;
      it=0;

      while( (bound=((sqre+sqim)<4.0))&&(it<64) )
      {
       im=2.0*re*im+cim;
       re=sqre-sqim+cre;
       sqre=re*re;
       sqim=im*im;

       it++;
      }

      if (!bound)
        ospmap[ospindex]=255-it;
      else
        ospmap[ospindex]=255;

      xc+=deltax;
      ospindex++;
    }

    yc+=deltay;
      }

      /* explore outer space */

      ospindex=(YRES>>1)*XRES;
      yc=deltay*0.5;

      /* okay, here it goes */

      for (yi=(YRES>>1); yi<YRES; yi++)
      {
    xc=-2.0;
    sqyc=yc*yc;

    for (xi=0; xi<XRES; xi++)
    {
      sqim=sqyc;
      sqre=xc*xc;
      re=xc;
      im=yc;
      it=0;
      jinx=ospindex;

      while ( (bound=((sqre+sqim)<4.0)) &&
              (empty=((hitval=ospmap[jinx])==0)) && (it<10000)
        )
      {
        orb[it]=jinx;
        ospmap[jinx]=255;

        im=2.0*re*im+cim;
        re=sqre-sqim+cre;
        if (im<0)
        {
          im=-im;
          re=-re;
        }
        sqre=re*re;
        sqim=im*im;

        jinx=( ((int) ((im+2.0)*idy+0.5)) <<XLOG)+(re+2.0)*idx+0.5;
        it++;
      }
      if (it!=0)
      {
        if ((!empty)&&(hitval!=255))
          do
          {
        if (--hitval==0)
          hitval=254;
        ospmap[orb[--it]]=hitval;
          }
          while (it>0);
        else
        if (!bound)
        {
          hitval=255;
          do
          {
        if (--hitval==0)
          hitval=254;
        ospmap[orb[--it]]=hitval;
          }
          while (it>0);
        }
      }

      xc+=deltax;
      ospindex++;
    }

    yc+=deltay;
      }

      pro=0;
      con=YRES*XRES-1;
      do
      {
    ospmap[pro++]=ospmap[con--];
      }
      while(pro<=con);

      lrectwrite(0, 0, XRES-1, YRES-1, (unsigned long *) ospmap);

    }
    while (1);

    return 0;                /* easy on the compiler */
}