An acceleration structure allows to skip the vast majority of intersection tests. For large numbers of primitives this results in dramatic speedups.
While there are many types of acceleration structures, and other options for BVHs, too, a BVH with Axis Aligned Bounding Boxes and binary split seems to be most commonly used now.
Code: Select all
DisableDebugger
ScreenWidth.i = 600 : ScreenHeight.i = 400
#MaxDist = 1000000
Structure xyz
x.f : y.f : z.f
EndStructure
Structure Ray
Origin.xyz
Vec.xyz
length.f
EndStructure
Structure triangle
P.xyz[3]
EndStructure
Global NewList Triangle.triangle()
Structure aabb ; Axis aligned bounding box
min.f[3]
max.f[3]
EndStructure
Structure BVH
min.f[3]
max.f[3]
Content.i
numPrimitives.i
EndStructure
Structure BVHPrimitive
AABB.aabb
Center.f[3]
EndStructure
Global Dim BVHPrimitive.BVHPrimitive(0)
Global Dim BVHPrimIdx(0)
Global nextEmptyNode.i
Global *BVH.BVH
Procedure.f DotProduct(*Vector1.xyz, *Vector2.xyz)
ProcedureReturn *Vector1\x * *Vector2\x + *Vector1\y * *Vector2\y + *Vector1\z * *Vector2\z
EndProcedure
Procedure CrossProduct(*Vector1.xyz, *Vector2.xyz, *ResulVec.xyz)
*ResulVec\x = *Vector1\y * *Vector2\z - *Vector1\z * *Vector2\y
*ResulVec\y = *Vector1\z * *Vector2\x - *Vector1\x * *Vector2\z
*ResulVec\z = *Vector1\x * *Vector2\y - *Vector1\y * *Vector2\x
EndProcedure
Procedure Normalize(*Vector.xyz)
Protected VectorLength.f
VectorLength = Sqr(*Vector\x * *Vector\x + *Vector\y * *Vector\y + *Vector\z * *Vector\z)
*Vector\x = *Vector\x / VectorLength
*Vector\y = *Vector\y / VectorLength
*Vector\z = *Vector\z / VectorLength
EndProcedure
Procedure AddTriangle(*p0.xyz, *p1.xyz, *p2.xyz)
AddElement(Triangle())
Triangle()\P[0]\x = *p0\x
Triangle()\P[0]\y = *p0\y
Triangle()\P[0]\z = *p0\z
Triangle()\P[1]\x = *p1\x
Triangle()\P[1]\y = *p1\y
Triangle()\P[1]\z = *p1\z
Triangle()\P[2]\x = *p2\x
Triangle()\P[2]\y = *p2\y
Triangle()\P[2]\z = *p2\z
EndProcedure
Procedure IntersectTriangle(*Ray.Ray, *Triangle.triangle)
Protected TriVec0.xyz, TriVec1.xyz, TriNormal.xyz, d.f, Rel.xyz, Vec0.xyz, Vec1.xyz, Vec2.xyz, duv.f, u.f, v.f, rLength.f
; We need the vectors p1->p0 and p2->p0 of the triangle
TriVec0\x = *Triangle\P[0]\x - *Triangle\P[1]\x
TriVec0\y = *Triangle\P[0]\y - *Triangle\P[1]\y
TriVec0\z = *Triangle\P[0]\z - *Triangle\P[1]\z
TriVec1\x = *Triangle\P[0]\x - *Triangle\P[2]\x
TriVec1\y = *Triangle\P[0]\y - *Triangle\P[2]\y
TriVec1\z = *Triangle\P[0]\z - *Triangle\P[2]\z
CrossProduct(TriVec0, TriVec1, TriNormal) ; ...and it's normal
d = DotProduct(TriNormal, *Ray\Vec)
If d <> 0 ; ray is not parallel to triangle
Rel\x = *Triangle\P[0]\x - *Ray\Origin\x
Rel\y = *Triangle\P[0]\y - *Ray\Origin\y
Rel\z = *Triangle\P[0]\z - *Ray\Origin\z
DistanceToIntersection = DotProduct(Rel, TriNormal) / d
If DistanceToIntersection > 0.0 And DistanceToIntersection < *Ray\length; Intersection is in front of the ray origin
; Compute the intersection coordinates (u,v) relative to the triangle (so that 0,0 | 1,0 | 0,1 defines the triangle)
CrossProduct(TriVec0, TriVec1, Vec0)
duv = DotProduct(*Ray\Vec, Vec0)
CrossProduct(TriVec1, *Ray\Vec, Vec1)
u = DotProduct(Vec1, Rel) / duv
CrossProduct(*Ray\Vec, TriVec0, Vec2)
v = DotProduct(Vec2, Rel) / duv
If u >= 0.0 And v >= 0.0 And u + v <= 1.0 ; inside the triangle
*Ray\length = DistanceToIntersection
EndIf
EndIf
EndIf
EndProcedure
Procedure.f AABoxIntersect(*AABox.aabb, *Ray.Ray)
Protected n.l, intersNear.f, intersFar.f, Inters1.f, Inters2.f
intersNear = 0.0
Inters1 = (*AABox\min[0] - *Ray\Origin\x) / *Ray\Vec\x
Inters2 = (*AABox\max[0] - *Ray\Origin\x) / *Ray\Vec\x
If Inters2 < Inters1
Swap Inters1, Inters2
EndIf
intersFar = Inters2
If Inters1 > intersNear
intersNear = Inters1
EndIf
Inters1 = (*AABox\min[1] - *Ray\Origin\y) / *Ray\Vec\y
Inters2 = (*AABox\max[1] - *Ray\Origin\y) / *Ray\Vec\y
If Inters2 < Inters1
Swap Inters1, Inters2
EndIf
If Inters1 > intersNear
intersNear = Inters1
EndIf
If Inters2 < intersFar
intersFar = Inters2
EndIf
Inters1 = (*AABox\min[2] - *Ray\Origin\z) / *Ray\Vec\z
Inters2 = (*AABox\max[2] - *Ray\Origin\z) / *Ray\Vec\z
If Inters2 < Inters1
Swap Inters1, Inters2
EndIf
If Inters1 > intersNear
intersNear = Inters1
EndIf
If Inters2 < intersFar
intersFar = Inters2
EndIf
If intersFar > intersNear
ProcedureReturn intersNear
EndIf
ProcedureReturn #MaxDist ; no hit
EndProcedure
; Recursively intersect the ray with the BVH nodes, to find the triangles to test for intersection
Procedure BVHTraversal(*Ray.Ray, *Node.BVH)
Protected DistanceToAABB.f
If *Node\numPrimitives = 0 ; Internal node, containing more nodes
; Does the ray hit any of the nodes?
DistanceToAABB = AABoxIntersect(*Node\Content, *Ray)
If *Ray\length > DistanceToAABB
BVHTraversal(*Ray, *Node\Content)
EndIf
DistanceToAABB = AABoxIntersect(*Node\Content + SizeOf(BVH), *Ray)
If *Ray\length > DistanceToAABB
BVHTraversal(*Ray, *Node\Content + SizeOf(BVH))
EndIf
Else ; Leaf, containing triangles
For i = 0 To *Node\numPrimitives - 1
IntersectTriangle(*Ray, SelectElement(Triangle(), BVHPrimIdx(*Node\Content + i)))
Next i
EndIf
EndProcedure
; Find the smallest Axis Aligned Bounding Box containing the triangles for the node
Procedure BVHFitAABB(*Node.BVH, FirstPrimitive.i, LastPrimitive.i)
Protected i.i, Axis.i, *BVHPrim.BVHPrimitive
*Node\min[0] = #MaxDist : *Node\min[1] = #MaxDist : *Node\min[2] = #MaxDist
*Node\max[0] = -#MaxDist : *Node\max[1] = -#MaxDist : *Node\max[2] = -#MaxDist
For i = FirstPrimitive To LastPrimitive
*BVHPrim = @BVHPrimitive(BVHPrimIdx(i))
For Axis = 0 To 2
If *BVHPrim\AABB\min[Axis] < *Node\min[Axis]
*Node\min[Axis] = *BVHPrim\AABB\min[Axis]
EndIf
If *BVHPrim\AABB\max[Axis] > *Node\max[Axis]
*Node\max[Axis] = *BVHPrim\AABB\max[Axis]
EndIf
Next
Next
EndProcedure
; Recursively split the node into two new nodes, sort the triangles into those nodes and compute an Axis Aligned Bounding Box around each group,
; until it is no longer beneficial. Then store an index to the first triangle and the number of triangles in the node (turning it into a leaf)
Procedure BVHCreate(*Node.BVH, first, last)
Protected Axis.i, maxLength.f, Center.f, i.i, Pivot.i, *lowNode.BVH, *highNode.BVH
If last - first > 3 ; Are there enough triangles in the mode that splitting it is worth it?
; Find the longest axis of this node...
maxLength = *Node\max[0] - *Node\min[0]
If *Node\max[1] - *Node\min[1] > maxLength
Axis = 1
maxLength = *Node\max[1] - *Node\min[1]
EndIf
If *Node\max[2] - *Node\min[2] > maxLength
Axis = 2
EndIf
; ... do a Center split
Center = (*Node\min[Axis] + *Node\max[Axis]) * 0.5
; ... and sort the primitives into the two new nodes
Pivot = first
For i = first To last
; instead of sorting the primitives directly, which would be slow and cause other problems, let's use an array of indices
If BVHPrimitive(BVHPrimIdx(i))\Center[Axis] < Center
Swap BVHPrimIdx(i), BVHPrimIdx(Pivot)
Pivot + 1
EndIf
Next
If Pivot > first ; The triangles are not bunched up, so splitting makes sense
; Create the two new nodes...
*lowNode = nextEmptyNode
*highNode = nextEmptyNode + SizeOf(BVH)
nextEmptyNode = nextEmptyNode + 2 * SizeOf(BVH)
*Node\Content = *lowNode ; We only need a pointer to the low node, because the high node is always next in memory
; ...fit the AABB of the node around the triangles
BVHFitAABB(*lowNode, first, Pivot - 1)
BVHFitAABB(*highNode, Pivot, last)
; ...and recursivly test them for further splits
BVHCreate(*lowNode, first, Pivot - 1)
BVHCreate(*highNode, Pivot, last)
ProcedureReturn
EndIf
EndIf
; No further splitting, so this becomes a leaf containing triangles
*Node\Content = first
*Node\numPrimitives = 1 + last - first
EndProcedure
Procedure BVHSetupTriangles()
Protected numTriangles, i.i, corner.i
numTriangles = ListSize(Triangle())
ReDim BVHPrimitive(numTriangles - 1)
ReDim BVHPrimIdx(numTriangles - 1)
; Compute an AABB and it's center for every triangle. Helps with BVH creation and makes it independant of the type of primitive used.
ForEach Triangle()
BVHPrimitive(i)\AABB\min[0] = #MaxDist : BVHPrimitive(i)\AABB\min[1] = #MaxDist : BVHPrimitive(i)\AABB\min[2] = #MaxDist
BVHPrimitive(i)\AABB\max[0] = -#MaxDist : BVHPrimitive(i)\AABB\max[1] = -#MaxDist : BVHPrimitive(i)\AABB\max[2] = -#MaxDist
For corner = 0 To 2
If Triangle()\P[corner]\x < BVHPrimitive(i)\AABB\min[0]
BVHPrimitive(i)\AABB\min[0] = Triangle()\P[corner]\x
EndIf
If Triangle()\P[corner]\x > BVHPrimitive(i)\AABB\max[0]
BVHPrimitive(i)\AABB\max[0] = Triangle()\P[corner]\x
EndIf
If Triangle()\P[corner]\y < BVHPrimitive(i)\AABB\min[1]
BVHPrimitive(i)\AABB\min[1] = Triangle()\P[corner]\y
EndIf
If Triangle()\P[corner]\y > BVHPrimitive(i)\AABB\max[1]
BVHPrimitive(i)\AABB\max[1] = Triangle()\P[corner]\y
EndIf
If Triangle()\P[corner]\z < BVHPrimitive(i)\AABB\min[2]
BVHPrimitive(i)\AABB\min[2] = Triangle()\P[corner]\z
EndIf
If Triangle()\P[corner]\z > BVHPrimitive(i)\AABB\max[2]
BVHPrimitive(i)\AABB\max[2] = Triangle()\P[corner]\z
EndIf
Next
BVHPrimitive(i)\Center[0] = (BVHPrimitive(i)\AABB\min[0] + BVHPrimitive(i)\AABB\max[0]) / 2
BVHPrimitive(i)\Center[1] = (BVHPrimitive(i)\AABB\min[1] + BVHPrimitive(i)\AABB\max[1]) / 2
BVHPrimitive(i)\Center[2] = (BVHPrimitive(i)\AABB\min[2] + BVHPrimitive(i)\AABB\max[2]) / 2
BVHPrimIdx(i) = i ; we need to sort the primitives into the BVH nodes. This index array will be used for that purpose
i + 1
Next
*BVH = AllocateMemory(2 * numTriangles * SizeOf(BVH)) ; this is the maximum size the BVH could have
BVHFitAABB(*BVH, 0, numTriangles - 1)
nextEmptyNode = *BVH + SizeOf(BVH)
BVHCreate(*BVH, 0, numTriangles - 1)
EndProcedure
Procedure SceneSetup()
Protected Size1.i, Size2.i, i.i, j.i, Angle1.f, Angle2.f, s1.f, c1.f, s2.f, c2.f
#Subdivs1 = 100
#Subdivs2 = 1000
Size1 = 200
Size2 = 60
Protected Dim p.xyz(#Subdivs1 - 1, #Subdivs2 - 1)
For j = 0 To #Subdivs2 - 1
Angle1 = j * 2 * #PI / #Subdivs2
s1 = Sin(Angle1)
c1 = Cos(Angle1)
For i = 0 To #Subdivs1 - 1
Angle2 = i * 2 * #PI / #Subdivs1
s2 = Sin(Angle2 + Angle1)
c2 = Cos(Angle2 + Angle1)
P(i, j)\x = s1 * (Size1 + s2 * Size2)
P(i, j)\y = c1 * (Size1 + s2 * Size2)
P(i, j)\z = 800 + c2 * Size2
Next
Next
For j = 0 To #Subdivs2 - 1
For i = 0 To #Subdivs1 - 1
AddTriangle(p(i, j), p((i + 1 + #Subdivs1) % #Subdivs1, j), p(i, (j + 1 + #Subdivs2) % #Subdivs2))
AddTriangle(p((i + 1 + #Subdivs1) % #Subdivs1, j), p((i + 1 + #Subdivs1) % #Subdivs1, (j + 1 + #Subdivs2) % #Subdivs2), p(i, (j + 1 + #Subdivs2) % #Subdivs2))
Next
Next
EndProcedure
If OpenWindow(0, 0, 0, ScreenWidth, ScreenHeight, "Ray Tracing Tutorial 7: BVH", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
If CreateImage(0, ScreenWidth, ScreenHeight)
If StartDrawing(ImageOutput(0))
SceneSetup()
BVHSetupTriangles()
Define ScreenX.i, ScreenY.i, Ray.Ray, RenderStartTime.i
RenderStartTime = ElapsedMilliseconds()
For ScreenY = 0 To ScreenHeight - 1
For ScreenX = 0 To ScreenWidth - 1
Ray\Origin\x = 0
Ray\Origin\y = 0
Ray\Origin\z = 0
Ray\Vec\x = 0.5 + ScreenX - ScreenWidth / 2
Ray\Vec\y = 0.5 + ScreenHeight / 2 - ScreenY
Ray\Vec\z = ScreenWidth
Normalize(Ray\Vec)
Ray\length = #MaxDist
BVHTraversal(Ray, *BVH)
If Ray\length < #MaxDist
Plot(ScreenX, ScreenY, RGB(2 * (880 - Ray\length), 0, 0))
EndIf
Next
Next
StopDrawing()
SetWindowTitle(0, Str(ListSize(Triangle())) + " Triangles / Render Time: " + Str(ElapsedMilliseconds() - RenderStartTime) + " ms")
EndIf
EndIf
ImageGadget(0, 0, 0, 0, 0, ImageID(0))
Repeat
Event = WaitWindowEvent()
Until Event = #PB_Event_CloseWindow
EndIf
End
For rather uniform scenes this works very well, but for more 'messy' scenes performance will suffer.
If you replace '#Subdivs1 = 100' with '#Subdivs1 = 3' in line 277, the number of triangle is down to 3%, but the render time more than doubles.
That is because long skinny triangles create large, mostly empty bounding boxes that overlap each other and don't help a lot.
But there are options to improve the creation of the BVH to deal with problems like that. (SAH, splitting triangles...)
