Raytracer in Julia

raytracer.jl

const delta = sqrt(eps(Float64))

immutable Vec
    x::Float64
    y::Float64
    z::Float64
end
+(a::Vec, b::Vec) = Vec(a.x+b.x, a.y+b.y, a.z+b.z)
-(a::Vec, b::Vec) = Vec(a.x-b.x, a.y-b.y, a.z-b.z)
*(a::Float64, b::Vec) = Vec(a*b.x, a*b.y, a*b.z)
*(a::Int64, b::Vec) = Vec(a*b.x, a*b.y, a*b.z)
*(a::Vec, b::Float64) = *(b,a)
dot(a::Vec, b::Vec) = (a.x*b.x + a.y*b.y + a.z*b.z)
unitize(a::Vec) = (1. / sqrt(dot(a, a)) * a)

type Ray
    orig::Vec
    dir::Vec
end

type Hit
    lambda::Float64
    normal::Vec
end

abstract Scene

immutable Sphere <: Scene
    center::Vec
    radius::Float64
end

function ray_sphere(s::Sphere, ray::Ray)
    v = s.center - ray.orig
    b = dot(v, ray.dir)
    disc = b*b - dot(v, v) + s.radius*s.radius
    if disc >= 0
        d = sqrt(disc)
        t2 = b + d
        if t2 >= 0
            t1 = b - d
            return t1 > 0 ? t1 : t2
        end
    end
    return Inf
end

function intersect(s::Sphere, i::Hit, ray::Ray)
    l = ray_sphere(s, ray)
    if l >= i.lambda
        return i
    else
        n = ray.orig + l * ray.dir - s.center
        return Hit(l, unitize(n))
    end
end

immutable Group <: Scene
    bound::Sphere
    objs::Array{Scene}
end

Group(b::Sphere) = Group(b, Scene[])

function intersect(g::Group, i::Hit, ray::Ray)
    l = ray_sphere(g.bound, ray)
    if l >= i.lambda
        return i
    else
        for j in g.objs
            i = intersect(j, i, ray)
        end
        return i
    end
end

function ray_trace(light::Vec, ray::Ray, scene::Scene)
    i = intersect(scene, Hit(Inf, Vec(0.,0.,0.)), ray)
    if i.lambda == Inf
        return 0
    end
    o = ray.orig + i.lambda*ray.dir + delta*i.normal
    g = dot(i.normal, light)
    if g >= 0
        return 0
    end
    sray = Ray(o, -1*light)
    si = intersect(scene, Hit(Inf, Vec(0.,0.,0.)), sray)
    return si.lambda == Inf ? -g : 0
end

function create(level, c::Vec, r)
    sphere = Sphere(c, r)
    if level == 1
        return sphere
    end
    group = Group(Sphere(c, 3*r))
    push!(group.objs, sphere)
    rn = 3*r/sqrt(12)
    for dz = -1:2:1
        for dx = -1:2:1
            c2 = c + Vec(dx*rn, rn, dz*rn)
            push!(group.objs, create(level-1, c2, r/2))
        end
    end
    return group
end


function Raytracer(levels, n, ss)
    scene = create(levels, Vec(0., -1., 0.), 1.)
    light = unitize(Vec(-1., -3., 2.))
    f = open("output.pgm", "w")
    write(f,string("P5\n",n," ",n,"\n",255,"\n"))
    # output = zeros(Int8, n, n)
    for y in (n-1):-1:0
        #println(y)
        for x in 0:1:(n-1)
            g = 0.
            for dx in 0:1:(ss-1)
                for dy in 0:1:(ss-1)
                    d = Vec(x+dx*1./ss-n/2., y+dy*1./ss-n/2., n*1.0)
                    ray = Ray(Vec(0., 0., -4.0), unitize(d))
                    g += ray_trace(light, ray, scene);
                end
            end
            # output[y+1,x+1] = int8(0.5 + 255.0 * g / (ss*ss))
            write(f, int8(0.5 + 255. * g / (ss*ss)))
        end
    end
    close(f)
end


# Warm up
Raytracer(2, 64, 4)
# Profile
using Profile
@sprofile Raytracer(int(ARGS[1]), int(ARGS[2]), int(ARGS[3]))
sprofile_flat()
# Time
tic()
Raytracer(int(ARGS[1]), int(ARGS[2]), int(ARGS[3]))
toc()