Search code examples
juliaraytracing

I have written a path tracer using julia programming language but i think it is slow


I have changed my post and posted the whole of my code! Could someone tell me how can I optimize it?

import Base: *, +, -, /, ^
using Images
const Π = convert(Float64, π)

#define vector
mutable struct Vec3
    x::Float64
    y::Float64
    z::Float64
end

function +(u::Vec3, v::Vec3)
    Vec3(u.x+v.x, u.y+v.y, u.z+v.z)
end

function -(u::Vec3, v::Vec3)
    Vec3(u.x-v.x, u.y-v.y, u.z-v.z)
end

function /(u::Vec3, v::Float64)
    Vec3(u.x/v, u.y/v, u.z/v)
end

function *(u, v::Vec3)
    if typeof(u) == Float64
        Vec3(u*v.x, u*v.y, u*v.z)
    elseif typeof(u) == Vec3
        Vec3(u.x*v.x, u.y*v.y, u.z*v.z)
    end
end

function ^(u::Vec3, v::Float64)
    Vec3(u.x^v, u.y^v, u.z^v)
end

function dot(u::Vec3, v::Vec3)
    u.x*v.x + u.y*v.y + u.z*v.z
end

function normalize(u::Vec3)
    u/sqrt(dot(u,u))
end

function cross(u::Vec3, v::Vec3)
    Vec3(u.y*v.z - v.y*u.z, u.z*v.x - v.z*u.x, u.x*v.y - v.x*u.y)
end

function gamma(u::Vec3)
    Vec3(u.x^(1/2.2), u.y^(1/2.2), u.z^(1/2.2))
end

function clamp(u::Vec3)
    u.x = u.x <= 1 ? u.x : 1
    u.y = u.y <= 1 ? u.y : 1
    u.z = u.z <= 1 ? u.z : 1
    u
end

#define ray
struct Ray
    s::Vec3
    d::Vec3
end

#define planes
struct xyRect
    z; x1; x2; y1; y2::Float64
    normal; emittance; reflectance::Vec3
    isLight::Bool
end

struct xzRect
    y; x1; x2; z1; z2::Float64
    normal; emittance; reflectance::Vec3
    isLight::Bool
end

struct yzRect
    x; y1; y2; z1; z2::Float64
    normal; emittance; reflectance::Vec3
    isLight::Bool
end

#define sphere
mutable struct Sphere
    radius::Float64
    center; normal; emittance; reflectance::Vec3
    isLight::Bool
end

#define empty object
struct Empty
    normal; emittance; reflectance::Vec3
end

#define surfaces
Surfaces = Union{xyRect, xzRect, yzRect, Sphere}

#define intersection function
function intersect(surface::Surfaces, ray::Ray)
if typeof(surface) == xyRect
    t = (surface.z - ray.s.z)/ray.d.z
    if surface.x1 < ray.s.x + t*ray.d.x < surface.x2 && surface.y1 < ray.s.y + t*ray.d.y < surface.y2 && t > 0
        t
    else
        Inf
    end

elseif typeof(surface) == xzRect
    t = (surface.y - ray.s.y)/ray.d.y
    if surface.x1 < ray.s.x + t*ray.d.x < surface.x2 && surface.z1 < ray.s.z + t*ray.d.z < surface.z2 && t > 0
        t
    else
        Inf
    end

elseif typeof(surface) == yzRect
    t = (surface.x - ray.s.x)/ray.d.x
    if surface.y1 < ray.s.y + t*ray.d.y < surface.y2 && surface.z1 < ray.s.z + t*ray.d.z < surface.z2 && t > 0
        t
    else
        Inf
    end

elseif typeof(surface) == Sphere
    a = dot(ray.d, ray.d)
    b = 2dot(ray.d, ray.s - surface.center)
    c = dot(ray.s - surface.center, ray.s - surface.center) - surface.radius*surface.radius
    Δ = b*b - 4*a*c
    if Δ > 0
        Δ = sqrt(Δ)
        t1 = 0.5(-b-Δ)/a
        t2 = 0.5(-b+Δ)/a
        if t1 > 0
            surface.normal = normalize(ray.s + t1*ray.d - surface.center)
            t1
        elseif t2 > 0
            surface.normal = normalize(ray.s + t2*ray.d - surface.center)
            t2
        else
            Inf
        end
    else
        Inf
    end
end
end
#define nearest function
function nearest(surfaces::Array{Surfaces, 1}, ray::Ray, tMin::Float64)
    hitSurface = Empty(Vec3(0,0,0), Vec3(0,0,0), Vec3(0,0,0))
    for surface in surfaces
        t = intersect(surface, ray)
        if t < tMin
            tMin = t
            hitSurface = surface
        end
    end
    tMin, hitSurface
end
#cosine weighted sampling of hemisphere
function hemiRand(n::Vec3)
    ξ1 = rand()
    ξ2 = rand()
    x = cos(2π*ξ2)*sqrt(ξ1)
    y = sin(2π*ξ2)*sqrt(ξ1)
    z = sqrt(1-ξ1)
    r = normalize(Vec3(2rand()-1, 2rand()-1, 2rand()-1))
    b = cross(n,r)
    t = cross(n,b)
    Vec3(x*t.x + y*b.x + z*n.x, x*t.y + y*b.y + z*n.y, x*t.z + y*b.z + z*n.z)
end
#trace the path
function trace(surfaces::Array{Surfaces, 1}, ray::Ray, depth::Int64, maxDepth::Int64)
    if depth >= maxDepth
        return Vec3(0,0,0)
    end
    t, material = nearest(surfaces, ray, Inf)
    if typeof(material) == Empty
        return Vec3(0,0,0)
    end
    if material.isLight == true
        return material.emittance
    end
    ρ = material.reflectance
    BRDF = ρ/Π
    n = material.normal
    R = hemiRand(n)
    In = trace(surfaces, Ray(ray.s + t*ray.d, R), depth+1, maxDepth)
    return Π*BRDF*In
end
#define camera
struct Camera
    eye; v_up; N::Vec3
    fov; aspect; distance::Float64
end    
#render function
function render(surfaces::Array{Surfaces,1},camera::Camera,xRes::Int64,yRes::Int64,numSamples::Int64,maxDepth::Int64)
    n = normalize(camera.N)
    e = camera.eye
    c = e - camera.distance*n
    θ = camera.fov*(π/180)
    H = 2*camera.distance*tan(θ/2)
    W = H*camera.aspect
    u = normalize(cross(camera.v_up,n))
    v = cross(n,u)
    img = zeros(3, xRes, yRes)
    pixHeight = H/yRes
    pixWidth = W/xRes
    L = c - 0.5*W*u - 0.5*H*v
    for i=1:xRes
        for j=1:yRes
            cl = Vec3(0,0,0)
            for s=1:numSamples
                pt = L + (i-rand())*pixWidth*u + (yRes-j+rand())*pixHeight*v
                cl = cl + trace(surfaces, Ray(e, pt-e), 0, maxDepth)
            end
            cl = gamma(clamp(cl/convert(Float64, numSamples)))
            img[:,j,i] = [cl.x, cl.y, cl.z]
        end
    end
    img
end
#the scene
p1 = xzRect(1.,0.,1.,-1.,0.,Vec3(0,-1,0),Vec3(0,0,0),Vec3(0.75,0.75,0.75),false)
p2 = xzRect(0.,0.,1.,-1.,0.,Vec3(0,1,0),Vec3(0,0,0),Vec3(0.75,0.75,0.75),false)
p3 = xyRect(-1.,0.,1.,0.,1.,Vec3(0,0,1),Vec3(0,0,0),Vec3(0.75,0.75,0.75),false)
p4 = yzRect(0.,0.,1.,-1.,0.,Vec3(1,0,0),Vec3(0,0,0),Vec3(0.75,0.25,0.25),false)
p5 = yzRect(1.,0.,1.,-1.,0.,Vec3(-1,0,0),Vec3(0,0,0),Vec3(0.25,0.25,0.75),false)
p6 = xzRect(0.999,0.35,0.65,-0.65,-0.35,Vec3(0,-1,0),Vec3(18,18,18),Vec3(0,0,0),true)
s1 = Sphere(0.15,Vec3(0.3,0.15,-0.6),Vec3(0,0,0),Vec3(0,0,0),Vec3(0.75,0.75,0.75),false)
surfs = Surfaces[p1,p2,p3,p4,p5,p6,s1]
cam = Camera(Vec3(0.5,0.5,2),Vec3(0,1,0),Vec3(0,0,1),28.07,1,2)
@time image = render(surfs, cam, 400, 400, 1, 4);
colorview(RGB, image)

I need to know why my code is bad and slow. I am a beginner programmer and I don't have enough experience. My path tracer scene contains 7 objects, its maximum depth is 4, and it takes more than 2 seconds to generate an image of size 400*400. I think It shouldn't be that slow because my cpu is core i7 4770. Sorry for changing my post.


Solution

  • To start with,

    struct yzRect
    x; y1; y2; z1; z2::Float64
    normal; emittance; reflectance::Vec3
    isLight::Bool
    end
    

    ends up only applying the type to the last variable for each line:

    julia> fieldtypes(yzRect)
    (Any, Any, Any, Any, Float64, Any, Any, Vec3, Bool)
    

    so Julia will basically not know about any types in your structs which slows things down.

    Also, your Vec3 should really be immutable and you then just create new instances of it when you want to "modify" the Vector.

    There might be many more issues but those are two standing out.

    Reading through https://docs.julialang.org/en/v1/manual/performance-tips/index.html and applying the guidelines in there is strongly recommended when analyzing performance.