Code is getting long, I know, and I'm sorry. Pretty soon I'll get a problem with post length.
Code: Select all
; Experiments In Ray Tracing 3
; more complex scene in realtime (about 50 FPS on default at the start with 4 core 4GHz Skylake)
; PB 5.42, done on Windows7 64bit
; Hades, 2016
CompilerIf #PB_Compiler_Thread = 0
MessageRequester( "Sorry...", "Please, set Compiler/Compiler Options... to threadsafe before compiling" )
End
CompilerEndIf
EnableExplicit
DisableDebugger
#TileSize = 32
#MaxDist = 1000000
#FpError = 0.001
Structure Vec3
x.f : y.f : z.f
EndStructure
Structure Matrix
m.f[16]
EndStructure
Structure Color
r.f : g.f : b.f
EndStructure
Structure Material
DiffuseColor.Color
DiffuseIntensity.f
SpecularIntensity.f
Glossiness.f
AmbientIntensity.f
Reflectance.f
Transmittance.f
IndexOfRefraction.f
Shader.i
EndStructure
Global NewList Material.Material()
Structure Light ; directional light, far away, like the sun
Direction.Vec3
EndStructure
Global Light.Light
Structure Ray
Origin.Vec3
Direction.Vec3
rDirection.Vec3
length.f
PrimitiveHit.i
Weight.f
Col.Color
inside.i
EndStructure
Global PrimRayOffX.Vec3
Global PrimRayOffY.Vec3
Structure Sphere
Position.Vec3
Radius.f
*Material.Material
EndStructure
Global NewList Sphere.Sphere()
Global Dim SphereX.f(0)
Global Dim SphereY.f(0)
Global Dim SphereZ.f(0)
Global Dim SphereRadius.f(0)
Global Dim *SphereMaterial.Material(0)
Structure Camera
Pos.Vec3
Matrix.Matrix
VScrTL.Vec3
VScrR.Vec3
VScrB.Vec3
PixOffsX.Vec3
PixOffsY.Vec3
EndStructure
Global Camera.Camera
Structure RayPattern
px.i
py.i
Buffer.i
EndStructure
Global Dim RayPattern.RayPattern(#TileSize * #TileSize - 1)
Structure Job
x.i
y.i
EndStructure
Global Dim Job.Job(0)
Global numJobs
Global JobIdx
Global mutJob
Global mutFrameSync
mutJob = CreateMutex()
mutFrameSync = CreateMutex()
Structure AABB ; Axis aligned bounding box
min.f[3]
max.f[3]
EndStructure
Structure AABoxHit
Near.f
Far.f
EndStructure
Structure BVH
min.f[3]
max.f[3]
Content.i
EndStructure
Structure TravStack
*Node.BVH
*Ray.Ray
*LaRay.Ray
EndStructure
Structure BVHPrimitive
AABB.AABB
Centre.f[3]
EndStructure
Global Dim BVHPrimitive.BVHPrimitive(0)
Global Dim BVHPrimIdx(0)
Global FreeNode.i
Structure Display
hdc.i
hrc.i
WindowNr.i
DisplList.i
pBuffer.i
BufferWidth.i
BufferHeight.i
BufferPitch.i
RenderWidth.i
RenderHeight.i
VisibleWidth.i
VisibleHeight.i
RenderQuality.i
HalfSize.i
EndStructure
Global Display.Display
Structure Main
Running.i
Time.i
InfoTime.i
NumThreads.i
*BVH.BVH
Material_Default.i
Material_Background.i
Shader_Default.i
MouseLook.i
MouseX.i
MouseY.i
MovFwd.f
MovSide.f
RotX.f
RotY.f
EndStructure
Global Main.Main
Global NewList Info.s()
Global Dim KeyDown(255)
Global JinJangSin.f, JinJangCos.f
Macro mColorSet(Color, cr, cg, cb)
Color\r = cr
Color\g = cg
Color\b = cb
EndMacro
Macro mVectorSet(Vector, VX, VY, VZ)
Vector\x = VX
Vector\y = VY
Vector\z = VZ
EndMacro
Macro mNormalize(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.Ray, Depth.i, pTravStack.i)
Procedure MatrixIdentity(*Matrix.Matrix)
*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 MatrixVec3Rotation(*Result.Vec3, *Vector.Vec3, *Matrix.Matrix)
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
*Result\z = *Matrix\m[8] * *Vector\x + *Matrix\m[9] * *Vector\y + *Matrix\m[10] * *Vector\z
*Result\x = TempX
*Result\y = TempY
EndProcedure
Procedure MatrixRotateX(*Matrix.Matrix, 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.Matrix, 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 SphereCreate(x.f, y.f, z.f, radius.f, *Material.Material)
AddElement(Sphere())
mVectorSet(Sphere()\Position, x, y, z)
Sphere()\radius = radius
Sphere()\Material = *Material
EndProcedure
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())
mColorSet(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 - (Material()\Reflectance + Material()\Transmittance)
Material()\Shader = Shader
ProcedureReturn @Material()
EndProcedure
Procedure.i PrimaryRayPatternInit(px.i, py.i, StepSize.i, i.i) ; Organizes the primary Rays into squares, and the squares into squares etc.,
Protected x.i, y.i ; because thightly packed rays are more likely to follow the same
For y = 0 To 1 ; path, and hit the same geometry. Helps to increase performance.
For x = 0 To 1
If StepSize = 1
RayPattern(i)\px = px + x
RayPattern(i)\py = py + y
RayPattern(i)\Buffer = RayPattern(i)\px * 4 + RayPattern(i)\py * Display\BufferPitch
i + 1
Else
i = PrimaryRayPatternInit(px + x * StepSize, py + y * StepSize, StepSize / 2, i)
EndIf
Next
Next
ProcedureReturn i
EndProcedure
Procedure RTCameraUpdate()
mVectorSet(Camera\VScrTL, 0.5 -0.5 * Display\RenderWidth, 0.5 + 0.5 * Display\RenderHeight, Display\RenderWidth)
mVectorSet(Camera\VScrB, 0.5 -0.5 * Display\RenderWidth, 0.5 -0.5 * Display\RenderHeight, Display\RenderWidth)
mVectorSet(Camera\VScrR, 0.5 + 0.5 * Display\RenderWidth, 0.5 + 0.5 * Display\RenderHeight, Display\RenderWidth)
MatrixVec3Rotation(Camera\VScrTL, Camera\VScrTL, Camera\Matrix)
MatrixVec3Rotation(Camera\VScrB, Camera\VScrB, Camera\Matrix)
MatrixVec3Rotation(Camera\VScrR, Camera\VScrR, Camera\Matrix)
mVectorSet(Camera\PixOffsX, (Camera\VScrR\x - Camera\VScrTL\x) / Display\RenderWidth, (Camera\VScrR\y - Camera\VScrTL\y) / Display\RenderWidth, (Camera\VScrR\z - Camera\VScrTL\z) / Display\RenderWidth)
mVectorSet(Camera\PixOffsY, (Camera\VScrB\x - Camera\VScrTL\x) / Display\RenderHeight, (Camera\VScrB\y - Camera\VScrTL\y) / Display\RenderHeight, (Camera\VScrB\z - Camera\VScrTL\z) / Display\RenderHeight)
EndProcedure
Procedure DisplayBufferInit(VisibleWidth.i, VisibleHeight.i)
Protected pExtensions.i, Extensions.s
Static Tex.i
If Display\pBuffer : FreeMemory(Display\pBuffer) : EndIf
Display\VisibleWidth = VisibleWidth : Display\VisibleHeight = VisibleHeight
If Display\HalfSize
Display\RenderWidth = VisibleWidth / 2 : Display\RenderHeight = VisibleHeight / 2
Else
Display\RenderWidth = VisibleWidth : Display\RenderHeight = VisibleHeight
EndIf
Display\BufferWidth = Display\RenderWidth : Display\BufferHeight = Display\RenderHeight
Display\BufferPitch = Display\BufferWidth * 4
Display\pBuffer = AllocateMemory(Display\BufferPitch * Display\BufferHeight)
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)
If Display\DisplList <> 0 : glDeleteLists_(Display\DisplList, 1) : EndIf
Display\DisplList = glGenLists_(1)
glNewList_(Display\DisplList, #GL_COMPILE)
glBegin_(#GL_QUADS)
glTexCoord2f_(0.0, 0.0)
glVertex2i_(0,0)
glTexCoord2f_(0.0, 1.0)
glVertex2i_(0, Display\VisibleHeight)
glTexCoord2f_(1.0, 1.0)
glVertex2i_(Display\VisibleWidth, Display\VisibleHeight)
glTexCoord2f_(1.0, 0.0)
glVertex2i_(Display\VisibleWidth, 0)
glEnd_()
glEndList_()
EndProcedure
Procedure glWindowOpen(Width.i, Height.i, Name.s = "Experiments In Ray Tracing 3", flags.i = #PB_Window_SystemMenu | #PB_Window_MinimizeGadget | #PB_Window_ScreenCentered)
If WindowID(0)
FreeGadget(0) : CloseWindow(0)
EndIf
OpenWindow(0, 0, 0, Width, Height, Name, flags)
OpenGLGadget(0, 0, 0, Width, Height, #PB_OpenGL_Keyboard | #PB_OpenGL_NoFlipSynchronization | #PB_OpenGL_NoDepthBuffer)
SetActiveGadget(0)
SetGadgetAttribute(0, #PB_OpenGL_FlipBuffers, #True)
glOrtho_(0.0, Width, Height, 0.0, 0.0, 1.0)
DisplayBufferInit(Width, Height)
PrimaryRayPatternInit(0, 0, #TileSize / 2, 0)
EndProcedure
Procedure DisplayBufferShow(pBuffer.i)
glEnable_(#GL_TEXTURE_2D)
glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, Display\BufferWidth, Display\BufferHeight, 0, #GL_BGRA_EXT, #GL_UNSIGNED_BYTE, pBuffer)
glCallList_(Display\DisplList)
SetGadgetAttribute(0, #PB_OpenGL_FlipBuffers, #True)
EndProcedure
Procedure Intersect_Spheres(*Ray.Ray, SphIdx.i)
Protected MaxD.f, t1.f, idx1, *SphereX.float, *SphereY.float, *SphereZ.float, *SphereRadius.float
*SphereX.float = @SphereX(SphIdx)
*SphereY.float = @SphereY(SphIdx)
*SphereZ.float = @SphereZ(SphIdx)
*SphereRadius.float = @SphereRadius(SphIdx)
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]
!movups xmm5, [eax]
!subps xmm0, xmm5 ; RelPos\x = *Ray\Origin\x - *Sphere\Position\x
!mov ecx, [p.p_SphereY]
!movups xmm6, [ecx]
!subps xmm1, xmm6 ; RelPos\y = *Ray\Origin\y - *Sphere\Position\y
!mov eax, [p.p_SphereZ]
!movups xmm7, [eax]
!subps xmm2, xmm7 ; 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
!movaps xmm2, xmm0 ; t = #MaxDist for d <= 0 or t1 <= 0
!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
If t1 < *Ray\length
*Ray\length = t1
*Ray\PrimitiveHit = SphIdx + idx1
EndIf
EndProcedure
Procedure AABoxIntersect(*AABox.AABB, *Ray.Ray, *AABoxHit.AABoxHit)
!MOV edx, [p.p_AABox]
!MOV ecx, [p.p_Ray]
!xorps xmm6, xmm6
!movups xmm0, [edx]
!movups xmm1, [edx + 12]
!movups xmm2, [ecx]
!subps xmm0, xmm2
!subps xmm1, xmm2
!movups xmm3, [ecx + 24]
!mulps xmm0, xmm3
!mulps xmm1, xmm3
!movaps xmm4, xmm0
!movaps xmm7, xmm1
!minps xmm4, xmm1
!maxps xmm7, xmm0
!maxps xmm6, xmm4
!movaps xmm0, xmm6
!shufps xmm0, xmm6, 01000110b
!maxps xmm6, xmm0
!movaps xmm2, xmm6
!shufps xmm6, xmm6, 10110001b
!maxps xmm6, xmm2
!movaps xmm3, xmm7
!shufps xmm3, xmm7, 01000110b
!minps xmm7, xmm3
!movaps xmm5, xmm7
!shufps xmm7, xmm7, 10110001b
!minps xmm7, xmm5
!mov edx, [p.p_AABoxHit]
!movss [edx], xmm6
!movss [edx + 4], xmm7
EndProcedure
Procedure BVHTraversal(*Ray.Ray, pTravStack.i)
Protected *StackPtr.integer, *Node.BVH, rRayDir.Vec3, AABoxHit.AABoxHit, near.f, pNode.i, SpheresIdx.i
*StackPtr = pTravStack
AABoxIntersect(Main\BVH, *Ray, @AABoxHit)
If AABoxHit\far > AABoxHit\near ; the ray is hitting the root node of the BVH
*Node = Main\BVH
Repeat
If *Node\Content > 0 ; this means this is not a leaf node
near = #MaxDist
AABoxIntersect(*Node\Content + SizeOf(BVH), *Ray, @AABoxHit)
If AABoxHit\far > AABoxHit\near And *Ray\length > AABoxHit\near ; if we have a hit, we need to check what's inside...
*StackPtr + SizeOf(*StackPtr) ; ...so push it onto the stack.
*StackPtr\i = *Node\Content + SizeOf(BVH) ; (But the big deal here is, if we miss, we have just discarded
near = AABoxHit\near ; about half of the primitives in this node, with a rather
EndIf ; cheap test, and that's where the speed comes from)
AABoxIntersect(*Node\Content, *Ray, @AABoxHit)
If AABoxHit\far > AABoxHit\near And *Ray\length > AABoxHit\near
If near < AABoxHit\near ; making sure to always check the closest node first,
pNode = *StackPtr\i ; so when we get an intersection with a primitive,
*StackPtr\i = *Node\Content ; we don't have to check any nodes that are further away,
*StackPtr + SizeOf(*StackPtr) ; witch makes a huge difference when looking at the scene from the wrong angle
*StackPtr\i = pNode
Else
*StackPtr + SizeOf(*StackPtr)
*StackPtr\i = *Node\Content
EndIf
EndIf
Else ; this is a leaf node containing the actual primitives (spheres in this case),...
SpheresIdx = -*Node\Content
Intersect_Spheres(*Ray, SpheresIdx) ; ...so check if the ray does hit something
EndIf
*Node = *StackPtr\i ; using stack based traversal vs recursive should have helped with the speed,
*StackPtr - SizeOf(*StackPtr) ; but strangely didn't. But it made traversing nodes in the right order a bit easier.
Until *StackPtr < pTravStack
EndIf
EndProcedure
Procedure BVHPacketTraversal(pPrimRays.i, pTravStack.i)
Protected i, *StackPtr.TravStack, *Node.BVH, *Ray.Ray, *LaRay.Ray, AABoxHit.AABoxHit, *AABox.AABB, dx1.f, dy1.f, dz1.f, dx2.f, dy2.f, dz2.f
For i = 0 To (#TileSize * #TileSize / 64) - 1
*Ray = pPrimRays + i * 64 * SizeOf(Ray)
*LaRay = *Ray + 63 * SizeOf(Ray)
*StackPtr = pTravStack
*Node = Main\BVH
Repeat
Repeat
AABoxIntersect(*Node, *Ray, @AABoxHit)
If AABoxHit\far > AABoxHit\near And *Ray\length > AABoxHit\near
If *Node\Content > 0
While *LaRay > *Ray
AABoxIntersect(*Node, *LaRay, @AABoxHit)
If AABoxHit\far > AABoxHit\near And *LaRay\length > AABoxHit\near
Break
EndIf
*LaRay - SizeOf(Ray)
Wend
*AABox = *Node\Content
dx1 = (*AABox\min[0] - *Ray\Origin\x) * *Ray\rDirection\x
dy1 = (*AABox\min[1] - *Ray\Origin\y) * *Ray\rDirection\y
dz1 = (*AABox\min[2] - *Ray\Origin\z) * *Ray\rDirection\z
*AABox = *Node\Content + SizeOf(BVH)
dx2 = (*AABox\min[0] - *Ray\Origin\x) * *Ray\rDirection\x
dy2 = (*AABox\min[1] - *Ray\Origin\y) * *Ray\rDirection\y
dz2 = (*AABox\min[2] - *Ray\Origin\z) * *Ray\rDirection\z
If dx1 * dx1 + dy1 * dy1 + dz1 * dz1 < dx2 * dx2 + dy2 * dy2 + dz2 * dz2
*StackPtr + SizeOf(TravStack)
*StackPtr\Node = *Node\Content + SizeOf(BVH)
*StackPtr\Ray = *Ray
*StackPtr\LaRay = *LaRay
*StackPtr + SizeOf(TravStack)
*StackPtr\Node = *Node\Content
*StackPtr\Ray = *Ray
*StackPtr\LaRay = *LaRay
Else
*StackPtr + SizeOf(TravStack)
*StackPtr\Node = *Node\Content
*StackPtr\Ray = *Ray
*StackPtr\LaRay = *LaRay
*StackPtr + SizeOf(TravStack)
*StackPtr\Node = *Node\Content + SizeOf(BVH)
*StackPtr\Ray = *Ray
*StackPtr\LaRay = *LaRay
EndIf
Break
Else
Intersect_Spheres(*Ray, -*Node\Content)
EndIf
EndIf
*Ray + SizeOf(Ray)
Until *Ray > *LaRay
*Node = *StackPtr\Node
*Ray = *StackPtr\Ray
*LaRay = *StackPtr\LaRay
*StackPtr - SizeOf(TravStack)
Until *StackPtr < pTravStack
Next
EndProcedure
; The idea behind the SAH algorithm is to split a node into the two bounding boxes with the lowest Surface Area, because it increases the chance of the ray missing, finding empty spots. And empty spots don't need any tests :)
Procedure.f BVHComputeSAH(pNode.l, rArea.f, left.l, right.l, Axis.l, SplitPos.f)
Protected n.l, i.l, AABB1.AABB, AABB2.AABB, ns1.l, ns2.l, *Prim.BVHPrimitive
Protected s1.Vec3, s2.Vec3, A1.f, A2.f, SAH.f
AABB1\min[0] = #MaxDist : AABB1\min[1] = #MaxDist : AABB1\min[2] = #MaxDist
AABB1\max[0] = -#MaxDist : AABB1\max[1] = -#MaxDist : AABB1\max[2] = -#MaxDist
AABB2\min[0] = #MaxDist : AABB2\min[1] = #MaxDist : AABB2\min[2] = #MaxDist
AABB2\max[0] = -#MaxDist : AABB2\max[1] = -#MaxDist : AABB2\max[2] = -#MaxDist
For n = left To right ; fit the bounding box around the primitives. For calculating the SAH value only.
*Prim = @BVHPrimitive(BVHPrimIdx(n))
If *Prim\Centre[Axis] < SplitPos
ns1 + 1
For i = 0 To 2
If *Prim\AABB\min[i] < AABB1\min[i]
AABB1\min[i] = *Prim\AABB\min[i]
EndIf
If *Prim\AABB\max[i] > AABB1\max[i]
AABB1\max[i] = *Prim\AABB\max[i]
EndIf
Next
Else
ns2 + 1
For i = 0 To 2
If *Prim\AABB\min[i] < AABB2\min[i]
AABB2\min[i] = *Prim\AABB\min[i]
EndIf
If *Prim\AABB\max[i] > AABB2\max[i]
AABB2\max[i] = *Prim\AABB\max[i]
EndIf
Next
EndIf
Next
; compute the SAH value...
If ns1 = 0 ; ...based on the number of primitives in the nodes...
ns1 = 1000 ; (no primitives in the node is bad, so set some bad values to avoid that case)
A1 = 1000000000.0
Else
s1\x = AABB1\max[0] - AABB1\min[0] ; ...the surface area of nodes...
s1\y = AABB1\max[1] - AABB1\min[1]
s1\z = AABB1\max[2] - AABB1\min[2]
A1 = s1\x * s1\y + s1\x * s1\z + s1\z * s1\y
EndIf
If ns2 = 0
ns2 = 1000
A2 = 1000000000.0
Else
s2\x = AABB2\max[0] - AABB2\min[0]
s2\y = AABB2\max[1] - AABB2\min[1]
s2\z = AABB2\max[2] - AABB2\min[2]
A2 = s2\x * s2\y + s2\x * s2\z + s2\z * s2\y
EndIf
SAH = (A1 * ns1 + A2 * ns2) * rArea ; ...and the area of the parent node
ProcedureReturn SAH
EndProcedure
Procedure BVHSplit(pNode, left, right, depth, UseSAH.i = #True)
Protected *Node.BVH, Axis, AxL.f, Centre.f, n, i
Protected *Prim.BVHPrimitive, PivotIdx, ok, *lNode.BVH, *rNode.BVH, Split.f, s.Vec3, rArea.f, SAH.f, SAH1.f
*Node = pNode
; if there are primitives in this node...
If right - left > 0 ; right is the last primitive, left the first in this node. The naming might seem strange, but is a result of the spatial sorting performed here
AxL = *Node\max[0] - *Node\min[0] ; ...find the longest axis of this node...
If *Node\max[1] - *Node\min[1] > AxL
Axis = 1
AxL = *Node\max[1] - *Node\min[1]
EndIf
If *Node\max[2] - *Node\min[2] > AxL
Axis = 2
AxL = *Node\max[2] - *Node\min[2]
EndIf
Centre = (*Node\min[Axis] + *Node\max[Axis]) * 0.5 ; ... and do a centre split
If UseSAH ; using the SAH will significantly improve render times in most cases, but will dramatically slow down BVH creation
s\x = *Node\max[0] - *Node\min[0]
s\y = *Node\max[1] - *Node\min[1]
s\z = *Node\max[2] - *Node\min[2]
rArea = 1.0 / (s\x * s\y + s\x * s\z + s\z * s\y)
SAH = BVHComputeSAH(pNode, rArea, left, right, Axis, Centre)
If (right - left) > 30 ; (value might need adjusting) if there is a significant amount of primitives in this node...
For i = 0 To 2 ; ...along each Axis..
For n = 1 To 9 ; ...make 9 splits...
Split = *Node\min[i] + ((*Node\max[i] - *Node\min[i]) * 0.1 * n)
SAH1 = BVHComputeSAH(pNode, rArea, left, right, i, Split) ; ...and check for the best SAH value
If SAH1 < SAH
SAH = SAH1
Axis = i
Centre = Split
EndIf
Next
Next
Else
For i = 0 To 2 ; ...otherwise do a complete sweep on all axis...
For n = left To right ; ...for all primitives
*Prim = @BVHPrimitive(BVHPrimIdx(n))
Split = *Prim\Centre[i]
SAH1 = BVHComputeSAH(pNode, rArea, left, right, i, Split)
If SAH1 < SAH
SAH = SAH1
Axis = i
Centre = Split
EndIf
Next
Next
EndIf
EndIf
PivotIdx = left ; ... and sort the primitives into left/right of the split
For n = left To right
*Prim = @BVHPrimitive(BVHPrimIdx(n))
If *Prim\Centre[Axis] < Centre
Swap BVHPrimIdx(n), BVHPrimIdx(PivotIdx)
PivotIdx + 1
EndIf
Next
If PivotIdx = left
PivotIdx = left + 1 + (right - left) / 2
EndIf
If PivotIdx > right
PivotIdx = left + 1 + (right - left) / 2
EndIf
If right - left >= 4 ; making sure there are no more than 4 spheres per node (because 4 intersections are computed at the same time)
*lNode = FreeNode ; ...create the two new nodes...
*rNode = *lNode + SizeOf(BVH)
FreeNode = FreeNode + 2 * SizeOf(BVH)
*Node\Content = *lNode
*lNode\min[0] = #MaxDist : *lNode\min[1] = #MaxDist : *lNode\min[2] = #MaxDist
*lNode\max[0] = -#MaxDist : *lNode\max[1] = -#MaxDist : *lNode\max[2] = -#MaxDist
For n = left To PivotIdx - 1 ; fit the bounding box of the node around the primitives in the node
*Prim = @BVHPrimitive(BVHPrimIdx(n))
For i = 0 To 2
If *Prim\AABB\min[i] < *lNode\min[i]
*lNode\min[i] = *Prim\AABB\min[i]
EndIf
If *Prim\AABB\max[i] > *lNode\max[i]
*lNode\max[i] = *Prim\AABB\max[i]
EndIf
Next
Next
*rNode\min[0] = #MaxDist : *rNode\min[1] = #MaxDist : *rNode\min[2] = #MaxDist
*rNode\max[0] = -#MaxDist : *rNode\max[1] = -#MaxDist : *rNode\max[2] = -#MaxDist
For n = PivotIdx To right
*Prim = @BVHPrimitive(BVHPrimIdx(n))
For i = 0 To 2
If *Prim\AABB\min[i] < *rNode\min[i]
*rNode\min[i] = *Prim\AABB\min[i]
EndIf
If *Prim\AABB\max[i] > *rNode\max[i]
*rNode\max[i] = *Prim\AABB\max[i]
EndIf
Next
Next
BVHSplit(*lNode, left, PivotIdx - 1, Depth + 1, UseSAH) ; ...and recursivly test them for further splits
BVHSplit(*rNode, PivotIdx, right, Depth + 1, UseSAH)
ProcedureReturn
EndIf
EndIf
*Node\Content = -left ; using a nagative number here to mark this as a leaf
EndProcedure
Procedure BVHBuild(*BVH.BVH, numPrimitives.i)
Protected n, i
*BVH\min[0] = #MaxDist : *BVH\min[1] = #MaxDist : *BVH\min[2] = #MaxDist
*BVH\max[0] = -#MaxDist : *BVH\max[1] = -#MaxDist : *BVH\max[2] = -#MaxDist
For n = 0 To numPrimitives - 1
BVHPrimIdx(n) = n ; we need to sort the primitives into the BVH nodes. This index array will be used for that purpose
For i = 0 To 2 ; Fit the root nodes AABB tightly around the scene
If BVHPrimitive(n)\AABB\min[i] < *BVH\min[i]
*BVH\min[i] = BVHPrimitive(n)\AABB\min[i]
EndIf
If BVHPrimitive(n)\AABB\max[i] > *BVH\max[i]
*BVH\max[i] = BVHPrimitive(n)\AABB\max[i]
EndIf
Next
Next
*BVH\Content = *BVH + SizeOf(BVH)
FreeNode = *BVH + SizeOf(BVH)
BVHSplit(*BVH, 0, numPrimitives - 1, 0)
EndProcedure
Procedure BVHSetupSpheres() ; Everything related to the type of primitive is moved out of the other BVH procedures and into this one
Protected n, numSpheres.i, *BVH.BVH ; (except for the intersect call in traversal), so for example switching to triangles is much easier now.
numSpheres = ListSize(Sphere())
*BVH = AllocateMemory(2 * numSpheres * SizeOf(BVH)) ; this is the maximum size the BVH could have
ReDim BVHPrimitive(numSpheres - 1)
ReDim BVHPrimIdx(numSpheres - 1)
ReDim SphereX(numSpheres + 6)
ReDim SphereY(numSpheres + 6)
ReDim SphereZ(numSpheres + 6)
ReDim SphereRadius(numSpheres + 6)
ReDim *SphereMaterial(numSpheres + 6)
ForEach Sphere() ; precompute a bounding box around every primitive, to ease the creation of the BVH
BVHPrimitive(n)\AABB\min[0] = Sphere()\Position\x - Sphere()\radius
BVHPrimitive(n)\AABB\min[1] = Sphere()\Position\y - Sphere()\radius
BVHPrimitive(n)\AABB\min[2] = Sphere()\Position\z - Sphere()\radius
BVHPrimitive(n)\AABB\max[0] = Sphere()\Position\x + Sphere()\radius
BVHPrimitive(n)\AABB\max[1] = Sphere()\Position\y + Sphere()\radius
BVHPrimitive(n)\AABB\max[2] = Sphere()\Position\z + Sphere()\radius
BVHPrimitive(n)\Centre[0] = Sphere()\Position\x ; for a triangle one would use the centre of the bounding box
BVHPrimitive(n)\Centre[1] = Sphere()\Position\y
BVHPrimitive(n)\Centre[2] = Sphere()\Position\z
n + 1
Next
BVHBuild(*BVH, numSpheres)
For n = 0 To numSpheres - 1
SelectElement(Sphere(), BVHPrimIdx(n))
SphereX(n) = Sphere()\Position\x
SphereY(n) = Sphere()\Position\y
SphereZ(n) = Sphere()\Position\z
SphereRadius(n) = Sphere()\Radius
*SphereMaterial(n) = Sphere()\Material
Next n
FirstElement(Sphere())
For n = numSpheres To numSpheres + 6
SphereX(n) = Sphere()\Position\x
SphereY(n) = Sphere()\Position\y
SphereZ(n) = Sphere()\Position\z
SphereRadius(n) = Sphere()\Radius
*SphereMaterial(n) = Sphere()\Material
NextElement(Sphere())
Next
ProcedureReturn *BVH
EndProcedure
Procedure.i Shadow(*HitPosition.Vec3, *Normal.Vec3, *LightDirection.Vec3, LightDistance.f, pTravStack.i)
Protected ShadowRay.Ray
If Display\RenderQuality > 1
mVectorSet(ShadowRay\Origin, *HitPosition\x, *HitPosition\y, *HitPosition\z)
mVectorSet(ShadowRay\Direction, *LightDirection\x, *LightDirection\y, *LightDirection\z)
mVectorSet(ShadowRay\rDirection, 1.0 / ShadowRay\Direction\x, 1.0 / ShadowRay\Direction\y, 1.0 / ShadowRay\Direction\z)
ShadowRay\length = LightDistance
ShadowRay\PrimitiveHit = -1
BVHTraversal(ShadowRay, pTravStack)
If ShadowRay\PrimitiveHit >= 0
ProcedureReturn #True
EndIf
EndIf
ProcedureReturn #False
EndProcedure
Procedure GetSphereNormal(SphIdx.i, *Ray.Ray, *HitPosition.Vec3, *Normal.Vec3)
Protected rLength.f
mVectorSet(*Normal, *HitPosition\x - SphereX(SphIdx), *HitPosition\y - SphereY(SphIdx), *HitPosition\z - SphereZ(SphIdx))
mNormalize(*Normal, rLength.f)
If *Normal\x * *Ray\Direction\x + *Normal\y * *Ray\Direction\y + *Normal\z * *Ray\Direction\z > 0.0 ; a hit from the inside?
mVectorSet(*Normal, -*Normal\x, -*Normal\y, -*Normal\z) ; so turn the normal around
EndIf
EndProcedure
Procedure GetLightSpecularColor(*Ray.Ray, *Light.Light, *Normal.Vec3, *LightDirection.Vec3, LightNormalAngle.f, *LightSpecularColor.Color)
Protected Angle.f, ReflectedLightVector.Vec3, Intensity.f, *Material.Material
*Material = *SphereMaterial(*Ray\PrimitiveHit)
Angle = 2.0 * LightNormalAngle
mVectorSet(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)
mColorSet(*LightSpecularColor, Intensity, Intensity, Intensity)
EndIf
EndProcedure
Procedure GetLight(*Ray.Ray, *HitPosition.Vec3, *Normal.Vec3, *LightDiffuseColor.Color, *LightSpecularColor.Color, pTravStack.i)
Protected LightNormalAngle.f, LightDistance.f, ShadowRayOrigin.Vec3, *Material.Material
*Material = *SphereMaterial(*Ray\PrimitiveHit)
mColorSet(*LightDiffuseColor, 0.0, 0.0, 0.0)
mColorSet(*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
mColorSet(*LightDiffuseColor, *Material\AmbientIntensity, *Material\AmbientIntensity, *Material\AmbientIntensity)
If LightNormalAngle > *Material\AmbientIntensity
mVectorSet(ShadowRayOrigin, *HitPosition\x + #FpError * *Normal\x, *HitPosition\y + #FpError * *Normal\y, *HitPosition\z + #FpError * *Normal\z)
If Shadow(ShadowRayOrigin, *Normal, Light\Direction, LightDistance, pTravStack) = #False
mColorSet(*LightDiffuseColor, LightNormalAngle, LightNormalAngle, LightNormalAngle)
If *SphereMaterial(*Ray\PrimitiveHit)\SpecularIntensity > 0.0
GetLightSpecularColor(*Ray, Light, *Normal, Light\Direction, LightNormalAngle, *LightSpecularColor)
EndIf
EndIf
EndIf
EndProcedure
Procedure ReflectTransmit(*Ray.Ray, *HitPosition.Vec3, *Material.Material, *Normal.Vec3, Depth.i, pTravStack.i)
Protected Angle.f, RayNormalAngle.f, Fresnel.f, rRay.Ray, OutsideIoR.f, RelativeIoR.f, SinT2.f, v.f
If Depth < 10 And Display\RenderQuality > 1
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)
rRay\Weight = *Ray\Weight * (*Material\Reflectance + (1 - *Material\Reflectance) * Fresnel)
If rRay\Weight > 0.003
mVectorSet(rRay\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)
mVectorSet(rRay\Direction, *Ray\Direction\x - Angle * *Normal\x, *Ray\Direction\y - Angle * *Normal\y, *Ray\Direction\z - Angle * *Normal\z)
mVectorSet(rRay\rDirection, 1.0 / rRay\Direction\x, 1.0 / rRay\Direction\y, 1.0 / rRay\Direction\z)
rRay\length = #MaxDist
rRay\PrimitiveHit = -1
rRay\Inside = *Ray\Inside
BVHTraversal(rRay, pTravStack)
Shade(rRay, Depth + 1, pTravStack)
mColorSet(*Ray\Col, *Ray\Col\r + rRay\Col\r, *Ray\Col\g + rRay\Col\g, *Ray\Col\b + rRay\Col\b)
EndIf
EndIf
If *Material\Transmittance > 0.0
rRay\Weight = *Ray\Weight * *Material\Transmittance * (1 - Fresnel)
If rRay\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
rRay\Inside = #False
Else
RelativeIoR = OutsideIoR / *Material\IndexOfRefraction
rRay\Inside = #True
EndIf
; we have to make sure the origin of the ray is below the surface
mVectorSet(rRay\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)
mVectorSet(rRay\Direction, RelativeIoR * *Ray\Direction\x - v * *Normal\x, RelativeIoR * *Ray\Direction\y - v * *Normal\y, RelativeIoR * *Ray\Direction\z - v * *Normal\z)
mVectorSet(rRay\rDirection, 1.0 / rRay\Direction\x, 1.0 / rRay\Direction\y, 1.0 / rRay\Direction\z)
rRay\length = #MaxDist
rRay\PrimitiveHit = -1
BVHTraversal(rRay, pTravStack)
Shade(rRay, Depth + 1, pTravStack)
mColorSet(*Ray\Col, *Ray\Col\r + rRay\Col\r, *Ray\Col\g + rRay\Col\g, *Ray\Col\b + rRay\Col\b)
EndIf
EndIf
EndIf
*Ray\Weight * (1 - (*Material\Reflectance + *Material\Transmittance)) * (1 - Fresnel)
EndIf
EndProcedure
Procedure Shader_Sky(*Ray.Ray, *Material.Material, Depth.i, pTravStack.i)
Protected Angle.f, Brightness.f
Angle = 0.1 + *Ray\Direction\y * 0.9
If Angle < 0.0
Angle = 0.0
EndIf
Angle = 1.0 - Angle
Angle = Angle * Angle
Brightness = 0.2 + 0.7 * Angle * Angle
*Ray\Col\r + *Ray\Weight * Brightness
*Ray\Col\g + *Ray\Weight * Brightness
*Ray\Col\b + *Ray\Weight * 0.9
EndProcedure
Procedure Shader_Floor(*Ray.Ray, *Material.Material, Depth.i, pTravStack.i)
Protected HitPosition.Vec3, Normal.Vec3, LightDiffuseColor.Color, LightSpecularColor.Color, Brightness.f, SphIdx.i
mVectorSet(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, pTravStack)
If HitPosition\x * HitPosition\x + HitPosition\z * HitPosition\z < 110000
Brightness = 0.3
*Ray\Col\r + *Ray\Weight * (Brightness * *Material\DiffuseIntensity * LightDiffuseColor\r + LightSpecularColor\r)
*Ray\Col\g + *Ray\Weight * (Brightness * *Material\DiffuseIntensity * LightDiffuseColor\g + LightSpecularColor\g)
*Ray\Col\b + *Ray\Weight * (Brightness * *Material\DiffuseIntensity * LightDiffuseColor\b + LightSpecularColor\b)
Else
*Ray\Col\r + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\r * LightDiffuseColor\r + LightSpecularColor\r)
*Ray\Col\g + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\g * LightDiffuseColor\g + LightSpecularColor\g)
*Ray\Col\b + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\b * LightDiffuseColor\b + LightSpecularColor\b)
EndIf
EndProcedure
Procedure Shader_JinJang(*Ray.Ray, *Material.Material, Depth.i, pTravStack.i)
Protected HitPosition.Vec3, Normal.Vec3, LightDiffuseColor.Color, LightSpecularColor.Color, Brightness.f
Protected Angle.f, s.f, c.f, rx.f, ry.f, dx.f, dy.f, d.f, SphIdx.i
mVectorSet(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, pTravStack)
SphIdx = *Ray\PrimitiveHit
dx = HitPosition\x - SphereX(SphIdx) : dy = HitPosition\y - SphereY(SphIdx)
rx = JinJangCos * dx + JinJangSin * dy : ry = JinJangCos * dy - JinJangSin * dx
d = Sqr(rx * rx + ry * ry)
Brightness = 1.0
If d < 42
Brightness = 0.2
If d < 40
If ry > 0
Brightness = 1.0
EndIf
rx - 20
ry * ry
d = Sqr(rx * rx + ry)
If d < 20
Brightness = 0.2
If d > 6
Brightness = 1.0
EndIf
EndIf
rx + 40
d = Sqr(rx * rx + ry)
If d < 20
Brightness = 1.0
If d > 6
Brightness = 0.2
EndIf
EndIf
EndIf
EndIf
*Ray\Col\r + *Ray\Weight * (Brightness * *Material\DiffuseIntensity * LightDiffuseColor\r + LightSpecularColor\r)
*Ray\Col\g + *Ray\Weight * (Brightness * *Material\DiffuseIntensity * LightDiffuseColor\g + LightSpecularColor\g)
*Ray\Col\b + *Ray\Weight * (Brightness * *Material\DiffuseIntensity * LightDiffuseColor\b + LightSpecularColor\b)
EndProcedure
Procedure Shader_Default(*Ray.Ray, *Material.Material, Depth.i, pTravStack.i)
Protected HitPosition.Vec3, Normal.Vec3, LightDiffuseColor.Color, LightSpecularColor.Color
mVectorSet(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, pTravStack)
ReflectTransmit(*Ray, HitPosition, *Material, Normal, Depth, pTravStack)
*Ray\Col\r + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\r * LightDiffuseColor\r + LightSpecularColor\r)
*Ray\Col\g + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\g * LightDiffuseColor\g + LightSpecularColor\g)
*Ray\Col\b + *Ray\Weight * (*Material\DiffuseIntensity * *Material\DiffuseColor\b * LightDiffuseColor\b + LightSpecularColor\b)
EndProcedure
Procedure Shade(*Ray.Ray, Depth.i, pTravStack.i)
Protected *Material.Material
If *Ray\PrimitiveHit >= 0
*Material = *SphereMaterial(*Ray\PrimitiveHit)
Else
*Material = Main\Material_Background
EndIf
CallFunctionFast(*Material\Shader, *Ray, *Material, Depth, pTravStack)
EndProcedure
Procedure PrimaryRaysCreate(TileX.i, TileY.i, pPrimRays.i)
Protected *Ray.Ray, TP.Vec3, i.i, rLength.f
*Ray = pPrimRays
mVectorSet(TP, Camera\VScrTL\x + TileX * Camera\PixOffsX\x + TileY * Camera\PixOffsY\x, Camera\VScrTL\y + TileX * Camera\PixOffsX\y + TileY * Camera\PixOffsY\y, Camera\VScrTL\z + TileX * Camera\PixOffsX\z + TileY * Camera\PixOffsY\z)
For i = 0 To #TileSize * #TileSize - 1
mVectorSet(*Ray\Origin, Camera\Pos\x, Camera\Pos\y, Camera\Pos\z)
mVectorSet(*Ray\Direction, TP\x + RayPattern(i)\px * Camera\PixOffsX\x + RayPattern(i)\py * Camera\PixOffsY\x, TP\y + RayPattern(i)\px * Camera\PixOffsX\y + RayPattern(i)\py * Camera\PixOffsY\y, TP\z + RayPattern(i)\px * Camera\PixOffsX\z + RayPattern(i)\py * Camera\PixOffsY\z)
mNormalize(*Ray\Direction, rLength)
mVectorSet(*Ray\rDirection, 1.0 / *Ray\Direction\x, 1.0 / *Ray\Direction\y, 1.0 / *Ray\Direction\z)
*Ray\PrimitiveHit = -1
*Ray\length = #MaxDist
mColorSet(*Ray\Col, 0, 0, 0)
*Ray\Weight = 1.0
*Ray\inside = 0
*Ray + SizeOf(Ray)
Next
EndProcedure
Procedure RenderTile(TileX.i, TileY.i, pPrimRays.i, pTravStack.i)
Protected BufferStart.i, *Buffer.Long, *Material.Material, i.i, *Ray.Ray, r.f, g.f, b.f
PrimaryRaysCreate(TileX, TileY, pPrimRays)
BVHPacketTraversal(pPrimRays, pTravStack)
*Ray = pPrimRays
BufferStart = Display\pBuffer + TileX * 4 + TileY * Display\BufferPitch
If Display\RenderQuality > 0
For i = 0 To #TileSize * #TileSize - 1
*Buffer = BufferStart + RayPattern(i)\Buffer
Shade(*Ray, 1, pTravStack)
If *Ray\Col\r > 1.0 : *Ray\Col\r = 1.0 : EndIf
If *Ray\Col\g > 1.0 : *Ray\Col\g = 1.0 : EndIf
If *Ray\Col\b > 1.0 : *Ray\Col\b = 1.0 : EndIf
*Buffer\l = RGB(*Ray\Col\b * 255, *Ray\Col\g * 255, *Ray\Col\r * 255)
*Ray + SizeOf(Ray)
Next
Else
For i = 0 To #TileSize * #TileSize - 1
*Buffer = BufferStart + RayPattern(i)\Buffer
b = 10.0 / *Ray\length
g = 50.0 / *Ray\length
r = 500.0 / *Ray\length
If b > 1.0 : b = 1.0 : EndIf
If g > 1.0 : g = 1.0 : EndIf
If r > 1.0 : r = 1.0 : EndIf
*Buffer\l = RGB(b * 255, g * 255, r * 255)
*Ray + SizeOf(Ray)
Next
EndIf
EndProcedure
Procedure ThreadsDoWork(ThreadIdx)
Protected iJob, pTravStack.i, pPrimRays.i
pTravStack = AllocateMemory(ListSize(Sphere()) * SizeOf(TravStack))
pPrimRays = AllocateMemory(#TileSize * #TileSize * SizeOf(Ray))
While Main\Running
iJob = -1
LockMutex(mutJob)
If numJobs > JobIdx
If JobIdx + 1 = numJobs
LockMutex(mutFrameSync)
EndIf
iJob = JobIdx
JobIdx + 1
EndIf
UnlockMutex(mutJob)
If iJob > -1
RenderTile(Job(iJob)\x, Job(iJob)\y, pPrimRays, pTravStack)
If iJob + 1 = numJobs
UnlockMutex(mutFrameSync)
EndIf
Else
Delay(0)
EndIf
Wend
EndProcedure
Procedure ThreadsCreate()
Protected i
Main\NumThreads = CountCPUs(#PB_System_ProcessCPUs)
For i = 0 To Main\NumThreads - 1
CreateThread(@ThreadsDoWork(), i)
Next
EndProcedure
Procedure ThreadsSetUpJobs()
Protected x, y, i
numJobs = -1
ReDim Job((Display\RenderWidth * Display\RenderHeight) / (#TileSize * #TileSize))
For y = 0 To Display\RenderHeight - 1 Step #TileSize
For x = 0 To Display\RenderWidth - 1 Step #TileSize
Job(i)\x = x
Job(i)\y = y
i + 1
Next x
Next y
JobIdx = i
numJobs = i
EndProcedure
Procedure SceneMaterialLibrary(r.f, g.f, b.f, var.f, size.i)
Protected i.i
For i = 0 To size - 1
RTMaterialCreate(r + var * Cos(9 * i / size), g + var * Sin(5 * i / size), b + var * Cos(7 * i / size))
Next i
ProcedureReturn ListSize(Material()) - size
EndProcedure
Procedure SceneTorus(px.f, py.f, pz.f)
Protected i.i, n.i, *Mat0.Material, *Mat1.Material, Angle1.f, Angle2.f, s1.f, c1.f, s2.f, c2.f, P.Vec3
*Mat0 = RTMaterialCreate(0.7, 0.1, 0.05)
For i = 0 To 320 - 1
Angle1 = i / 320 * #PI * 2
s1 = Sin(Angle1)
c1 = Cos(Angle1)
For n = 0 To 3
Angle2 = n / 4 * #PI * 2
s2 = Sin(Angle2 + Angle1 * 8)
c2 = Cos(Angle2 + Angle1 * 8)
P\x = s1 * (100 + s2 * 30)
P\y = c1 * (100 + s2 * 30)
P\z = c2 * 30
*Mat1 = RTMaterialCreate(0.4 + 0.3 * c1, 0.4 + 0.3 * s1, 0.4 - 0.3 * c1)
SphereCreate(px + P\x, py + P\y , pz + P\z, 6, *Mat1)
Next n
Next i
EndProcedure
Procedure ScenePillar(px.f, py.f, pz.f, MatIdx.i)
Protected i.i, n.i, Angle1.f, Angle2.f, s1.f, c1.f, mr.i
mr = Random(1023)
For i = 0 To 320 - 1
Angle1 = i / 320 * #PI * 2
If i & 3 = 2
s1 = Sin(Angle1)
c1 = Cos(Angle1)
SphereCreate(px + s1 * 17, py , pz + c1 * 17, 3, SelectElement(Material(), MatIdx + (mr + i) & 1023))
SphereCreate(px + s1 * 17, py + 160 , pz + c1 * 17, 3, SelectElement(Material(), MatIdx + (mr + i + 300) & 1023))
EndIf
For n = 0 To 3
Angle2 = n / 4 * #PI * 2
s1 = Sin(Angle2 + Angle1 * 8)
c1 = Cos(Angle2 + Angle1 * 8)
SphereCreate(px + s1 * 15, py + 0.5 * i , pz + c1 * 15, 3, SelectElement(Material(), MatIdx + (mr + i * n * 8 + n * 350) & 1023))
Next n
Next i
EndProcedure
Procedure SceneCult(px.f, py.f, pz.f)
Protected i.i, Angle.f, *mat.Material, MatIdx.i
MatIdx = SceneMaterialLibrary(0.8, 0.8, 0.8, 0.02, 1024)
*mat = RTMaterialCreate(0.3, 0.3, 0.3)
*mat\Reflectance = 0.7
For i = 0 To 7
Angle = (i + 0.5) / 8 * #PI * 2
ScenePillar(px + 300 * Sin(Angle), py, pz + 300 * Cos(Angle), MatIdx)
SphereCreate(px + 300 * Sin(Angle), py + 185 , pz + 300 * Cos(Angle), 30, *mat)
Next i
SceneTorus(px, py + 130, pz)
*mat = RTMaterialCreate(0.9, 0.9, 0.9, @Shader_JinJang())
SphereCreate(px, py + 130, pz, 64, *mat)
*mat = RTMaterialCreate(0.0, 0.2, 0.4)
*mat\Reflectance = 0.01
*mat\Transmittance = 0.95
*mat\IndexOfRefraction = 1.5
*mat\SpecularIntensity = 0.0
SphereCreate(px, py + 40, pz - 200, 10, *mat)
EndProcedure
Procedure SceneCreate()
Protected i.i, rLength.f, MatIdx.i, *mat.Material
*mat = RTMaterialCreate(0.3, 0.8, 0.1, @Shader_Floor())
SphereCreate(0, -100000, 0, 100000, *mat)
SceneCult(0, 0, 0)
mVectorSet(Light\Direction, 1, 1, -1)
mNormalize(Light\Direction, rLength)
mVectorSet(Camera\Pos, 0, 30, -700)
Main\RotX = 0.15
Main\BVH = BVHSetupSpheres()
EndProcedure
Procedure Init(Width.i, Height.i)
Display\RenderQuality = 2
glWindowOpen(Width, Height)
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())
Main\MouseX = DesktopMouseX() : Main\MouseY = DesktopMouseY()
MatrixIdentity(Camera\Matrix)
SceneCreate()
AddElement(Info()) : Info() = "'escape' to exit"
AddElement(Info()) : Info() = "+/- to change window size"
AddElement(Info()) : Info() = "WSAD or numpad to move"
AddElement(Info()) : Info() = "left click or 'return' to toggle mouse look"
AddElement(Info()) : Info() = "H to toggle half resolution"
AddElement(Info()) : Info() = "Q to cycle render quality"
AddElement(Info()) : Info() = Str(ListSize(Sphere())) + " Spheres"
ThreadsSetUpJobs()
Main\Running = #True
ThreadsCreate()
EndProcedure
Procedure Move(Angle.f)
Protected s.f, c.f, Dir.Vec3, Angle1.f
s = Sin(Main\RotY)
c = Cos(Main\RotY)
mVectorSet(Dir, -s * Main\MovFwd - c * Main\MovSide, 0 , c * Main\MovFwd - s * Main\MovSide)
mVectorSet(Camera\Pos, Camera\Pos\x + Dir\x, Camera\Pos\y + Dir\y, Camera\Pos\z + Dir\z)
MatrixIdentity(Camera\Matrix)
MatrixRotateX(Camera\Matrix, Main\RotX)
MatrixRotateY(Camera\Matrix, Main\RotY)
Angle1 = -0.0002 * Main\Time
JinJangSin = Sin(Angle1) : JinJangCos = Cos(Angle1)
EndProcedure
Procedure StartFrame(Angle.f)
LockMutex(mutFrameSync)
UnlockMutex(mutFrameSync)
DisplayBufferShow(Display\pBuffer)
Move(Angle)
RTCameraUpdate()
JobIdx = 0
EndProcedure
Procedure CheckEvents(AvTime.f)
Protected Event.i, key, size.i, ToggleML.i
Main\movFwd = 0
Main\MovSide = 0
Repeat
Event = WindowEvent()
Select Event
Case #PB_Event_CloseWindow
Main\Running = #False
EndSelect
If EventType() = #PB_EventType_LeftButtonDown
ToggleML = 1
EndIf
key = GetGadgetAttribute(0, #PB_OpenGL_Key)
If EventType() = #PB_EventType_KeyDown
If key < 256 : KeyDown(key) = 1 : EndIf
EndIf
If EventType() = #PB_EventType_KeyUp
If key < 256 : KeyDown(key) = 0 : EndIf
Select key
Case #PB_Shortcut_Escape
Main\Running = #False
Case #PB_Shortcut_Add
size = #TileSize
Case #PB_Shortcut_Subtract
size = -#TileSize
Case #PB_Shortcut_Q
Display\RenderQuality - 1
If Display\RenderQuality < 0
Display\RenderQuality = 2
EndIf
Case #PB_Shortcut_H
Display\HalfSize = 1 - Display\HalfSize
Delay(100)
glWindowOpen(Display\VisibleWidth, Display\VisibleHeight)
ThreadsSetUpJobs()
Case #PB_Shortcut_Left
Main\InfoTime = ElapsedMilliseconds()
If Not PreviousElement(Info())
LastElement(Info())
EndIf
Case #PB_Shortcut_Right
Main\InfoTime = ElapsedMilliseconds()
If Not NextElement(Info())
FirstElement(Info())
EndIf
EndSelect
EndIf
If KeyDown(#PB_Shortcut_W) Or KeyDown(#PB_Shortcut_Pad8)
Main\MovFwd = 0.1 * AvTime
EndIf
If KeyDown(#PB_Shortcut_S) Or KeyDown(#PB_Shortcut_Pad5)
Main\MovFwd = -0.1 * AvTime
EndIf
If KeyDown(#PB_Shortcut_A) Or KeyDown(#PB_Shortcut_Pad4)
Main\MovSide = 0.1 * AvTime
EndIf
If KeyDown(#PB_Shortcut_D) Or KeyDown(#PB_Shortcut_Pad6)
Main\MovSide = -0.1 * AvTime
EndIf
If KeyDown(#PB_Shortcut_Return)
ToggleML = 1
EndIf
Until Event = 0
If ToggleML
Main\MouseLook = 1 - Main\MouseLook
If Main\MouseLook
SetGadgetAttribute(0, #PB_OpenGL_Cursor, #PB_Cursor_Invisible)
CompilerIf #PB_Compiler_OS = #PB_OS_Windows
SetCursorPos_(WindowX(0) + WindowWidth(0) / 2, WindowY(0) + WindowHeight(0) / 2)
CompilerEndIf
Main\MouseX = DesktopMouseX() : Main\MouseY = DesktopMouseY()
Else
SetGadgetAttribute(0, #PB_OpenGL_Cursor, #PB_Cursor_Default)
EndIf
EndIf
If Main\MouseLook
Main\RotY - 0.002 * (DesktopMouseX() - Main\MouseX)
Main\RotX - 0.002 * (DesktopMouseY() - Main\MouseY)
CompilerIf #PB_Compiler_OS = #PB_OS_Windows
SetCursorPos_(WindowX(0) + WindowWidth(0) / 2, WindowY(0) + WindowHeight(0) / 2)
CompilerEndIf
Main\MouseX = DesktopMouseX() : Main\MouseY = DesktopMouseY()
If Main\RotX > 0.45 * #PI
Main\RotX = 0.45 * #PI
EndIf
If Main\RotX < -0.45 * #PI
Main\RotX = -0.45 * #PI
EndIf
EndIf
If size <> 0
Delay(100)
If WindowX(0) - 3 * size > 0 And WindowY(0) - 2 * size > 0
If Display\VisibleWidth + 6 * size >= #TileSize * 2 And Display\VisibleHeight + 4 * size >= #TileSize * 2
glWindowOpen(Display\VisibleWidth + 6 * size, Display\VisibleHeight + 4 * size)
ThreadsSetUpJobs()
If Main\MouseLook
SetGadgetAttribute(0, #PB_OpenGL_Cursor, #PB_Cursor_Invisible)
CompilerIf #PB_Compiler_OS = #PB_OS_Windows
SetCursorPos_(WindowX(0) + WindowWidth(0) / 2, WindowY(0) + WindowHeight(0) / 2)
CompilerEndIf
EndIf
EndIf
EndIf
EndIf
EndProcedure
Procedure Main()
Protected time.i, AvTime.f, NewAvTime.f, Event.i, Angle.f, InfoTime.i
Init(640, 384)
Main\InfoTime = ElapsedMilliseconds()
While Main\Running
If JobIdx >= numJobs
CheckEvents(AvTime)
If ElapsedMilliseconds() - Main\InfoTime > 4000
Main\InfoTime = ElapsedMilliseconds()
If Not NextElement(Info())
FirstElement(Info())
EndIf
EndIf
StartFrame(Angle)
NewAvTime = ElapsedMilliseconds() - Main\Time
If Abs(NewAvTime - AvTime) > NewAvTime * 0.2
AvTime = NewAvTime
Else
AvTime = AvTime * 0.97 + NewAvTime * 0.03
EndIf
Angle + 0.0002 * AvTime
Main\Time = ElapsedMilliseconds()
SetWindowTitle(Display\WindowNr, "FPS: " + Right(" " + Str(1000 / AvTime), 4) + " / <crsr> " + info())
EndIf
Delay(0)
Wend
Delay(100) ; Wait for threads to finish
EndProcedure
Main()




