mode 1
autosave
'==================================
' Gabriel Gambetta's ZX Spectrum Raytracer
' from <https://gabrielgambetta.com/zx-raytracer.html>
' Ported to PicoMiteVGA by Vegipete
'
mode 2 : CLS
Dim ctab(15):For i = 0 To 15: Read ctab(i) : Next
Data RGB(BLACK),RGB(BLUE),RGB(GREEN),RGB(CYAN),RGB(RED)
Data RGB(MAGENTA),RGB(YELLOW),RGB(WHITE),RGB(MYRTLE)
Data RGB(COBALT) ,RGB(MIDGREEN),RGB(CERULEAN)
Data RGB(RUST), RGB(FUCHSIA),RGB(BROWN),RGB(LILAC)

' ===== Initialize scene =====
LX = -1/SQR(2)
LY =  1/SQR(2)
LZ =  0
AI = 0.075
DI = 1 - AI
RESTORE sphere
READ NS
DIM SPH(NS,5)
FOR S = 1 TO NS
 FOR I = 1 TO 5
   READ SPH(S, I)
 NEXT I
NEXT S

' ===== Initialize 8x8 Bayer matrix =====
DIM HBay(64)
RESTORE Bayer
FOR I = 1 TO 64
 READ HBay(I):  HBay(I) = HBay(I) / 64
NEXT I

DIM Cary(64), Lary(64), Aary(8)

timer = 0
DZ = 1
FOR X = 0 TO 312 STEP 8
 FOR Y = 0 TO 175 STEP 8

   ' --- Fast approximate path:
   ' if all 4 corners of this 8x8 block black, ignore the block ---
   DX = (X - 128) / 256 : DY = (Y - 88) / 256
   TraceRay
   CTL = PC : LTL = PL
   
   DX = (X - 128 + 7) / 256
   TraceRay
   CTR = PC : LTR = PL
   
   DY = (Y - 88 + 7) / 256
   TraceRay
   CBR = PC : LBR = PL
   
   DX = (X - 128) / 256
   TraceRay
   CBL = PC : LBL = PL
   
   IF CTL = 0 AND CTR = 0 AND CBR = 0 AND CBL = 0 THEN continue for '(Y)

   ' --- For each 8x8 block, collect the pixel colors and their counts ---
   CI = 1
   PL = 0

   FOR U = X TO X+7
     DX = (U - 128) / 256

     FOR V = Y TO Y+7
       select case CI
         case  1
           PC = CTL:  PL = LTL
         case  8
           PC = CBL:  PL = LBL
         case 57
           PC = CTR:  PL = LTR
         case 64
           PC = CBR:  PL = LBR
         case else
           DY = (V - 88) / 256
           TraceRay
       end select

       Aary(PC+1) = Aary(PC+1) + 1
       Cary(CI) = PC
       Lary(CI) = PL
       CI = CI + 1

     NEXT V
   NEXT U

   ' --- Find most frequent color in this 8x8 block ---
   MFC = 0
   FOR C = 2 TO 8
     IF Aary(C) > MFC THEN  MFC = Aary(C):  MFI = C
   NEXT C
   FC = MFI - 1
   
   ' --- PLOT the non-zero pixels if they're below the dithering threshold --
   CI = 1
   FOR U = X TO X+7
     FOR V = Y TO Y+7
       IF Cary(CI) > 0 AND HBay(CI) <= Lary(CI) THEN pixel U,220-V,ctab(FC)
       CI = CI + 1
     NEXT V
   NEXT U
 NEXT Y
 
 print @(10,225) timer "      ";
NEXT X

end


' ===== TraceRay =====
' Params: (DX, DY, DZ): ray direction
' Returns: PC: pixel color; PL: pixel light intensity
' Optimizations: ray origin hardcoded to (0, 0, 0); (TMIN, TMAX) hardcoded to (0, +inf)
' Clobbers: A, CS, B, C, D, IX, IY, IZ, MT, NL, NX, NY, NZ, S, T
sub TraceRay
 local H = 0
 MT = 1E10
 A = 2*(DX*DX + DY*DY + DZ*DZ)

 FOR S = 1 TO NS
   B = 2*(DX*SPH(S,1) + DY*SPH(S,2) + DZ*SPH(S,3))
   C = SPH(S,1)*SPH(S,1) + SPH(S,2)*SPH(S,2) + SPH(S,3)*SPH(S,3) - SPH(S,4)

   D = B*B - 2*A*C
   IF D < 0 THEN continue for
   D = SQR(D)

   T = (B + D) / A
   IF T > 0 AND T < MT THEN  PC = SPH(S, 5) :  MT = T :  CS = S
   T = (B - D) / A
   IF T > 0 AND T < MT THEN  PC = SPH(S, 5) :  MT = T :  CS = S
 NEXT S

 IF MT = 1E10 THEN  PC = 0: exit sub

 IX = DX*MT:  IY = DY*MT:  IZ = DZ*MT
 NX = IX - SPH(CS, 1):  NY = IY - SPH(CS, 2):  NZ = IZ - SPH(CS, 3)
 PL = AI

 CheckShadow(H)
 IF H = 1 THEN exit sub

 NL = (NX*LX + NY*LY + NZ*LZ)
 IF NL > 0 THEN  PL = PL + DI * NL / SQR(NX*NX + NY*NY + NZ*NZ)

end sub


' ----- Specialized TraceRay for shadow checks -----
' Params: (IX, IY, IZ): ray start; (LX, LY, LZ): ray direction (directional light vector)
' Returns: H = 1 if the ray intersects any sphere, H = 0 otherwise
' Optimizations: (TMIN, TMAX) hardcoded to (epsilon, +inf)
sub CheckShadow(H)
 A = 2*(LX*LX + LY*LY + LZ*LZ)

 FOR S = 1 TO NS
   CX = IX - SPH(S,1):  CY = IY - SPH(S,2):  CZ = IZ - SPH(S,3)
   B = -2*(CX*LX + CY*LY + CZ*LZ)
   C = (CX*CX + CY*CY + CZ*CZ) - SPH(S, 4)
 
   D = B*B - 2*A*C
   IF D < 0 THEN continue for
   D = SQR(D)
 
   T = (B + D) / A
   IF T > 0.01 THEN  H = 1: exit sub
   T = (B - D) / A
   IF T > 0.01 THEN  H = 1: exit sub
 
 NEXT S
 H = 0
end sub

Bayer:
DATA  0, 32,  8, 40,  2, 34, 10, 42
DATA 48, 16, 56, 24, 50, 18, 58, 26
DATA 12, 44,  4, 36, 14, 46,  6, 38
DATA 60, 28, 52, 20, 62, 30, 54, 22
DATA  3, 35, 11, 43,  1, 33,  9, 41
DATA 51, 19, 59, 27, 49, 17, 57, 25
DATA 15, 47,  7, 39, 13, 45,  5, 37
DATA 63, 31, 55, 23, 61, 29, 53, 21

' ----- Sphere data -----
' Sphere count, followed by (CX, CY, CZ, SR*SR, SC)
sphere:
DATA 4
DATA  0, -1, 4, 1, 2
DATA  2,  0, 4, 1, 1
DATA -2,  0, 4, 1, 4
DATA 0, -5001, 0, 5000^2, 6

mode 1
autosave
'==================================
' Gabriel Gambetta's ZX Spectrum Raytracer
' from <https://gabrielgambetta.com/zx-raytracer.html>
' Ported to PicoMiteVGA by Vegipete
'   Simple changes involved removing line numbers, let, goto, gosub,
'   altering program flow/structure to suit MMBasic and adjusting
'   variable names.
'   Deeper changes involved removing code to handle the ZX Spectrum
'   single-colour-per-8x8-block limitations.
'
MODE 2 : CLS
Dim ctab(15):For i = 0 To 15: Read ctab(i) : Next
Data RGB(BLACK),RGB(BLUE),RGB(GREEN),RGB(CYAN),RGB(RED)
Data RGB(MAGENTA),RGB(YELLOW),RGB(WHITE),RGB(MYRTLE)
Data RGB(COBALT) ,RGB(MIDGREEN),RGB(CERULEAN)
Data RGB(RUST), RGB(FUCHSIA),RGB(BROWN),RGB(LILAC)

Dim DX,DY,DZ  ' ray direction
Dim PC,PL     ' pixel color, intensity

' ===== Initialize scene =====
LX =  1/Sqr(2)  ' direction to light source from 0,0,0
LY =  1/Sqr(2)
LZ =  0

AI = 0.075
DI = 1 - AI
Restore sphere
Read NS
Dim SPH(NS,5)
For S = 1 To NS
 For I = 1 To 5
   Read SPH(S, I)
 Next I
Next S

' ===== Initialize 8x8 Bayer matrix =====
Dim HBay(64)
Restore Bayer
For I = 1 To 64
 Read HBay(I) : HBay(I) = HBay(I) / 64
Next I

'DIM Cary(64), Lary(64), Aary(16)

Timer = 0
DZ = 1
For X = 0 To 312 Step 8
 For Y = 0 To 175 Step 8

   ' --- Fast approximate path:
   ' if all 4 corners of this 8x8 block black, ignore the block ---
   DX = (X - 128) / 256 : DY = (Y - 88) / 256
   TraceRay
   CTL = PC : LTL = PL

   DX = (X - 128 + 7) / 256
   TraceRay
   CTR = PC : LTR = PL

   DY = (Y - 88 + 7) / 256
   TraceRay
   CBR = PC : LBR = PL

   DX = (X - 128) / 256
   TraceRay
   CBL = PC : LBL = PL

   If CTL = 0 And CTR = 0 And CBR = 0 And CBL = 0 Then Continue For '(Y)

   ' --- For each 8x8 block, collect the pixel colors and their counts ---
   CI = 1
   PL = 0

   For U = X To X+7
     DX = (U - 128) / 256

     For V = Y To Y+7
       Select Case CI
         Case  1
           PC = CTL : PL = LTL
         Case  8
           PC = CBL : PL = LBL
         Case 57
           PC = CTR : PL = LTR
         Case 64
           PC = CBR : PL = LBR
         Case Else
           DY = (V - 88) / 256
           TraceRay
       End Select

       ' --- PLOT the non-zero pixels if above the dithering threshold --
       If PC > 0 And HBay(CI) <= PL Then Pixel U,220-V,ctab(PC)

       Inc CI
     Next V
   Next U
 Next Y

 Print @(10,225) Str$(Timer/1000,7,3);
Next X

End

' ===== TraceRay =====
' Params: (DX, DY, DZ): ray direction
' Returns: PC: pixel color; PL: pixel light intensity
' Optimizations: ray origin hardcoded to (0, 0, 0);
'                (TMIN, TMAX) hardcoded to (0, +inf)
' Clobbers: A, CS, B, C, D, IX, IY, IZ, Tmax, NL, NX, NY, NZ, S, T
Sub TraceRay
 Local S, H = 0
 Tmax = 1E10
 A = 2*(DX*DX + DY*DY + DZ*DZ)

 For S = 1 To NS
   B = 2*(DX*SPH(S,1) + DY*SPH(S,2) + DZ*SPH(S,3))
   C = SPH(S,1)^2 + SPH(S,2)^2 + SPH(S,3)^2 - SPH(S,4)

   D = B*B - 2*A*C
   If D < 0 Then Continue For
   D = Sqr(D)

   T = (B + D) / A
   If T > 0 And T < Tmax Then  PC = SPH(S, 5) : Tmax = T : CS = S
   T = (B - D) / A
   If T > 0 And T < Tmax Then  PC = SPH(S, 5) : Tmax = T : CS = S
 Next S

 If Tmax = 1E10 Then  PC = 0: Exit Sub

 IX = DX*Tmax : IY = DY*Tmax : IZ = DZ*Tmax
 NX = IX - SPH(CS, 1) : NY = IY - SPH(CS, 2) : NZ = IZ - SPH(CS, 3)
 PL = AI

 CheckShadow(H)
 If H = 1 Then Exit Sub

 NL = (NX*LX + NY*LY + NZ*LZ)
 If NL > 0 Then  PL = PL + DI * NL / Sqr(NX*NX + NY*NY + NZ*NZ)

End Sub

' ----- Specialized TraceRay for shadow checks -----
' Params: (IX, IY, IZ): ray start;
'         (LX, LY, LZ): ray direction (directional light vector)
' Returns: H = 1 if the ray intersects any sphere, H = 0 otherwise
' Optimizations: (TMIN, TMAX) hardcoded to (epsilon, +inf)
Sub CheckShadow(H)
 A = 2*(LX*LX + LY*LY + LZ*LZ)

 For S = 1 To NS
   CX = IX - SPH(S,1) : CY = IY - SPH(S,2) : CZ = IZ - SPH(S,3)
   B = -2*(CX*LX + CY*LY + CZ*LZ)
   C = CX*CX + CY*CY + CZ*CZ - SPH(S,4)

   D = B*B - 2*A*C
   If D < 0 Then Continue For
   D = Sqr(D)

   T = (B + D) / A
   If T > 0.01 Then  H = 1 : Exit Sub
   T = (B - D) / A
   If T > 0.01 Then  H = 1 : Exit Sub

 Next S
 H = 0
End Sub

' 8x8 Dithering Block
Bayer:
Data  0, 32,  8, 40,  2, 34, 10, 42
Data 48, 16, 56, 24, 50, 18, 58, 26
Data 12, 44,  4, 36, 14, 46,  6, 38
Data 60, 28, 52, 20, 62, 30, 54, 22
Data  3, 35, 11, 43,  1, 33,  9, 41
Data 51, 19, 59, 27, 49, 17, 57, 25
Data 15, 47,  7, 39, 13, 45,  5, 37
Data 63, 31, 55, 23, 61, 29, 53, 21

' ----- Sphere data -----
' Sphere count, followed by (CX, CY, CZ, SR*SR, SC)
' center x, center y, center z, radius squared, colour
sphere:
Data 4
Data  0, -1, 4, 1, 2
Data  2,  0, 4, 1, 1
Data -2,  0, 4, 1, 4
Data 0, -5001, 0, 5000^2, 6