Ray Traced Spheres

[SimonSideburns]

I am a member of a BASIC Programming on the ZX Spectrum group on Facebook and the members share their ideas, programs, effects, etc.

One of our chaps, Daniel A Nagy has written a short program to render ray-traced spheres above a checker board surface, but with the attribute clash that is prevalent on the ZX Spectrum this was in Monochrome only. Another chap, Paul Dunn has created a modern version of ZX Basic with additional keywords, etc. that runs on Windows and his environment doesn't suffer with colour clash, so he adapted the code to create a more colourful display.

I have taken this code and created a couple of versions of it for the BBC Micro. One for MODE 0 (still in monochrome) and the other for MODE 2 in 8 glorious colours. The upshot of this is that when run (albeit slowly) in an emulator in B/W mode it gives a quite passable greyscale rendered image. The coloured version isn't too bad either.

I thought I'd show the code here, include a disk image with the two programs on it, and include a few screenshots for your delight.

The beauty of the code is that the number of spheres can be adjusted along with their positions within the render, and the rendered image will reflect (no pun intended) the changes made. I'm yet to fully study and work out how this is done, but it does seem to be fairly simple.

In BeebEm it is possible to speed up by around 9x by selecting the 1 FPS option on the Speed menu, but even then you're looking at well over 10 minutes to draw the full screen.

I'm sure there could be tweaks made, and what about the possibility of writing this in machine code instead? Something way above my current skill set I would think. One option would be to utilise a second processor to see what sort of speed increases that would give. I don't have one at the moment so feel free to give it a go and report back your findings.

Here's the MODE 0 screen:

First coloured MODE 2 screen:

First MODE 2 screen in B/W:

Second version in MODE 2 (better shadow contrast?):

Second version in MODE 2 in B/W:

Here's a copy of the code:
5MODE2:VDU 23;8202;0;0;0;
10DATA 4,1,5,2,6,3
12FOR n%=1 TO 6:READ c%:VDU 19,n%,c%,0,0,0:NEXT n%
15DATA 0,6,1,7,5
17DIM cl%(5):FOR n%=1 TO 5:READ cl%(n%):NEXT n%
20s%=4:DIM c(s%,3):DIM r(s%):DIM q(s%)
30FOR k=1 TO s%:READ c(k,1),c(k,2),c(k,3),r:r(k)=r:q(k)=r*r:NEXT k
40DATA -0.8,-1,3.2,0.7
45DATA 0,-0.45,2,0.3
50DATA 1.2,-0.7,2.5,0.5
55DATA 0.4,-1,4,0.4
60FOR i=0 TO 1024 STEP 4:FOR j=0 TO 1280 STEP 8
70x=0.3:y=-0.5:z=0:ba%=3:dx=j-640:dy=512-i:dz=1300:dd=dx*dx+dy*dy+dz*dz
80GOSUB 100:NEXT j:NEXT i
90END
100n=(y>=0 OR dy<=0):IF NOT n THEN s=-y/dy
110FOR k=1 TO s%
120px=c(k,1)-x:py=c(k,2)-y:pz=c(k,3)-z
130pp=px*px+py*py+pz*pz
140sc=px*dx+py*dy+pz*dz
150IF sc<=0 THEN GOTO 200
160bb=sc*sc/dd
170aa=q(k)-pp+bb
180IF aa<=0 THEN GOTO 200
190sc=(SQR bb-SQR aa)/SQR dd:IF sc<s OR n<0 THEN n=k:s=sc
200NEXT k
210IF n<0 THEN GCOL 0,1+(dy*dy/dd)*7:PLOT 69,j,i:n=0:RETURN
220dx=dx*s:dy=dy*s:dz=dz*s:dd=dd*s*s
230x=x+dx:y=y+dy:z=z+dz
240IF n=0 THEN GOTO 300
250nx=x-c(n,1):ny=y-c(n,2):nz=z-c(n,3)
260nn=nx*nx+ny*ny+nz*nz
270l=2*(dx*nx+dy*ny+dz*nz)/nn
280dx=dx-nx*l:dy=dy-ny*l:dz=dz-nz*l
290GOTO 100
300FOR k=1 TO s%
310u=c(k,1)-x:v=c(k,3)-z: IF u*u+v*v<=q(k) THEN ba%=1
320NEXT k
330IF (x-INT x>.5)=(z-INT z>.5) THEN GCOL 0,cl%(ba%) ELSE GCOL 0,cl%(ba%+1)
335PLOT 69,j,i
340RETURN

Here's the floppy image:

Small Type-Ins.dsd

(7.5 KiB) Downloaded 11 times

[danielj]

Very cool! Just a cursory glance suggests replacing variables with integers might speed things up?

d.

[flaxcottage]

If you copy the code into BBC BASIC for Windows it takes less than 2 seconds to render.

The supplied download does not work on a real BBC. All the files report 'Bad program' when loaded, except C.Spheres. The download is a DSD file but IMPORTS as a SSD file into a Datacentre.

I tried the SPHERES program on my Master with the Pi Co-pro but get a Division by zero error at line 100.

[Phlamethrower]

Nice!

One easy optimisation: Move line 130 down to 165 (or just get rid of the pp variable completely, since it's only used once)

[AndyF]

This is excellent. I did the 130 to 165 tweak.

Ran it on BeebEm in 'B mode but using the GFX extension (confirmed active by PAGE being at 7168 instead of 6400) , unsure if it actually made any difference as I ran it at a higher speed anyway but to my mind for these type of demo's I prefer to have it there with DFS (1770) and nothing else. Seems a sensible ROM setup for such things.

[kieranhj]

Nice! I do love me some ray-traced spheres above a checkerboard floor. I wrote a ray-tracer from scratch for my Masters degree back in the 90's, all in C++ mind. Always wondered how much work it would be to write a simpler version in 6502. (Still terribly slow undoubtedly with lots of painful 8-bit maths.)

If you are interested in learning more, this classic text by Andrew Glassner (my original copy is still on my book shelf!) is highly recommended: https://www.glassner.com/portfolio/an-i ... y-tracing/

[DutchAcorn]

Nice! It looks almost exactly like a program my brother wrote for the Arc. He included rotation of the spheres. Saved every image after tracing and made a nice smooth animation from the saved screens. I'll see if I can find the code...

[SimonSideburns]

Not sure what happened with the disk image. I created a blank floppy and saved the files into that disk from within the BeebEm emulator.

I have tried loading a few of the files from within BeebEm and it seems fine to load them.

[lurkio]

This is a great little program!



Here's the BASIC listing, refactored to remove GOTOs, GOSUBs and line numbers (for no real reason):

Code: Select all

NEW
AUTO

MODE2:PROCi

FORI%=0TO1024STEP4:FORJ%=0TO1280STEP8

x=0.3:y=-0.5:z=0:ba%=3:dx=J%-640:dy=512-I%:dz=1300:dd=dx*dx+dy*dy+dz*dz

REPEAT

n=(y>=0 OR dy<=0):IF NOT n s=-y/dy

FORK%=1TOS%
px=c(K%,1)-x:py=c(K%,2)-y:pz=c(K%,3)-z
sc=px*dx+py*dy+pz*dz
IF sc<=0 NEXT ELSE bb=sc*sc/dd:aa=q(K%)-(px*px+py*py+pz*pz)+bb:IF aa<=0 NEXT ELSE sc=(SQR bb-SQR aa)/SQR dd:IF sc<s OR n<0 n=K%:s=sc:NEXT ELSE NEXT

IF n<0 GCOL 0,1+(dy*dy/dd)*7:PLOT 69,J%,I%:n=0:UNTIL1:NEXT,:END

dx=dx*s:dy=dy*s:dz=dz*s:dd=dd*s*s
x=x+dx:y=y+dy:z=z+dz

IF n<>0 nx=x-c(n,1):ny=y-c(n,2):nz=z-c(n,3):nn=nx*nx+ny*ny+nz*nz:l=2*(dx*nx+dy*ny+dz*nz)/nn:dx=dx-nx*l:dy=dy-ny*l:dz=dz-nz*l:UNTIL0

FORK%=1TOS%
u=c(K%,1)-x:v=c(K%,3)-z: IF u*u+v*v<=q(K%) ba%=1
NEXT

IF (x-INT x>.5)=(z-INT z>.5) GCOL 0,cl%(ba%) ELSE GCOL 0,cl%(ba%+1)
PLOT 69,J%,I%

UNTIL1:NEXT,
END

DEFPROCi:VDU 23;8202;0;0;0;
DATA 4,1,5,2,6,3
FORN%=1TO6:READ C%:VDU19,N%,C%,0,0,0:NEXT
DATA 0,6,1,7,5
DIM cl%(5):FORN%=1TO5:READ cl%(N%):NEXT
S%=4:DIM c(S%,3):DIM r(S%):DIM q(S%)
FOR K%=1 TO S%:READ c(K%,1),c(K%,2),c(K%,3),r:r(K%)=r:q(K%)=r*r:NEXT
DATA -0.8,-1,3.2,0.7
DATA 0,-0.45,2,0.3
DATA 1.2,-0.7,2.5,0.5
DATA 0.4,-1,4,0.4
ENDPROC

Here's a link to run the program online in JSBeeb (with a CPU acceleration-factor that doesn't actually seem to speed things up very much, for some reason):

• http://bbcmicro.co.uk/jsbeeb/play.php?c ... 9A05E10DFF

And here's a .SSD disc-image with the program on it, plus a compressed version of the program (compressed using the PRES Advanced BASIC Editor ROMs):

• SPHERE.ssd

  DFS .SSD disc-image of ray-traced spheres program

  (200 KiB) Downloaded 8 times

[sbadger]

This is really nice, nice work.
Begging to be done though.....here's a NULA version.

I just added it to Lurkio's disc img. The file is called N.SPHERE

[Lardo Boffin]

It still amazes me what can be achieved with such a relatively small amount of code! Incredible.

[SimonSideburns]

I thought some of you more artistic people would have increased the number of spheres, changed their positions, etc. and rendered a new scene or two. Maybe due to the time it takes to generate each scene that soon would get boring.

I love the nULA version with increased grey shades. Very nice, but how about making the checkerboard in two contrasting colours (say, white and brown) and the sky a fetching set of blues? That would really make it look authentic. I don't have a nULA so someone else would have to do that bit.

Good to have given you all something else to play with though.

I did create a desktop background for my windows laptop at 1600x900 resolution which you are more than welcome to a copy of. This one has 5 spheres, and the colouring is how it was changed to be used in the SpecBAS software (by Paul Dunn) that I mentioned in the original post:

[marcusjambler]

Ray tracing space hopper

[marcusjambler]

and another one....

[kieranhj]

I'm a bit rusty on my ray tracing maths, so I might have missed something, but here is how the program works.

Code: Select all

\\ MODE 2 turn off cursor
5MODE2:VDU 23;8202;0;0;0;

\\ Set colour palette 0-7 from darkest to brightest 
10DATA 4,1,5,2,6,3
12FOR n%=1 TO 6:READ c%:VDU 19,n%,c%,0,0,0:NEXT n%

\\ Checkerfloor colours inc. shadow flag
15DATA 0,6,1,7,5
17DIM cl%(5):FOR n%=1 TO 5:READ cl%(n%):NEXT n%

\\ 4x spheres each with centre c(n)(x,y,z), radius: r(n), compute radius squared: q(n)
20s%=4:DIM c(s%,3):DIM r(s%):DIM q(s%)

\\ Read data for spheres x,y,z,radius
30FOR k=1 TO s%:READ c(k,1),c(k,2),c(k,3),r:r(k)=r:q(k)=r*r:NEXT k
40DATA -0.8,-1,3.2,0.7
45DATA 0,-0.45,2,0.3
50DATA 1.2,-0.7,2.5,0.5
55DATA 0.4,-1,4,0.4

\\ Loop over y pixels, loop over x pixels
60FOR i=0 TO 1024 STEP 4:FOR j=0 TO 1280 STEP 8

\\ Set ray origin to camera position x,y,z = (0.3,-0.5,0)
\\ Camera always looks down +z axis
\\ Compute ray direction vector d through screen pixel (j,i) at distance 1300
\\ FOV ~60 degrees
\\ Compute dd = length of vector d

70x=0.3:y=-0.5:z=0:ba%=3:dx=j-640:dy=512-i:dz=1300:dd=dx*dx+dy*dy+dz*dz

\\ Compute pixel colour and plot, loop to next pixel
80GOSUB 100:NEXT j:NEXT i

\\ End
90END

\\ Fixed infinite floor plane at (0,0,0) parallel to x,z axes
\\ We seem to have -ve y pointing upwards?
\\ If camera is below the plane then ray will always hit it
\\ If ray points downwards then it will always intersect somewhere
\\ If intersection then set object n=0 and compute intersection distance s = -y/dy
100n=(y>=0 OR dy<=0):IF NOT n THEN s=-y/dy

\\ Loop over all spheres
110FOR k=1 TO s%
\\ Calculate vector p from ray origin to centre of sphere
120px=c(k,1)-x:py=c(k,2)-y:pz=c(k,3)-z
\\ Calculate pp = length of p
130pp=px*px+py*py+pz*pz
\\ Calculate sc = p dot product d = length of projection of p onto ray direction d
140sc=px*dx+py*dy+pz*dz
\\ If sc <= 0 then ray points away from sphere so no intersection - goto next sphere
150IF sc<=0 THEN GOTO 200
\\ Compute bb = distance along ray to where centre of sphere is projected
160bb=sc*sc/dd
\\ Compute aa = radius squared - squared distance to sphere centre + bb
170aa=q(k)-pp+bb
\\ If aa <= 0 then ray misses sphere - goto next sphere
180IF aa<=0 THEN GOTO 200
\\ Calculate sc = distance along ray that intersects sphere
\\ If sc < current closest distance s then set object n = sphere k
190sc=(SQR bb-SQR aa)/SQR dd:IF sc<s OR n<0 THEN n=k:s=sc
\\ Next sphere
200NEXT k

\\ If n FALSE then ray has not collided with floor or sphere so compute colour for sky, plot and return for next pixel
\\ Intensity = squared length of ray y component / squared length of ray
\\ Simulates a directional light source coming directly from above
210IF n<0 THEN GCOL 0,1+(dy*dy/dd)*7:PLOT 69,j,i:n=0:RETURN

\\ Shorten ray direction vector d to length of intersection distance
\\ Compute dd = new length of d
220dx=dx*s:dy=dy*s:dz=dz*s:dd=dd*s*s

\\ Calculate new ray origin as intersection point = old origin + direction ray
230x=x+dx:y=y+dy:z=z+dz

\\ If object hit was floor compute colour
240IF n=0 THEN GOTO 300

\\ Object was a sphere
\\ Compute vector n = sphere normal at intersection point = intersection point - sphere centre
250nx=x-c(n,1):ny=y-c(n,2):nz=z-c(n,3)
\\ Compute nn = length of n
260nn=nx*nx+ny*ny+nz*nz
\\ Compute incident value l = 2 * d dot product n / length of n^2
270l=2*(dx*nx+dy*ny+dz*nz)/nn
\\ Spheres have perfect reflection
\\ Compute new ray direction d = old ray direction reflected around normal = d - n * incident value
280dx=dx-nx*l:dy=dy-ny*l:dz=dz-nz*l
\\ Trace new ray with origin on sphere
290GOTO 100
\\ This will repeat forever until either the floor or the sky is hit!

\\ Compute floor colour
\\ For each sphere
300FOR k=1 TO s%
\\ Project vector from intersection point to centre of sphere onto 2D plane x,z
\\ If length of that vector is less than radius of sphere then we are in shadow
\\ Simulates a direction light source coming from directly above
310u=c(k,1)-x:v=c(k,3)-z: IF u*u+v*v<=q(k) THEN ba%=1
320NEXT k

\\ Ray hit floor
\\ Round off fraction from x,z coordinates of intersection point
\\ If both fractions are > 0.5 or both <= 0.5 then choose colour 0+shadow flag
\\ else choose colour 1+shadow flag - gives checkerboard of size 0.5x0.5
330IF (x-INT x>.5)=(z-INT z>.5) THEN GCOL 0,cl%(ba%) ELSE GCOL 0,cl%(ba%+1)
\\ Plot pixel
335PLOT 69,j,i
\\ Return for next pixel
340RETURN

[jgharston]

I've always wondered what program generated this image, and it was probably something similar.

[RobC]

Great thread - must get this running on a Pi co-pro

I don't have a nULA so someone else would have to do that bit.

If you want to develop for VideoNuLA without using the actual hardware, the B-Em and b2 emulators support it.

I also still have some boards available for those who want the hardware

[kieranhj]

Great thread - must get this running on a Pi co-pro

Selecting 64 MHz 6502 copro in b-em makes a big difference to the speed. Was hoping someone would have chance to try on 300+MHz Pi.

After pondering an assembly implementation we concluded BASIC + copro is probably the way to go due to all the floating point calcs.

[SimonSideburns]

I'm a bit rusty on my ray tracing maths, so I might have missed something, but here is how the program works.

Code: Select all

\\ MODE 2 turn off cursor
5MODE2:VDU 23;8202;0;0;0;

\\ Set colour palette 0-7 from darkest to brightest 
10DATA 4,1,5,2,6,3
12FOR n%=1 TO 6:READ c%:VDU 19,n%,c%,0,0,0:NEXT n%

\\ Checkerfloor colours inc. shadow flag
15DATA 0,6,1,7,5
17DIM cl%(5):FOR n%=1 TO 5:READ cl%(n%):NEXT n%

\\ 4x spheres each with centre c(n)(x,y,z), radius: r(n), compute radius squared: q(n)
20s%=4:DIM c(s%,3):DIM r(s%):DIM q(s%)

\\ Read data for spheres x,y,z,radius
30FOR k=1 TO s%:READ c(k,1),c(k,2),c(k,3),r:r(k)=r:q(k)=r*r:NEXT k
40DATA -0.8,-1,3.2,0.7
45DATA 0,-0.45,2,0.3
50DATA 1.2,-0.7,2.5,0.5
55DATA 0.4,-1,4,0.4

\\ Loop over y pixels, loop over x pixels
60FOR i=0 TO 1024 STEP 4:FOR j=0 TO 1280 STEP 8

\\ Set ray origin to camera position x,y,z = (0.3,-0.5,0)
\\ Camera always looks down +z axis
\\ Compute ray direction vector d through screen pixel (j,i) at distance 1300
\\ FOV ~60 degrees
\\ Compute dd = length of vector d

70x=0.3:y=-0.5:z=0:ba%=3:dx=j-640:dy=512-i:dz=1300:dd=dx*dx+dy*dy+dz*dz

\\ Compute pixel colour and plot, loop to next pixel
80GOSUB 100:NEXT j:NEXT i

\\ End
90END

\\ Fixed infinite floor plane at (0,0,0) parallel to x,z axes
\\ We seem to have -ve y pointing upwards?
\\ If camera is below the plane then ray will always hit it
\\ If ray points downwards then it will always intersect somewhere
\\ If intersection then set object n=0 and compute intersection distance s = -y/dy
100n=(y>=0 OR dy<=0):IF NOT n THEN s=-y/dy

\\ Loop over all spheres
110FOR k=1 TO s%
\\ Calculate vector p from ray origin to centre of sphere
120px=c(k,1)-x:py=c(k,2)-y:pz=c(k,3)-z
\\ Calculate pp = length of p
130pp=px*px+py*py+pz*pz
\\ Calculate sc = p dot product d = length of projection of p onto ray direction d
140sc=px*dx+py*dy+pz*dz
\\ If sc <= 0 then ray points away from sphere so no intersection - goto next sphere
150IF sc<=0 THEN GOTO 200
\\ Compute bb = distance along ray to where centre of sphere is projected
160bb=sc*sc/dd
\\ Compute aa = radius squared - squared distance to sphere centre + bb
170aa=q(k)-pp+bb
\\ If aa <= 0 then ray misses sphere - goto next sphere
180IF aa<=0 THEN GOTO 200
\\ Calculate sc = distance along ray that intersects sphere
\\ If sc < current closest distance s then set object n = sphere k
190sc=(SQR bb-SQR aa)/SQR dd:IF sc<s OR n<0 THEN n=k:s=sc
\\ Next sphere
200NEXT k

\\ If n FALSE then ray has not collided with floor or sphere so compute colour for sky, plot and return for next pixel
\\ Intensity = squared length of ray y component / squared length of ray
\\ Simulates a directional light source coming directly from above
210IF n<0 THEN GCOL 0,1+(dy*dy/dd)*7:PLOT 69,j,i:n=0:RETURN

\\ Shorten ray direction vector d to length of intersection distance
\\ Compute dd = new length of d
220dx=dx*s:dy=dy*s:dz=dz*s:dd=dd*s*s

\\ Calculate new ray origin as intersection point = old origin + direction ray
230x=x+dx:y=y+dy:z=z+dz

\\ If object hit was floor compute colour
240IF n=0 THEN GOTO 300

\\ Object was a sphere
\\ Compute vector n = sphere normal at intersection point = intersection point - sphere centre
250nx=x-c(n,1):ny=y-c(n,2):nz=z-c(n,3)
\\ Compute nn = length of n
260nn=nx*nx+ny*ny+nz*nz
\\ Compute incident value l = 2 * d dot product n / length of n^2
270l=2*(dx*nx+dy*ny+dz*nz)/nn
\\ Spheres have perfect reflection
\\ Compute new ray direction d = old ray direction reflected around normal = d - n * incident value
280dx=dx-nx*l:dy=dy-ny*l:dz=dz-nz*l
\\ Trace new ray with origin on sphere
290GOTO 100
\\ This will repeat forever until either the floor or the sky is hit!

\\ Compute floor colour
\\ For each sphere
300FOR k=1 TO s%
\\ Project vector from intersection point to centre of sphere onto 2D plane x,z
\\ If length of that vector is less than radius of sphere then we are in shadow
\\ Simulates a direction light source coming from directly above
310u=c(k,1)-x:v=c(k,3)-z: IF u*u+v*v<=q(k) THEN ba%=1
320NEXT k

\\ Ray hit floor
\\ Round off fraction from x,z coordinates of intersection point
\\ If both fractions are > 0.5 or both <= 0.5 then choose colour 0+shadow flag
\\ else choose colour 1+shadow flag - gives checkerboard of size 0.5x0.5
330IF (x-INT x>.5)=(z-INT z>.5) THEN GCOL 0,cl%(ba%) ELSE GCOL 0,cl%(ba%+1)
\\ Plot pixel
335PLOT 69,j,i
\\ Return for next pixel
340RETURN

Thanks for that explanation. It does make sense. Am I right in thinking that if you were able to place a sphere directly in front of the 'view point' and one behind where you are standing in exactly the right point that the ray would bounce indefinitely? In that case the program would never end.

In that case am I right in thinking it would be better to have an extra flag that would count the number of reflections and once a certain number was reached then say that enough is enough?

[marcusjambler]

If you want to develop for VideoNuLA without using the actual hardware, the B-Em and b2 emulators support it.

Need help setting B-em up for VideoNula - I cant see any settings in the drop downs

Marcus

[RobC]

Need help setting B-em up for VideoNula - I cant see any settings in the drop downs

I don't think there are any dropdowns - it just recognises the extra NuLA register addresses at FE22 & FE23.

I'm using v2.2 but I'm not sure which branches the NuLA support is in.

[DutchAcorn]

Here's the version my brother made. The files on the disk:
RAY_TRACE Does the ray tracing and saves the screens
MERGESCR merges the screendumps into an "ALLSCREENS" sprite file.
SHOWTRACE shows an animation based on the sprites.

Make sure you have at least 800K of sprite size configured.
The ALLSCREENS file is on the disk.

It is written for an ARM2...

Edit: replaced attachment

[IanS]

DEADDEAD is generally not a good sign, looks like you haven't let it finish copying.

[DutchAcorn]

DEADDEAD is generally not a good sign, looks like you haven't let it finish copying.
ray.PNG

That's what you get if you try to do things too fast.. Attachment replaced...

[marcusjambler]

not sure which branches the NuLA support is in.

Thanks Rob I found the pre release in the github repo.

And a quick Nula....

The code does run way quicker on the Pi-copro

Marcus

[RobC]

And a quick Nula....

Looks great. Have you got the code you used to generate it? I've been thinking of putting together a little animation...

[marcusjambler]

Looks great. Have you got the code you used to generate it?

I used the code uploaded by sbadger further up the thread.
And changed the palette setup loop for blue/green from 17 to 15 for the above image

Here's one set to '16'