PureBasic Forums - English

Code: Select all

; Experiments In Ray Tracing 2 
; some pretty spheres in realtime (around 70 - 80 FPS on 4GHz Skylake) 
; PB 5.42, done on Windows7 64bit
; Hades, 2016

CompilerIf #PB_Compiler_Thread = 0
  MessageRequester( "Sorry...", "Please, set the compiler to threadsafe before compiling" )
  End
CompilerEndIf

EnableExplicit
DisableDebugger

#TileSize = 16

#RenderWidth = 640  ; RenderWidth/Height needs to be a multiple of #TileSize, otherwise you'll crash!
#RenderHeight = 384 ; Window size and render size are independent, but choosing one as the multiple of the other gives better results.
#WindowWidth = 640
#WindowHeight = 384

#MaxDist = 1000000
#FpError = 0.001

Structure stVector
  x.f
  y.f
  z.f
EndStructure

Structure stMatrix
  m.f[16]
EndStructure  

Structure stColor
  r.f
  g.f
  b.f
EndStructure

Structure stMaterial
  DiffuseColor.stColor
  DiffuseIntensity.f
  SpecularIntensity.f
  Glossiness.f
  AmbientIntensity.f
  Reflectance.f
  Transmittance.f
  IndexOfRefraction.f
  Shader.i
EndStructure
Global NewList Material.stMaterial()

Structure stLight   ; directional light, far away, like the sun
  Direction.stVector
EndStructure
Global Light.stLight

Structure stRay
  Origin.stVector
  Direction.stVector
  length.f
  PrimitiveHit.i
  Inside.i
  Weight.f
EndStructure

Structure stSphere
  Position.stVector
  radius.f
  *Material.stMaterial
EndStructure
Global NewList Sphere.stSphere()
Global Dim SphereX.f(3)
Global Dim SphereY.f(3)
Global Dim SphereZ.f(3)
Global Dim SphereRadius.f(3)
Global Dim *SphereMaterial.stMaterial(3)

Structure stRenderBuffer
  pAddr.i[2]  ; I am using 2 render buffers, so while the last threads are still working on the...
  Width.i     ; ...current frame, the one's that are done can already start with the next one,...
  Height.i
EndStructure
Global RenderBuffer.stRenderBuffer

Structure stCamera
  Pos.stVector
  Matrix.stMatrix
  VScrTL.stVector
  VScrR.stVector
  VScrB.stVector
EndStructure
Global Dim Camera.stCamera(1) ; ...but that means i need two versions for stuff like this, too.

Structure stJob
  x.i
  y.i
  Frame.i
EndStructure
Global NumThreads
Global Dim Job.stJob(0)
Global numJobs
Global JobIdx
Global mutJob
Global mutFrameSync
mutJob = CreateMutex()
mutFrameSync = CreateMutex()

Structure stAABB ; Axis aligned bounding box
  min.f[3]
  max.f[3]
EndStructure

Structure stAABoxHit
  Near.f
  Far.f
EndStructure  

Structure stBVH
  min.f[3]
  max.f[3]
  Content.i
EndStructure

Structure stBVHPrimitive
  AABB.stAABB
  Centre.f[3]
EndStructure  
Global Dim BVHPrimitive.stBVHPrimitive(0)
Global Dim BVHPrimIdx(0)
Global FreeNode.i


Structure stMain
  Running.i
  Frame.i
  DisplList.i
  Material_Default.i
  Shader_Default.i
  Material_Background.i
EndStructure
Global Main.stMain


Macro ColorSet(Color, cr, cg, cb)
  Color\r = cr
  Color\g = cg
  Color\b = cb
EndMacro

Macro VectorSet(Vector, VX, VY, VZ)
  Vector\x = VX
  Vector\y = VY
  Vector\z = VZ
EndMacro

Macro Normalize(Vector, rLength)
  rLength =  1.0 / Sqr(Vector\x * Vector\x + Vector\y * Vector\y + Vector\z * Vector\z)
  Vector\x * rLength
  Vector\y * rLength
  Vector\z * rLength
EndMacro  

Declare Shade(*Ray.stRay, *Color.stColor, Depth.i, ShadowRayFlag.i = #False)

Procedure.i RTMaterialCreate(r.f = 1.0, g.f = 1.0, b.f = 1.0, Shader.i = 0)
  If Shader = 0
    Shader = Main\Shader_Default
  EndIf
  AddElement(Material())
  ColorSet(Material()\DiffuseColor, r, g, b)
  Material()\SpecularIntensity = 0.3
  Material()\Glossiness = 5
  Material()\AmbientIntensity = 0.1
  Material()\Reflectance = 0.0
  Material()\Transmittance = 0.0
  Material()\IndexOfRefraction = 1.0
  Material()\DiffuseIntensity = 1.0 - 1.0 * (Material()\Reflectance + Material()\Transmittance)
  Material()\Shader = Shader
  ProcedureReturn @Material()
EndProcedure

Procedure SphereCreate(x.f, y.f, z.f, radius.f, *Material.stMaterial)
  AddElement(Sphere())
  VectorSet(Sphere()\Position, x, y, z)
  Sphere()\radius = radius
  Sphere()\Material = *Material
EndProcedure

Procedure DisplayBufferInit(BufferWidth.i, BufferHeight.i, VisibleWidth.i, VisibleHeight.i)
  Static Tex.i
  glShadeModel_(#GL_FLAT)
  glEnable_(#GL_TEXTURE_2D)
  glBindTexture_(#GL_TEXTURE_2D, @Tex)
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR)
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR)
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_CLAMP)
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_CLAMP)
  Main\DisplList = glGenLists_(1)
  glNewList_(Main\DisplList, #GL_COMPILE)
  glBegin_(#GL_QUADS)
  glTexCoord2f_(0.0, 0.0)
  glVertex2i_(0,0)
  glTexCoord2f_(0.0, 1.0)
  glVertex2i_(0, VisibleHeight)
  glTexCoord2f_(1.0, 1.0)
  glVertex2i_(VisibleWidth, VisibleHeight)
  glTexCoord2f_(1.0, 0.0)
  glVertex2i_(VisibleWidth, 0)
  glEnd_()
  glEndList_()
EndProcedure

Procedure glWindowOpen(Width.i, Height.i, Name.s = "Experiments In Ray Tracing 2", flags.i = #PB_Window_SystemMenu | #PB_Window_MinimizeGadget | #PB_Window_ScreenCentered)
  OpenWindow(0, 0, 0, Width, Height, Name, flags)
  OpenGLGadget(0, 0, 0, Width, Height, #PB_OpenGL_NoFlipSynchronization | #PB_OpenGL_NoDepthBuffer)
  SetGadgetAttribute(0, #PB_OpenGL_FlipBuffers, #True)
  glOrtho_(0.0, Width, Height, 0.0, 0.0, 1.0)
EndProcedure

Procedure DisplayBufferShow(pBuffer.i)
  glEnable_(#GL_TEXTURE_2D)
  glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, #RenderWidth, #RenderHeight, 0, #GL_BGRA_EXT, #GL_UNSIGNED_BYTE, pBuffer)
  glCallList_(Main\DisplList)
  SetGadgetAttribute(0, #PB_OpenGL_FlipBuffers, #True)
EndProcedure

Procedure Intersect_Spheres(*Ray.stRay)
  Protected MaxD.f, t1.f, t2.f, idx1, idx2, *SphereX.float, *SphereY.float, *SphereZ.float, *SphereRadius.float
  *SphereX.float = @SphereX(0)
  *SphereY.float = @SphereY(0)
  *SphereZ.float = @SphereZ(0)
  *SphereRadius.float = @SphereRadius(0)
  MaxD = #MaxDist
  !mov edx, [p.p_Ray]
  !movss xmm0, [edx]        ; *Ray\Origin\x
  !shufps xmm0, xmm0, 0
  !movss xmm1, [edx + 4]    ; *Ray\Origin\y
  !shufps xmm1, xmm1, 0
  !movss xmm2, [edx + 8]    ; *Ray\Origin\z
  !shufps xmm2, xmm2, 0

  !mov eax, [p.p_SphereX]
  !subps xmm0, [eax]              ; RelPos\x = *Ray\Origin\x - *Sphere\Position\x
  !mov ecx, [p.p_SphereY]
  !subps xmm1, [ecx]              ; RelPos\y = *Ray\Origin\y - *Sphere\Position\y
  !mov eax, [p.p_SphereZ]
  !subps xmm2, [eax]              ; RelPos\z = *Ray\Origin\z - *Sphere\Position\z
  !movss xmm3, [edx + 12]   ; *Ray\Direction\x
  !shufps xmm3, xmm3, 0
  !movss xmm4, [edx + 16]   ; *Ray\Direction\y
  !shufps xmm4, xmm4, 0
  !movss xmm5, [edx + 20]   ; *Ray\Direction\z
  !shufps xmm5, xmm5, 0  
  ;b = *Ray\Direction\x * RelPos\x + *Ray\Direction\y * RelPos\y + *Ray\Direction\z * RelPos\z
  !mulps xmm3, xmm0               ; *Ray\Direction\x * RelPos\x  
  !mulps xmm4, xmm1               ; *Ray\Direction\y * RelPos\y  
  !mulps xmm5, xmm2               ; *Ray\Direction\z * RelPos\z
  
  !addps xmm3, xmm4 
  !addps xmm3, xmm5
  ;c = (RelPos\x * RelPos\x + RelPos\y * RelPos\y + RelPos\z * RelPos\z) - *Sphere\radius * *Sphere\radius
  !mov edx, [p.p_SphereRadius]
  !movups xmm7, [edx]
  !mulps xmm7, xmm7               ; *Sphere\radius * *Sphere\radius
  
  !mulps xmm0, xmm0               ; RelPos\x * RelPos\x  
  !mulps xmm1, xmm1               ; RelPos\y * RelPos\y  
  !mulps xmm2, xmm2               ; RelPos\z * RelPos\z
  
  !addps xmm0, xmm1
  !addps xmm0, xmm2
  !subps xmm0, xmm7
  
  ;d = b * b - c
  !movaps xmm4, xmm3
  !mulps xmm4, xmm4
  !subps xmm4, xmm0
  ; setting everything below 0.0 to 0.0, to not cause trouble with sqrt
  !xorps xmm6, xmm6               ; all 0.0
  !maxps xmm4, xmm6
  
  !sqrtps xmm4, xmm4              ; d = sqr(d)
  
  ; b = -b
  !pcmpeqw xmm7, xmm7             ; all bits to 1
  !pslld xmm7, 31                 ; only sign bits set
  !xorps xmm3, xmm7               ; flip sign bits
  
  ; There are two intersections with a sphere possible: t1 = b - d and t2 = b + d
  !movaps xmm1, xmm3
  !subps xmm3, xmm4               ; t1 = b - d
  !addps xmm1, xmm4               ; t2 = b + d
  ; set t1 to #MaxDist if d <= 0.0 or t1 <= 0.0, so no hit will be detected, to avoid conditonals
  !movss xmm7, [p.v_MaxD]         ; load #MaxDist
  !shufps xmm7, xmm7, 0
  !cmpps xmm4, xmm6, 2            ; create mask d <= 0
  !movaps xmm0, xmm3
  !cmpps xmm0, xmm6, 2            ; create mask t1 <= 0
  !orps xmm0, xmm4                ; combine
  
  ;!blendvps xmm3, xmm7, xmm0      ; t = #MaxDist for d <= 0 or t1 <= 0
  !movaps xmm2, xmm0
  !movaps xmm5, xmm3
  !andps xmm7, xmm2
  !andnps xmm2, xmm5
  !orps xmm7, xmm2
  !movaps xmm3, xmm7

  ; find closest hit
  !movaps xmm5, xmm3
  !movaps xmm0, xmm3
  !shufps xmm0, xmm5, 01001110b   ; shuffle to 1032
  !minps xmm5, xmm0
  !movaps xmm2, xmm5
  !shufps xmm5, xmm5, 10110001b   ; shuffle to 2301
  !minps xmm5, xmm2               ; now smallest t1 is in all positions
  !movss [p.v_t1], xmm5           ; store closest hit
  !cmpeqps xmm5, xmm3             ; set mask for smallest t1 in original
  !movmskps eax, xmm5             ; move sign to rax
  !bsf eax, eax                   ; turn sign into index
  !mov [p.v_idx1], eax            ; store index of closest hit
  
  ; same for t2 now:
  !movaps xmm0, xmm1
  !cmpps xmm0, xmm6, 2
  !orps xmm0, xmm4
  ;!blendvps xmm1, xmm7, xmm0
  !movaps xmm2, xmm0
  !movaps xmm5, xmm1
  !andps xmm7, xmm2
  !andnps xmm2, xmm5
  !orps xmm7, xmm2
  !movaps xmm1, xmm7
  
  !movaps xmm5, xmm1
  !movaps xmm0, xmm1
  !shufps xmm0, xmm5, 01001110b
  !minps xmm5, xmm0
  !movaps xmm2, xmm5
  !shufps xmm5, xmm5, 10110001b
  !minps xmm5, xmm2
  !movss [p.v_t2], xmm5
  !cmpeqps xmm5, xmm1
  !movmskps eax, xmm5
  !bsf eax, eax
  !mov [p.v_idx2], eax
  If t1 < *Ray\length
    *Ray\length = t1
    *Ray\PrimitiveHit = idx1
  EndIf  
  If t2 < *Ray\length
    *Ray\length = t2
    *Ray\PrimitiveHit = idx2
    *Ray\inside = 1
  EndIf  
EndProcedure

Procedure Intersect_Spheres1(*Ray.stRay) ; only here for reference
  Protected RelPos.stVector, i, b.f, c.f, d.f, t1.f, t2.f
  FirstElement(Sphere())
  For i = 0 To ListSize(Sphere()) - 1
    VectorSet(RelPos, *Ray\Origin\x - SphereX(i), *Ray\Origin\y - SphereY(i), *Ray\Origin\z - SphereZ(i))
    b = *Ray\Direction\x * RelPos\x + *Ray\Direction\y * RelPos\y + *Ray\Direction\z * RelPos\z
    c = (RelPos\x * RelPos\x + RelPos\y * RelPos\y + RelPos\z * RelPos\z) - SphereRadius(i) * SphereRadius(i)
    d = b * b - c
    If d > 0.0
      d = Sqr(d)
      t1 = (-b - d)
      If t1 > 0.0 
        If t1 < *Ray\length 
          *Ray\length = t1                  
          *Ray\PrimitiveHit = i
        EndIf
      Else  
        t2 = (-b + d)
        If t2 > 0.0 And t2 < *Ray\length ; a hit from inside the sphere
          *Ray\length = t2
          *Ray\PrimitiveHit = i
          *Ray\inside = 1
        EndIf  
      EndIf
    EndIf
  Next
EndProcedure

Procedure.i Shadow(*HitPosition.stVector, *Normal.stVector, *LightDirection.stVector, LightDistance.f)
  Protected ShadowRay.stRay, Col.stColor
  VectorSet(ShadowRay\Origin, *HitPosition\x, *HitPosition\y, *HitPosition\z)
  VectorSet(ShadowRay\Direction, *LightDirection\x, *LightDirection\y, *LightDirection\z)
  ShadowRay\length = LightDistance
  ShadowRay\PrimitiveHit = -1
  
  Intersect_Spheres(ShadowRay)
  Shade(ShadowRay, Col, 1, #True)
  If ShadowRay\PrimitiveHit >= 0
    ProcedureReturn #True
  EndIf
  ProcedureReturn #False
EndProcedure

Procedure GetSphereNormal(SphIdx.i, *Ray.stRay, *HitPosition.stVector, *Normal.stVector)
  Protected rLength.f
  VectorSet(*Normal, *HitPosition\x - SphereX(SphIdx), *HitPosition\y - SphereY(SphIdx), *HitPosition\z - SphereZ(SphIdx))
  Normalize(*Normal, rLength.f)
  If *Normal\x * *Ray\Direction\x + *Normal\y * *Ray\Direction\y + *Normal\z * *Ray\Direction\z > 0.0
    VectorSet(*Normal, -*Normal\x, -*Normal\y, -*Normal\z)
  EndIf
EndProcedure

Procedure GetLightSpecularColor(*Ray.stRay, *Light.stLight, *Normal.stVector, *LightDirection.stVector, LightNormalAngle.f, *LightSpecularColor.stColor)
  Protected Angle.f, ReflectedLightVector.stVector, Intensity.f, *Material.stMaterial
  *Material = *SphereMaterial(*Ray\PrimitiveHit)
  Angle = 2.0 * LightNormalAngle
  VectorSet(ReflectedLightVector, *LightDirection\x - Angle * *Normal\x, *LightDirection\y - Angle * *Normal\y, *LightDirection\z - Angle * *Normal\z)
  Angle = *Ray\Direction\x * ReflectedLightVector\x + *Ray\Direction\y * ReflectedLightVector\y + *Ray\Direction\z * ReflectedLightVector\z
  If Angle > 0.0
    Intensity = *Material\SpecularIntensity * Pow(Angle, *Material\Glossiness)
    ColorSet(*LightSpecularColor, Intensity, Intensity, Intensity)
  EndIf
EndProcedure

Procedure GetLight(*Ray.stRay, *HitPosition.stVector, *Normal.stVector, *LightDiffuseColor.stColor, *LightSpecularColor.stColor)
  Protected LightNormalAngle.f, LightDistance.f, ShadowRayOrigin.stVector, *Material.stMaterial
  *Material = *SphereMaterial(*Ray\PrimitiveHit)
  ColorSet(*LightDiffuseColor, 0.0, 0.0, 0.0)
  ColorSet(*LightSpecularColor, 0.0, 0.0, 0.0)
  LightDistance = #MaxDist
  LightNormalAngle = *Normal\x * Light\Direction\x + *Normal\y * Light\Direction\y + *Normal\z * Light\Direction\z
  ColorSet(*LightDiffuseColor, *Material\AmbientIntensity, *Material\AmbientIntensity, *Material\AmbientIntensity)
  If LightNormalAngle > *Material\AmbientIntensity
    VectorSet(ShadowRayOrigin, *HitPosition\x + #FpError * *Normal\x, *HitPosition\y + #FpError * *Normal\y, *HitPosition\z + #FpError * *Normal\z)
    If Shadow(ShadowRayOrigin, *Normal, Light\Direction, LightDistance) = #False
      ColorSet(*LightDiffuseColor, LightNormalAngle, LightNormalAngle, LightNormalAngle)
      If *SphereMaterial(*Ray\PrimitiveHit)\SpecularIntensity > 0.0
        GetLightSpecularColor(*Ray, Light, *Normal, Light\Direction, LightNormalAngle, *LightSpecularColor)
      EndIf  
    EndIf
  EndIf
EndProcedure

Procedure Shader_Holes(*Ray.stRay, *Material.stMaterial, *Color.stColor, Depth.i, ShadowRayFlag.i)
  Protected HitPosition.stVector, Normal.stVector, LightDiffuseColor.stColor, LightSpecularColor.stColor, x, y, z, Rel.stVector, d.f, TransmittedRay.stRay
  VectorSet(HitPosition, *Ray\Origin\x + *Ray\length * *Ray\Direction\x, *Ray\Origin\y + *Ray\length * *Ray\Direction\y, *Ray\Origin\z + *Ray\length * *Ray\Direction\z)
  GetSphereNormal(*Ray\PrimitiveHit, *Ray, HitPosition, Normal)
  GetLight(*Ray, HitPosition, Normal, LightDiffuseColor, LightSpecularColor)
  VectorSet(Rel, HitPosition\x * 0.1, HitPosition\y * 0.1, HitPosition\z * 0.1)
  x = Rel\x : y = Rel\y : z = Rel\z
  VectorSet(Rel, Rel\x - x, Rel\y - y, Rel\z - z)
  d = 10 * Sqr(Rel\x * Rel\x + Rel\y * Rel\y + Rel\z * Rel\z)
  If d > 4.8
    If ShadowRayFlag
      ProcedureReturn
    EndIf
    If *Ray\inside
      *Color\r + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\g * LightDiffuseColor\r + LightSpecularColor\r)
      *Color\g + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\r * LightDiffuseColor\g + LightSpecularColor\g)
      *Color\b + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\b * LightDiffuseColor\b + LightSpecularColor\b)
    Else  
      *Color\r + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\r * LightDiffuseColor\r + LightSpecularColor\r)
      *Color\g + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\g * LightDiffuseColor\g + LightSpecularColor\g)
      *Color\b + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\b * LightDiffuseColor\b + LightSpecularColor\b)
    EndIf
  Else
    If ShadowRayFlag
      VectorSet(HitPosition, HitPosition\x - #FpError * Normal\x, HitPosition\y - #FpError * Normal\y, HitPosition\z - #FpError * Normal\z)
      VectorSet(Normal, -Normal\x, -Normal\y, -Normal\z)
      *Ray\PrimitiveHit = -1
      If Shadow(HitPosition, Normal, Light\Direction, #MaxDist)
        *Ray\PrimitiveHit = 0
      EndIf  
      ProcedureReturn
    EndIf  
    VectorSet(TransmittedRay\Origin, HitPosition\x - #FpError * Normal\x, HitPosition\y - #FpError * Normal\y, HitPosition\z - #FpError * Normal\z)
    VectorSet(TransmittedRay\Direction, *Ray\Direction\x, *Ray\Direction\y, *Ray\Direction\z)
    TransmittedRay\length = #MaxDist
    TransmittedRay\PrimitiveHit = -1
    TransmittedRay\Weight = *Ray\Weight
    Intersect_Spheres(TransmittedRay)
    Shade(TransmittedRay, *Color, Depth + 1)
  EndIf  
EndProcedure  

Procedure Shader_Tiles(*Ray.stRay, *Material.stMaterial, *Color.stColor, Depth.i, ShadowRayFlag.i)
  Protected HitPosition.stVector, Normal.stVector, LightDiffuseColor.stColor, LightSpecularColor.stColor
  Protected Tile.i, TileX.i, TileY.i, Brightness.f
  If ShadowRayFlag
    ProcedureReturn
  EndIf  
  VectorSet(HitPosition, *Ray\Origin\x + *Ray\length * *Ray\Direction\x, *Ray\Origin\y + *Ray\length * *Ray\Direction\y, *Ray\Origin\z + *Ray\length * *Ray\Direction\z)
  GetSphereNormal(*Ray\PrimitiveHit, *Ray, HitPosition, Normal)
  GetLight(*Ray, HitPosition, Normal, LightDiffuseColor, LightSpecularColor)
  TileX = HitPosition\x
  TileY = HitPosition\z
  Brightness = 0.6
  TileX & 30
  TileY & 30
  If TileX And TileY
    Brightness + 0.4
    *Color\r + *Ray\Weight * (*Material\DiffuseIntensity * Brightness * LightDiffuseColor\r + LightSpecularColor\r)
    *Color\g + *Ray\Weight * (*Material\DiffuseIntensity * Brightness * LightDiffuseColor\g + LightSpecularColor\g)
    *Color\b + *Ray\Weight * (*Material\DiffuseIntensity * Brightness * LightDiffuseColor\b + LightSpecularColor\b)
  Else
    *Color\r + *Ray\Weight * (*Material\DiffuseIntensity * Brightness * LightDiffuseColor\r)
    *Color\g + *Ray\Weight * (*Material\DiffuseIntensity * Brightness * LightDiffuseColor\g)
    *Color\b + *Ray\Weight * (*Material\DiffuseIntensity * Brightness * LightDiffuseColor\b)
  EndIf
EndProcedure

Procedure Shader_Sky(*Ray.stRay, *Material.stMaterial, *Color.stColor, Depth.i, ShadowRayFlag.i)
  Protected Angle.f, Brightness.f
  If ShadowRayFlag
    ProcedureReturn
  EndIf  
  Angle = 0.15 + *Ray\Direction\y * 0.85
  If Angle < 0.0
    Angle = 0.0
  EndIf
  Angle = 1.0 - Angle
  Angle = Angle * Angle
  Brightness = 0.2 + 0.7 * Angle * Angle
  *Color\r + *Ray\Weight * Brightness
  *Color\g + *Ray\Weight * Brightness
  *Color\b + *Ray\Weight * 0.9
EndProcedure

Procedure Shader_Default(*Ray.stRay, *Material.stMaterial, *Color.stColor, Depth.i, ShadowRayFlag.i)
  Protected HitPosition.stVector, Normal.stVector, LightDiffuseColor.stColor, LightSpecularColor.stColor
  Protected ReflectedRay.stRay, Angle.f, TransmittedRay.stRay, OutsideIoR.f, RelativeIoR.f, SinT2.f, v.f, Weight.f, RayNormalAngle.f, Fresnel.f
  If ShadowRayFlag
    ProcedureReturn
  EndIf  
  VectorSet(HitPosition, *Ray\Origin\x + *Ray\length * *Ray\Direction\x, *Ray\Origin\y + *Ray\length * *Ray\Direction\y, *Ray\Origin\z + *Ray\length * *Ray\Direction\z)
  GetSphereNormal(*Ray\PrimitiveHit, *Ray, HitPosition, Normal)
  GetLight(*Ray, HitPosition, Normal, LightDiffuseColor, LightSpecularColor)
  If Depth < 10
    If *Material\Reflectance > 0.0
      RayNormalAngle = (*Ray\Direction\x * Normal\x + *Ray\Direction\y * Normal\y + *Ray\Direction\z * Normal\z)
      Fresnel = Pow((1.0 + RayNormalAngle), 5.0)
      ReflectedRay\Weight = *Ray\Weight * (*Material\Reflectance + (1 - *Material\Reflectance) * Fresnel)
      If ReflectedRay\Weight > 0.003
        VectorSet(ReflectedRay\Origin, HitPosition\x + #FpError * Normal\x, HitPosition\y + #FpError * Normal\y, HitPosition\z + #FpError * Normal\z)
        Angle = 2.0 * (*Ray\Direction\x * Normal\x + *Ray\Direction\y * Normal\y + *Ray\Direction\z * Normal\z)
        VectorSet(ReflectedRay\Direction, *Ray\Direction\x - Angle * Normal\x, *Ray\Direction\y - Angle * Normal\y, *Ray\Direction\z - Angle * Normal\z)
        ReflectedRay\length = #MaxDist
        ReflectedRay\PrimitiveHit = -1
        ReflectedRay\Inside = *Ray\Inside
        Intersect_Spheres(ReflectedRay)
        Shade(ReflectedRay, *Color, Depth + 1)
      EndIf
    EndIf
    If *Material\Transmittance > 0.0
      TransmittedRay\Weight = *Ray\Weight * *Material\Transmittance * (1 - Fresnel)
      If TransmittedRay\Weight > 0.003
        ; we assume that outside of the refractive object is air or vacuum (index of refraction = 1.0)
        OutsideIoR = 1.0
        If *Ray\Inside
          RelativeIoR = *Material\IndexOfRefraction / OutsideIoR
          TransmittedRay\Inside = #False  
        Else
          RelativeIoR = OutsideIoR / *Material\IndexOfRefraction
          TransmittedRay\Inside = #True
        EndIf
        ; we have to make sure the origin of the ray is below the surface
        VectorSet(TransmittedRay\Origin, HitPosition\x - #FpError * Normal\x, HitPosition\y - #FpError * Normal\y, HitPosition\z - #FpError * Normal\z)
        Angle = *Ray\Direction\x * Normal\x + *Ray\Direction\y * Normal\y + *Ray\Direction\z * Normal\z
        SinT2 = RelativeIoR * RelativeIoR * (1.0 - Angle * Angle)
        If SinT2 <= 1.0 ; if there is no total internal reflection:
          v = RelativeIoR * Angle + Sqr(1.0 - SinT2)
          VectorSet(TransmittedRay\Direction, RelativeIoR * *Ray\Direction\x - v * Normal\x, RelativeIoR * *Ray\Direction\y - v * Normal\y, RelativeIoR * *Ray\Direction\z - v * Normal\z)
          TransmittedRay\length = #MaxDist
          TransmittedRay\PrimitiveHit = -1
          Intersect_Spheres(TransmittedRay)
          Shade(TransmittedRay, *Color, Depth + 1)
        EndIf
      EndIf
    EndIf
  EndIf
  
  Weight = *Ray\Weight * (1 - (*Material\Reflectance + *Material\Transmittance)) * (1 - Fresnel)
  *Color\r + Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\r * LightDiffuseColor\r + LightSpecularColor\r)
  *Color\g + Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\g * LightDiffuseColor\g + LightSpecularColor\g)
  *Color\b + Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\b * LightDiffuseColor\b + LightSpecularColor\b)
EndProcedure  

Procedure Shade(*Ray.stRay, *Color.stColor, Depth.i, ShadowRayFlag.i = #False)
  Protected *Material.stMaterial
  If *Ray\PrimitiveHit >= 0
    *Material = *SphereMaterial(*Ray\PrimitiveHit)
  Else
    *Material = Main\Material_Background
  EndIf
  CallFunctionFast(*Material\Shader, *Ray, *Material, *Color, Depth, ShadowRayFlag.i)
EndProcedure

Procedure RenderTile(TileX.i, TileY.i, Frame.i)
  Protected *Buffer.Long, *Material.stMaterial, x.i, y.i, i.i, rLength.f, Ray.stRay, Col.stColor, OffX.stVector, OffY.stVector, RayLS.stVector
  Protected HitPosition.stVector, Normal.stVector, LightNormalAngle.f, Angle.f, ReflectedLightVector.stVector, SpecIntensity.f, Brightness.f
  OffX\x = (Camera(Frame)\VScrR\x - Camera(Frame)\VScrTL\x) / #RenderWidth
  OffX\y = (Camera(Frame)\VScrR\y - Camera(Frame)\VScrTL\y) / #RenderWidth
  OffX\z = (Camera(Frame)\VScrR\z - Camera(Frame)\VScrTL\z) / #RenderWidth
  OffY\x = (Camera(Frame)\VScrB\x - Camera(Frame)\VScrTL\x) / #RenderHeight
  OffY\y = (Camera(Frame)\VScrB\y - Camera(Frame)\VScrTL\y) / #RenderHeight
  OffY\z = (Camera(Frame)\VScrB\z - Camera(Frame)\VScrTL\z) / #RenderHeight
  VectorSet(Ray\Origin, Camera(Frame)\Pos\x, Camera(Frame)\Pos\y, Camera(Frame)\Pos\z)
  For y = TileY To TileY + #TileSize - 1 
    *Buffer = RenderBuffer\pAddr[Frame] + TileX * 4 + y * #RenderWidth * 4
    VectorSet(RayLS, 0.5 + Camera(Frame)\VScrTL\x + y * OffY\x, 0.5 + Camera(Frame)\VScrTL\y + y * OffY\y, Camera(Frame)\VScrTL\z + y * OffY\z)
    For x = TileX To TileX + #TileSize - 1
      VectorSet(Ray\Direction, RayLS\x + x * OffX\x, RayLS\y + x * OffX\y, RayLS\z + x * OffX\z)
      Normalize(Ray\Direction, rLength)
      Ray\PrimitiveHit = -1
      Ray\length = #MaxDist
      Ray\Weight = 1.0
      Ray\Inside = 0
      ColorSet(Col, 0, 0, 0)
      Intersect_Spheres(Ray)
      Shade(Ray, Col, 1)
      *Buffer\l = RGB(Col\b * 255, Col\g * 255, Col\r * 255)
      *Buffer + 4
    Next
  Next
EndProcedure  

Procedure WorkerThread(ThreadIdx)
  Protected iJob, pTravStack.i
  While Main\Running
    iJob = -1
    LockMutex(mutJob)
    If JobIdx < numJobs
      If JobIdx + 1 = numJobs
        LockMutex(mutFrameSync)
      EndIf  
      iJob = JobIdx
      JobIdx + 1    
    EndIf  
    UnlockMutex(mutJob)
    If iJob > -1
      RenderTile(Job(iJob)\x, Job(iJob)\y, Job(iJob)\Frame)
      Job(iJob)\Frame = 1 - Job(iJob)\Frame
      If iJob + 1 = numJobs
        UnlockMutex(mutFrameSync)
      EndIf  
    Else  
      Delay(0)
    EndIf  
  Wend
EndProcedure

Procedure SetUpJobs(RenderWidth.i, RenderHeight.i)
  Protected x, y, i
  ReDim Job((RenderHeight * RenderWidth) / (#TileSize * #TileSize))
  For y = 0 To RenderHeight - 1 Step #TileSize
    For x = 0 To RenderWidth - 1 Step #TileSize
      Job(i)\x = x
      Job(i)\y = y
      i + 1
    Next x  
  Next y  
  numJobs = i
  JobIdx = i
EndProcedure

Procedure CreateThreads()
  Protected i
  NumThreads = CountCPUs(#PB_System_ProcessCPUs)
  For i = 0 To NumThreads - 1
    CreateThread(@WorkerThread(), i)
  Next
EndProcedure

Procedure SceneCreate()
  Protected *Material.stMaterial, i, x, y, z, rLength.f
  *Material = RTMaterialCreate(0.0, 0.0, 0.0, @Shader_Tiles())
  *Material\SpecularIntensity = 0.3
  *Material\Glossiness = 20
  SphereCreate(0, -5025, 0, 5000, *Material)
  
  *Material = RTMaterialCreate(0.7, 0.2, 0.1, @Shader_Holes())
  *Material\SpecularIntensity = 0.2
  *Material\Glossiness = 5
  SphereCreate(-40, -6, 40, 20, *Material)
  
  *Material = RTMaterialCreate(0.0, 0.2, 0.4)
  *Material\Reflectance = 0.01
  *Material\Transmittance = 0.95
  *Material\IndexOfRefraction = 1.5
  *Material\SpecularIntensity = 0.0
  SphereCreate(14, 0, 0, 25, *Material)
  
  *Material = RTMaterialCreate(0.0, 0.2, 0.4)
  *Material\Reflectance = 0.80
  *Material\SpecularIntensity = 0.0
  SphereCreate(5, 5, 50, 15, *Material)
  
  VectorSet(Light\Direction, 1, 1, -1)
  Normalize(Light\Direction, rLength)
  
  FirstElement(Sphere())
  For i = 0 To ListSize(Sphere()) - 1       ; we need to turn the primitive data
    SphereX(i) = Sphere()\Position\x        ; into a SIMD (SSE) friendly format
    SphereY(i) = Sphere()\Position\y        ; Lists and multithreading don't play nice, anyway
    SphereZ(i) = Sphere()\Position\z
    SphereRadius(i) = Sphere()\radius
    *SphereMaterial(i) = Sphere()\Material
    NextElement(Sphere())
  Next i
EndProcedure

Procedure MatrixIdentity(*Matrix.stMatrix)
  *Matrix\m[0] = 1.0 : *Matrix\m[1] = 0.0 : *Matrix\m[2] = 0.0 : *Matrix\m[3] = 0.0
  *Matrix\m[4] = 0.0 : *Matrix\m[5] = 1.0 : *Matrix\m[6] = 0.0 : *Matrix\m[7] = 0.0
  *Matrix\m[8] = 0.0 : *Matrix\m[9] = 0.0 : *Matrix\m[10] = 1.0 : *Matrix\m[11] = 0.0
  *Matrix\m[12] = 0.0 : *Matrix\m[13] = 0.0 : *Matrix\m[14] = 0.0 : *Matrix\m[15] = 1.0
EndProcedure

Procedure MatrixRotateX(*Matrix.stMatrix, Angle.f)
  Protected Help.f, tsin.f, tcos.f
  tsin = Sin(Angle) 
  tcos = Cos(Angle) 
  
  Help = *Matrix\m[4] * tcos + *Matrix\m[8] * tsin
  *Matrix\m[8] = *Matrix\m[4] * -tsin + *Matrix\m[8] * tcos
  *Matrix\m[4] = Help 
  
  Help = *Matrix\m[5] * tcos + *Matrix\m[9] * tsin
  *Matrix\m[9] = *Matrix\m[5] * -tsin + *Matrix\m[9] * tcos
  *Matrix\m[5] = Help 
  
  Help = *Matrix\m[6] * tcos + *Matrix\m[10] * tsin
  *Matrix\m[10] = *Matrix\m[6] * -tsin + *Matrix\m[10] * tcos
  *Matrix\m[6] = Help 
EndProcedure

Procedure MatrixRotateY(*Matrix.stMatrix, Angle.f)
  Protected Help.f, tsin.f, tcos.f
  tsin = Sin(Angle) 
  tcos = Cos(Angle) 
  
  Help = *Matrix\m[0] * tcos - *Matrix\m[8] * tsin
  *Matrix\m[8] = *Matrix\m[0] * tsin + *Matrix\m[8] * tcos
  *Matrix\m[0] = Help 
  
  Help = *Matrix\m[1] * tcos - *Matrix\m[9] * tsin
  *Matrix\m[9] = *Matrix\m[1] * tsin + *Matrix\m[9] * tcos
  *Matrix\m[1] = Help 
  
  Help = *Matrix\m[2] * tcos - *Matrix\m[10] * tsin
  *Matrix\m[10] = *Matrix\m[2] * tsin + *Matrix\m[10] * tcos
  *Matrix\m[2] = Help 
EndProcedure

Procedure MatrixVectorRotation(*out.stVector, *Vector.stVector, *Matrix.stMatrix)
  Protected TempX.f, TempY.f
  TempX = *Matrix\m[0] * *Vector\x + *Matrix\m[1] * *Vector\y + *Matrix\m[2] * *Vector\z 
  TempY = *Matrix\m[4] * *Vector\x + *Matrix\m[5] * *Vector\y + *Matrix\m[6] * *Vector\z
  *out\z = *Matrix\m[8] * *Vector\x + *Matrix\m[9] * *Vector\y + *Matrix\m[10] * *Vector\z
  *out\x = TempX
  *out\y = TempY
EndProcedure 

Procedure Animate(Angle.f)
  Protected Help.f, tsin.f, tcos.f
  tsin = Sin(Angle) 
  tcos = Cos(Angle) 
  Camera(Main\Frame)\Pos\x = -150 * tsin
  Camera(Main\Frame)\Pos\y = 30
  Camera(Main\Frame)\Pos\z = -150 * tcos
  MatrixIdentity(Camera(Main\Frame)\Matrix)
  MatrixRotateX(Camera(Main\Frame)\Matrix, -0.3)
  MatrixRotateY(Camera(Main\Frame)\Matrix, -Angle)
  VectorSet(Camera(Main\Frame)\VScrTL, -0.5 * #RenderWidth, 0.5 * #RenderHeight, #RenderWidth)          ; instead of rotating every single ray,
  MatrixVectorRotation(Camera(Main\Frame)\VScrTL, Camera(Main\Frame)\VScrTL, Camera(Main\Frame)\Matrix) ; I'm using the camera matrix to rotate
  VectorSet(Camera(Main\Frame)\VScrB, -0.5 * #RenderWidth, -0.5 * #RenderHeight, #RenderWidth)          ; a virtual screen, and interpolate the
  MatrixVectorRotation(Camera(Main\Frame)\VScrB, Camera(Main\Frame)\VScrB, Camera(Main\Frame)\Matrix)   ; rays over it, when creating them.
  VectorSet(Camera(Main\Frame)\VScrR, 0.5 * #RenderWidth, 0.5 * #RenderHeight, #RenderWidth)            ; (just to gain some speed)
  MatrixVectorRotation(Camera(Main\Frame)\VScrR, Camera(Main\Frame)\VScrR, Camera(Main\Frame)\Matrix)
  SphereY(3) = -11 + Abs(15 * Sin(Angle * 5)) ; Bounce
EndProcedure

Procedure Init()
  glWindowOpen(#WindowWidth, #WindowHeight, "Experiments In Raytracing 2")
  RenderBuffer\pAddr[0] = AllocateMemory(#RenderWidth * #RenderHeight * 4)
  RenderBuffer\pAddr[1] = AllocateMemory(#RenderWidth * #RenderHeight * 4)
  RenderBuffer\Width = #RenderWidth
  RenderBuffer\Height = #RenderHeight
  DisplayBufferInit(#RenderWidth, #RenderHeight, #WindowWidth, #WindowHeight)
  Main\Shader_Default = @Shader_Default()
  Main\Material_Default = RTMaterialCreate(1.0, 1.0, 1.0)
  Main\Material_Background = RTMaterialCreate(0.0, 0.0, 0.0, @Shader_Sky())
  SceneCreate()
  SetUpJobs(#RenderWidth, #RenderHeight)
  Main\Running = #True
  CreateThreads()
EndProcedure

Procedure Main()
  Protected time.i, AvTime.f, NewTime.f, Event.i, Frame.i, Angle.f, Speed.f
  Init()
  While Main\Running
    If JobIdx >= numJobs
      LockMutex(mutFrameSync)
      UnlockMutex(mutFrameSync)
      DisplayBufferShow(RenderBuffer\pAddr[Main\Frame])
      JobIdx = 0
      Main\Frame = 1 - Main\Frame
      NewTime = ElapsedMilliseconds() - time
      time = ElapsedMilliseconds()
      If Abs(NewTime - AvTime) > NewTime * 0.3
        AvTime = NewTime
      Else  
        AvTime = AvTime * 0.95 + NewTime * 0.05
      EndIf  
      speed = 0.0003 * AvTime
      If speed > 0.01
        speed = 0.01
      EndIf  
      Animate(Angle) : Angle + speed
      SetWindowTitle(0, "Experiments In Raytracing 2 / FPS: " + Str(1000 / AvTime))
    EndIf  
    Delay(0)
    Event = WindowEvent()
    If Event = #PB_Event_CloseWindow
      Main\Running = #False
    EndIf  
  Wend 
  Delay(10)
EndProcedure

Main()