IFS is what Michael Barnsley calls his Iterated Function Systems, and IFS is another part of fractals (remember those Mandelbrot all-nighters ?). Barnsley did IFS a few years ago and has gone on to bigger and better stuff now. He explains how IFS can be done by anyone with a personal computer using a “collage” approach. Ready for a collage education? Read on.
What’s the big idea ? – data compression! with the IFS algorithm and just a few well chosen coefficients, you can represent a very complicated picture that otherwise would require many kilobytes to describe. In other words ‘one picture is worth a thousand words’ but an IFS only needs a dozen or so coefficients.
The concept is that a bigger picture is composed of self-similar smaller pictures that just happen to cover the bigger picture nicely. For example, a checkerboard is just a big square that can be neatly covered with 64 smaller squares. Actually we only need 32, since half are blank. You get the idea. Now it can be more complicated, since we can also use rectangles, which can be considered as self-similar to squares where we have changed the scales along the x and/or y axes. Not only that, but we are also allowed to tilt and rotate the little pieces to help us get a collage with better coverage. Not too useful for squares, but real handy for more organic swirly ‘natural’ shapes like Barnsley’s popular fern patterns.
Once you have the big outline, and then get a covering of it with a collage of smaller self-similar shapes, you can work out the coefficients. A convenient way is to use three point affine transformations to develop sets of simultaneous equations. Solving these three at a time produces the coefficients needed by the IFS. You can then assign probability values to each transformation set based on how dark you want it to be in the final picture.
With the coefficients and probabilities in hand, the drudge work begins, such as calculating and plotting 10,000 points where the coordinates of each point depends on the coordinates of the previous point, (with a little randomness thrown in). This tedious task is the kind of thing computers do well. It takes about 10 minutes for the QL to do 10000 points. As you watch and wait you will see a nebulous scatter of points take on shape and substance.
I’ve written a SuperBasic program to do 2-dimensional IFS, copying some of the better examples from the literature and working out a few on my own. I’ve also programmmed IFS in Small-C.
Barnsley has gone on to 3-D and color, and real time HDTV encoding/decoding etc. with dedicated microchips. One of Barnsley’s pet IFS pictures is a fernlike image, a complicated structure that is encoded with only a few coefficients. If you’d like to try it on your computer the data given by Barnsley for a black spleenwort fern is:
| Translation | Rotation | Scale | Probability |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 0 | .16 | .005 | 0 |
| 1.6 | -2.5 | -2.5 | .85 |
| .85 | .8 | 0 | 1.6 |
| 49 | 49 | .3 | .34 |
| .0975 | 0 | .44 | 120 |
| -50 | .3 | .37 | .0975 |
We can convert these into a simpler form by using trigonometric functions to put the rotations into their equivalent x,y shifts. After doing that we end up with a table (from BYTE Jan ’88) such as:
| W | a | b | c | d | e | f | p |
|---|---|---|---|---|---|---|---|
| 1 | 0 | 0 | 0 | .16 | 0 | 0 | .01 |
| 2 | .2 | -.26 | .23 | .22 | 0 | 1.6 | .07 |
| 3 | -.15 | .28 | .26 | .24 | 0 | .44 | .07 |
| 4 | .85 | .04 | -.04 | .85 | 0 | 1.6 | .85 |
These values are used in the SuperBasic program ‘IFS4_bas’ and can be entered at the prompts in Small-C IFS4_c. The algorithm is said to be ‘robust’, and small changes in the values only change the final result slightly. Enjoy the show by changing one or more of them and allow about 10 minutes (or slightly less using the Small-C version) for the program to do its thing.
If there is interest we have coefficients for a few other figures, such as snowflakes, leaves, curlicues, galaxies, pentagons, dragon curves, etc.
References:
- BYTE January and April 1988
- Algorithms May/June 1990
- Scientific American August 1990
- Electronics & Wireless World August 1990
There are a series of books out by Robert L. Devaney, some at the popular level, and others for hard-core mathematicians. He also made some beautiful educational videos. I recommend the books and videos to those who are curious about things like fractals, chaos, Mandelbrot and Julia sets, strange attractors, etc.
SuperBASIC
100 REMark IFS4_bas Iterative Function System
110 REMark H. L. Schaaf June 3, 1991
120 n$="fern" :REMark change to suit
130 REMark set up to do fern
140 RESTORE
150 WMON : MODE 4:WINDOW #1,512,256,0,0
160 PAPER #1,4: INK #1,0:CLS#0
170 PRINT "probabilities",
180 REMark m is the number of transformations
190 m=4 : REMark change to suit, fern only uses 4
200 DIM a(m), b(m), c(m), d(m), e(m), f(m), p(m)
210 DATA 0, 0, 0, .16, 0, 0, 1E-2
220 DATA .2, -.26, .23, .22, 0, 1.6, 7E-2
230 DATA -.15, .28, .26, .24, 0, .44, 7E-2
240 DATA .85, 4E-2, -4E-2, .85, 0, 1.6, .85
250 pt=0 :REMark total of probabilities
260 FOR j=1 TO m
270 READ a(j), b(j), c(j), d(j), e(j), f(j), p(j)
280 pt = pt + p(j)
290 END FOR j
300 FOR j = 1 TO m: p(j)= p(j)/pt :PRINT p(j),:END FOR j
310 pt=0
320 FOR i = 1 TO m
330 pt=pt+p(i) : p(i) = pt
340 PRINT\ a(i),b(i),c(i),d(i),e(i),f(i),p(i)
350 END FOR i
360 x=0 : y=0 : maxy=-1E6 : maxx=-1E6 : minx=1E6 : miny=1E6
370 FOR n = 1 TO 100
380 AT 5,0: PRINT n,
390 pk = RND
400 PRINT pk,
410 k = 1
420 REPeat loopa
430 IF pk <= p(k) THEN EXIT loopa
440 k = k+1
450 END REPeat loopa
460 xnxt=a(k)*x+b(k)*y+e(k)
470 ynxt=c(k)*x+d(k)*y+f(k)
480 x=xnxt : y=ynxt
490 IF y > maxy THEN maxy = y
500 IF x > maxx THEN maxx = x
510 IF x < minx THEN minx = x
520 IF y < miny THEN miny = y
530 PRINT x,y\minx,maxx,miny,maxy
540 END FOR n
550 CLS
560 yscale = 1.1* ABS(maxy-miny)
570 xscale = 1.1* ABS(maxx-minx)
580 PRINT "X & Y & P scales = ",xscale,yscale,
590 IF xscale>yscale : pscale=xscale :ELSE pscale = yscale
600 PRINT pscale\
610 xoffset = minx: yoffset = miny
620 xfill = ABS(maxx-minx) : xmarg=(1.3*pscale)-xfill
630 PRINT "X range & margin ",xfill,xmarg\
640 yfill = ABS(maxy-miny) : ymarg=pscale-yfill
650 PRINT "Y range & margin ",yfill,ymarg\
660 xoffset = minx -.5*xmarg
670 yoffset = miny -.5*ymarg
680 PRINT "Scale = ";pscale;" lower left (x,y)= ";
xoffset, yoffset
690 CLS#0: PRINT#0;"any key to begin 5000 iterations";
" = +/- 5 minutes"
700 PAUSE
710 CLS
720 SCALE pscale,xoffset,yoffset
730 iters = 5000 :REMark 5000 iterations change to suit
740 FOR n = 1 TO iters
750 pk = RND
760 k = 1
770 REPeat loopb
780 IF pk<=p(k) THEN EXIT loopb
790 k = k+1
800 END REPeat loopb
810 xnxt=a(k)*x+b(k)*y+e(k)
820 ynxt=c(k)*x+d(k)*y+f(k)
830 x = xnxt : y = ynxt
840 POINT x,y
850 END FOR n
860 BEEP 30000,2 :REMark lets you know it's done
870 SBYTES "ram1_"&n$&"_pix",131072,32768
880 REMark saves it 'clean'
890 PRINT #0\"Clean screen in ram1_"&n$&"_pix"
900 REMark you can save it to flp, etc.
C
/* IFS4b_c H.L.Schaaf June 3, 1991
* Iterated Function Systems
*/
#include <stdio_h>
int a[21], b[21], c[21], d[21], e[21], f[21], p[21],
ps[21], pk[3], fpw[3], fpx[3], fpy[3], fpz[3],
rnd[3], ra[3], rd[3], rm[3], rr[3], det[3], tmp[3],
w[3], x[3], y[3], z[3], xs[3], ys[3],
maxx[3], maxy[3], minx[3], miny[3],
yscale[3], xscale[3], pscale[3],
xoffset[3], yoffset[3], i, j, k, n, t, iter;
char *ow, ch, str[16];
main() {
initrand();
prompt("choose number of transforms"); choose_n();
prompt("enter array values"); selarray();
prompt("determine raw probabilities"); rawprob();
prompt("equalize probabilities, sum = 1"); sumprob();
prompt("get ranges of graphics values"); sizeup();
prompt("see IFS develop, 1/2 minute for K points");
forshow(x,y);
at(ow,0,0); printf("any key to exit"); getchar();
exit();
}
getsit() {
printf("?"); gets(str);
*tmp = atof(str); return tmp;
}
selarray() {
csize(ow,0,0); cls(ow);
printf("\n%7c%14c%14c%14c%14c%14c\n"
,'a','b','c','d','e','f');
for ( j = 0; j < n; j++ ) {
k = 2 * j;
at(ow,( 5 + k ), 0 );
printf("%u",j);
}
for ( j = 0; j < n; ++j ) {
k = 2 * j;
at(ow,( 5 + k),(7+(0*14))); getsit(); fputa(a,j,tmp);
at(ow,(5 + k),6); fpica(a,j,tmp);
printf("%s \n",ftoa(tmp,str));
at(ow,( 5 + k),(7+(1*14))); getsit(); fputa(b,j,tmp);
at(ow,(5 + k),20); fpica(b,j,tmp);
printf("%s \n",ftoa(tmp,str));
at(ow,( 5 + k),(7+(2*14))); getsit(); fputa(c,j,tmp);
at(ow,(5 + k),34); fpica(c,j,tmp);
printf("%s \n",ftoa(tmp,str));
at(ow,( 5 + k),(7+(3*14))); getsit(); fputa(d,j,tmp);
at(ow,(5 + k),48); fpica(d,j,tmp);
printf("%s \n",ftoa(tmp,str));
at(ow,( 5 + k),(7+(4*14))); getsit(); fputa(e,j,tmp);
at(ow,(5 + k),62); fpica(e,j,tmp);
printf("%s \n",ftoa(tmp,str));
at(ow,( 5 + k),(7+(5*14))); getsit(); fputa(f,j,tmp);
at(ow,(5 + k),76); fpica(f,j,tmp);
printf("%s \n",ftoa(tmp,str));
}
printf("\nIs this OK ? [y/n]");
if((ch = getchar()) == 'n') { selarray(); }
else { csize(ow,1,0); }
}
askmore() {
at(ow,24,0);
printf("%uK want another K ?",iter); ch = getchar();
if(ch != 'n') { at(ow,24,0);
printf(" "); doakay();
}
}
forshow(x,y) int x[3], y[3]; {
scale(ow,pscale,xoffset,yoffset); ink(ow,6); cls(ow);
point(ow,maxx,maxy); point(ow,maxx,miny);
point(ow,minx,maxy); point(ow,minx,miny);
ink(ow,4); iter = 0; doakay();
}
doakay() {
for ( i = 0 ; i < 1024 ; i++) {
iterate(x,y); point(ow,x,y);
}
iter ++; askmore();
}
iterate(x,y) int x[3], y[3]; {
k = selectk(); fpica(a,k,fpw); fpica(b,k,fpx);
fpica(c,k,fpy); fpica(d,k,fpz); *xs = fmult(fpw,x);
*tmp = fmult(fpx,y); *xs = fadd(xs,tmp);
fpica(e,k,fpw); *xs = fadd(xs,fpw);
*ys = fmult(fpy,x); *tmp = fmult(fpz,y);
*ys = fadd(ys,tmp); fpica(f,k,fpw);
*ys = fadd(ys,fpw); fmove(x,xs);
fmove(y,ys); return (x,y);
}
selectk() {
fmove(pk,rand());
for(k=0; k<n; k++) { fpica(ps,k,fpw);
*fpz = fsub(pk,fpw); if(int(fpz) < 0) break;
}
return (k);
}
sizeup() {
trace("BEGINNING OF SIZEUP");
*x = float(0); *y = float(0); *xs = float(0);
*ys = float(0); *maxy = float(-1000);
*maxx = float(-1000); *minx = float(1000);
*miny = float(1000); cls(stdout);
trace("LOOPING INTO 200 ITERATIONS");
for(i=0; i<200; i++) { iterate(x,y);
if (fcmp(y,maxy) >0) { fmove(maxy,y); }
if (fcmp(x,maxx) >0) { fmove(maxx,x); }
if (fcmp(y,miny) <0) { fmove(miny,y); }
if (fcmp(x,minx) <0) { fmove(minx,x); }
at(stdout,10,31); printf("%u \n",i);
at(stdout,10,8);
printf("X min = %s ", ftoa(minx,str));
at(stdout,10,36);
printf("X max = %s \n",ftoa(maxx,str));
at(stdout,16,25);
printf("Y min = %s ", ftoa(miny,str));
at(stdout,4,25);
printf("Y max = %s \n",ftoa(maxy,str));
}
trace("BEGINNING GRAPHICS SCALING");
*w = atof("1.2"); /* 10% margin */
*tmp = fsub(maxy,miny); *tmp = fabs(tmp);
*yscale = fmult(tmp,w); *tmp = fsub(maxx,minx);
*tmp = fabs(tmp); *xscale = fmult(tmp,w);
if(fcmp(xscale,yscale) >0 ) fmove(pscale,xscale);
else fmove(pscale,yscale);
fmove(xoffset,minx); fmove(yoffset,miny);
*fpx = fsub(maxx,minx); *fpx = fabs(fpx);
*fpy = fsub(maxy,miny); *fpy = fabs(fpy);
*w = atof("1.3"); *fpw = fmult(pscale,w);
*fpw = fsub(fpw,fpx); *fpz = fsub(pscale,fpy);
*w = atof(".5"); *fpw = fmult(w,fpw);
*fpz = fmult(w,fpz); *xoffset = fsub(minx,fpw);
*yoffset = fsub(miny,fpz);
trace("END OF SCALING\n\n");
trace("Scaling factors based on 200 iterations:");
printf("pscale = %s \n",ftoa(pscale,str));
printf("xoffset = %s \n",ftoa(xoffset,str));
printf("yoffset = %s \n",ftoa(yoffset,str));
printf("x = %s \n",ftoa(x,str));
printf("y = %s \n",ftoa(y,str));
scale(stdout,pscale,xoffset,yoffset);
trace("Ready to do a IFS\n");
return(pscale,xoffset,yoffset);
}
sumprob() {
z[0] = float(0);
for( i = 0; i < n; i++) { fpica(p,i,fpw);
*tmp = fdiv(fpw,pk); *z = fadd(z,tmp);
printf("cumulative probability at %u is %s \n"
,i,ftoa(z,str));
fputa(ps,i,z); fputa(p,i,tmp);
}
printf("sum of probabilities = %s\n",ftoa(z,str));
}
rawprob() {
*pk = float(0);
for( i=0; i<n; i++) { fpica(a,i,fpw);
fpica(d,i,fpx); fpica(b,i,fpy);
fpica(c,i,fpz); *det = fmult(fpw,fpx);
*tmp = fmult(fpy,fpz); *det = fsub(det,tmp);
*det = fabs(det);
if( *det == float(0) ) { *det = atof("+.01");
trace("ZERO DETERMINANT !");
}
*pk = fadd(pk,det); fputa(p,i,det);
printf("det(%u) = %s\n",i,ftoa(det,str));
}
printf("Sum of raw probabilities = %s \n",ftoa(pk,str));
printf("\nChange probability value(s) ? [y/n]");
while((ch = getchar()) != 'n') {
iter = 99;
while( iter <0 || iter >(n-1)) {
printf(" Which transform 0 to %u ? ",n-1);
fgets(ch,1,stdin);
iter = atoi(ch);
printf(" probability for %u ",iter);
getsit();
fputa(p,iter,tmp);
printf(" = %12s \n",ftoa(tmp,str));
}
printf("\n Review of probabilities \n\n");
*pk = float(0);
for( i = 0; i < n; i++) {
fpica(p,i,fpw);
printf(" transform %u has probability of %12s\n"
,i,ftoa(fpw,str));
*pk = fadd(pk,fpw);
}
printf("\nChange probability value(s) ? [y/n]");
}
}
fputa(array,element,value) int array[21], element,
value[3]; { int i; element *= 3;
for(i = 0; i < 3 ; i++) { array[i+element] = value[i];
}
}
fpica(array,element,value) int array[21], element,
value[3]; { int i; element *= 3;
for(i = 0; i < 3 ; i++) { value[i] = array[i+element];
}
}
rand() {
*ra = fmult(ra,rm); *ra = fdiv(ra,rd); t = int(ra);
*rr = float(t); *rnd = fsub(ra,rr);
*ra = fmult(rnd,rd); return rnd ;
}
choose_n() {
cls(stdout); n = 0;
while( n<2 || n >7 ) {
trace("Choose a number between 2 and 6\n");
fgets(ch,1,stdin); n = atoi(ch);
printf("%2u was chosen \n",n);
}
return(n);
}
prompt(str) char str[]; {
printf("\nTouch [ SPACE BAR ] to %s\n",str);
getchar();
}
trace(str) char str[]; {
printf("\ntrace %s\n",str);
}
initrand() {
*rm = float(125); *rd = atof("2796203");
*ra = atof("100001");
for( i = 0; i<21; i++) {
a[i]= b[i]= c[i]= d[i]= e[i]= f[i]= p[i]= ps[i]= 0;
}
}
_console() {
mode(4); fopen(ow,"con","w"); window(ow,512,256,0,0);
paper(ow,0); ink(ow,4); border(ow,1,2); cls(ow);
}
/* JUNE 3 1991 10pm * */