Julia Set (shader)

[Richard Russell]

Here's a shader-based program to plot the Julia Set, showing that using the GPU from BBC BASIC need not be complicated. Indeed the shader code is virtually identical to the BBC BASIC code you would use to plot it natively. It will run in BBC BASIC for Windows or BBC BASIC for SDL 2.0.

To emphasise its simplicity I've omitted the code needed to allow real-time panning-and-zooming. This can easily be transplanted from the supplied example program mandel.bbc into this one if you wish.

Code: Select all

      MODE 9
      INSTALL @lib$ + "shaderlib"
      DIM Vertex$(10), Fragment$(20)
      PROC_readshader(Vertex$())
      PROC_readshader(Fragment$())

      PROC_shaderinit(oVertex%, oFragment%)
      PROC_compileshader(oVertex%, Vertex$(), "Vertex")
      PROC_compileshader(oFragment%, Fragment$(), "Fragment")
      oProgram% = FN_useshaders(oVertex%, oFragment%)

      REPEAT : PROC_render : WAIT 2 : UNTIL FALSE
      END

      REM Vertex shader:
      DATA "attribute vec4 vPosition;"
      DATA "void main()"
      DATA "{"
      DATA "gl_Position = vPosition;"
      DATA "}"
      DATA ""

      REM Fragment shader:
      DATA "void main()"
      DATA "{"
      DATA "float v, ci = 0.5213, cr = -0.5125;"
      DATA "float zr = gl_FragCoord.x / 639.0 * 3.0 - 1.5;"
      DATA "float zi = gl_FragCoord.y / 511.0 * 2.0 - 1.0;"
      DATA "for( int i = 0; i <= 510; i++ )"
      DATA " {"
      DATA " v = float(i) / 512.0;"
      DATA " float zn = (zr + zi) * (zr - zi) + cr;"
      DATA " zi = 2.0 * zr * zi + ci; zr = zn;"
      DATA " if ((zr * zr + zi * zi) > 3.0) break;"
      DATA " }"
      DATA "gl_FragColor = vec4(sqrt(v), v, v*v, 1.0 );"
      DATA "}"
      DATA ""

julia.jpg

[Richard Russell]

For comparison, here's a native version. It will run in BB4W, BBCSDL and Matrix Brandy but probably more than 1000 times slower than the shader version! To run it in Acorn's ARM BASIC 5 change 'EXIT FOR' to 'GOTO 200'.

Code: Select all

   10 MODE 9
   20 FOR K% = 0 TO 15
   30   PROCpalette(K%, K%/15)
   40 NEXT
   50 PROCjulia(639, 511)
   60 END
   70
   80 DEF PROCjulia(W%, H%)
   90 LOCAL I%, X%, Y%, zr, zi, zn, ci, cr
  100 ci = 0.5213 : cr = -0.5125
  110 FOR Y% = 0 TO H%
  120   FOR X% = 0 TO W%
  130     zr = X% / W% * 3 - 1.5
  140     zi = Y% / H% * 2 - 1.0
  150     FOR I% = 0 TO 510
  160       zn = (zr + zi) * (zr - zi) + cr
  170       zi = 2 * zr * zi + ci : zr = zn
  180       IF zr * zr + zi * zi > 3 EXIT FOR
  190     NEXT
  200     GCOL I% DIV 32
  210     LINE X%*2, Y%*2, X%*2, Y%*2
  220   NEXT X%
  230 NEXT Y%
  240 ENDPROC
  250
  260 DEF PROCpalette(K%, v)
  270 COLOUR K%, 255*SQR(v), 255*v, 255*v*v
  280 ENDPROC

[Richard Russell]

Note the similarity between the BBC BASIC code for calculating the Julia Set:

Code: Select all

      FOR I% = 0 TO 510
        zn = (zr + zi) * (zr - zi) + cr
        zi = 2 * zr * zi + ci : zr = zn
        IF zr * zr + zi * zi > 3 EXIT FOR
      NEXT

and the shader code to do the same thing:

Code: Select all

      for (int i = 0; i <= 510; i++)
        {
        float zn = (zr + zi) * (zr - zi) + cr;
        zi = 2.0 * zr * zi + ci; zr = zn;
        if ((zr * zr + zi * zi) > 3.0) break;
        }

[DDRM]

Hi Richard,

That's very nice! I'd point out, however, that there isn't "A" Julia set, there's one for each point in the complex plane (here you've chosen the point ci = 0.5213 : cr = -0.5125), though the interesting ones lie on the fringes of the Mandelbrot set, say in the range Cr = -2 to 2 and Ci =-1.25 to 1.25.

Is there a way (using shaderlib) to pass these parameters to the shader in real time, rather than compiling them as hard coded values? Obviously I could modify the relevant strings in the Fragment$ array, and then recompile, but that's not really an optimal solution.

I note that shaderlib doesn't appear to be documented in either the included or online manual (though it's existence is noted). I had a search on the wiki, and found some old pages of yours from 2015, but they don't use shaderlib... Looking at the code of the library, I can't see where it is possible to pass values into the shader [pipline] when it is called.

You probably don't remember (it was in 2016 and a fringe interest even then!) that I wrote a Julia viewer using assembler, which allowed the user to move the mouse about in the window (which represents the space x(real) = -2 to 2, y(imaginary) = -1.25 to 1.25 ). Your shader code should be able to achieve much the same, but in a much easier-to-write and read form.

My assembler code is below: run and move your mouse about over the window to see the Julia set corresponding to the underlying point (in the space, not the set!) - that's the kind of thing I'd like to be able to achieve.

Best wishes,

D

Code: Select all

      MODE 8
      xres=@vdu%!208
      yres=@vdu%!212

      minx=-2.0
      maxx=2.0
      xrange=maxx-minx
      xstep=xrange/xres
      miny=-1.25
      maxy=1.25
      yrange=maxy-miny
      ystep=yrange/yres

      REM Reserve space for the assembly routine itself, making sure it's in its own 2K block to avoid cache thrashing.
      size% = 2048
      DIM code% NOTEND AND 2047, code% size%-1
      REM These are just dummy variables, while setting up the assembler
      nr=3.5
      ni=1.5
      cr=0.29
      ci=0.01
      mag%=0
      temp=0.0
      col%=TRUE
      fixed%=FALSE

      PROCAssJulia

      REM Set up a DIBsection for screen output, so we can write data to it directly

      DIM BITMAPINFOHEADER{Size%, Width%, Height%, Planes{l&,h&}, BitCount{l&,h&}, \
      \                    Compression%, SizeImage%, XPelsPerMeter%, YPelsPerMeter%, \
      \                    ClrUsed%, ClrImportant%}

      DIM bmi{Header{} = BITMAPINFOHEADER{}, Palette%(255)}

      bmi.Header.Size% = DIM(BITMAPINFOHEADER{})
      bmi.Header.Width% = @vdu%!208
      bmi.Header.Height% = @vdu%!212
      bmi.Header.Planes.l& = 1
      bmi.Header.BitCount.l& = 8

      REM We've made an 8 bit per pixel bitmap definition: define a  palette to go with it

      SYS "CreateDIBSection", @memhdc%, bmi{}, 0, ^bits%, 0, 0 TO hbitmap%
      IF hbitmap% = 0 ERROR 100, "Couldn't create DIBSection"

      SYS "SelectObject", @memhdc%, hbitmap% TO oldhbm%
      SYS "DeleteObject", oldhbm%

      PROCSetColour(col%)

      CLS
      bytesperpixel% = bmi.Header.BitCount.l& DIV 8
      bytesperline% = ((bmi.Header.Width% * bytesperpixel%) + 3) AND -4

      REPEAT
        MOUSE x%,y%,z%
        IF NOT fixed% THEN
          cr=minx+x%*xstep/2
          ci=miny+y%*ystep/2
          PRINT TAB(0,0);cr,ci
        ELSE
          IF z%>0 THEN
            tx%=x%
            ty%=y%
            REPEAT MOUSE x%,y%,z% UNTIL z%=0
            minx+=(tx%-x%)*xstep/2
            maxx+=(tx%-x%)*xstep/2
            miny+=(ty%-y%)*ystep/2
            maxy+=(ty%-y%)*ystep/2
          ENDIF
        ENDIF
        CALL code%
        SYS "InvalidateRect", @hwnd%, 0, 0
        q$=INKEY$(0)
        CASE q$ OF
          WHEN "c","C":col%=NOT col%:PROCSetColour(col%)
          WHEN "f","F": fixed%=NOT fixed%
          WHEN CHR$(140):
            fixed%=TRUE
            minx+=xrange/4
            maxx-=xrange/4
            xrange=maxx-minx
            xstep=xrange/xres
            miny+=yrange/4
            maxy-=yrange/4
            yrange=maxy-miny
            ystep=yrange/yres
          WHEN CHR$(141):
            fixed%=TRUE
            minx-=xrange/4
            maxx+=xrange/4
            xrange=maxx-minx
            xstep=xrange/xres
            miny-=yrange/4
            maxy+=yrange/4
            yrange=maxy-miny
            ystep=yrange/yres
        ENDCASE
  
      UNTIL q$="q" OR q$="Q"
      QUIT
      :
      DEFPROCSetColour(col%)
      LOCAL i%,r%,g%,b%
      IF col% THEN
        FOR i% = 0 TO 255
          r% = (i% MOD 16)*16
          g% = (i% MOD 64)*4
          b% = 255-i%/2
          bmi.Palette%(i%) = b% + (g% << 8) + (r% << 16)
        NEXT
      ELSE
        FOR i% = 0 TO 255
          bmi.Palette%(i%) = i%  + (i% << 8) + (i% << 16)
        NEXT
      ENDIF
      SYS "SetDIBColorTable", @memhdc%, 0, 256, ^bmi.Palette%(0)

      ENDPROC
      :
      DEFPROCAssJulia
      LOCAL opt%
      FOR opt%=0 TO 2 STEP 2
        P%=code%
        [
        OPT opt%
        mov esi,[^yres]
        mov edx,[^bits%]
        fld tbyte [^miny]
        fstp tbyte [^ni]
        .yloop
        mov ecx,[^xres]
        fld tbyte [^minx]
        fstp tbyte [^nr]
        .xloop
        mov eax,0
        fld tbyte [^ci]
        fld tbyte [^cr]
        fld tbyte [^nr]
        fst st3
        fmul st0,st0
        fstp st4      ;After this we have cr,ci,nr and nr2 in ST0-3, and 2 pushes
        fld tbyte [^ni]
        fst st5
        fmul st0,st0
        fstp st6      ;After this we have cr,ci,nr, nr2,ni and ni2 in ST0-5, and 2 pushes
        .lpt
        ;OK now we should be able to calculate the new values of nr, ni
        fld st2
        fadd st0,st0
        fmul st0,st5
        fadd st0,st2
        fstp st5     ;we should be in the same place, but with the new value of ni
        fld st3
        fsub st0,st6
        fadd st0,st1
        fstp st3     ;After this we now have cr,ci, the new nr, the old nr2, the new ni, and the old ni2 in ST0-5, and 2 pushes
        fld st2
        fmul st0,st0
        fstp st4    ;now with new nr2, still 2 pushes
        fld st4
        fmul st0,st0
        fst st6    ;now with new ni2, now 3 pushes
        fadd st0,st4  ;gives magnitude squared
        fistp dword [^mag%]  ;store it as an integer, so it can be read into a standard register and compared, Back to 2 pushes
        inc eax
        mov ebx,[^mag%]
        cmp ebx,4
        ja esc
        cmp eax,255
        jb lpt
        .esc
        ;Now we have the problem of popping the stack twice
        fstp tbyte [^temp]     ;should pop cr into temp
        fstp tbyte [^temp]     ;should pop ci into temp
        mov byte [edx],al
        add edx,[^bytesperpixel%]
        fld tbyte [^nr]
        fld tbyte [^xstep]
        faddp st1,st0
        fstp tbyte [^nr]
        dec ecx
        jnz near xloop
        fld tbyte [^ni]
        fld tbyte [^ystep]
        faddp st1,st0
        fstp tbyte [^ni]
        dec esi
        jnz near yloop
        ret
        ]
      NEXT opt%
      ENDPROC

[Richard Russell]

there isn't "A" Julia set, there's one for each point in the complex plane (here you've chosen the point ci = 0.5213 : cr = -0.5125)

I know, but it's not a point that I chose but rather the point which is (probably) most commonly used when illustrating the Julia Set. For example it's used in the very first example illustrated at the relevant Wikipedia page and also in this post from the Facebook BASIC Programming Language group (which is where my code, and the specific palette used, originated).

Is there a way (using shaderlib) to pass these parameters to the shader in real time, rather than compiling them as hard coded values?

I stated in my post that "I've omitted the code needed to allow real-time panning-and-zooming. This can easily be transplanted from the supplied example program mandel.bbc into this one if you wish" and of course any other values you might want to pass at run-time can use exactly the same technique. Something else most of the shader examples do is pass the window dimensions at run time so that if you re-size (or maximize) it the rendering will automatically adjust to suit.

I note that shaderlib doesn't appear to be documented in either the included or online manual (though it's existence is noted).

All shaderlib does is to allow you to compile and run the shaders, and trigger the rendering cycle. It doesn't have any rôle in what actually happens at run-time, either in respect of the shader code itself or the transfer of parameters from BASIC into the shader code. As such, in my judgement, its use is best 'documented' by example, and there are several such examples supplied: fluid.bbc, mandel.bbc, raytrace.bbc, seascape,bbc, slitscan.bbc and voronoi.bbc.

Whenever BBC BASIC interacts with a complex external library, one cannot expect the supplied documentation to do anything other than scratch the surface of what is possible. So it is with shader code, but I hope the Julia Set example, having been stripped of all the ancillary stuff such as run-time controls, shows that at heart it's not complicated.

[DDRM]

Hi Richard,

OK, thanks. I looked at mandel.bbc, in my BB4W examples folder, but it's the old assembler version. I've now looked in the BBC_SDL examples folder, and I see that's the shader-based version. I'll study that.

Best wishes,

D

Edit,having spent some time with the code: Ah, yes, it's beginning to come back to me! I spent quite a bit of time on shader code (in C) a few years ago, but not using it has led it to fade a bit in my memory...

I see you've included @hwnd% as a parameter for your calls to glGetUniformLocation and glUniform1fv, even though the specification doesn't seem to require it, and taking it out doesn't seem to harm the program - is there a reason for it?

[DDRM]

OK, had a bit of a play, and got it working: see code box if interested! I've made only small changes, basically to add two variables to the shader code, and read the mouse position within the window and sent those coordinates to the routine.

Any idea why you need to pass the values in the form of a structure? I tried assigning it (using FN_f4) to an integer and passing that, but it crashed!

Best wishes,

D

Code: Select all

      MODE 9
      REM MODE 9 is 640 x 512 pixels, or 1280 x 1024 graphics units
      xres%=1280
      yres%=1024
      INSTALL @lib$ + "shaderlib"
      DIM Vertex$(10), Fragment$(20),Float{0%,1%}
      PROC_readshader(Vertex$())
      PROC_readshader(Fragment$())

      PROC_shaderinit(oVertex%, oFragment%)
      PROC_compileshader(oVertex%, Vertex$(), "Vertex")
      PROC_compileshader(oFragment%, Fragment$(), "Fragment")
      oProgram% = FN_useshaders(oVertex%, oFragment%)

      REM Here we store the locations of our shader variables ("uniforms")
      REM in BBC BASIC variables (pReal and pImag)
      SYS `glGetUniformLocation`, oProgram%, "cr" TO pReal%
      SYS `glGetUniformLocation`, oProgram%, "ci" TO pImag%
[img]D:\Basic programmes\Things from forum[/img]
      REPEAT :
        REM We get the mouse coordinates, and check it's in the window area
        REM (That way we know the calculated start point will lie within (-2 to 2,-1.25 to 1.25) )
        REM Though I think the shader actually only shows the area (-1.5 to 1.5, -1 to 1)?
        MOUSE px%,py%,pz%
        IF px%>0 AND px%<xres% AND py%>0 AND py%<yres% THEN
          REM Scale the mouse X coordinate into the required space. Not sure why it needs to go into a structure...
          Float.0% = FN_f4(4*px%/xres% - 3)
          SYS `glUniform1fv`, pReal%, 1, Float{}
          REM Ditto for Y coordinate
          Float.0% = FN_f4(2.5*py%/yres% - 1.25)
          SYS `glUniform1fv`, pImag%, 1, Float{}
          PRINT TAB(0,0),4*px%/1280 - 2,4*py%/1024 - 1.25
        ENDIF
  
        PROC_render
        WAIT 2
      UNTIL FALSE
      END

      REM Vertex shader:
      DATA "attribute vec4 vPosition;"
      DATA "void main()"
      DATA "{"
      DATA "gl_Position = vPosition;"
      DATA "}"
      DATA ""

      REM Fragment shader:
      REM We'll define two parameters, representing
      REM the real and imaginary parts of the start location
      DATA "uniform float cr;"
      DATA "uniform float ci;"

      DATA "void main()"
      DATA "{"
      DATA "float v;"
      DATA "float zr = gl_FragCoord.x  * 3.0 / 639.0 - 1.5;"
      DATA "float zi = gl_FragCoord.y * 2.0 / 511.0  - 1.0;"
      DATA "for( int i = 0; i <= 510; i++ )"
      DATA " {"
      DATA " v = float(i) / 512.0;"
      DATA " float zn = (zr + zi) * (zr - zi) + cr;"
      DATA " zi = 2.0 * zr * zi + ci; zr = zn;"
      DATA " if ((zr * zr + zi * zi) > 3.0) break;"
      DATA " }"
      DATA "gl_FragColor = vec4(sqrt(v), v, v*v, 1.0 );"
      DATA "}"
      DATA ""

[Richard Russell]

I see you've included @hwnd% as a parameter for your calls to glGetUniformLocation and glUniform1fv, even though the specification doesn't seem to require it, and taking it out doesn't seem to harm the program - is there a reason for it?

Do you mean @memhdc%? I'm finding that it is necessary here (BBC BASIC for SDL 2.0, Windows 11); if I delete it mandel.bbc gives a black screen.

More generally, you can only conclude that "taking it out doesn't seem to harm the program" if you've tried it on every supported platform, which I'm pretty confident you won't have!

Even I struggle to keep a sample of every platform working for testing purposes, my last 32-bit x86 Android device bit the dust some while ago so I'm unable to test on that platform now.

[DDRM]

You're right, I mean @memhdc% - I'd deleted them all by then, and was too lazy to check back to your original code...

You're also right that I haven't tested it on any other platform (just BB4W on a Windows machine). However, the OpenGL manual pages (on khronos.org) don't suggest that there is a parameter available for it:

GLint glGetUniformLocation( GLuint program, const GLchar *name);
void glUniform1fv( GLint location, GLsizei count, const GLfloat *value);

...and it works for me (here on this setup)...

Interesting that the coordinates show up on the screen grab, but not "in real life": presumably they are being overwritten too fast, or something. I was rather surprised it worked at all (since it isn't going into the OpenGL environment).

Best wishes,

D

Edit: I wonder if it would be interesting to use the mouse scroll wheel (or even the buttons, changing by a power of 2) to adjust the maximum depth of search, and the mapping to the colour palette: a lot of the area is all in dark colours, and it would be nice to be able to bring them out. The bits inside the Mandelbrot set will always be white!

[Richard Russell]

Any idea why you need to pass the values in the form of a structure?

I doubt that you do, the shader shouldn't be able to tell whether it is a structure or not. But structures are guaranteed to be DWORD-aligned (i.e. at a multiple-of-four address) which regular variables aren't, so it's entirely possible that what you are seeing is an alignment dependence, not a structure dependence.

You could try forcing your integer variable to be aligned (most straightforwardly by using one of the static variables A%-Z%) to see if that works. If it does it lends weight to the alignment theory. Acorn's ARM BASIC 5 aligns even regular variables I think, because the original ARM required it, but x86 and modern ARMs don't so my BASICs don't bother.

Very occasionally you may find an API function mandating a 16-byte alignment (it's not common, but it happens) and there's no way of forcing that alignment in BBC BASIC other than doing it explicitly by address manipulation, e.g. an AND -16 operation. That can be a pain.