Sinclair BASIC Mandelbrot smooth colouring

[simon.snake]

Hi.

I am on a facebook group called BASIC Programming on the ZX Spectrum.

For fun, I wrote a program to generate images of the Mandelbrot.  The program is pretty slow (but the programming environment can speed it up substantially) and it was never written to be highly optimised.  Due to a limitation on the original ZX Spectrum hardware, only two colours are allowed in every 8x8 character on the screen.


Rather than try to add colouring and suffering from the aforementioned attribute clash, I have implemented two colouring methods:

The first is simple plain monochrome alternating iteration bands.

The second is white, checkerboard black/white, black and finally checkerboard white/black.  It then cycles around these 4.

Here are a couple of samples (of the same location) but with the different colouring options:

EDIT: Sorry, I didn't know facebook image links expire after a certain time, so the link has been removed!

Now, I'd like to experiment with smooth colouring.  It will still use only two colours - in this case white and black (or any 2 of the 8 primary colours of the machine), but I'd like to be able to plot only some of the pixels in order to get some sort of dithering.  In my mind I am looking for a value between 0 and 1, which when compared with a pseudo-random value from the machines RNG will be used to decide whether to plot the point or not, and that way I will achieve some sort of gradient.

I've read up on it, and the closest I've got to something at my level of understanding is the following site: linas.org/art-gallery/escape/smooth.html but the problem I have is knowing when to perform the code.  Is it calculated every iteration or only at the end?  I'm not very familiar with the notation used (I only ever got to GCSE C grade and didn't progress it further).

ZX BASIC has a natural log function (LN) and I'm hoping this is what is required.  In my code I have separate variables for the real and imaginary values of z, etc. so how do I convert what is being described on the site I have linked to (which looks more like FractInt formula parser code) into Sinclair BASIC?

I don't know if it helps but I'm happy sharing the code I have written.  It's in a format that is used by a programming environment for Windows called BasinC:

Code:

10 REM Hires Mandelbrot
20 REM 2016 Simon Ferre'
30 BORDER 7: PAPER 7: INK 0: CLS 
40 LET dither=1
50 LET iter=64
60 LET width=4
70 LET height=3
80 LET xstart=width/2-width
90 LET ystart=height/2-height
100 LET xstep=width/255
110 LET ystep=height/175
120 LET res=16
130 PRINT #0;AT 1,0;"{vi}z{vn}oom {vi}i{vn}ters {vi}d{vn}ither"
140 FOR y=0 TO 175 STEP res
150 IF res=1 THEN PRINT #0;AT 1,28;INT (y/1.75);"%";
160 FOR x=0 TO 255 STEP res
170 PRINT #0;AT 0,0;"x,y=(";x;",";y;")  "
180 REM {vi}Iterate{vn}
190 IF POINT (x,y) THEN GO TO 450
200 LET i=0: LET it=1: REM Iter count and Iter type
210 LET zr=xstart+x*xstep: LET JuliaR=zr
220 LET zi=ystart+y*ystep: LET JuliaI=zi
230 REM {vi}Inner iterations{vn}
240 IF k$="" THEN LET k$=INKEY$: IF k$="" THEN GO TO 280
250 IF k$="d" THEN INPUT "Dither (y/n)? "; LINE k$: LET dither=0+(1 AND k$="y"): CLS : GO TO 120
260 IF k$="z" THEN GO TO 480
270 IF k$="i" THEN INPUT ("Old iters=";iter;", new value=");iter: CLS : GO TO 120
280 LET zrsqr=zr*zr
290 LET zisqr=zi*zi
300 IF zisqr+zrsqr>4 OR i>=iter THEN GO TO 360
310 LET zi=zr*zi*2+JuliaI
320 LET zr=zrsqr-zisqr+JuliaR
330 LET i=i+1: LET it=it+1: IF it=5 THEN LET it=1
340 GO TO 230
350 REM {vi}Escaped or max iter{vn}
360 IF dither THEN GO TO 390
370 IF i/2=INT (i/2) THEN PLOT x,y
380 GO TO 450
390 IF i>=iter THEN LET it=1
400 IF it=1 THEN GO TO 450
410 IF it=3 THEN PLOT x,y: GO TO 450
420 LET p=0: LET xy=x+y: IF xy/2=INT (xy/2) THEN LET p=1
430 IF it=2 THEN PLOT INVERSE NOT p;x,y
440 IF it=4 THEN PLOT INVERSE p;x,y
450 NEXT x
460 NEXT y
470 LET res=res/2: IF res>=1 THEN GO TO 140
480 LET xl=0: LET xr=255: LET yb=0: LET yt=175: LET ratio=175/255
490 PRINT #0;AT 0,0;"{vi}zZ{vn}oom, {vi}Enter{vn}=zoom in": PRINT #0;AT 1,0;"{vi}QAOP{vn} move box, {vi}Delete{vn}=zoom out"
500 PLOT OVER 1;xl,yb: DRAW OVER 1;xr-xl,0: DRAW OVER 1;0,yt-yb: DRAW OVER 1;-(xr-xl),0: DRAW OVER 1;0,-(yt-yb): PLOT OVER 1;xl,yb
510 LET k$=INKEY$: IF k$="" THEN GO TO 510
520 PLOT OVER 1;xl,yb: DRAW OVER 1;xr-xl,0: DRAW OVER 1;0,yt-yb: DRAW OVER 1;-(xr-xl),0: DRAW OVER 1;0,-(yt-yb): PLOT OVER 1;xl,yb
530 IF k$="z" AND xr-xl>16 THEN LET xl=xl+2: LET xr=xr-2: LET yb=yb+2*ratio: LET yt=yt-2*ratio: GO TO 500
540 IF k$<>"Z" THEN GO TO 600
550 IF xl>1 THEN LET xl=xl-2: LET xr=xr+2: LET yb=yb-2*ratio: LET yt=yt+2*ratio
560 IF xl<0 THEN LET o=0-xl: LET xl=0: IF xr+o<=255 THEN LET xr=xr+o
570 IF xr>255 THEN LET o=xr-255: LET xr=255: IF xl-o>=0 THEN LET xl=xl-o
580 IF yb<0 THEN LET o=0-yb: LET yb=0: IF yt+o<=175 THEN LET yt=yt+o
590 IF yt>175 THEN LET o=yt-175: LET yt=175: IF yb-o>=0 THEN LET yb=yb-o
600 IF k$="o" AND xl>2 THEN LET xl=xl-2: LET xr=xr-2
610 IF k$="p" AND xr<253 THEN LET xl=xl+2: LET xr=xr+2
620 IF k$="q" AND yt<173 THEN LET yb=yb+2: LET yt=yt+2
630 IF k$="a" AND yb>2 THEN LET yb=yb-2: LET yt=yt-2
640 IF k$<>CHR$ 13 THEN GO TO 710
650 REM set new offsets and redraw
660 LET width=(xr-xl)*xstep
670 LET height=(yt-yb)*ystep
680 LET xstart=xstart+xl*xstep
690 LET ystart=ystart+yb*ystep
700 CLS : LET k$="": GO TO 100
710 IF k$<>CHR$ 12 THEN GO TO 730
720 REM Zoom out
730 IF k$="0" THEN STOP 
740 GO TO 500

In the code {vi} and z{vn} are directives that specify normal or inverse video (and can safely be ignored).  The code between lines 180 and 340 is performing the maths for each point.

Simon

[quaz0r]

it's just a calculation you do with the last value of Z after youve escaped, to compute a fractional value to add to your iteration count, which is what allows you to create the smooth coloring.  i guess the main notation in the formula on that site is |Z|, which means the magnitude of Z, which is .  and yeah the log is the natural log.

In my mind I am looking for a value between 0 and 1

so, this part is the fractional value between 0 and 1 then:

1 - (log (log |Z|) / log 2)

with the formula usually stated as that plus N, or the iteration you escaped at.

[simon.snake]

1 - (log (log |Z|) / log 2)

with the formula usually stated as that plus N, or the iteration you escaped at.

Ok, I think I get this.  I entered what I believe to be the right code and it started to generate a screen, but then threw an error.  It would appear that on the ZX Spectrum the number passed to LN (natural log) has to be positive.  If I use instead 1-(log(abs(log(abs(|z|)))) / log 2) is it equivalent?

I'll give it a go anyway as I've got nothing to lose.  I will report back with any useful images I get.

[quaz0r]

taking the magnitude of a complex number always results in a positive value, as you are squaring the real and imaginary components and adding them together.   if log|Z| is resulting in a negative, maybe your bailout isnt high enough ?  personally i use the smoothing formula posted in this thread.  its a little more complicated and im not sure what the difference between the different formulas are, in fact there are several smoothing formulas out there i think..

[simon.snake]

I think the issue with LN is that it was producing a negative value.  Not sure if it is meant to but the ABS have at least allowed it to run.

Here's the location in the original post but with the smoothing applied:

EDIT: Sorry, I didn't know facebook image links expire after a certain time, so the link has been removed!

And the initial starting location of the Mandelbrot (as a side effect of my code I now have shading for the inside colouring too):

EDIT: Sorry, I didn't know facebook image links expire after a certain time, so the link has been removed!

It's not as smooth as I was hoping, but that's either due to the low resolution, limitations within the floating point routines within Sinclair BASIC or maybe something else that can be resolved.

For comparison, here is the start image (that the smooth one is based on) with the usual black/white stripes:

EDIT: Sorry, I didn't know facebook image links expire after a certain time, so the link has been removed!

And here's a zoom into the spike:

EDIT: Sorry, I didn't know facebook image links expire after a certain time, so the link has been removed!

As you can see, it's not that smooth really.

One last thing, here's my updated code (it has become a bit of a mess but I'll attempt to tidy it up once it works):

Code:

10 REM {vi}Hires Mandelbrot{vn}
20 REM {vi}*2016 Simon Ferre'{vn}
30 BORDER 7: PAPER 7: INK 0: CLS : LET k$=""
40 LET dither=0
45 LET smooth=0
50 LET iter=64
60 LET width=4
70 LET height=3
80 LET xstart=width/2-width
90 LET ystart=height/2-height
100 LET xstep=width/255
110 LET ystep=height/175
120 LET res=1
130 PRINT #0;AT 1,0;"{vi}z{vn}oom {vi}i{vn}ters {vi}d{vn}ither {vi}s{vn}mooth"
140 FOR y=0 TO 175 STEP res
150 IF res=1 THEN PRINT #0;AT 1,28;INT (y/1.75);"%";
160 FOR x=0 TO 255 STEP res
170 PRINT #0;AT 0,0;"x,y=(";x;",";y;")  "
180 REM {vi}Iterate{vn}
190 IF POINT (x,y) THEN GO TO 450
200 LET i=0: LET it=1: REM Iter count and Iter type
210 LET zr=xstart+x*xstep: LET JuliaR=zr
220 LET zi=ystart+y*ystep: LET JuliaI=zi
230 REM {vi}Inner iterations{vn}
240 IF k$="" THEN LET k$=INKEY$: IF k$="" THEN GO TO 280
250 IF k$="d" THEN INPUT "Dither (y/n)? "; LINE k$: LET dither=0+(1 AND k$="y"): CLS : GO TO 120
255 IF k$="s" THEN INPUT "Smooth (y/n)? "; LINE k$: LET smooth=0+(1 AND k$="y"): CLS : GO TO 120
260 IF k$="z" THEN GO TO 480
270 IF k$="i" THEN INPUT ("Old iters=";iter;", new value=");iter: CLS : GO TO 120
280 LET zrsqr=zr*zr
290 LET zisqr=zi*zi
300 IF zisqr+zrsqr>4 OR i>=iter THEN GO TO 360
310 LET zi=zr*zi*2+JuliaI
320 LET zr=zrsqr-zisqr+JuliaR
330 LET i=i+1: LET it=it+1: IF it=5 THEN LET it=1
340 GO TO 230
350 REM {vi}Escaped or max iter{vn}
360 IF dither THEN GO TO 390
364 IF NOT smooth THEN GO TO 370
365 LET s=(LN ABS (LN ABS (SQR (ABS (zisqr+zrsqr)))/LN 2)): IF i/2=INT (i/2) THEN LET s=1-s
367 IF s>RND THEN PLOT x,y
368 GO TO 450
370 IF i/2=INT (i/2) THEN PLOT x,y
380 GO TO 450
390 IF i>=iter THEN LET it=1
400 IF it=1 THEN GO TO 450
410 IF it=3 THEN PLOT x,y: GO TO 450
420 LET p=0: LET xy=x+y: IF xy/2=INT (xy/2) THEN LET p=1
430 IF it=2 THEN PLOT INVERSE NOT p;x,y
440 IF it=4 THEN PLOT INVERSE p;x,y
450 NEXT x
460 NEXT y
470 LET res=res/2: IF res>=1 THEN GO TO 140
480 LET xl=0: LET xr=255: LET yb=0: LET yt=175: LET ratio=175/255
490 PRINT #0;AT 0,0;"{vi}zZ{vn}oom, {vi}Enter{vn}=zoom in": PRINT #0;AT 1,0;"{vi}QAOP{vn} move box, {vi}Delete{vn}=zoom out"
500 PLOT OVER 1;xl,yb: DRAW OVER 1;xr-xl,0: DRAW OVER 1;0,yt-yb: DRAW OVER 1;-(xr-xl),0: DRAW OVER 1;0,-(yt-yb): PLOT OVER 1;xl,yb
510 LET k$=INKEY$: IF k$="" THEN GO TO 510
520 PLOT OVER 1;xl,yb: DRAW OVER 1;xr-xl,0: DRAW OVER 1;0,yt-yb: DRAW OVER 1;-(xr-xl),0: DRAW OVER 1;0,-(yt-yb): PLOT OVER 1;xl,yb
530 IF k$="z" AND xr-xl>16 THEN LET xl=xl+2: LET xr=xr-2: LET yb=yb+2*ratio: LET yt=yt-2*ratio: GO TO 500
540 IF k$<>"Z" THEN GO TO 600
550 IF xl>1 THEN LET xl=xl-2: LET xr=xr+2: LET yb=yb-2*ratio: LET yt=yt+2*ratio
560 IF xl<0 THEN LET o=0-xl: LET xl=0: IF xr+o<=255 THEN LET xr=xr+o
570 IF xr>255 THEN LET o=xr-255: LET xr=255: IF xl-o>=0 THEN LET xl=xl-o
580 IF yb<0 THEN LET o=0-yb: LET yb=0: IF yt+o<=175 THEN LET yt=yt+o
590 IF yt>175 THEN LET o=yt-175: LET yt=175: IF yb-o>=0 THEN LET yb=yb-o
600 IF k$="o" AND xl>2 THEN LET xl=xl-2: LET xr=xr-2
610 IF k$="p" AND xr<253 THEN LET xl=xl+2: LET xr=xr+2
620 IF k$="q" AND yt<173 THEN LET yb=yb+2: LET yt=yt+2
630 IF k$="a" AND yb>2 THEN LET yb=yb-2: LET yt=yt-2
640 IF k$<>CHR$ 13 THEN GO TO 710
650 REM set new offsets and redraw
660 LET width=(xr-xl)*xstep
670 LET height=(yt-yb)*ystep
680 LET xstart=xstart+xl*xstep
690 LET ystart=ystart+yb*ystep
700 CLS : LET k$="": GO TO 100
710 IF k$<>CHR$ 12 THEN GO TO 730
720 REM Zoom out
730 IF k$="0" THEN STOP 
740 GO TO 500

[quaz0r]

i think the smoothing works better with a higher bailout, so maybe see if >100000 or something looks better than >4 ?  though i dont know how buttery smooth dithered monotone is really going to get  

what would probably work best for you actually is to use the distance estimate.  iterate the derivative of Z alongside (but before) Z:

Z₀ = 0
Z'₀ = 0
Z'_n+1 = 2Z_n'Z_n + 1
Z_n+1 = Z_n² + C

then when you escape the distance estimate is , which you can then use to index into a palette or something similar to escape time coloring, or check it against a threshold value to color either black or white, or something like that.

[DarkBeam]

Home computers still kick ... Great 

[simon.snake]

i think the smoothing works better with a higher bailout, so maybe see if >100000 or something looks better than >4 ?  though i dont know how buttery smooth dithered monotone is really going to get  

Yes, I do agree.  I tried to increase the bailout but that didn't look any better.  I think my issue is the sharp turn from white to black at the edge of some of the iteration bands.  What's strange is the other side looks fine.  Once I implement my zoom out code I'll get a better view of the whole of the first iteration band - it may offer clues.

The first image is with a bailout of 3, the second is with a bailout of 16.  A bailout of 100 or even 1000 just gives me back a black/white striped image with no dithering.

EDIT: Sorry, I didn't know facebook image links expire after a certain time, so the link has been removed!

Am I supposed to be changing anything in the equation with all the LNs when I change the bailout?

[hermann]

A Sinclair ZX81 was my first computer! At a time I was still in University. About the mid of the eighties.
Hermann

[simon.snake]

Researching further, I'm getting errors because the output of the LN command on the ZX Spectrum is returning negative values.  Whether it is meant to or not is something I have no idea about.

Subsequently, the LN function on the ZX Spectrum does not allow one to use a negative value.  This is (from what I've read on the subject) the usual case for LN.

Now I am at a loss on what to do.  Whatever tweaks I make just don't seem to make the images look better than those I have already posted.

I'm not going to give up though.  It's keeping my grey matter active.

[quaz0r]

if you look at the graph you can see the range of inputs and outputs to the log function.

i think if you use a higher bailout you probably need to use the extended formula i linked to earlier which calculates the log(p) part instead of using log(2) or whatever.

also i still think you might want to check out distance estimation since you are working with a b&w system anyway    you can get nice b&w renders of the boundary like this

[simon.snake]

I'm pretty happy with it now.

The following code:

smooth iter = iter + 1 + ( log(log(bailout)-log(log(cabs(z))) )/log(2)

Was mentioned as a smoothing method for use in ultrafractal in the following site:  https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/MandelbrotSetExterior

So I tweaked my code to include this additional log(log(bailout)) (even with its missing closing bracket) and the following image was calculated:

EDIT: Sorry, I didn't know facebook image links expire after a certain time, so the link has been removed!

That's exactly how I envisioned the end result looking.