ray tracing tutorial
Re: ray tracing tutorial
pf shadoko, thanks for the reply.
I love using macros, pointers, assembly... but I simply can't tell what makes stuff simpler or more confusing for the regular PB user.
Have you looked at the short lived 'experiments in ray tracing' series I've done here?
I thought I've kept that simple at the time, but it seemes like I've lost everybody pretty quickly. At least I haven't seen anybody doing something with it.
That's my goal, to get a few people here to get really deep into ray tracing, so we can help, inspire and push each other to something awesome.
But I don't really know how to get there, and maybe it's unrealistic and ray tracing is just an obsession of mine.
Any thoughts/ideas?
I love using macros, pointers, assembly... but I simply can't tell what makes stuff simpler or more confusing for the regular PB user.
Have you looked at the short lived 'experiments in ray tracing' series I've done here?
I thought I've kept that simple at the time, but it seemes like I've lost everybody pretty quickly. At least I haven't seen anybody doing something with it.
That's my goal, to get a few people here to get really deep into ray tracing, so we can help, inspire and push each other to something awesome.
But I don't really know how to get there, and maybe it's unrealistic and ray tracing is just an obsession of mine.
Any thoughts/ideas?
- Michael Vogel
- Addict
- Posts: 2677
- Joined: Thu Feb 09, 2006 11:27 pm
- Contact:
Re: ray tracing tutorial
First of all, quite interesting stuff, Hades - thanks for your great code!
Not sure, if macros could make your tutorial easier to understand, I've transformed your first example which looks like this now:
Not sure, if macros could make your tutorial easier to understand, I've transformed your first example which looks like this now:
Code: Select all
; Define
Structure Ray
x.f
y.f
z.f
EndStructure
Macro SetRay(v,x_,y_,z_)
v\x=x_ : v\y=y_ : v\z=z_
EndMacro
Macro SubRay(p,p1,p2)
p\x=p1\x-p2\x
p\y=p1\y-p2\y
p\z=p1\z-p2\z
EndMacro
Macro AddRay(p,p1,p2)
p\x=p1\x+p2\x
p\y=p1\y+p2\y
p\z=p1\z+p2\z
EndMacro
Macro DivRay(p1,v)
p1\x/v : p1\y/v : p1\z/v
EndMacro
Macro MulRay(p1,v)
p1\x*v : p1\y*v : p1\z*v
EndMacro
Macro RaySqr(v)
(v\x*v\x + v\y*v\y + v\z*v\z)
EndMacro
Macro RayLen(v)
Sqr(RaySqr(v))
EndMacro
; EndDefine
ScreenWidth = 600 : ScreenHeight = 400
If OpenWindow(0, 0, 0, ScreenWidth, ScreenHeight, "Ray Tracing Tutorial 1", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
If CreateImage(0, ScreenWidth, ScreenHeight)
If StartDrawing(ImageOutput(0))
ScreenZ = ScreenWidth ; Distance between the (virtual) eye and screen = width of screen for a 45 degrees horizontal field of view
; A sphere straight ahead and as high as the screen, but slightly further away
SphereRadius.f = ScreenHeight / 2
Sphere.Ray
SetRay(Sphere,0,0,ScreenWidth * 1.2)
; Create a ray from the eye to the screen.
RayOrigin.Ray
RelPos.Ray
RayVec.Ray
; For every pixel on the screen
For ScreenY = 0 To ScreenHeight - 1
For ScreenX = 0 To ScreenWidth - 1
; Eye position
SetRay(RayOrigin,0,0,0)
; Shift 0,0 to the center of the screen and making sure that up is positive
SetRay(RayVec,ScreenX-ScreenWidth/2,ScreenHeight/2-ScreenY,ScreenZ)
; Rays have a starting point and a direction (vector). When doing math with vectors the result will
; be in multiples of the length of the vector, so we need to make sure the lenght is 1, or it gets messy.
; This is called normalizing:
RayLength.f=RayLen(RayVec)
DivRay(RayVec,RayLength)
; Let's check if the ray hits the sphere.
SubRay(RelPos,RayOrigin,Sphere)
a.f=RayVec\x*RelPos\x + RayVec\y*RelPos\y + RayVec\z*RelPos\z
a=a*a-(RaySqr(RelPos)-SphereRadius*SphereRadius)
If a > 0.0 ; Ray vector intersects sphere. (Not testing if the sphere is in front or behind the ray origin for now)
Plot(Screenx, Screeny, RGB(200,0,0))
EndIf
Next
Next
StopDrawing()
EndIf
EndIf
ImageGadget(0, 0,0, 0,0, ImageID(0))
Repeat
Event = WaitWindowEvent()
Until Event = #PB_Event_CloseWindow
EndIf
Re: ray tracing tutorial
Thanks for the nice words and for the 'transformation'.
Stuff like that can improve readability or confuse the heck out of people. And I have a hard time judging that, as I know what I'm coding (well, most of the time anyways).
So for now I will stick to the simple implementation.
Stuff like that can improve readability or confuse the heck out of people. And I have a hard time judging that, as I know what I'm coding (well, most of the time anyways).
So for now I will stick to the simple implementation.
Re: ray tracing tutorial
Refraction.
I am rather unhappy with this one as it looks crappy and I can't really explain the math.
Making it look better would require much more code for scene and material setup, and probably make it more confusing.
So I just post it as is, for now, to get it out of the way.
I am rather unhappy with this one as it looks crappy and I can't really explain the math.
Making it look better would require much more code for scene and material setup, and probably make it more confusing.
So I just post it as is, for now, to get it out of the way.
Code: Select all
ScreenWidth = 600 : ScreenHeight = 400
#Epsilon = 0.001
Structure xyz
x.f
y.f
z.f
EndStructure
Structure sphere
radius.f
x.f
y.f
z.f
EndStructure
Global NewList Sphere.sphere()
Global LightVec.xyz
Procedure.f DotProduct(*Vector1.xyz, *Vector2.xyz)
ProcedureReturn *Vector1\x * *Vector2\x + *Vector1\y * *Vector2\y + *Vector1\z * *Vector2\z
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.f IntersectSphere(*RayOrigin.xyz, *RayVec.xyz, *Sphere.sphere)
Protected RelPos.xyz, b.f, c.f, d.f, DistanceToIntersection.f
RelPos\x = *RayOrigin\x - *Sphere\x
RelPos\y = *RayOrigin\y - *Sphere\y
RelPos\z = *RayOrigin\z - *Sphere\z
b = *RayVec\x * RelPos\x + *RayVec\y * RelPos\y + *RayVec\z * RelPos\z
c = (RelPos\x * RelPos\x + RelPos\y * RelPos\y + RelPos\z * RelPos\z) - *Sphere\radius * *Sphere\radius
d = b * b - c
If d > 0.0
d = Sqr(d)
DistanceToIntersection = -b - d
If DistanceToIntersection > 0.0
ProcedureReturn DistanceToIntersection
Else
; Intersection from the inside of the sphere?
DistanceToIntersection = -b + d
If DistanceToIntersection > 0.0
ProcedureReturn DistanceToIntersection
EndIf
EndIf
EndIf
ProcedureReturn -1.0 ; no hit
EndProcedure
Procedure AddSphere(radius.f, x.f, y.f, z.f)
AddElement(Sphere())
Sphere()\radius = radius
Sphere()\x = x
Sphere()\y = y
Sphere()\z = z
EndProcedure
Procedure.f TraceRay(*RayOrigin.xyz, *RayVec.xyz, RayInsideSolid.i, RecursionDepth.i)
Protected ClosestIntersection.f, SphereHit.i, Brightness.f, A.f, Intersection.xyz, SurfaceNormal.xyz, ReflectedRayVec.xyz, RelativeIoR.f, SinT2.f, v.f, TransmittedRayVec.xyz
ClosestIntersection = 1000000.0
SphereHit = -1
ForEach Sphere()
DistanceToIntersection.f = IntersectSphere(*RayOrigin, *RayVec, Sphere())
If DistanceToIntersection > 0.0 And DistanceToIntersection < ClosestIntersection
ClosestIntersection = DistanceToIntersection
SphereHit = ListIndex(Sphere())
EndIf
Next
If SphereHit >= 0
SelectElement(Sphere(), SphereHit)
Intersection\x = *RayOrigin\x + *RayVec\x * ClosestIntersection
Intersection\y = *RayOrigin\y + *RayVec\y * ClosestIntersection
Intersection\z = *RayOrigin\z + *RayVec\z * ClosestIntersection
SurfaceNormal\x = Intersection\x - Sphere()\x
SurfaceNormal\y = Intersection\y - Sphere()\y
SurfaceNormal\z = Intersection\z - Sphere()\z
Normalize(SurfaceNormal)
Brightness = DotProduct(SurfaceNormal, LightVec)
If Brightness > 0.0
ForEach Sphere()
If ListIndex(Sphere()) <> SphereHit
If IntersectSphere(Intersection, LightVec, Sphere()) > 0.0
Brightness = 0.0 ; Shadow
Break
EndIf
EndIf
Next
Else
Brightness = 0.0
EndIf
If SphereHit = 1 And ((Int((1000 + Intersection\x) * 0.01) + Int(1000 + Intersection\z * 0.01)) & 1)
Brightness * 0.4
EndIf
If RecursionDepth < 10
If SphereHit = 0
; IoR: Index of Refraction
OutsideIoR.f = 1.0 ; Let's assume that outside of the refractive object is air or vacuum (index of refraction = 1.0)
InsideIoR.f = 1.5 ; Good value for glass
If RayInsideSolid
RelativeIoR = InsideIoR / OutsideIoR
RayInsideSolid = #False
; making sure the normal points inward when the ray comes from the inside of the sphere
SurfaceNormal\x = -SurfaceNormal\x
SurfaceNormal\y = -SurfaceNormal\y
SurfaceNormal\z = -SurfaceNormal\z
Else
RelativeIoR = OutsideIoR / InsideIoR
RayInsideSolid = #True
EndIf
A = DotProduct(*RayVec, SurfaceNormal)
SinT2 = RelativeIoR * RelativeIoR * (1.0 - A * A)
If SinT2 <= 1.0 ; if there is no total internal reflection:
v = RelativeIoR * A + Sqr(1.0 - SinT2)
TransmittedRayVec\x = RelativeIoR * *RayVec\x - v * SurfaceNormal\x
TransmittedRayVec\y = RelativeIoR * *RayVec\y - v * SurfaceNormal\y
TransmittedRayVec\z = RelativeIoR * *RayVec\z - v * SurfaceNormal\z
; Due to rounding errors the intersection may end up on the wrong side of the surface. By pushing it a short distance along the normal we can fix that.
Intersection\x - SurfaceNormal\x * #Epsilon
Intersection\y - SurfaceNormal\y * #Epsilon
Intersection\z - SurfaceNormal\z * #Epsilon
Brightness = TraceRay(Intersection, TransmittedRayVec, RayInsideSolid, RecursionDepth + 1)
EndIf
EndIf
EndIf
EndIf
ProcedureReturn Brightness
EndProcedure
If OpenWindow(0, 0, 0, ScreenWidth, ScreenHeight, "Ray Tracing Tutorial 5: Refraction", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
If CreateImage(0, ScreenWidth, ScreenHeight)
If StartDrawing(ImageOutput(0))
ScreenZ = ScreenWidth
AddSphere(ScreenHeight / 2, 0, 0, ScreenWidth * 1.2)
AddSphere(ScreenHeight * 50, 0, -ScreenHeight * 50.5, ScreenWidth * 1.2)
AddSphere(ScreenHeight / 10, -ScreenHeight / 4, ScreenHeight / 4, ScreenWidth * 1.2 - ScreenHeight / 2)
LightVec\x = 1.0
LightVec\y = 1.0
LightVec\z = -1.0
Normalize(LightVec)
For ScreenY = 0 To ScreenHeight - 1
For ScreenX = 0 To ScreenWidth - 1
Define RayOrigin.xyz, RayVec.xyz
RayOrigin\x = 0
RayOrigin\y = 0
RayOrigin\z = 0
RayVec\x = ScreenX - ScreenWidth / 2
RayVec\y = ScreenHeight / 2 - ScreenY
RayVec\z = ScreenZ
Normalize(RayVec)
Brightness.f = TraceRay(RayOrigin, RayVec, #False, 1)
Plot(ScreenX, ScreenY, RGB(Brightness * 200, 0 ,0))
Next
Next
StopDrawing()
EndIf
EndIf
ImageGadget(0, 0, 0, 0, 0, ImageID(0))
Repeat
Event = WaitWindowEvent()
Until Event = #PB_Event_CloseWindow
EndIf
End
Re: ray tracing tutorial
Impressive!
Thank you!
Thank you!
Re: ray tracing tutorial
Refined code.
- pf shadoko
- Enthusiast
- Posts: 291
- Joined: Thu Jul 09, 2015 9:07 am
Re: ray tracing tutorial
I saw all your codes (the most impressive ones of the forum), I suspected that it was a choice if you don't use macros.
the refraction is also interesting with a low value (ex: 1.02) it simulates a hollow sphere
it will also be interesting to manage the opacity of the glass.
(reduce the transparency depending on the glass thickness through which it passes)
you're not interested in classical 3D?
the refraction is also interesting with a low value (ex: 1.02) it simulates a hollow sphere
it will also be interesting to manage the opacity of the glass.
(reduce the transparency depending on the glass thickness through which it passes)
you're not interested in classical 3D?
Re: ray tracing tutorial
Thanks guys!
@pf shadoko
I'd like to implement a material system (and maybe a camera) before going any further down the quality/realism road. It's already quite messy as it is.
But that stuff will make the code much more complex, so my thought was to take a step back, and do tutorials about triangle intersection and acceleration structures first, as those would be much more useful, even outside of rendering.
About 'classical 3D' ... I've done software rasterizing in the 80's and 90's and quite a lot with Directx and OpenGL after that, but I am no game developer, so using an engine like Ogre gets me nowhere.
There are a lot of big ideas in my head that right now would work only with rasterizing. But I am currently having a hard time getting anything done, so just sticking to the thing i like most, ray tracing, and doing a few small programs with it is the only real option for me. I hope I can ease myself back into doing some serious stuff, but I'll just have to wait and see.
Seeing someone here doing the stuff that I'd like to do would be the next best thing, which is part of the reason for this tutorial.
@pf shadoko
I'd like to implement a material system (and maybe a camera) before going any further down the quality/realism road. It's already quite messy as it is.
But that stuff will make the code much more complex, so my thought was to take a step back, and do tutorials about triangle intersection and acceleration structures first, as those would be much more useful, even outside of rendering.
About 'classical 3D' ... I've done software rasterizing in the 80's and 90's and quite a lot with Directx and OpenGL after that, but I am no game developer, so using an engine like Ogre gets me nowhere.
There are a lot of big ideas in my head that right now would work only with rasterizing. But I am currently having a hard time getting anything done, so just sticking to the thing i like most, ray tracing, and doing a few small programs with it is the only real option for me. I hope I can ease myself back into doing some serious stuff, but I'll just have to wait and see.
Seeing someone here doing the stuff that I'd like to do would be the next best thing, which is part of the reason for this tutorial.
Re: ray tracing tutorial
The triangle intersection method I've chose here is relative easy to understand and very flexible.
It's main downside is, it suffers from high rounding errors, which can result in jagged edges. Switching to doubles can fix that, but I plan to show a different method without that drawback in the future.
Also it is rather slow, mainly because of how I have implemented it. But when going for clarity and flexibility speed isn't really on the table anymore.
Experiment with lines 60+ for alternate results.
It's main downside is, it suffers from high rounding errors, which can result in jagged edges. Switching to doubles can fix that, but I plan to show a different method without that drawback in the future.
Also it is rather slow, mainly because of how I have implemented it. But when going for clarity and flexibility speed isn't really on the table anymore.
Experiment with lines 60+ for alternate results.
Code: Select all
ScreenWidth = 600 : ScreenHeight = 400
Structure xyz
x.f
y.f
z.f
EndStructure
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 IntersectTriangle(*RayOrigin.xyz, *RayVec.xyz, *TriP0.xyz, *TriP1.xyz, *TriP2.xyz)
Protected TriVec0.xyz, TriVec1.xyz, TriNormal.xyz, d.f, Rel.xyz, Vec0.xyz, Vec1.xyz, Vec2.xyz, duv.f, u.f, v.f
; We need the vectors p1->p0 and p2->p0 of the triangle
TriVec0\x = *TriP0\x - *TriP1\x
TriVec0\y = *TriP0\y - *TriP1\y
TriVec0\z = *TriP0\z - *TriP1\z
TriVec1\x = *TriP0\x - *TriP2\x
TriVec1\y = *TriP0\y - *TriP2\y
TriVec1\z = *TriP0\z - *TriP2\z
; ...and it's normal
CrossProduct(TriVec0, TriVec1, TriNormal) ; Sadly we can't do TriNormal = CrossProduct(*TriVec0, *TriVec1)
d = DotProduct(TriNormal, *RayVec)
If d <> 0 ; ray is not parallel to triangle (d>0 for backface culling)
; Compute triangle position relative to ray origin
Rel\x = *TriP0\x - *RayOrigin\x
Rel\y = *TriP0\y - *RayOrigin\y
Rel\z = *TriP0\z - *RayOrigin\z
DistanceToIntersection = DotProduct(Rel, TriNormal) / d
If DistanceToIntersection > 0.0 ; Intersection (with the plane the triangle is on) 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(*RayVec, Vec0)
CrossProduct(TriVec1, *RayVec, Vec1)
u = DotProduct(Vec1, Rel) / duv
CrossProduct(*RayVec, TriVec0, Vec2)
v = DotProduct(Vec2, Rel) / duv
; Now let's determin if the intersection is inside the
If u >= 0.0 And v >= 0.0 And u + v <= 1.0 ; triangle
;If u >= 0.0 And v >= 0.0 And u <= 1.0 And v <= 1.0 ; parallelogram
;If Sqr(u*u + v*v) <= 1.0 ; ellipse
;If u >= 0.0 And v >= 0.0 And Sqr(u*u + v*v) <= 1.0 ; ellipse sector
;If u + v >= 1.0 And Sqr(u*u + v*v) <= 1.0 ; ellipse segment
;If Abs(u * 10 - Int(u * 10)) < 0.1 Or Abs(v * 10 - Int(v * 10)) < 0.1; whatever you can come up with
ProcedureReturn DistanceToIntersection
EndIf
EndIf
EndIf
ProcedureReturn -1.0 ; no hit
EndProcedure
If OpenWindow(0, 0, 0, ScreenWidth, ScreenHeight, "Ray Tracing Tutorial 6: Triangle Intersection", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
If CreateImage(0, ScreenWidth, ScreenHeight)
If StartDrawing(ImageOutput(0))
ScreenZ = ScreenWidth
; define the 3 corners of a triangle
Dist.f = ScreenWidth * 2
Size.f = ScreenHeight / 2
Define TriP0.xyz, TriP1.xyz, TriP2.xyz
TriP0\x = -Size
TriP0\y = -Size
Trip0\z = Dist
TriP1\x = Size
TriP1\y = 0 ; -Size ; for axis alignment
Trip1\z = Dist
TriP2\x = -Size
TriP2\y = Size
Trip2\z = Dist
For ScreenY = 0 To ScreenHeight - 1
For ScreenX = 0 To ScreenWidth - 1
Define RayOrigin.xyz, RayVec.xyz
RayOrigin\x = 0
RayOrigin\y = 0
RayOrigin\z = 0
RayVec\x = 0.5 + ScreenX - ScreenWidth / 2 ; added 0.5 to send the rays through the center of the pixel,
RayVec\y = 0.5 + ScreenHeight / 2 - ScreenY ; to help with jagged edges caused by rounding errors
RayVec\z = ScreenZ
Normalize(RayVec)
DistanceToIntersection = IntersectTriangle(RayOrigin, RayVec, TriP0, TriP1, TriP2)
If DistanceToIntersection > 0.0
Plot(ScreenX, ScreenY, RGB(200, 0 ,0))
EndIf
Next
Next
StopDrawing()
EndIf
EndIf
ImageGadget(0, 0, 0, 0, 0, ImageID(0))
Repeat
Event = WaitWindowEvent()
Until Event = #PB_Event_CloseWindow
EndIf
End
Re: ray tracing tutorial
BVH (Bounding Volume Hierarchy)
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.
Here is the most simple implementation I was able to come up with.
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...)
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.
Here is the most simple implementation I was able to come up with.
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
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...)