Interactive Exploration of Ray Tracing in Julia
Summary
Introduction
What is ray tracing?
Ray tracing is a rendering technique that can produce incredibly realistic lighting effects. Essentially, an algorithm can trace the path of light, and then simulate the way that the light interacts with the virtual objects it ultimately hits in the computer-generated world. Source
Running This Notebook
This blog post was made entirely in a Jupyter Notebook so you can download it (link to GitHub at top of post) and run it for yourself. I chose to write in Julia as opposed to Python for a number of reasons:
- Julia is compiled and runs as almost as fast as C
- Interactive and can be used in a Jupyter Notebook (the Ju in Jupyter stands for Julia after all)
- Built for parallel computing, GPU computing, and shared memory computing
- Supports mathematical symbols in code
- Rather similar to use for someone used to
numpy - I'm learning Julia right now and the best way to learn a new programming language is to use it!
I want to try using BibTeX in my blog posts (default is APA). Here is a test...
Test BibTeX citation (Rew & Davis, 1990)
Setup
{:width="70%"}
We need to model
- Light propagation
- Timestepping: useful for clouds/smoke, or in our case, an atmosphere with changing refractive index
- Event driven: traditional approach. only recompute on interaction with a surface
- Reflection
- Refraction
mutable struct Ray
v̅::Vector{Float64} # velocity vector
x̅::Vector{Float64} # position vector
# c::RGB # colour
# the alpha channel can deal with reflection/refraction splitting into two rays (recursive step)
end
Let's now define a FlatSurface. Other interesting easy to model surfaces are spheres. We can also model partial reflection and refraction.
abstract type Surface end
"""
The equation for a plane ax+by+cz+d=0
is given by the normal n̂=(a,b,c) and
a point x̅=(x,y,z)
"""
struct FlatSurface <: Surface
n̂::Vector{Float64} # normal vector
x̅::Vector{Float64} # position
# n::Float64 # refractive index
# the reflectivity and transmittivity can be calculated by the refractive index
end
struct SphereSurface <: Surface
r::Float64 # radius
x̅::Vector{Float64} # centre point
#n::Float64 # refractive index
end
And work out how to draw these surfaces
using PyPlot
function drawSurface(s::FlatSurface,bbox::Array=[0,-10,-1,10,10,10])
# since we are drawing on x-z plane, ⟂ vector is ±(-nz,0,nx)
x_perp = [-s.n̂[3],0,s.n̂[1]]
x_perp ./= norm(x_perp)
bbox_r = norm([bbox[4]-bbox[1],bbox[5]-bbox[2],bbox[6]-bbox[3]])
p₁ = s.x̅ .- bbox_r.*x_perp
p₂ = s.x̅ .+ bbox_r.*x_perp
plt.plot([p₁[1],p₂[1]],[p₁[3],p₂[3]],"-k")
plt.axis("equal")
bbox = [bbox[1],bbox[4],bbox[3],bbox[6]]
plt.axis(bbox)
end
function drawSurface(s::SphereSurface,bbox::Array=[-1,-10,-1,10,10,10])
#y = 0 # drawing on a x-z plane
plt.gcf().gca().add_artist(plt.Circle((s.x̅[1],s.x̅[3]), s.r, fill=false))
plt.axis("equal")
bbox = [bbox[1],bbox[4],bbox[3],bbox[6]]
plt.axis(bbox)
end
Now I create a reflect function that modifies the ray's direction in place
Important: The ! symbol means the function can operate on the arguments in-place. As opposed to working on a copy.
function reflect!(ray::Ray,n̂::Vector{Float64})
ray.v̅ .= ray.v̅ .- 2*dot(ray.v̅,n̂).*n̂
#x̅ = ⟶.x̅ .+ 0.001.*⟶.x̅
#return Ray(v̅, ⟶.x̅)
end
And now we have two intersect functions that work in place
Tip: This is an example of Julia’s multiple dispatch feature. It’s really powerful!
using LinearAlgebra
"""
check if ray is travelling towards a surface, and if so, reflect it
"""
function intersect!(ray::Ray,s::FlatSurface)
if dot(ray.v̅,s.n̂) < 0 && dot(ray.x̅,s.n̂) ≤ dot(s.x̅,s.n̂)
reflect!(ray,s.n̂)
# and refract
end
end
function intersect!(ray::Ray,s::SphereSurface)
r̅ = ray.x̅-s.x̅
if dot(ray.v̅,r̅) < 0 && norm(r̅) ≤ s.r
reflect!(ray,r̅/norm(r̅))
# and refract
end
end
work out how to propagate the rays
function propagate!(ray::Ray, s::Surface, Δt::Float64)
intersect!(ray,s)
ray.x̅ .+= ray.v̅*Δt # divide by refractive index
end
s = FlatSurface([0,0,1],[0,0,0])
rays = [ Ray([0.5+0.2*i,0,-1],[i,0,10]) for i ∈ -1:0.5:1 ]
Δt = 0.05 # s
T = 15 # s
drawSurface(s)
for dt in 0:Δt:T
for ray in rays
#quiver(ray.x̅[1],ray.x̅[3],ray.v̅[1],ray.v̅[3],headaxislength=0,headlength=0,width=1e-4)
scatter(ray.x̅[1],ray.x̅[3],s=5e-2)
propagate!(ray,s,Δt)
end
end
for ray in rays
quiver(ray.x̅[1],ray.x̅[3],ray.v̅[1],ray.v̅[3],width=4e-3)
end
s = SphereSurface(3,[4,0,2])
rays = [ Ray([0.55+0.2*i,0,-1],[i,0,8]) for i ∈ -1:0.5:1 ]
Δt = 0.05 # s
T = 6 # s
drawSurface(s)
for dt in 0:Δt:T
for ray in rays
#quiver(ray.x̅[1],ray.x̅[3],ray.v̅[1],ray.v̅[3],headaxislength=0,headlength=0,width=1e-4)
scatter(ray.x̅[1],ray.x̅[3],s=5e-2)
propagate!(ray,s,Δt)
end
end
for ray in rays
quiver(ray.x̅[1],ray.x̅[3],ray.v̅[1],ray.v̅[3],width=4e-3)
end
that took 4 seconds...
Warning: None of this was optimised for parallel computing!
- Now think about plotting this in 3D.
- Add transparent materials with different refractive indicies.
- Deal with coloured spheres (ray will have a RGBA value) and
- could be possible to do a BRDF.
$~$
This is a fun side project to get me understanding how to write Julia code!
using Plots
using LinearAlgebra
struct Ray
v̅::Vector{Float64} # velocity vector
x̅::Vector{Float64} # position vector
#ζ::Vector{Float64} # polarisation vector
n::Float64 # refractive index of the medium the ray is travelling in
# c::RGB # colour
# the alpha channel can deal with reflection/refraction splitting into two rays (recursive step)
end
struct Camera
f::Float64 # focal length in mm
pixels::Vector{Int64} # rows, columns in a single frame
dim::Vector{Float64} # physical dimension of total pixels in mm
x̅::Vector{Float64} # position in m
point::Vector{Float64} # direction where camera is pointing
rot::Float64 # rotation around pointing direction in rad
# or a quaternion?
end
#cam = Camera(100,[135,240],[135,240],[0,0,10],[0,0,-1],0)
cam = Camera(100,[10,16],[135,240],[0,0,10],[0,0,-1],0)
using LinearAlgebra
ray_queue = []
mid = cam.pixels./2
for i in 1:cam.pixels[1], j in 1:cam.pixels[2]
v̅ = normalize([ (i-mid[1])/cam.f,(j-mid[2])/cam.f,cam.point[3] ])
push!(ray_queue,Ray(v̅,cam.x̅,1.))
end
ray_queue
using Plots
RGBA(12,10,8,.9)
const DEPTH = 5
abstract type Object end
abstract type Sphere <: Object end
"""
The equation for a plane ax+by+cz+d=0 is given by the normal n̂=(a,b,c) and
a point x̅=(x,y,z). A mirror reflects all light rays.
"""
struct Mirror <: Object
n̂::Vector{Float64} # normal vector
x̅::Vector{Float64} # position
end
"""
I am modelling these as dielectric spheres which will reflect and refract light.
"""
struct ClearSphere <: Sphere
r::Float64 # radius
x̅::Vector{Float64} # centre point
n::Float64 # refractive index
end
struct SolidSphere <: Sphere
r::Float # radius
x̅::Vector{Float64} # centre point
c::RGB # colour
end
struct MirrorSphere <: Sphere
r::Float # radius
x̅::Vector{Float64} # centre point
end
struct SkyBox <: Sphere
r::Float # radius
x̅::Vector{Float64} # centre point
# define a colour gradient somehow
c₁::RGB # colour 1
c₂::RBG # colour 2
end
struct Miss end
struct Intersection{T<:Object}
object::T
distance::Float64
point::Vector{Float64}
end
function intersect_dist(photon::Photon,s::Mirror)
D = -dot(photon.x̅ - s.x̅,s.n̂)
D /= dot(photon.v̅,s.n̂)
return
end
function intersect(photon::Photon,s::Mirror,ϵ::Float64=3e-3)
x̅_intersect = v̅ + D*x̅_1
end
function intersect(photon::Photon,s::Sphere,ϵ::Float64=3e-3)
r̅ = photon.x̅ - s.x̅
a = dot(photon.v̅,photon.v̅)
b = 2*dot(photon.v̅,r̅)
c = dot(r̅,r̅) - s.r
Δ = b^2-4a*c
n_roots = (Δ < 0) ? 0 : 2
if Δ == 0 n_roots = 1 end
if n_roots == 0
return Miss
else
dist = (-b - √Δ)/(2a)
return Intersection(s,dist,)
end
end
function plot_object!(p, s::Mirror)
# old_xlims = xlims(p)
# old_ylims = ylims(p)
perp = [s.n̂[3], -s.n̂[1]] # perpendicular vector for plotting in x-z axis
a = wall.x̅ .+ 20.*perp
b = wall.x̅ .- 20.*perp
line = [a, b]
plot!(p, first.(line), last.(line), label="Wall")
# xlims!(p, old_xlims)
# xlims!(p, old_xlims)
end
function plot_object!(p,s::Union{SolidSphere,ClearSphere})
# old_xlims = xlims(p)
# old_ylims = ylims(p)
# xlims!(p, old_xlims)
# xlims!(p, old_xlims)
end
function drawSurface(s::FlatSurface,bbox::Array=[0,-10,-1,10,10,10])
# since we are drawing on x-z plane, ⟂ vector is ±(-nz,0,nx)
x_perp = [-s.n̂[3],0,s.n̂[1]]
x_perp ./= norm(x_perp)
bbox_r = norm([bbox[4]-bbox[1],bbox[5]-bbox[2],bbox[6]-bbox[3]])
p₁ = s.x̅ .- bbox_r.*x_perp
p₂ = s.x̅ .+ bbox_r.*x_perp
plt.plot([p₁[1],p₂[1]],[p₁[3],p₂[3]],"-k")
plt.axis("equal")
bbox = [bbox[1],bbox[4],bbox[3],bbox[6]]
plt.axis(bbox)
end
Meta.parse("x+y") |> typeof
- Rasterisation
- Ray Tracing
- Ray Marching
@time begin
using Plots
using LinearAlgebra
using Test
using ColorVectorSpace
using Images, FileIO
end
struct Ray
v̅::Vector{Float64} # velocity (direction) vector
x̅::Vector{Float64} # position vector
c::RGB # colour
n::Float64 # refractive index of the medium the ray is travelling in
#ζ::Vector{Float64} # polarisation vector
# the alpha channel can deal with reflection/refraction splitting into two rays (recursive step)
end
#=
struct Photon
"Position vector"
p::Vector{Float64}
"Direction vector"
l::Vector{Float64}
"Color associated with the photon"
c::RGB
ior::Real
end
=#
abstract type Object end
struct Miss end
struct Intersection{T<:Object}
obj::T # intersection object
dist::Float64 # distance
x̅::Vector{Float64} # position vector
end
Base.isless(a::Miss, b::Miss) = false
Base.isless(a::Miss, b::Intersection) = false
Base.isless(a::Intersection, b::Miss) = true
Base.isless(a::Intersection, b::Intersection) = a.dist < b.dist
reflect(v̅::Vector{Float64}, n̂::Vector{Float64})::Vector{Float64} = v̅ - 2*dot(v̅,n̂)*n̂
#=
normalize(ℓ₁ - 2 * dot(ℓ₁, n̂) * n̂)
function reflect!(ray::Ray,n̂::Vector{Float64})
ray.v̅ .= ray.v̅ .- 2*dot(ray.v̅,n̂).*n̂
#x̅ = ⟶.x̅ .+ 0.001.*⟶.x̅
#return Ray(v̅, ⟶.x̅)
end
=#
function refract_copy(
ℓ₁::Vector, n̂::Vector,
old_ior, new_ior
)
r = old_ior / new_ior
n̂_oriented = if -dot(ℓ₁, n̂) < 0
-n̂
else
n̂
end
c = -dot(ℓ₁, n̂_oriented)
if abs(c) > 0.999
ℓ₁
else
f = 1 - r^2 * (1 - c^2)
if f < 0
ℓ₁
else
normalize(r * ℓ₁ + (r*c - sqrt(f)) * n̂_oriented)
end
end
end
function refract(v̅::Vector{Float64},n̂::Vector{Float64},n_old::Float64,n_new::Float64)::Vector{Float64}
r = n_old / n_new
#=
v_n̂ = dot(v̅,n̂) ≥ 0 ? -n̂ : n̂
v_n = -dot(v̅,v_n̂)
=#
v_n = -dot(v̅,n̂)
s = sign(v_n)
if abs(v_n) > 0.999 # angle of incidence essentially zero
return v̅
else
Δ = 1 - r^2 * (1 - v_n^2)
if Δ < 0
#return v̅ # no refraction (this should be total internal reflection)
return reflect(v̅,s*n̂)
else
return r*v̅ + (r*s*v_n - √Δ)*s*n̂
end
end
end
#args = (normalize([0,1,0.1]),normalize([0,0,1]), 1.5, 1.)
#refract2(args...) == refract(args...)
#refract2(args...), refract(args...)
struct Surface
r::Float64 # Reflectivity
t::Float64 # Transmission
c::RGBA # Color
n::Float64 # index of refraction
end
abstract type Object end
#abstract type Sphere <: Object end
"""
The equation for a plane ax+by+cz+d=0 is given by the normal n̂=(a,b,c) and
a point x̅=(x,y,z). A mirror reflects all light rays.
"""
struct Mirror <: Object
n̂::Vector{Float64} # normal vector
x̅::Vector{Float64} # position
end
"""
I am modelling these as dielectric spheres which will reflect and refract light.
"""
struct Sphere <: Object
r::Float64 # radius
x̅::Vector{Float64} # centre point
s::Surface # properties of the sphere
#n::Float64 # refractive index
end
struct SkyBox <: Object
r::Float64 # radius
x̅::Vector{Float64} # centre point
# define a colour gradient somehow
c::Function # colour function
end
#=
function intersection(photon::Photon, sphere::S; ϵ=1e-3) where {S <: Union{Skybox,Sphere}}
a = dot(photon.v̅,photon.v̅)
b = 2*dot(photon.v̅, photon.x̅ - sphere.x̅)
c = dot(photon.x̅ - sphere.x̅, photon.x̅ - sphere.x̅) - sphere.r^2
Δ = b^2 - 4*a*c
=#
function intersect(ray::Ray,s::S,ϵ::Float64=3e-3) where {S <: Union{SkyBox,Sphere}}
r̅ = ray.x̅ - s.x̅
a = dot(ray.v̅,ray.v̅)
b = 2*dot(ray.v̅,r̅)
c = dot(r̅,r̅) - s.r^2
Δ = b^2-4a*c
n_roots = (Δ < 0) ? 0 : 2
if Δ == 0 n_roots = 1 end
if n_roots == 0
return Miss()
else
dist = (-b - √Δ)/(2a)
return Intersection(s,dist,ray.x̅ + dist*ray.v̅)
end
end
struct Camera <: Object
f::Float64 # focal length in mm
resolution::Tuple{Int64,Int64} # rows, columns in a single frame
aperture_width::Float64 # physical dimension of total pixels in mm
x̅::Vector{Float64} # position in m
#point::Vector{Float64} # direction where camera is pointing
#rot::Float64 # rotation around pointing direction in rad
# or a quaternion?
end
function init_rays(cam::Camera)
# Physical size of the aperture/image/grid
aspect_ratio = cam.resolution[1] / cam.resolution[2]
dim = (
cam.aperture_width,
cam.aperture_width / aspect_ratio
)
# The x, y coordinates of every pixel in our image grid
# relative to the image center
xs = LinRange(-0.5* dim[1], 0.5 * dim[1], cam.resolution[1])
ys = LinRange(0.5* dim[2], -0.5 * dim[2], cam.resolution[2])
pixel_positions = [[x, y, cam.f] for y in ys, x in xs]
directions = normalize.(pixel_positions)
Ray.( directions,[cam.x̅], [zero(RGB)], [1.0])
end
test_cam = Camera(-1,(20,10),9,[0,0,100])
test_rays = init_rays(test_cam)
function closest_hit(ray::Ray, objects::Vector{<:Object})
hits = intersect.([ray], objects)
return minimum(hits)
end
# +(::RGB{FixedPointNumbers.Normed{UInt8,8}}, ::RGB{Float64})
#isa(zero(RGB),RGB{Float64})
function gradient_skybox_color(x̅, skybox)
r = skybox.r
c = zero(RGB)
#@show norm(x̅),x̅
if x̅[1] < r && x̅[1] > -r
c += RGB((x̅[1]+r)/(2.0*r), 0, 0)
end
if x̅[2] < r && x̅[2] > -r
c += RGB(0,0,(x̅[2]+r)/(2.0*r))
end
if x̅[3] < r && x̅[3] > -r
c += RGB(0,(x̅[3]+r)/(2.0*r), 0)
end
return c
end
function interact(ray::Ray, hit::Intersection{SkyBox}, ::Any, ::Any)
ray_colour = hit.obj.c(hit.x̅, hit.obj)
return Ray(hit.x̅, ray.v̅, ray_colour, ray.n)
end
interact(ray::Ray, ::Miss, ::Any, ::Any) = ray
function interact(ray::Ray,hit::Intersection{Sphere},::Any,::Any)
return Ray(ray.x̅, reflect(ray.v̅,normalize(ray.x̅-hit.x̅)),ray.c,ray.n)
end
#reflect(v̅::Vector{Float64}, n̂::Vector{Float64})::Vector{Float64} = v̅ - 2*dot(v̅,n̂)*n̂
function step_ray(ray::Ray, objs::Vector{T}, num_intersects) where {T <: Object}
if num_intersects == 0
ray
else
hit = closest_hit(ray, objs)
interact(ray, hit, num_intersects, objs)
end
end
extract_colors(rays) = map(ray -> ray.c, rays)
extract_colors(test_rays)
function ray_trace(objs::Vector{T}, cam::Camera, num_intersects::Int64 = 4) where {T <: Object}
rays = init_rays(cam)
new_rays = step_ray.(rays, [objs],[num_intersects])
extract_colors(new_rays)
end
struct Sphere <: Object
r::Float64 # radius
x̅::Vector{Float64} # centre point
s::Surface # properties of the sphere
#n::Float64 # refractive index
end
Sphere(20., [0.,0.,-25.], Surface(1.0, 0.0, RGBA(1,1,1,0.0), 1.5))
sky = SkyBox(1000, [0.0, 0.0, 0.0], gradient_skybox_color)
main_scene = [
sky,
Sphere(20., [0.,0.,-25.], Surface(1.0, 0.0, RGBA(1,1,1,0.0), 1.5)),
Sphere(20., [0.,50.,-100.], Surface(0.0, 1.0, RGBA(0,0,0,0.0), 1.0)),
Sphere(20., [-50.,0.,-25.], Surface(0.0, 0.0, RGBA(0, .3, .8, 1), 1.0)),
Sphere(20., [30., 25., -60.], Surface(0.0, 0.75, RGBA(1,0,0,0.25), 1.5)),
Sphere(20., [50., 0., -25.], Surface(0.5, 0.0, RGBA(.1,.9,.1,0.5), 1.5)),
Sphere(20., [-30., 25., -60.], Surface(0.5, 0.5, RGBA(1,1,1,0), 1.5)),
]
sky.r
function gradient_skybox_color2(x̅, skybox)
r = skybox.r
c = RGB(0.5,0.5,0.5)
return c
end
basic_camera = Camera(-5,(300,200), 160, [0,20,100])
#test_cam = Camera(-10,(20,10),9,[0,0,100])
@time let
scene = [sky]
ray_trace(main_scene, basic_camera)
end
x = Array{T,1} where T<:Object
isa(main_scene,x)
Pkg.add("ImageMagick")
using Images, FileIO
earth = load(download("https://upload.wikimedia.org/wikipedia/commons/thumb/8/8f/Whole_world_-_land_and_oceans_12000.jpg/1280px-Whole_world_-_land_and_oceans_12000.jpg"))
function get_index_rational(A, x, y)
a, b = size(A)
u = clamp(floor(Int, x * (a-1)) + 1, 1, a)
v = clamp(floor(Int, y * (b-1)) + 1, 1, b)
A[u,v]
end
function image_skybox(img)
f = function(x̅, skybox)
lon = atan(-x̅[1], x̅[3])
lat = -atan(x̅[2], norm(x̅[[1,3]]))
get_index_rational(img, (lat/(π)) + .5, (lon/2π) + .5)
end
SkyBox(1000, [0.0, 0.0, 0.0], f)
end
escher_sphere = Sphere(10,[0,0,0], Surface(1.0, 0.0, RGBA(1,1,1,0.0), 1.5))
struct Sphere <: Object
r::Float64 # radius
x̅::Vector{Float64} # centre point
s::Surface # properties of the sphere
#n::Float64 # refractive index
end
escher_cam = Camera(-3,(300,300),30,[0,0,30])
#Camera(-5,(30,20), 160, [0,20,100])
earth_skybox = image_skybox(earth)
let
scene = [earth_skybox, escher_sphere]
ray_trace(scene, escher_cam, 4)
end
escher_sphere = Sphere(00,[0,0,0], Surface(1.0, 0.0, RGBA(1,1,1,0.0), 1.5))
scene = [earth_skybox, escher_sphere]
ray_trace(scene, escher_cam, 4)