Author Topic: Raytracing  (Read 4554 times)

ZXDunny

  • Guest
Raytracing
« on: September 30, 2018, 06:45:55 PM »
I know it's probably been done before, but this is kinda fun.

A chap on one the BASIC facebook groups I'm involved in wrote a small Sinclair BASIC raytracer - monochrome (naturally) and took 8 hours to complete.

It pretty much pasted into SpecBAS with only one line needing to be altered (and a bug in my INT implementation to be fixed) before it worked - and a bonus was that it took less than 1 second to produce the same image :)

Get the SpecBAS update here: https://drive.google.com/open?id=0B6gXsz43xir_eV9JbWhUZ1Rxc0k

Here's the SpecBAS code:

Code: [Select]
10 PAPER 0: INK 15: CLS:
 palette 16,255,255,255:
 palette 32,0,192,255:
 palette 255,0,0,192:
 rainbow 16 to 32:
 rainbow 32 to 255
20 read spheres:t=msecs:
 DIM c(spheres,3),r(spheres),q(spheres),cl(4):
 w=scrw/2,h=scrh/2,s=0:
 cl(1)=6,cl(2)=1,
 cl(3)=cl(1)+8,cl(4)=cl(2)+8
30 FOR k=1 TO spheres:
    READ c(k,1),c(k,2),c(k,3),r:
    r(k)=r,q(k)=r*r:
 NEXT k
40 data 6:
 DATA -0.3,-0.8,3,0.6
50 DATA 0.9,-1.4,3.5,0.35:
 data 0.7,-0.45,2.5,0.4:
 data -0.5,-0.3,1.5,0.15:
 DATA 1.0,-0.2,1.5,0.1:
 DATA -0.1,-0.2,1.25,0.2
60 FOR i=1 TO scrh:
    FOR j=0 TO scrw-1
70       x=0.3,y=-0.5,z=0,ba=3:
       dx=j-w,dy=h-i,dz=(scrh/480)*600,
       dd=dx*dx+dy*dy+dz*dz
80       n=-(y>=0 OR dy<=0):
       IF n=0 THEN
          s=-y/dy
90       FOR k=1 TO spheres:
          px=c(k,1)-x,py=c(k,2)-y,pz=c(k,3)-z,
          pp=px*px+py*py+pz*pz,
          sc=px*dx+py*dy+pz*dz:
          IF sc<=0 THEN
             GO TO 120
100          bb=sc*sc/dd,
          aa=q(k)-pp+bb:
          IF aa<=0 THEN
             GO TO 120
110          sc=(SQR bb-SQR aa)/SQR dd:
          IF sc<s OR n<0 THEN
             n=k,s=sc
120       NEXT k
130       IF n<0 THEN
          plot ink 16+(dy*dy/dd)*240;j,scrh-i:
          go to 200
140       dx=dx*s,dy=dy*s,dz=dz*s,dd=dd*s*s,
       x+=dx,y+=dy,z+=dz:
       IF n<>0 THEN
           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:
          GO TO 80
160       FOR k=1 TO spheres:
          u=c(k,1)-x,
          v=c(k,3)-z:
          IF u*u+v*v<=q(k) THEN
             ba=1:
             go to 180
170       NEXT k
180       IF (x-INT x>.5)=(z-INT z>.5) THEN
          ik=cl(ba)
       else
          ik=cl(ba+1)
190       plot ink ik;j,scrh-i       
200    NEXT j:
 NEXT i
210 print at 0,0;transparent 1;ink 0;"Time: ";(msecs-t)/1000;" secs":
 pause 0:

It looks better in the SpecBAS editor, I promise.

And finally, here's the image (I changed the output up to 1920x1024 so it takes almost 30 seconds to make this image):



Had a nice fun day with this one - it's not something I personally would have done in Sinclair BASIC but I'm glad someone did!

Mike Lobanovsky

  • Guest
Re: Raytracing
« Reply #1 on: October 01, 2018, 01:39:49 PM »
Impressive! 8)

Exactly what was the problem with INT? Was it the need to round it rather than truncate or vice versa?

ZXDunny

  • Guest
Re: Raytracing
« Reply #2 on: October 01, 2018, 02:39:01 PM »
Exactly what was the problem with INT? Was it the need to round it rather than truncate or vice versa?

On the Sinclair Spectrum, INT rounds down to negative infinity - 2.8 == 2, -3.6 == -4. In SpecBAS, I was rounding towards zero instead which caused the grid at the bottom (a simple XOR pattern) to fail to draw for any x < 0. So a bonus is that I'm now a little more compatible with the original 80s BASIC.

Mike Lobanovsky

  • Guest
Re: Raytracing
« Reply #3 on: October 01, 2018, 06:56:12 PM »
... 2.8 == 2, -3.6 == -4.

That is, it rounds towards the nearest even whole number. IIRC VB and PB also behave in that manner.

ZXDunny

  • Guest
Re: Raytracing
« Reply #4 on: October 01, 2018, 07:06:32 PM »
... 2.8 == 2, -3.6 == -4.

That is, it rounds towards the nearest even whole number. IIRC VB and PB also behave in that manner.

Almost, as -3.1 will also round to -4.

Mike Lobanovsky

  • Guest
Re: Raytracing
« Reply #5 on: October 02, 2018, 06:38:21 AM »
Almost, as -3.1 will also round to -4.

That's what I'd expect but what will +3.1 round to, +3 or +4?

ZXDunny

  • Guest
Re: Raytracing
« Reply #6 on: October 02, 2018, 08:42:11 AM »
Almost, as -3.1 will also round to -4.

That's what I'd expect but what will +3.1 round to, +3 or +4?

+3.1 rounds to 3 on a real Spectrum.

Mike Lobanovsky

  • Guest
Re: Raytracing
« Reply #7 on: October 02, 2018, 06:37:43 PM »
Ah yes,

Then what Sinclair INT rightfully does is truncation down towards negative infinity. Traditional BASIC truncation towards zero is called FIX.

ZXDunny

  • Guest
Re: Raytracing
« Reply #8 on: October 02, 2018, 06:40:16 PM »
Ah yes,

Then what Sinclair INT rightfully does is truncation down towards negative infinity. Traditional BASIC truncation towards zero is called FIX.

Indeed. I spent quite a lot of time in the z80 ROM disassembly figuring that one out... Surprising how much of a difference the rounding can make especially in something as simple as an XOR operation...

But anyway, it works now and I'm a fraction of a percent closer to the original BASIC so it's a win :)

Mike Lobanovsky

  • Guest
Re: Raytracing
« Reply #9 on: October 02, 2018, 06:53:01 PM »
Great!

The reason I asked was because I remembered once I posted an FBSL BASIC script (ported from VB6) that drew varying ocean sunset views programmatically. You tried to replicate it but to no avail. IIRC the script was using VB6 Int() on floating-point arrays heavily, and now we know why SpecBAS failed. :)

ZXDunny

  • Guest
Re: Raytracing
« Reply #10 on: October 02, 2018, 09:20:00 PM »
Great!

The reason I asked was because I remembered once I posted an FBSL BASIC script (ported from VB6) that drew varying ocean sunset views programmatically. You tried to replicate it but to no avail. IIRC the script was using VB6 Int() on floating-point arrays heavily, and now we know why SpecBAS failed. :)

HOLY CARP! You're right!

I even added 32bit graphics to SpecBAS for that reason. Damn, I need to see if I can dig the code out.

Mike Lobanovsky

  • Guest
Re: Raytracing
« Reply #11 on: October 02, 2018, 09:54:07 PM »
I don't think I still have it in FBSL BASIC but I remember that I got the original VB6 script from the pouet.net demo scene. If you can't find it at your place, I can search for the original there.

Mike Lobanovsky

  • Guest
Re: Raytracing
« Reply #12 on: October 02, 2018, 11:45:50 PM »
OK Paul,

I failed to find the VB6 original because the demo dates back to 2010 and the very few links that used to lead to the VB6 source code rather than its 4K executable are now broken with error code 404.

But I've found the FBSL BASIC equivalent (rather slow, I must admit). :) It renders a non-persistent image directly to the on-screen window, and you shouldn't try to drag the window, or change focus, or change the window z-order before the calc completes, and drawing stops, and the overall time is displayed in the accompanying console. Otherwise the app will just freeze. The whole process takes approx. 30 secs to complete on my 3.6GHz PC. (it's a single-threaded demo, after all)

Whenever you see Floor() in the code below, it stands for VB6/Sinclair/Spectrum INT() in its truncate-towards-minus-infinity sense. Regretfully, the var names have been shortened/obfuscated, and comments, removed to yield an equivalent 4K-only tiny FBSL executable:

Code: [Select]
// VB6 code  (c)2010 Mikle           http://www.fbsl.net/phpbb2
// FBSL port (c)2014 Mike Lobanovsky http://www.fbsl.net/phpbb2

#AppType Console
#Include <Include/Windows.inc>

Dim %NZ[511, 511], WB[1024, 384 To 768], WX[1023, 384 To 767], WY[1023, 384 To 767]
Dim %Col[1023, 767], %CC[128, 8]
Dim FC, SX, SY, $PS * 64

SetWindowLong(ME, GWL_STYLE, &H6000000)
Resize(ME, 0, 0, 1024, 768)
Center(ME): Show(ME)

Begin Events
  Select Case CBMSG
    Case WM_NCHITTEST
      Return HTCAPTION
    Case WM_COMMAND
      If CBWPARAM = 2 Then PostMessage(ME, WM_CLOSE, 0, 0)
    Case WM_PAINT
      InvalidateRect(ME, NULL, FALSE)
      Render(BeginPaint(ME, PS)): EndPaint(ME, PS): Return 0
  End Select
End Events

Sub Render(hDC)
  Initialize()
  Sky()
  Colorize()
  Water()
  Air(hDC)
End Sub

Sub Initialize()
  Dim x, y, d = 64, d2 = 128, gtc = GetTickCount()
 
  Print "Running Initialize() ";
  Randomize
  Do
  For y = 0 To 511 Step d2
  For x = 0 To 511 Step d2
  NZ[(x + d) BAnd 511, y] = (NZ[x, y] + NZ[(x + d2) BAnd 511, y]) * 0.5 + d * (Rnd() - 0.5)
  NZ[x, (y + d) BAnd 511] = (NZ[x, y] + NZ[x, (y + d2) BAnd 511]) * 0.5 + d * (Rnd() - 0.5)
  NZ[(x + d) BAnd 511, (y + d) BAnd 511] = (NZ[x, y] + NZ[(x + d2) BAnd 511, (y + d2) BAnd 511] + NZ[x, (y + d2) BAnd 511] + NZ[(x + d2) BAnd 511, y]) * 0.25 + d * (Rnd() - 0.5)
  Next
  Next
  If d = 1 Then Exit Do
  d = d \ 2: d2 = d + d
  Loop
  Print GetTickCount() - gtc, " msec"
End Sub

Sub Air(hDC)
  Dim x, y, c, k1, k2, s, gtc = GetTickCount()

  Print "Running Air() ";
  For y = 0 To 767
    k1 = (1 - Abs(383.5 - y) / 384) ^ 5
    For x = 0 To 1023
      If y = SY Then
        k2 = 0.25
      Else
        k2 = ATn((x - SX) / (y - SY)) / M_TWOPI + 0.25
      End If
      If y - SY < 0 Then k2 = k2 + 0.5
      k2 = BN(k2 * 512, 0) * 0.03
      k2 = 0.2 - k2 ^ 2: If k2 < 0 Then k2 = 0
      s = 30 / SqR((x - SX) ^ 2 + (y - SY) ^ 2)
      If s > 1 Then s = 1
      c = Lerp(&HFFFFFF, FC, k2 * (1 - s))
      SetPixelV(hDC, x, y, Lerp(c, Col[x, y], k1))
    Next
  Next
  Print GetTickCount() - gtc, " msec"
End Sub

Sub Water()
  Dim x, y, x1, y1, k, kx, sx1, sy1, sx2, sy2, gtc = GetTickCount()
 
  Print "Running Water() ";
  For y = 767 DownTo 384
    k = (y - 383) * 0.5: kx = (900 - y) / 580
    For x = 1023 DownTo 0
      sy1 = 64000 / (y - 380)
      sx1 = (x - 511.5) * sy1 * 0.002
      sy2 = sy1 * 0.34 - sx1 * 0.71
      sx2 = sx1 * 0.34 + sy1 * 0.71
      sy1 = sy2 * 0.34 - sx2 * 0.21
      sx1 = sx2 * 0.34 + sy2 * 0.21
      WB[x, y] = BN(sx1, sy1) - BN(sx2, sy2)
      WX[x, y] = (WB[x + 1, y] - WB[x, y]) * k * kx
      WY[x, y] = (WB[x, y + 1] - WB[x, y]) * k
      x1 = Abs(x + WX[x, y])
      y1 = 768 - y + WY[x, y]
      If y1 < 0 Then
        y1 = 0
      ElseIf y1 > 383 Then
        y1 = 383
      End If
      Col[x, y] = Lerp(BC(x1 / 8 / 2, y1 / 48 / 2), &H251510, kx) ' water tint
    Next
  Next
  Print GetTickCount() - gtc, " msec"
End Sub

Sub Sky()
  Dim x, y, c1, c2, k, s, sx1, sy1, dy, gtc = GetTickCount()
 
  Print "Running Sky() ";
  SX = 100 + Rnd() * 824: SY = 192 + Rnd() * 157
  For y = 0 To 383
    sy1 = 100000 / (390 - y)
    For x = 0 To 1023
      sx1 = (x - 511.5) * sy1 * 0.0005
      k = BN(sx1, sy1) - BN(sx1 * 0.14 + sy1 * 0.21, sy1 * 0.14 - sx1 * 0.21)
      If k < -8 Then
        k = 0
      Else
        k = (k + 8) * 0.02 ' cloud density
      End If
      If k > 1 Then k = 1
      dy = y / 384
      FC = &H908000 + (SY + 500) * 0.2 ' haze tint
      c1 = Lerp(FC + 25, &H906050, dy)
      c2 = Lerp(&H807080, &HD0D0D0, dy)
      s = 30 / SqR((x - SX) ^ 2 + (y - SY) ^ 2) ' sun size
      If s > 1 Then s = 1
      c1 = Lerp(&HFFFFFF, c1, s)
      Col[x, y] = Lerp(c2, c1, k)
    Next
  Next
  Print GetTickCount() - gtc, " msec"
End Sub

Sub Colorize()
  Dim x, y, xx, yy, c, r, g, b, gtc = GetTickCount()
 
  Print "Running Colorize() ";
  For x = 0 To 127
  For y = 0 To 7
  Let(r, g, b) = 0
  For yy = 0 To 47
  For xx = 0 To 7
  c = Col[xx + x * 8, yy + y * 48]
  r = r + (c BAnd &HFF)
  g = g + (c BAnd &HFF00)
  b = b + ((c BAnd &HFF0000) >> 8)
  Next
  Next
  CC[x, y] = r \ 384 + ((g \ 384) BAnd &HFF00) + (((b \ 384) BAnd &HFF00) << 8)
  Next
  CC[x, 8] = CC[x, 7]
  Next
  Print GetTickCount() - gtc, " msec"
End Sub

Function BC(x, y)
  Dim ix = Floor(x), iy = Floor(y), SX = x - ix, SY = y - iy, c0, c1, c2, c3
  Dim ixy = (1 - SX) * (1 - SY), isxy = SX * (1 - SY), isyx = SY * (1 - SX), xy = SX * SY
 
  c0 = CC[ix BAnd 127, iy Mod 9]
  c1 = CC[(ix + 1) BAnd 127, iy Mod 9]
  c2 = CC[ix BAnd 127, (iy + 1) Mod 9]
  c3 = CC[(ix + 1) BAnd 127, (iy + 1) Mod 9]
 
  Return (c0 BAnd &HFF) * ixy + (c1 BAnd &HFF) * isxy + (c2 BAnd &HFF) * isyx + (c3 BAnd &HFF) * xy + _
  ((c0 BAnd &HFF00) * ixy + (c1 BAnd &HFF00) * isxy + (c2 BAnd &HFF00) * isyx + (c3 BAnd &HFF00) * xy BAnd &HFF00) + _
  ((c0 BAnd &HFF0000) * ixy + (c1 BAnd &HFF0000) * isxy + (c2 BAnd &HFF0000) * isyx + (c3 BAnd &HFF0000) * xy BAnd &HFF0000)
End Function

Function BN(x, y)
  Dim ix = Floor(x), iy = Floor(y), SX = x - ix, SY = y - iy, isx = 1 - SX, isy = 1 - SY, dx = (ix + 1) BAnd 511, dy = (iy + 1) BAnd 511
 
  ix = ix BAnd 511: iy = iy BAnd 511
  Return NZ[ix, iy] * isx * isy + NZ[dx, iy] * SX * isy + NZ[ix, dy] * isx * SY + NZ[dx, dy] * SX * SY
End Function

Function Lerp(c1, c2, k)
Return((c1 BAnd &HFF) * k + (c2 BAnd &HFF) * (1 - k)) BOr (((c1 BAnd &HFF00) * k + (c2 BAnd &HFF00) * (1 - k)) BAnd &HFF00) BOr (((c1 BAnd &HFF0000) * k + (c2 BAnd &HFF0000) * (1 - k)) BAnd &HFF0000)
End Function
« Last Edit: October 03, 2018, 12:00:54 AM by Mike Lobanovsky »